LU decomposition

The formula

The core identity:

A = L · U

L is lower triangular, U is upper triangular. Reading: everything Gaussian elimination did to reduce A to a triangle is stored in L; the triangle it produced is U.

In practice the rows get reordered for numerical safety, so what a library actually computes is:

P · A = L · U

P is a permutation matrix — a bookkeeping record of the row swaps. Reading: some row order of A factors cleanly; P remembers which one.

Once factored, A·x = b becomes two triangular sweeps:

solve  L·y = P·b   by forward substitution   (top to bottom)
solve  U·x = y     by back substitution      (bottom to top)

Reading: a triangular system needs no elimination at all — each unknown falls out in one line once the ones above (or below) it are known.

Two naming conventions and one refinement you'll see in library docs:

Doolittle:  L has 1's on its diagonal, U carries the pivots
Crout:      U has 1's on its diagonal, L carries the pivots
LDU:        A = L · D · U   — 1's on both diagonals, pivots pulled into diagonal D

Reading: same factorization, different choice of where to park the scaling. The determinant comes along free:

det(A) = ±( u₁₁ · u₂₂ · ⋯ · uₙₙ )    sign flips once per row swap in P

Where you meet it

Every time you type x = A\b or numpy.linalg.solve(A, b). Neither MATLAB nor NumPy "inverts" anything. For a general square system both hand the matrix to LAPACK, which runs LU with partial pivoting, then does the two triangular sweeps. When a solver throws a "matrix is singular to working precision" warning, LU just hit a pivot it couldn't live with. Knowing that turns a cryptic warning into a diagnosis: your model has a redundant or missing equation.
Many load cases on one structure. A stress group runs one airframe model against dozens of load cases — maneuver, gust, hard landing, thermal. The stiffness matrix doesn't change; only the right-hand side does. Factoring costs about 2n³/3 operations; each additional solve costs about 2n². For a 1,000-equation model that's roughly 6.7×10⁸ operations to factor and 2×10⁶ per load case afterward — the second load case costs about a third of a percent of the first. Factor once, sweep many. That asymmetry is the entire economics of production FEA and the reason load-case sweeps are cheap while remeshing is expensive.
SPICE, every timestep. Circuit simulation builds a nodal-analysis matrix that is sparse and unsymmetric — Cholesky doesn't apply. At every Newton iteration of every timestep, SPICE re-solves that system with sparse LU. When a transient run of a power-supply board crawls, the time is going into repeated LU factorizations of the circuit matrix.
At the review board, when someone says "just invert the matrix." Nobody credible inverts. Computing A⁻¹ explicitly costs about three times the factorization, loses accuracy, and destroys sparsity. If a briefing chart shows x = A⁻¹·b as the implementation rather than the notation, that's a finding. The correct answer to "how did you solve it" is "factored once, forward/back substitution per case."

How it works

LU is not a different algorithm from Gaussian elimination — it is Gaussian elimination, with the multipliers you'd normally discard filed into L. Eliminate column by column; each "subtract m times row i from row j" writes m into L, and what's left when you finish is U. Small worked case:

A = |  2   1   1 |        L = |  1   0   0 |     U = | 2   1   1 |
    |  4  -6   0 |    →       |  2   1   0 |         | 0  -8  -2 |
    | -2   7   2 |            | -1  -1   1 |         | 0   0   1 |

Multiply L·U and A comes back exactly. Row 2 of L reads "I subtracted 2× row 1"; that 4 in A became 2·2. And det(A) = 2·(−8)·1 = −16 with no extra work.

Pivoting is not optional. Plain LU fails on perfectly good matrices — A = |0 1; 1 0| is invertible but its first pivot is zero, and the factorization dies on the first division. The formal condition: an invertible matrix has an LU factorization without pivoting only when every leading principal minor is nonzero ^[1]. Worse than outright failure is quiet damage. Factor A = |1e-20 1; 1 1| with right-hand side (1, 2) in double precision without pivoting and you get x = (0, 1); the true answer is essentially (1, 1). The tiny pivot manufactures a multiplier of 10²⁰, and the real data drowns in it. Partial pivoting — swap up the largest available entry in the column before eliminating — fixes it for a cost that's only quadratic. The theoretical worst-case element growth under partial pivoting is 2ⁿ⁻¹, but matrices that achieve it are contrived; decades of production use say pivoted LU is stable in practice ^[1]^[3]. Every serious library pivots by default. If you're hand-rolling LU on an embedded target and skipped pivoting to save code space, you've built the failure above into your flight software.

A clean factorization proves nothing about conditioning. LU can succeed brilliantly on a nearly singular matrix and hand you a solution that's noise. The factorization answers "can I eliminate?" — the condition number answers "should I trust the result?" They are separate checks; do both.

Limits and the standard mistakes:

Symmetric positive definite? Use Cholesky. Half the operations, half the storage, no pivoting. LU on a stiffness matrix works but wastes the structure.
Sparse matrices reorder before they factor. Naive elimination order fills a sparse matrix with new nonzeros ("fill-in"). Sparse LU packages spend real effort choosing the elimination order first; that's why they beat dense solvers by orders of magnitude on FEA-sized problems.
The classic loop bug: calling solve(A, b) inside a loop over right-hand sides, refactoring the same matrix a thousand times. Factor once, keep L, U, P, and substitute per case. Every numerical library exposes this (scipy.linalg.lu_factor / lu_solve, MATLAB decomposition(A)); most of the slow scripts in any test group ignore it.

History

Elimination itself is ancient. Chapter Eight of the Chinese Nine Chapters on the Mathematical Art — a text referenced by 179 AD, with parts dating to roughly 150 BC — solves simultaneous linear systems by exactly the array manipulations we'd now call row reduction ^[2]^[6]. Carl Friedrich Gauss earned the modern name a different way: around 1810 he devised a compact notation for symmetric elimination that nineteenth-century hand computers adopted wholesale for the normal equations of least-squares work in geodesy and astronomy ^[2]^[6].

The hand computers kept tuning it. Myrick Doolittle, a computer at the United States Coast Survey, reorganized Gauss's process into tables built for speed with products and reciprocals — his name is still attached to the unit-lower-triangular convention ^[1]^[2]. In 1941 Prescott Crout of MIT published a rearrangement in the Transactions of the AIEE that accumulated sums of products, tailored to the desk calculating machines of the day; a leading calculator manufacturer publicized the method, and the unit-upper-triangular convention carries his name ^[2]^[5].

The decisive conceptual move — seeing elimination not as a procedure but as a factorization — came from an astronomer. Tadeusz Banachiewicz (1882–1954), director of the Kraków observatory, had invented his own column-oriented matrix algebra, the "Cracovian" calculus, in 1925 for orbit and least-squares computations ^[4]. In 1938 he published the decomposition view of elimination for general matrices: the matrix splits into triangular factors, and solving becomes substitution ^[1]^[2].

The idea then landed at exactly the right moment. In 1947, John von Neumann and Herman Goldstine, analyzing whether the new electronic computers could be trusted with elimination, opened their rounding-error study by describing the method as the combination of two tricks — first decompose A into a product of two triangular matrices ^[2]^[3]. A year later Alan Turing, working at the National Physical Laboratory, published "Rounding-off errors in matrix processes" in the Quarterly Journal of Mechanics and Applied Mathematics, writing Gaussian elimination explicitly as the LU factorization and introducing, in the same paper, the condition number of a matrix ^[2]^[3]. The influential numerical analysis textbooks of the 1960s reprised Turing's presentation ^[2], which is why elimination is taught today as A = L·U — and why what runs when your solver "solves" is a method with two thousand years of field history and a seventy-five-year-old proof of trust.

Related tools

/tools/beam-deflection
/tools/vibration-natural-freq
/tools/series-rlc-impedance

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →