The formula
Gaussian elimination doesn't have a single formula so much as a recipe that turns the system A·x = b into something you can read off by inspection. Three moves are legal, and they never change the solution:
Swap two rows. (reorder your equations)
Scale a row by a nonzero constant. (multiply an equation through)
Add a multiple of one row to another. (the actual elimination step)
The one that does the work is the elimination step. For each row i below the pivot row k:
multiplier m = A[i,k] / A[k,k]
row_i = row_i − m · row_k
Reading: pick the number sitting under the pivot, divide it by the pivot to get m, then subtract m times the pivot row so that the entry under the pivot becomes zero. Repeat down the column, then move to the next pivot.
When every entry below the diagonal is zero, you have an upper-triangular system:
[ u11 u12 u13 ] [x1] [c1]
[ 0 u22 u23 ] [x2] = [c2]
[ 0 0 u33 ] [x3] [c3]
Now back-substitute from the bottom: x3 = c3/u33, then x2 = (c2 − u23·x3)/u22, and so on up. Gauss-Jordan keeps going and clears the entries above the diagonal too, landing on the identity so the answer sits in the right-hand column with no back-substitution. It costs more arithmetic and is rarely worth it for solving one system — its real use is inverting a matrix or teaching.
Where you meet it
- Any resistive network you node-analyze by hand. Write Kirchhoff's current law at each node and you get one linear equation per node. A five-node circuit is a 5×5 system, and the way you solve it on paper — add this equation to that one to kill a term — is Gaussian elimination whether you call it that or not. SPICE does the same thing in its inner loop, just with pivoting and sparse bookkeeping.
- Calibration fits on the test bench. Fitting a load cell or thermocouple to a polynomial, or solving the normal equations of a least-squares fit, ends in a small dense linear system. The
A·x = bthat pops out of a curve fit is solved by elimination or its cousin, the Cholesky factorization. - Finite-element and structural solvers at a design review. The stiffness equation
K·u = Fis Gaussian elimination on a very large, very sparseK. When a colleague's homegrown solver returns garbage on one load case and fine on the next, the fault is almost always pivoting — see below. - Sensor fusion and state estimation. A Kalman update inverts a small covariance matrix every step. That inversion is elimination underneath, and it runs in a tight loop on flight hardware.
How it works
The idea is embarrassingly simple: subtract equations to make variables disappear. The trouble is entirely in the arithmetic, and it shows up the moment you use floating-point instead of exact fractions. Note also that the recipe as written assumes a square system with exactly one solution — as many independent equations as unknowns.
The pivot can't be zero — and shouldn't be small. The multiplier m = A[i,k]/A[k,k] divides by the pivot. If the pivot is zero the step is undefined. If it's merely tiny, you divide by a small number, get a huge multiplier, and subtract a huge row from a normal one. The normal row's real information gets swamped by rounding error in the huge one. Watch it happen:
[ 1e-20 1 ] [x1] [1] true answer: x1 ≈ 1, x2 ≈ 1
[ 1 1 ] [x2] = [2]
Eliminate straight down without thinking, using 1e-20 as the pivot, and in double precision you get x1 = 0 — dead wrong. The multiplier is 1e20, and 1 − 1e20 rounds to exactly −1e20; the "1" that carried x1's information is gone. Swap the two rows first so the pivot is 1, run the identical algorithm, and you get x1 = 1. Same math, same equations, right answer. That row swap is partial pivoting: before each elimination step, scan the column at and below the pivot, and swap up whichever row has the largest absolute value. It costs almost nothing, and it is why production solvers don't blow up on problems like the one above. One case a swap can't rescue: if every entry in the pivot column, at the pivot and below, is zero, there is nothing to swap up — the matrix is singular and the system has no unique solution. The algorithm hasn't failed; the problem has. This is the mistake people make — a hand-built or textbook-copied solver that eliminates in fixed order with no pivot search will run clean for months and then detonate the first time a small diagonal entry shows up.
Pivoting makes the algorithm stable — it can't make the problem good. Partial pivoting buys you backward stability: the answer you compute is the exact solution of a nearby problem. Whether that nearby answer is anywhere close to your answer depends on the matrix itself, measured by its condition number, cond(A). The rule of thumb is brutal and simple: expect to lose about log10(cond(A)) of your 16 double-precision digits, and no amount of pivoting gets them back. A matrix with cond(A) = 1e12 leaves you roughly four trustworthy digits even from a flawless solve. The warning sign inside the elimination is a pivot that stays tiny even after the swap — the biggest entry left in the column is still small, which means you're brushing up against singularity. This bites hardest on the calibration-fit use case above: forming the normal equations squares the condition number of the design matrix, and a polynomial fit on raw units (say, temperatures from 273 to 373 raised to the fourth and fifth powers) produces a Vandermonde-flavored matrix that is spectacularly ill-conditioned before you square anything. The fix is upstream of the solver — center and scale the data before forming the fit, or skip the normal equations entirely and use a QR-based least-squares solve. Check the condition number when the answer matters; a perfect elimination on a rotten matrix still returns rot.
Pivoting bounds the damage but doesn't erase it. The quantity that governs stability is the growth factor — how large the intermediate numbers get relative to the original entries. With partial pivoting the growth factor can, in the worst case, reach 2^(n−1) for an n×n matrix. Wilkinson's classic example makes it happen: ones on the diagonal, −1 below it, and a column of ones on the right. For a 6×6 that matrix pushes the largest value in the triangular factor to exactly 2^5 = 32, and partial pivoting can't stop it because the diagonal is already the biggest entry in every column, so no swap ever fires. In practice you almost never see this. On real matrices the growth factor typically creeps up like n^(2/3), not 2^n, which is why partial pivoting has been the default for seventy years despite the ugly worst case. If you genuinely can't tolerate the risk, complete pivoting searches the whole remaining submatrix for the largest entry and swaps both rows and columns; it's provably tamer but costs enough extra searching that it's reserved for the paranoid.
The cost is cubic. Elimination on a dense n×n system takes about (2n³ + 3n² − 5n)/6 multiplications and the same count of subtractions — roughly 2/3 · n³ operations for the forward sweep, plus an n² back-substitution. Double the size of the matrix and the work goes up eightfold. That n³ is why doubling mesh density in an FE model quietly octuples solve time, and why anyone with a large sparse K never runs plain dense elimination — they exploit the zeros.
Factor once, solve many. The elimination is really computing an L·U factorization of A (the multipliers you used are L; the triangle you're left with is U). If you need to solve A·x = b for several different right-hand-side vectors b — many load cases, many time steps — factor A once for 2/3·n³, then each new b is a cheap n² pair of triangular solves. Recomputing the full elimination per right-hand side is the second-most-common rookie mistake after skipping the pivot.
History
The method is old enough that "Gaussian" is a bit of a swindle. The earliest complete write-up is Chinese: Chapter Eight of the Nine Chapters on the Mathematical Art — Jiuzhang Suanshu — titled "Rectangular Arrays," lays out eighteen problems solved by eliminating unknowns between equations, using counting rods arranged in columns and even handling negative numbers along the way [1][2]. The text was compiled across the Han dynasty; parts trace to roughly the second century BCE, with the received version settling by around 179 CE [1][2]. That is essentially the algorithm on this page, columns instead of rows, two thousand years early.
Europe caught up slowly. Isaac Newton, around 1669–1670, complained that none of the algebra books he owned taught how to solve simultaneous equations, and wrote out a systematic elimination to fill the gap; Cambridge eventually printed those notes as the Arithmetica Universalis in 1707, well after he'd left the subject behind [1][3]. Carl Friedrich Gauss came to it in 1810 as a computational tool, working out a notation for symmetric elimination while grinding through the least-squares arithmetic of tracking asteroids — the kind of overdetermined fitting problem that gave the modern method its polish [1][3]. Wilhelm Jordan described the reduce-all-the-way variant in 1888, in a handbook on surveying, and B.-I. Clasen published something equivalent the same year, apparently independently — which is why the "reduce to the identity" form carries the compound name Gauss-Jordan [1][3]. The label "Gaussian elimination" for the plain version is younger than the algorithm by two millennia: it stuck only in the 1950s, out of some confusion over who had actually done what [1].
Related tools
- /tools/wheatstone-bridge — a bridge balance is a small linear system in the arm resistances; solving it by hand is elimination in miniature
- /tools/series-rlc-impedance — network analysis that scales up into the node-equation systems Gaussian elimination solves
- /tools/pi-attenuator — resistor-network design that reduces to solving simultaneous equations for the branch values
- /tools/t-attenuator — the T-network counterpart, same simultaneous-equation core
- /tools/strain-gauge-bridge — multi-gauge rosette reduction ends in a small linear system for the strain components