Gaussian Elimination and LU Decomposition¶

• Goal: solve linear equation

For simplicity we consider a square matrix $\mathbf{A} \in \mathbb{R}^{n \times n}$.

• History: the method is named after Carl Friedrich Gauss (1777–1855), although it stems from the notes of Isaac Newton. If fact, GE was known to Chinese mathematicians as early as 179 AD.

• A toy example.

What happens when we call A \ b to solve a linear equation?

Elementary operator matrix¶

• Elementary operator matrix is the identity matrix with the 0 in position $(j,k)$ replaced by $c$

• $\mathbf{E}_{jk}(c)$ is unit triangular, full rank, with inverse $\mathbf{E}_{jk}^{-1}(c) = \mathbf{E}_{jk}(-c)$.

• $\mathbf{E}_{jk}(c)$ left-multiplied to an $n \times m$ matrix $\mathbf{X}$ effectively replaces its $j$-th row $\mathbf{x}_{j\cdot}$ by $c \mathbf{x}_{k \cdot} + \mathbf{x}_{j \cdot}$

+ $2n$ flops (why?).

• Gaussian elimination applies a sequence of elementary operator matrices to transform the linear system $\mathbf{A} \mathbf{x} = \mathbf{b}$ to an upper triangular system

Column 1:

Column 2:

Gauss transformations¶

• For the $k$-th column,

• $\mathbf{M}_1, \ldots, \mathbf{M}_{n-1}$ are called the Gauss transformations.

• Gauss transformations $\mathbf{M}_k$ are unit triangular, full rank, with inverse

LU decomposition¶

• Gaussian elimination does $\mathbf{M}_{n-1} \cdots \mathbf{M}_1 \mathbf{A} = \mathbf{U}$.

• Let

• Thus we have the LU decomposition

where L is unit lower triangular and U is upper triangular. + LU decomposition exists if the principal sub-matrix $\mathbf{A}[1:k, 1:k]$ is non-singular for $k=1,\ldots,n-1$.

• If the LU decomposition exists and A is non-singular, then the LU decmpositon is unique and $$ \small \det(\mathbf{A}) = \det(\mathbf{L}) \det(\mathbf{U}) = \prod_{k=1}^n u_{kk}. $$
  • This is basically how Julia and R compute determinants.

• The whole LU algorithm is done in place, i.e., A is overwritten by L and U.

for k=1:n-1
    for i=k+1:n
        A[i, k] /= A[k, k]
        for j=k+1:n
            A[i, j] -= A[i, k] * A[k, j]
        end
    end
end

• The LU decomposition costs $2(n-1)^2 + 2(n-2)^2 + \cdots + 2 \cdot 1^2 \approx \frac 23 n^3$ flops.

• Actual implementations can differ: outer product LU (kij loop), block outer product LU (higher level-3 fraction), Crout's algorithm (jki loop).

• Given LU decomposition $\mathbf{A} = \mathbf{L} \mathbf{U}$, solving $(\mathbf{L} \mathbf{U}) \mathbf{x} = \mathbf{b}$ for one right hand side costs $2n^2$ flops:

  • One forward substitution ($n^2$ flops) to solve

  $$ \mathbf{L} \mathbf{y} = \mathbf{b} $$

  • One back substitution ($n^2$ flops) to solve

  $$ \mathbf{U} \mathbf{x} = \mathbf{y} $$

• Therefore to solve $\mathbf{A} \mathbf{x} = \mathbf{b}$ via LU, we need a total of $\frac 23 n^3 + 2n^2$ flops.

• If there are multiple right hand sides, LU only needs to be done once.

Matrix inversion¶

• For matrix inversion, there are $n$ right hand sides $\mathbf{e}_1, \ldots, \mathbf{e}_n$. Taking advantage of 0s reduces $2n^3$ flops to $\frac 43 n^3$ flops. So matrix inversion via LU costs $\frac 23 n^3 + \frac 43 n^3 = 2n^3$ flops.

• No inversion mentality:

  • Whenever you see matrix inverse, you should think in terms of solving linear equations. In short, no inv(A).

• We do not compute matrix inverse unless

  1. it is necessary to compute standard errors;
  2. number of right hand sides is much larger than $n$;
  3. $n$ is small.

Pivoting¶

• Let

Is there a solution to $\mathbf{A} \mathbf{x} = \mathbf{b}$ for an arbitrary $\mathbf{b}$? Does GE/LU work for $\mathbf{A}$?

• What if, during LU procedure, the pivot $a_{kk}$ is 0 or nearly 0 due to underflow?

Example: $$ \begin{split} 0.001x_1 + x_2 &= 1 \\ x_1 + x_2 &= 2 \end{split} $$ with solution

• With two significant digits, GE yields

or $$ \begin{split} 0.001x_1 + x_2 &= 1 \\ -1000 x_2 &= -1000 \end{split} , $$ yielding $$ \begin{bmatrix} x_1 \\ x_ 2 \end{bmatrix} = \begin{bmatrix} 0.0 \\ 1.0 \end{bmatrix} ! $$

• Solution: pivoting.

• With two significant digits, GE yields

or $$ \begin{split} x_1 + x_2 &= 2 \\ 1.0 x_2 &= 1.0 \end{split} $$

yielding $$ \begin{bmatrix} x_1 \\ x_ 2 \end{bmatrix} = \begin{bmatrix} 1.0 \\ 1.0 \end{bmatrix} . $$

• Partial pivoting: before zeroing the $k$-th column, move the row with $\max_{i=k}^n |a_{ik}|$ is into the $k$-th row (called GEPP).

• LU with partial pivoting yields

where $\mathbf{P}$ is a permutation matrix, $\mathbf{L}$ is unit lower triangular with $|\ell_{ij}| \le 1$, and $\mathbf{U}$ is upper triangular.

• Complete pivoting (GECP): Do both row and column interchanges so that the largest entry in the submatrix A[k:n, k:n] is permuted to the $(k,k)$-th entry. This yields the decomposition

where $|\ell_{ij}| \le 1$.

• GEPP is the most commonly used methods for solving general linear systems. GECP is more stable but costs more computation. Partial pivoting is stable most of times (which is partially because GEPP guarantees $|\ell_{ij}| \le 1$).

Implementation¶

• LAPACK: ?GETRF does $\mathbf{P} \mathbf{A} = \mathbf{L} \mathbf{U}$ (LU decomposition with partial pivoting) in place.

• R: solve() implicitly performs LU decomposition (wrapper of LAPACK routine DGESV). solve() allows specifying a single or multiple right hand sides. If none, it computes the matrix inverse. The matrix package contains lu() function that outputs L, U, and pivot.

• Julia:
  • LinearAlgebra.lu.
  • LinearAlgebra.lu!. In-place version. Input matrix gets overwritten with L and U.
  • Or call LAPACK wrapper function LAPACK.getrf! directly.
  • Other LU-related LAPACK wrapper functions: getrs, gesv, gesvx, trtri (inverse of triangular matrix), trtrs.

Further reading¶

• Sections II.5.2 and II.5.3, Computational Statistics by James Gentle (2010).

• Chapter 2, Applied Numerical Linear Algebra by James W. Demmel (1997).

• Chapter 3, Matrix Computation by Gene Golub and Charles Van Loan (2013).