\[ \newcommand{R}{\mathbb{R}} \newcommand{defarrow}{\quad \stackrel{\text{def}}{\Longleftrightarrow} \quad} \]

Gaussian elimination

Linear systems

Any linear system of equations \[ \begin{aligned} a_{11} x_1 + a_{12} x_2 + \cdots + a_{1n} x_n &= b_1\\ a_{21} x_1 + a_{22} x_2 + \cdots + a_{2n} x_n &= b_2\\ \vdots \qquad \vdots \qquad\quad \vdots \quad &\quad \vdots\\ a_{m1} x_1 + a_{m2} x_2 + \cdots + a_{mn} x_n &= b_m\\ \end{aligned} \] where \(a_{ij}, b_i \in \mathbb C\) for \(i = 1, \ldots,m\), \(j = 1, \ldots, n\), and \(x_1, \ldots, x_n\) are unknowns, can be written in matrix form: \[ \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n}\\ a_{21} & a_{22} & \cdots & a_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ a_{m1} & a_{m2} & \cdots & a_{mn}\\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \end{bmatrix}. \] Finding solutions \((x_1, \ldots, x_n\) of the linear system of equations is then equivalent to the following question: given a matrix \(A \in M_{m\times n}(\mathbb C)\) and a vector \(b \in \mathbb C^m\), is there a vector \(x \in \mathbb C^n\) such that \(Ax = b\)? ¹⁾
The answer depends on \(A\) and, for some \(A\), also on \(b\).

Ex.

\( \begin{bmatrix} 1 & 0 \\ 2 & 0 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \end{bmatrix} \quad\Longleftrightarrow\quad \begin{aligned} x_1 &= 1\\ 2x_1& = 1 \end{aligned} \qquad\text{ has no solution.}\)

\( \begin{bmatrix} 1 & 0 \\ 2 & 0 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \end{bmatrix} \quad\Longleftrightarrow\quad \begin{aligned} x_1 &= 1\\ 2x_1& = 2 \end{aligned} \qquad\text{ has infinitely many solutions: } x_1 = 1, x_2 \in \mathbb \R. \)

\( \begin{bmatrix} 1 & 1 \\ 2 & 0 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \end{bmatrix} \quad\Longleftrightarrow\quad \begin{aligned} x_1 + x_2 &= b_1\\ 2x_1 \:\:\quad &= b_2 \end{aligned} \qquad\text{ has a unique solution for any } b_1,b_2: x_1 = \frac{b_2}{2}, x_2 = b_1 - \frac{b_2}{2}. \)

Gaussian elimination and the row echelon form of a matrix

A matrix is in row echelon form if i) the left-most non-zero entry of each row (pivot) is strictly to the right of the left-most non-zero entry of any row above, and ii) all-zero rows are at the bottom of the matrix: \[ \begin{bmatrix} 1 & \cdots & & &\\ 0 & 2 & \cdots & &\\ 0 & 0 & 0 & 7 & \cdots\\ 0 & 0 & 0 & 0 & \cdots \end{bmatrix} \]

› Any linear system can be brought into row echelon form

Proof

If \(A = 0\) is the zero matrix, we are done.

Else, assume for simplicity that \(a_{11} \neq 0\) (If not rearrange the rows, or relabel the x_j:s, or both). To each row \((a_{i1}, \ldots, a_{in})\) below the first, add \(-\frac{a_{i1}}{a_{11}} \times \) the first row: \[ \begin{aligned} a_{11} x_1 + a_{12} x_2 + \cdots a_{1n} x_n &= b_1 \qquad -\frac{a_{i1}}{a_{11}}\\ a_{i1} x_1 + a_{i2} x_2 + \cdots a_{in} x_n &= b_2\\ \hline 0 + \bigg(a_{i2} - \frac{a_{i1}}{a_{11}} a_{12}\bigg) x_2 + \cdots \bigg(a_{in} - \frac{a_{i1} }{a_{11}}a_{1n} \bigg) x_n &= b_2 -\frac{a_{i1}}{a_{11}}b_1\\ \end{aligned} \] Then \(a_{i1} = 0\) for all \(i=2, \ldots, m\).

Now, either \(a_{ij} = 0\) for all \(i,j \geq 2\), or we can restart this procedure (looking at the matrix for indices \(i,j \geq 2\)). Since there are finitely many rows this procedure must eventually terminate, yielding a matrix \(\tilde A = (\tilde a_{ij})_{ij}\) in row echelon form.

This algorithm is called Gaussian elimination.

A neat trick is the following: write \(Ax = b\) as \(Ax = Ib\): \[ \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n}\\ a_{21} & a_{22} & \cdots & a_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ a_{m1} & a_{m2} & \cdots & a_{mn}\\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} = \begin{bmatrix} 1 & 0 & \cdots & 0\\ 0 & 1 & \cdots & 0\\ \vdots & 0 & \ddots & \ldots\\ 0 & & \cdots & 1 \end{bmatrix} \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \end{bmatrix}, \qquad I \in M_{m \times m}. \] Gaussian elimination affects only the matrices \(A\) and \(I\), not the vectors \(x\) and \(b\). We therefore introduce the \(m \times (n+m) \) augmented matrix \[ \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} & & 1 & 0 & \cdots & 0\\ a_{21} & a_{22} & \cdots & a_{2n} & & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots & & \vdots & 0 & \ddots & \ldots\\ a_{m1} & a_{m2} & \cdots & a_{mn} & & 0 & & \cdots & 1 \end{bmatrix}. \] The linear operations applied to \(A\) in Gaussian elimination are 'stored' in the augmented matrix.

Ex.
Solve \[ \underbrace{ \begin{bmatrix} 1 & 3 & 1\\ 2 & 2 & 0\\ 2 & 2 & -1 \end{bmatrix}}_{A} \underbrace{ \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}}_{x} = \underbrace{ \begin{bmatrix} 1 \\ 4 \\ 2 \end{bmatrix}}_{b} \quad\Longleftrightarrow\quad \begin{bmatrix} 1 & 3 & 1\\ 2 & 2 & 0\\ 2 & 2 & -1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 \\ 4 \\ 2 \end{bmatrix}. \] Now, performing the same row operations on the whole augmented matrix, \[ \begin{align*} \begin{bmatrix} 1 & 3 & 1 & & 1 & 0 & 0\\ 2 & 2 & 0 & & 0 & 1 & 0\\ 2 & 2 & -1 & & 0 & 0 & 1 \end{bmatrix} \quad&\Longleftrightarrow\quad \begin{bmatrix} 1 & 3 & 1 & & 1 & 0 & 0\\ 0 & -4 & -2 & & -2 & 1 & 0\\ 0 & -4 & -3 & & -2 & 0 & 1 \end{bmatrix}\Longleftrightarrow\quad \begin{bmatrix} 1 & 3 & 1 & & 1 & 0 & 0\\ 0 & -4 & -2 & & -2 & 1 & 0\\ 0 & 0 & -1 & & 0 & -1 & 1 \end{bmatrix}, \end{align*} \] we find the row echelon form \[ \underbrace{ \begin{bmatrix} 1 & 3 & 1\\ 0 & -4 & -2\\ 0 & 0 & -1 \end{bmatrix}}_{U: \text{ upper triangular}} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \underbrace{ \begin{bmatrix} 1 & 0 & 0\\ -2 & 1 & 0 \\ 0 & -1 & 1 \end{bmatrix}}_{\tilde L: \text{ lower triangular}} \begin{bmatrix} 1 \\ 4 \\ 2 \end{bmatrix}. \] Note that \(\tilde L\) describes the (linear) transformation applied to \(A\) to obtain \(U\): we have \(\tilde L A = U\).

LU-decompositions

As we shall see, if \(\tilde L A = U\), and if the inverse matrix \(\tilde L^{-1}\) exists, bringing \(\tilde L\) into row echelon row form yields \(\tilde L^{-1}\). Define \(L := \tilde L^{-1}\). Then \[ \tilde L A = U \quad \Longleftrightarrow \quad L \tilde L A = L U \quad \Longleftrightarrow \quad A = LU. \] This is the LU-decomposition of the matrix \(A\) (it does not always exist).²⁾

N.b. If the rows of \(A\) are not in correct order (for example, if \(a_{11} = 0\)), they can be rearranged by applying a permutation matrix \(P\) ³⁾: \[ \begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} \text{row } 1 \\ \text{row } 2 \\ \text{row } 3 \end{bmatrix} = \begin{bmatrix} \text{row } 2 \\ \text{row } 3 \\ \text{row } 1 \end{bmatrix}. \] This is known as an LUP-factorization: \(PA = LU\). For square matrices, an LUP-factorization always exists (but it is not necessarily unique).

Ex.
Gaussian elimination for \(\tilde L\) yields \[ \begin{bmatrix} 1 & 0 & 0 & & 1 & 0 & 0\\ -2 & 1 & 0 & & 0 & 1 & 0 \\ 0 & -1 & 1 & & 0 & 0 & 1 \end{bmatrix} \quad \Longleftrightarrow \quad \begin{bmatrix} 1 & 0 & 0 & & 1 & 0 & 0\\ 0 & 1 & 0 & & 2 & 1 & 0 \\ 0 & -1 & 1 & & 0 & 0 & 1 \end{bmatrix} \quad \Longleftrightarrow \quad \begin{bmatrix} 1 & 0 & 0 & & 1 & 0 & 0\\ 0 & 1 & 0 & & 2 & 1 & 0 \\ 0 & 0 & 1 & & 2 & 1 & 1 \end{bmatrix}, \] meaning that \[ \underbrace{ \begin{bmatrix} 1 & 0 & 0\\ 2 & 1 & 0 \\ 2 & 1 & 1 \end{bmatrix}}_{L = \tilde L^{-1}} \underbrace{ \begin{bmatrix} 1 & 0 & 0\\ -2 & 1 & 0 \\ 0 & -1 & 1 \end{bmatrix}}_{\tilde L} = \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix}. \] Hence, the LU-factorization of \(A\) is \[ \begin{bmatrix} 1 & 3 & 1\\ 2 & 2 & 0\\ 2 & 2 & -1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0\\ 2 & 1 & 0 \\ 2 & 1 & 1 \end{bmatrix} \begin{bmatrix} 1 & 3 & 1\\ 0 & -4 & -2\\ 0 & 0 & -1 \end{bmatrix}. \]

Gauss-Jordan elimination

Let \(\tilde L A = U\). The algorithm that applies Gaussian elimination (backwards) to the matrix \(U\) is called Gauss-Jordan elimination. It places zeros below and above each pivot; this is called the reduced row echelon form of the matrix \(A\). If \(A^{-1}\) exists, Gauss–Jordan elimination finds it.

Above, \(\tilde L^{-1}\) gave us \(A = \tilde L^{-1} U\) (the LU-factorization of \(A\)).
In the following, \(U^{-1}\) gives us \( U^{-1} \tilde L A = I\), so that \(U^{-1} \tilde L = U^{-1}L^{-1} = A^{-1} \) (the inverse of \(A\)).

Ex.
Still solving \[ \underbrace{ \begin{bmatrix} 1 & 3 & 1\\ 2 & 2 & 0\\ 2 & 2 & -1 \end{bmatrix}}_{A} \underbrace{ \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}}_{x} = \underbrace{ \begin{bmatrix} 1 \\ 4 \\ 2 \end{bmatrix}}_{b} \quad\Longleftrightarrow\quad \begin{bmatrix} 1 & 3 & 1\\ 2 & 2 & 0\\ 2 & 2 & -1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 \\ 4 \\ 2 \end{bmatrix}. \] Gauss: \[ \begin{bmatrix} 1 & 3 & 1 & & 1 & 0 & 0\\ 0 & -4 & -2 & & -2 & 1 & 0\\ 0 & 0 & -1 & & 0 & -1 & 1 \end{bmatrix}. \] Gauss–Jordan⁴⁾: \[ \begin{bmatrix} 1 & 3 & 1 & & 1 & 0 & 0\\ 0 & -4 & -2 & & -2 & 1 & 0\\ 0 & 0 & -1 & & 0 & -1 & 1 \end{bmatrix} \:\Leftrightarrow\: \begin{bmatrix} 1 & 3 & 0 & & 1 & -1 &1\\ 0 & -4 & 0 & & -2 & 3 & -2\\ 0 & 0 & -1 & & 0 & -1 & 1 \end{bmatrix} \:\Leftrightarrow\: \begin{bmatrix} 1 & 0 & 0 & & \textstyle{-\frac{1}{2}} & \textstyle{\frac{5}{4}} & \textstyle{-\frac{1}{2}}\\ 0 & -4 & 0 & & -2 & 3 & -2\\ 0 & 0 & -1 & & 0 & -1 & 1 \end{bmatrix} \] Dividing each row with its pivot gives us the inverse of \(A\): \[ \begin{bmatrix} 1 & 0 & 0 & & \textstyle{-\frac{1}{2}} & \textstyle{\frac{5}{4}} & \textstyle{-\frac{1}{2}}\\ 0 & 1 & 0 & & \textstyle{\frac{1}{2}} & -\textstyle{\frac{3}{4}} & \textstyle{\frac{1}{2}}\\ 0 & 0 & 1 & & 0 & 1 & -1 \end{bmatrix}, \] meaning that \[ A^{-1} = U^{-1} \tilde L = \begin{bmatrix} \textstyle{-\frac{1}{2}} & \textstyle{\frac{5}{4}} & \textstyle{-\frac{1}{2}}\\ \textstyle{\frac{1}{2}} & -\textstyle{\frac{3}{4}} & \textstyle{\frac{1}{2}}\\ 0 & 1 & -1 \end{bmatrix} \] is the linear transformation that brings \(A\) into the unit matrix. The solution to the original problem is \[ x = A^{-1}b \quad\text{ for any } b \in \mathbb R^3. \]

N.b. The product of the final pivots in the Gauss–Jordan elimination (the diagonal elements \( (1)(-4)(-1) = 4\) in the example above) is the determinant of \(A\). It is zero exactly if the Gauss(–Jordan) elimination produces a zero row. In this case the equation \(Ax = b\) either has no, or infinitely many, solutions (depending on \(b\)).

¹⁾

When \(A\) and \(b\) are real, one looks for \(x\) real.

²⁾

When it exists, it is unique if one requires the diagonal elements of \(L\) to be all ones.

³⁾

A permutation matrix has (a permutation of) the standard basis \(\{e_j\}_{j}\) as rows.

⁴⁾

If we start here with \(I\) to the right in the augmented matrix, we obtain \(U^{-1}\) instead of \(A^{-1}\).