Gaussian Elimination

The three legal row operations

Gaussian elimination only ever does three things to the augmented matrix [A | b]:

1. Swap two rows. The system has the same solution either way — equation order is meaningless.
2. Multiply a row by a nonzero scalar. Scaling 2x + 4y = 10 to x + 2y = 5 doesn't change which (x, y) satisfies it.
3. Add a multiple of one row to another. If row i and row j are both true, then row i + c·row j is also true.

Every step of elimination is a sequence of operation 3, occasionally interleaved with operation 1 for stability. The point is to use these moves to drive every entry below the diagonal to zero.

Forward sweep — getting to upper triangular

Process columns left to right. For column k:

1. Choose a pivot — the entry on the diagonal in column k. (With partial pivoting, swap in the row whose entry in column k has the largest absolute value.)
2. For every row i below the pivot, compute the multiplier m_ik = a_ik / a_kk.
3. Subtract m_ik times row k from row i. This zeroes a_ik and updates the rest of row i.
4. Apply the same multiplier to b_i: b_i ← b_i − m_ik·b_k.

After processing column k, every entry below the diagonal in that column is zero. After processing all n columns the augmented matrix has the shape:

[ u₁₁ u₁₂ u₁₃ ... u₁ₙ | c₁ ]
[  0  u₂₂ u₂₃ ... u₂ₙ | c₂ ]
[  0   0  u₃₃ ... u₃ₙ | c₃ ]
[  ...                  ... ]
[  0   0   0  ... uₙₙ | cₙ ]

Back-substitution — reading the answer

The bottom row is now a one-variable equation: u_nn·x_n = c_n, so x_n = c_n/u_nn. Substitute that into row n−1, which now has only one unknown left, and solve for x_n−1. Walk up the rows one at a time. The general formula:

xₖ = ( cₖ − Σⱼ₌ₖ₊₁ⁿ uₖⱼ · xⱼ ) / uₖₖ

Worked example — a 3×3 system

Solve:

2x +  y −  z =  8
−3x −  y + 2z = −11
−2x +  y + 2z =  −3

Augmented matrix:

[  2   1  −1 |  8 ]
[ −3  −1   2 | −11]
[ −2   1   2 | −3 ]

Step 1 — eliminate column 1. Pivot is 2 in row 1. Multipliers:

• m₂₁ = −3/2. Row 2 ← Row 2 − (−3/2)·Row 1 = Row 2 + 1.5·Row 1.
• m₃₁ = −2/2 = −1. Row 3 ← Row 3 + 1·Row 1.

[  2   1   −1   |  8 ]
[  0   1/2  1/2 |  1 ]
[  0   2    1   |  5 ]

Step 2 — eliminate column 2. Pivot is 1/2 in row 2. Multiplier m₃₂ = 2 / (1/2) = 4. Row 3 ← Row 3 − 4·Row 2.

[  2   1   −1  |  8 ]
[  0   1/2  1/2|  1 ]
[  0   0   −1  |  1 ]

Step 3 — back-substitute.

• From row 3: −1·z = 1, so z = −1.
• From row 2: (1/2)y + (1/2)(−1) = 1 → (1/2)y = 3/2 → y = 3.
• From row 1: 2x + 3 − (−1) = 8 → 2x = 4 → x = 2.

Solution: (x, y, z) = (2, 3, −1). Check: 2(2) + 3 − (−1) = 4 + 3 + 1 = 8. ✓

Pivoting — partial, complete, and rook

The single most important practical detail is pivoting. Naive elimination divides by whatever sits on the diagonal, which can be tiny or zero. Three strategies:

Strategy	What you swap	Cost per step	Stability
No pivoting	Nothing	0	Fails outright if a pivot is 0; numerically unstable in general
Partial pivoting	Swap rows so the largest |entry| in column k is on the diagonal	O(n) compares per step	Backward-stable in practice; almost always used
Complete pivoting	Search the whole remaining submatrix and swap both rows and columns	O((n−k)²) compares	Provably stable; rarely worth the constant factor
Rook pivoting	Search a row and column alternately for a local max	O(n) on average	Compromise — better than partial, cheaper than complete

LAPACK's dgesv, NumPy's linalg.solve, and MATLAB's backslash all use partial pivoting. You should too.

Gaussian elimination vs Gauss–Jordan vs LU

	Gaussian + back-sub	Gauss–Jordan	LU factorization
Final shape of A	Upper triangular U	Identity I	L (lower) and U (upper) stored separately
Forward-sweep cost	≈ n³/3	≈ n³/2	≈ n³/3 (same as Gaussian)
Back-substitution cost	≈ n²/2	0 (read off directly)	≈ n² per right-hand side
Best for one Ax = b	Yes — minimum work	No — wasted arithmetic	No — same cost as Gaussian without reuse
Best for many b's, same A	Have to redo each time	Have to redo each time	Yes — factor once, solve cheaply per b
Useful for inversion	No (need extra work)	Yes — apply to [A | I]	Yes — solve LU·X = I column by column
Output reusable	Just the answer x	Just the answer x	L, U, P — reusable for solve, det, refinement

Elimination as a sequence of matrix multiplications

Each elimination step is left-multiplication by an elementary matrix. Subtracting m·row k from row i corresponds to multiplying by I minus a single m·e_ie_k^T. After all n−1 columns, you have

Lₙ₋₁ ⋯ L₂ L₁ · A = U

and inverting the L_k's shows that A = (L₁⁻¹ ⋯ L_n−1⁻¹) · U = L · U where L collects all multipliers below its diagonal. That's the LU decomposition: a record of the elimination steps.

Where Gaussian elimination shows up

• Solving Ax = b — the canonical use, from circuit equations to finite-element stiffness systems.
• Computing the determinant — det(A) = ±(product of diagonal entries of U), where the sign comes from the parity of row swaps.
• Computing a matrix inverse — apply elimination to [A | I]; when the left side becomes I, the right side is A⁻¹.
• Finding rank, null space, column space — the row-echelon form exposes pivot columns (basis of column space) and free columns (null-space directions).
• Polynomial and spline interpolation — Vandermonde and tridiagonal systems are still solved by elimination, often with structured-matrix shortcuts that bring the cost down to O(n²) or O(n).
• Linear programming pivots — the simplex method is essentially Gaussian elimination with a cost-vector test driving which pivot to take next.

Common mistakes

• Skipping pivoting. Hand-rolled Gaussian elimination without row swaps is one of the classic numerical traps. A pivot of 10⁻¹⁵ produces multipliers around 10¹⁵ and destroys 15 digits of precision in a single step.
• Forgetting to update the right-hand side. Every operation on row k of A must also be applied to b_k. Working on the augmented matrix [A | b] enforces this automatically.
• Treating row swaps as an afterthought. The row-swap parity flips the sign of the determinant, and the swap order matters when you want LU = PA later. Track the permutation as you go.
• Mistaking a zero pivot for "the algorithm broke". A zero pivot means the column has a structural problem — either swap to fix it or accept that the matrix is singular. The algorithm is fine; the data is exposing rank deficiency.
• Using Gauss–Jordan when Gaussian + back-sub would do. Gauss–Jordan is taught first because it produces the identity directly, but it does roughly 50% more arithmetic. For Ax = b, stop at upper triangular and back-substitute.
• Inverting A to solve Ax = b. Computing A⁻¹b is slower and less accurate than direct elimination. The numerical-analysis maxim "never invert a matrix unless you really need A⁻¹" applies here.