Cholesky decomposition

The formula

The whole thing is one identity. If A is symmetric and positive definite:

A = L · Lᵀ

L is lower triangular with positive numbers on its diagonal. Reading: A is the "square" of a triangular matrix, the matrix version of writing a positive number as x².

The entries fall out one column at a time:

L_jj = √( A_jj − Σ(k=1..j-1) L_jk² )              (diagonal)
L_ij = ( A_ij − Σ(k=1..j-1) L_ik · L_jk ) / L_jj   (below diagonal, i > j)

Reading: each new entry is the original matrix entry minus the dot product of the parts of L you already computed, then scaled. Nothing to the right or below is touched — you march top-left to bottom-right and never look back.

Once you have L, solving A·x = b is two triangular sweeps:

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

There is a square-root-free cousin, the LDLᵀ form:

A = L · D · Lᵀ

Here L has 1's on the diagonal and D is a diagonal matrix. Reading: same idea, but the scaling is pulled out into D so you never call sqrt. That matters on hardware where a square root is expensive or where you want to handle a matrix that is only semi-definite.

Where you meet it

FEA stiffness solves. Assemble a finite-element model of a bracket or an airframe rib and you get K·u = F, where K is the global stiffness matrix. K is symmetric and — for a properly constrained structure — positive definite. Every mainstream solver (NASTRAN, ANSYS, Abaqus) factors K with Cholesky or its sparse variant before it ever computes a displacement. When the run dies with "negative pivot" or "matrix not positive definite," Cholesky is the thing that just failed, and it's telling you the structure has a rigid-body mode or a mechanism — a missing constraint, not a bad material card.
Kalman filter square roots. A Kalman filter carries a covariance matrix P that must stay symmetric positive definite. Roundoff slowly erodes that, and a filter with an indefinite P diverges silently. Square-root filters store the Cholesky factor of P instead of P itself and propagate the factor. You literally cannot represent a non-PD covariance that way, so the filter stays honest. This shows up in every INS/GPS box and missile seeker tracker in the valley.
Monte Carlo with correlated inputs. You have a covariance matrix Σ for a set of correlated variables — say the joint scatter on a strapdown IMU's three gyro biases — and you need to draw correlated random samples for a dispersion run. Factor Σ = L·Lᵀ, draw a vector z of independent unit-normal numbers, and L·z has exactly the covariance you wanted. One factorization, thousands of correlated draws.
Least-squares normal equations, on the bench. Fit a calibration curve or level a survey net and you form the normal equations AᵀA·x = Aᵀb. The matrix AᵀA is symmetric and, provided the columns of A are independent (a full-rank fit), positive definite — so Cholesky is the textbook solver ^[3]. If the fit is rank-deficient (a survey net with a datum defect, a calibration model with collinear terms), AᵀA is only semi-definite and the factorization breaks down. One more caution: forming AᵀA squares the condition number of the problem, so for ill-conditioned fits — high-order polynomials, near-collinear parameters — factor A directly with QR instead ^[3]. Still, this is exactly the problem Cholesky himself was solving with survey data.

How it works

The engine is the same as Gaussian elimination, but symmetry lets you throw away half the bookkeeping. Ordinary LU factorization spends about 2n³/3 floating-point operations. Cholesky spends about n³/3 — roughly twice as fast, and it stores only one triangle instead of two. On a million-DOF FEA model that factor of two is real wall-clock time.

The bigger win is that Cholesky needs no pivoting. General Gaussian elimination has to hunt for a large pivot each step to stay numerically stable, which wrecks the memory-access pattern and, on sparse matrices, scrambles the fill-in. Positive definiteness guarantees the diagonal you divide by is always positive and never dangerously small, so you march straight down the diagonal. That predictability is why sparse Cholesky is the backbone of structural solvers — you can plan the elimination order (the "symbolic factorization") before touching a single number.

Cholesky doubles as a positive-definiteness test, and this is the gotcha that bites people. The algorithm only works if every A_jj − Σ L_jk² under the square root stays positive. The moment it goes zero or negative, you're taking √(negative) and the factorization halts. That is not a bug — it's the matrix telling you it is not positive definite. In FEA that means an unconstrained model; in a Kalman filter it means your covariance has gone bad; in optimization it means your Hessian isn't at a minimum. Trap the failure and read it as a diagnostic, not a crash.

Limits of validity, stated plainly:

Symmetric only. If A isn't symmetric, Cholesky doesn't apply — reach for LU.
Positive definite only. Merely symmetric isn't enough. If your matrix is indefinite (some positive, some negative eigenvalues), use LDLᵀ with symmetric pivoting (Bunch–Kaufman), not plain Cholesky.
Watch the square roots. The classic LLᵀ form calls sqrt on the diagonal. On tight embedded targets, or when the matrix is only positive semi-definite, use the LDLᵀ variant to sidestep the roots and the near-zero divisions.

The mistake people make is feeding it a matrix that is almost positive definite — a covariance that's mathematically PD but has drifted slightly indefinite from roundoff — and then treating the failure as a numerical fluke to be patched by "adding a small number to the diagonal." Sometimes that jitter is legitimate regularization. Often it's papering over a real modeling error you should have caught.

Worked check you can run in your head. Take:

A = | 4   12  -16 |
    | 12  37  -43 |
    | -16 -43  98 |

The factor is:

L = |  2   0   0 |
    |  6   1   0 |
    | -8   5   3 |

Multiply L·Lᵀ and you get A back exactly — top-left 2·2 = 4, and the bottom-right (-8)² + 5² + 3² = 64 + 25 + 9 = 98. That last sum being positive is precisely why the factorization exists.

History

André-Louis Cholesky was not a mathematician by trade — he was a French army officer and a surveyor, and the method fell out of the day job. Born 15 October 1875 in Montguyon, in southwestern France, he came up through the École Polytechnique and the artillery school, then landed in the Geodesic Section of the Army's Geographic Service around 1905 ^[1]^[2]. His work was triangulating the ground — running survey nets across Tunisia, Algeria, Morocco, and Crete — and every large survey ends in the same numerical headache: a big pile of overdetermined measurements you have to reconcile by least squares.

Between roughly 1905 and 1910 Cholesky worked out a compact hand-computation scheme for solving the symmetric normal equations those least-squares problems produce ^[1]^[2]. It was a computer's method in the old sense of the word — a human-computer's method — laid out so a surveyor with a table of logarithms could grind through it without mistakes. He wrote it up in a manuscript but never published it in a journal ^[1]^[2].

He never got the chance. Cholesky was killed on 31 August 1918, from wounds taken in the field in the last months of the First World War ^[1]^[2]. His method surfaced afterward: Commandant Ernest Benoit, a fellow officer, published an account of it in 1924 in the Bulletin géodésique, which is how the wider world first read it ^[1]^[2].

For decades it stayed a geodesist's tool. It became a household name in numerical linear algebra only after the digital-computer era, when John Todd, Alan Turing, and others in the 1940s and 1950s catalogued the direct methods for solving linear systems and Cholesky's symmetric, pivot-free scheme turned out to be exactly the right fit for symmetric positive-definite matrices ^[1]. A survey trick from the age of theodolites is now running inside every finite-element and Kalman-filter solver on Redstone.

Related tools

/tools/vibration-natural-freq
/tools/beam-deflection
/tools/stress-strain

Sources

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