In general, computing the inverse of a matrix or solving a system of equations straight from the textbook formula costs on the order of \(n!\) operations, reaching millions of operations for a \(10\times 10\) matrix. This is not just a computational inconvenience; it also causes significant accumulation of floating-point errors in computers, leading to inaccuracies that can't be overlooked. Fortunately a lot of work can be reduced by a single \(O(n^3)\) procedure: Gaussian elimination, a.k.a. row reduction. As with all things in math, we're not satisfied with just knowing the how, but also the why. This article justifies Gaussian elimination from the ground up.

Throughout, \(\mathbb{K}\) is a field (think \(\mathbb{R}\) or \(\mathbb{C}\)), \(M_n(\mathbb{K})\) the \(n\times n\) matrices over it, and \(I\) the identity.

The HOW

Is this matrix invertible?

Take \(A=\begin{bmatrix} 2 & 1 & 1 \\ 4 & 3 & 3 \\ 8 & 7 & 9 \end{bmatrix}\). Subtract \(2\times\)row 1 from row 2 and \(4\times\)row 1 from row 3, then \(3\times\)row 2 from row 3:

\[ \begin{bmatrix} 2 & 1 & 1 \\ 4 & 3 & 3 \\ 8 & 7 & 9 \end{bmatrix} \longrightarrow \begin{bmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 3 & 5 \end{bmatrix} \longrightarrow \begin{bmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 2 \end{bmatrix}. \]

The result is upper triangular with non-zero diagonal entries. A square matrix is invertible if and only if row reduction leaves no zero on the diagonal.

Finding \(A^{-1}\)

To find the inverse matrix \(A^{-1}\), simply glue the identity onto the right of \(A\) to form the wide augmented matrix \([\,A\mid I\,]\), and apply row operations to the whole thing:

\[ \left[\begin{array}{ccc|ccc} 2 & 1 & 1 & 1 & 0 & 0 \\ 4 & 3 & 3 & 0 & 1 & 0 \\ 8 & 7 & 9 & 0 & 0 & 1 \end{array}\right] \longrightarrow \left[\begin{array}{ccc|ccc} 2 & 1 & 1 & 1 & 0 & 0 \\ 0 & 1 & 1 & -2 & 1 & 0 \\ 0 & 0 & 2 & 2 & -3 & 1 \end{array}\right]. \]

The left block is now upper triangular. We keep going until the left block is the identity. This extended sweep is called Gauss–Jordan elimination. Scale row 1 by \(\tfrac12\) and row 3 by \(\tfrac12\), then clear column 3 (via \(R_2\leftarrow R_2-R_3\), \(R_1\leftarrow R_1-\tfrac12 R_3\)) and finally column 2 (via \(R_1\leftarrow R_1-\tfrac12 R_2\)):

\[ \left[\begin{array}{ccc|ccc} 2 & 1 & 1 & 1 & 0 & 0 \\ 0 & 1 & 1 & -2 & 1 & 0 \\ 0 & 0 & 2 & 2 & -3 & 1 \end{array}\right] \longrightarrow \left[\begin{array}{ccc|ccc} 1 & 0 & 0 & \tfrac32 & -\tfrac12 & 0 \\ 0 & 1 & 0 & -3 & \tfrac52 & -\tfrac12 \\ 0 & 0 & 1 & 1 & -\tfrac32 & \tfrac12 \end{array}\right]. \]

The right block is \(A^{-1}\):

\[ A^{-1} = \begin{bmatrix} \tfrac32 & -\tfrac12 & 0 \\ -3 & \tfrac52 & -\tfrac12 \\ 1 & -\tfrac32 & \tfrac12 \end{bmatrix} \]

Here is a handwavy version of "the why": Every row operation is left-multiplication by some matrix \(E_i\), so reducing the left block to \(I\) means the accumulated product \(E_k\cdots E_1\) satisfies \(E_k\cdots E_1 A = I\), hence by the uniqueness of inverses, \(E_k\cdots E_1 = A^{-1}\). The right block started as \(I\) and underwent the very same operations, so it became \(E_k\cdots E_1 I = A^{-1}\).

Solving \(Ax=b\)

To solve \(x_1 + 2x_2 - x_3 = 2,\ \ 2x_1 + 3x_2 + x_3 = 11,\ \ x_1 - x_2 + 2x_3 = 5\), stack the coefficients and right-hand side into one augmented matrix \([\,A\mid b\,]\) and reduce (subtract \(2\times\)row 1 from row 2, row 1 from row 3, then \(3\times\)row 2 from row 3):

\[ \left[\begin{array}{ccc|c} 1 & 2 & -1 & 2 \\ 2 & 3 & 1 & 11 \\ 1 & -1 & 2 & 5 \end{array}\right] \longrightarrow \left[\begin{array}{ccc|c} 1 & 2 & -1 & 2 \\ 0 & -1 & 3 & 7 \\ 0 & 0 & -6 & -18 \end{array}\right]. \]

The system is now triangular, so read it off bottom-up — back-substitution:

\[ -6x_3 = -18 \implies x_3 = 3, \qquad -x_2 + 9 = 7 \implies x_2 = 2, \qquad x_1 + 4 - 3 = 2 \implies x_1 = 1. \]

So \(x=(1,2,3)\), and indeed \(Ax=(2,11,5)\). No inverse was computed: a few row subtractions plus three substitutions, total cost \(O(n^3)\).

What about the determinant?

Most courses introduce determinants as a way to test invertibility, compute inverses of matrices, and calculate eigenvalues, and use row reduction technique to reduce the complexity of calculating determinants. We're not going to talk about determinants yet because determinants deserve their own proper introduction: they are so so so much more than a complicated formula (the Leibniz formula) which generalizes the \(2\times 2\) and \(3\times 3\) formulae we all know by heart. Determinants emerge as a natural mathematical object and generalize concepts from linear algebra to commutative algebra with deep connections to algebraic geometry.

The WHY: every row operation is an elementary matrix

Everything below rests on one identity. Write the rows of \(A\) as \(r_1,\dots,r_n\). For any \(n\times n\) matrix \(E\), the \(i\)-th row of the product \(EA\) is

\[ (EA)_i = \Big[\, \textstyle\sum_{k=1}^{n} E_{ik}A_{k1}, \ \dots, \ \sum_{k=1}^{n} E_{ik}A_{kn} \,\Big] = \sum_{k=1}^{n} E_{ik}\,r_k, \]

Read the right-hand side carefully. Each \(E_{ik}\) is a scalar (the \((i,k)\) entry of \(E\)) and each \(r_k\) is a row vector (the \(k\)-th row of \(A\)), so every term \(E_{ik}\,r_k\) is a scalar–vector product, and the sum is a linear combination of \(A\)'s rows with the weights taken from row \(i\) of \(E\).

Each of row reduction's three moves is therefore realized by taking \(E\) to be the identity with a single tweak, and we define the corresponding elementary matrix to be exactly that tweak of \(I\). In each case the effect on \(A\) is read straight off the identity above:

  1. Swap \(S_{ij}\) (\(i\neq j\)): \(I\) with rows \(i\) and \(j\) exchanged, so its row \(i\) is \(e_j^{\mathsf T}\) and its row \(j\) is \(e_i^{\mathsf T}\). Then \((S_{ij}A)_i = r_j\), \((S_{ij}A)_j = r_i\), leaving all other rows fixed.
  2. Scale \(D_i(c)\) (\(c\neq 0\)): \(I\) with its \(i\)-th diagonal entry replaced by \(c\), so its row \(i\) is \(c\,e_i^{\mathsf T}\). Then \((D_i(c)A)_i = c\,r_i\), leaving all other rows fixed.
  3. Transvection \(T_{ij}(c)\) (\(i\neq j\)): \(I + c\,e_i e_j^{\mathsf T}\), i.e. \(I\) with an extra \(c\) in entry \((i,j)\). Then \((T_{ij}(c)A)_i = r_i + c\,r_j\).

Every elementary matrix is invertible

Lemma. Each elementary matrix is invertible, and its inverse is elementary of the same type:

\[ S_{ij}^{-1}=S_{ij}, \qquad D_i(c)^{-1}=D_i\!\left(\tfrac1c\right), \qquad T_{ij}(c)^{-1}=T_{ij}(-c). \]

Proof. Each identity exhibits a same-type operation that undoes the first — equivalently, by the correspondence above, a same-type matrix whose product with the original is \(I\). A swap repeated restores the rows: \(S_{ij}^2=I\). Scaling row \(i\) by \(c\) and then by \(1/c\) restores it: \(D_i(c)\,D_i(1/c)=I\). For the transvection, \(i\neq j\) gives \(e_j^{\mathsf T}e_i=0\), so

\[ \big(I+c\,e_ie_j^{\mathsf T}\big)\big(I-c\,e_ie_j^{\mathsf T}\big) = I - c^2\,e_i\underbrace{(e_j^{\mathsf T}e_i)}_{=\,0}e_j^{\mathsf T} = I. \qquad \blacksquare \]

Since the composition of invertible matrices is also invertible, \(E=E_k\cdots E_1\) is invertible.

Reduction preserves invertibility and yields \(A^{-1}\)

Each elementary matrix is invertible, and a composition of invertible linear maps is again invertible, so \(A\) is invertible \(\iff EA=E_k\cdots E_1 A\) is. Notice that \(EA=I \iff E=A^{-1}\) by the uniqueness of the inverse map, so the product of elementary matrices gives us the literal inverse of \(A\).

This is why the augmented sweep \([\,A\mid I\,] \to [\,I\mid A^{-1}\,]\) works: carrying \(I\) through the same operations records \(E_k\cdots E_1 I = A^{-1}\) in the right block.

Solving \(Ax=b\)

Since \(E\) is invertible, \(Ax=b\iff EAx=Eb\), hence

\[ \{x : Ax=b\} = \{x : EAx = Eb\}. \]

So each elementary row operation on the augmented matrix \([\,A\mid b\,]\) yields an equivalent system. Reduce to \([\,R\mid d\,]\) where \(R=EA\) is the reduced row echelon form of \(A\). This holds whether or not \(A\) is invertible! When \(A\) is not invertible, the solutions to \(Ax=b\) are either not unique or non-existent. To make sense of these facts, we need something called the null space.

Definition. The null space (or kernel) of \(A\) is \(N(A) = \{x\in\mathbb{K}^n : Ax=0\}\).

This is a subspace: if \(Ax=Ay=0\) and \(\lambda\in\mathbb{K}\) then \(A(x+\lambda y)=Ax+\lambda Ay = 0\). Furthermore, row reduction does not alter the null space of \(A\), since \(Ax=0 \iff (E_k\cdots E_1A)x = 0\); hence \(N(A)=N(R)\).

Suppose one particular solution \(x_p\) exists, \(Ax_p=b\). Then for any \(x\),

\[ Ax=b \iff A(x-x_p)=0 \iff x-x_p \in N(A), \]

so the full solution set is the coset \(x_p + N(A)\). This is the structure theorem for \(Ax=b\): whenever a solution exists, the solution set is one particular solution plus the entire null space. Since \(N(A)=\{0\}\) if and only if \(A\) is invertible, \(Ax=b\) has a unique solution if and only if \(A\) is invertible. That a singular \(A\) has \(N(A)\neq\{0\}\) follows from the rank–nullity theorem, fully explained in the next section. Intuitively, if the null space were trivial the map would be injective and hence, as an endomorphism on a finite-dimensional space, invertible.

Now let's consider the case where \(A\) is singular, i.e. \(N(A)\) is non-zero.

By the structure theorem the solution set is either empty or the coset \(x_p + N(A)\); now that \(N(A)\neq\{0\}\), the second case is no longer a single point. Which case occurs is decided purely by consistency: whether a particular solution exists at all. The reduced augmented matrix \([\,R\mid d\,]\) shows it directly.

Consistency criterion. The system \(Rx=d\) is equivalent to \(Ax=b\), so inspect the zero rows of \(R\). If some zero row of \(R\) is paired with \(d_i\neq 0\), that row reads \(0=d_i\), which is impossible: the system is inconsistent and \(Ax=b\) has no solution (equivalently \(b\notin\operatorname{im}A\)). Otherwise every zero row reads \(d_i=0\), a particular solution \(x_p\) exists, and the full solution set is the coset \(x_p+N(A)\), which is infinite over an infinite field.

An optional but important digression: the rank–nullity theorem

Regard \(A\) as the linear map \(T:\mathbb{K}^n\to\mathbb{K}^n\), \(T(x)=Ax\). Its null space \(N(A)=\ker T\) measures how far \(T\) is from being injective: \(T(x)=T(y)\iff T(x-y)=0\iff x-y\in N(A)\), so \(T\) is injective exactly when \(N(A)=\{0\}\), and a larger null space means more inputs are crushed onto the same output. An invertible map is in particular injective, so its null space is necessarily trivial. The rank–nullity theorem precisely describes the relationship between the size of the map's image \(\operatorname{im}T\) (the column space, whose dimension is the rank of \(A\)) and the size of its null space.

Matrices vs linear transformations made precise!

It is worth being careful about what we just did. A matrix \(A\) and a linear transformation \(T\) are not quite the same thing: the matrix is what \(T\) looks like once a basis is fixed. Choose a basis \(\{v_1,\dots,v_n\}\) of \(V\); then the \(i\)-th column of \(A\) holds the coordinates of \(T(v_i)\) in that basis — i.e. the unique scalars \(a_{1i},\dots,a_{ni}\) with \(T(v_i)=\sum_k a_{ki}v_k\). With the standard basis \(\{e_1,\dots,e_n\}\) of \(\mathbb{K}^n\) this is the familiar statement "the \(i\)-th column of \(A\) is \(T(e_i)\)." Multiplying \(A\) by a coordinate vector then just computes \(T\) on the corresponding vector, which is why \(T(x)=Ax\) and why the column space is \(\operatorname{im}T\).

The point of the abstraction is that vectors need not "look like" columns of numbers. A finite-dimensional vector space over a field \(\mathbb{K}\) is defined purely by axioms — an abelian group under addition equipped with a compatible scalar multiplication — and its elements can be far more abstract: polynomials of degree \(<n\), \(m\times n\) matrices, solutions of a linear differential equation, and so on. Fixing a basis is what turns any such space into \(\mathbb{K}^n\) and any linear map into a matrix. In general, column (or row) vectors filled with numbers are a representation of vectors with respect to a chosen basis, not the vectors themselves.

Theorem (rank–nullity). For a linear map \(T:V\to W\) with \(V\) finite-dimensional,

\[ \dim V = \dim\ker T + \dim\operatorname{im}T. \]

Proof. Choose a basis \(u_1,\dots,u_k\) of \(\ker T\) and extend it to a basis \(u_1,\dots,u_k,\,w_1,\dots,w_m\) of \(V\), so \(\dim\ker T=k\) and \(\dim V=k+m\). It suffices to show \(Tw_1,\dots,Tw_m\) is a basis of \(\operatorname{im}T\), for then \(\dim\operatorname{im}T=m\) and \(\dim V=k+m=\dim\ker T+\dim\operatorname{im}T\). They span: any element of \(\operatorname{im}T\) is \(T(v)\) for some \(v=\sum_i a_iu_i+\sum_j b_jw_j\), and since \(Tu_i=0\), \(T(v)=\sum_j b_j\,Tw_j\). They are independent: if \(\sum_j c_j\,Tw_j=0\) then \(T\!\big(\sum_j c_jw_j\big)=0\), so \(\sum_j c_jw_j\in\ker T=\operatorname{span}(u_1,\dots,u_k)\); writing it as \(\sum_i d_iu_i\) gives \(\sum_j c_jw_j-\sum_i d_iu_i=0\), and since the \(u\)'s and \(w\)'s together form a basis, every \(c_j=0\). \(\blacksquare\)

For \(A\in M_n(\mathbb{K})\) we have \(\dim V=n\), \(\ker T=N(A)\), and \(\dim\operatorname{im}T=\operatorname{rank}A\), so

\[ \operatorname{rank}A + \dim N(A) = n. \]

This settles the claim left open earlier. \(A\) is invertible iff \(T\) is bijective iff \(\operatorname{rank}A=n\) iff \(\dim N(A)=0\). A singular \(A\) has \(\operatorname{rank}A< n\), so \(\dim N(A)=n-\operatorname{rank}A>0\) and \(N(A)\neq\{0\}\), and each lost unit of rank is one extra parameter in the solution family whenever the system is consistent.

None of this is special to matrices. Rank–nullity is the first isomorphism theorem of abstract algebra specialized to linear algebra: the isomorphism \(V/\ker T\cong\operatorname{im}T\) prescisely describes the relationship between the size of the image of a group homomoprhism and the size of its kernel, where the group in question is a vector space, and the group homomorphism is a linear transformation in the realm of linear algebra.