Cholesky Factorization as a Positive-Definiteness Test

A symmetric matrix that is also positive definite has a factorization that is unique, cheap, and requires no decisions. The matrix A is written as L L^T with L lower triangular and positive diagonal, and the computation that produces L either succeeds or proves that no such L exists. That last property is unusual: most matrix algorithms produce an answer and leave the question of whether the answer means anything to a separate diagnostic.

The factorization is the standard way to solve the normal equations of a least-squares problem, the stiffness systems of implicit finite element codes, and the covariance systems of Gaussian models. Each of those applications has its own reason for caring about the definiteness verdict, and each pays for the factorization in a different currency: flops, storage, or the number of triangular solves.

The recurrence, read column by column

Equating entries of L L^T with those of A gives, for the diagonal and the entries below it,

a_jj = sum_{k<j} l_jk^2 + l_jj^2
a_ij = sum_{k<j} l_ik l_jk + l_ij l_jj     (i > j)

so each column of L is determined by the entries of A in that column and the columns already computed. The diagonal entry requires a square root; the off-diagonal entries require a division by that root. Nothing else is needed, and no comparison between candidate pivots appears anywhere in the recurrence.

#include <math.h>

/* In-place Cholesky of a symmetric matrix stored row-major, lower triangle
   only. On success the lower triangle of A holds L. Returns 0 on success,
   or the 0-based index of the column where the pivot was not positive. */
int cholesky(int n, double *A)
{
    for (int j = 0; j < n; j++) {
        double s = A[j * n + j];
        for (int k = 0; k < j; k++)
            s -= A[j * n + k] * A[j * n + k];
        if (!(s > 0.0))                 /* also catches NaN from earlier steps */
            return j;
        double ljj = sqrt(s);
        A[j * n + j] = ljj;
        for (int i = j + 1; i < n; i++) {
            double t = A[i * n + j];
            for (int k = 0; k < j; k++)
                t -= A[i * n + k] * A[j * n + k];
            A[i * n + j] = t / ljj;
        }
    }
    return 0;
}

The inner loops are dot products of length j, which is the same shape as the unblocked LU update but with half the terms. The total is n^3/3 multiplications and additions, one third of the 2n^3/3 that the general factorization costs, and the storage is n(n+1)/2 entries rather than n^2.

Why the pivot search disappears

For an LU factorization, the pivot at step k is whatever the elimination has left on the diagonal, and it can be zero or tiny even when the matrix is invertible. Cholesky is protected by a structural fact: if A is positive definite, every leading principal submatrix is positive definite as well, so the Schur complement after any number of steps is positive definite, and its diagonal entries are all strictly positive.

There is a stronger statement, and it is the one that removes the need for growth analysis. The diagonal entries of the Schur complement at step k are bounded above by the corresponding diagonal entries of A. Elimination cannot make a diagonal entry larger; it can only shrink it. The growth factor, which measures how much larger the intermediate quantities become than the input, is therefore exactly one, for every positive definite matrix, with no search and no permutation.

That gives the backward error directly. The computed factor satisfies L L^T = A + E with ||E||_2 of order n * u * ||A||_2, where u is the unit roundoff. The constant does not depend on any property of A beyond symmetry and definiteness, and in particular it does not depend on the condition number. The forward error on a solve still does, through the condition number, exactly as described in Condition Number, Residual, and Backward Error, but the factorization itself is as stable as a factorization can be.

Where the factor of two comes from

The saving is not a trick of implementation. Half the matrix is redundant because A is symmetric, and the recurrence exploits that by never touching the upper triangle. The flop count drops from 2n^3/3 to n^3/3 and the storage from n^2 to n(n+1)/2, which for large n is the difference between a matrix that fits in memory and one that does not.

The same structure helps at the level of the kernel. The blocked version of the algorithm computes a panel of nb columns, then updates the trailing submatrix with a symmetric rank-nb update, which is a matrix-matrix product against its own transpose. Only one triangle of the result needs to be written, so the update moves half as many bytes as the corresponding LU update and performs half as many flops. LAPACK’s potrf is organized around this split, and on machines where the matrix-matrix kernel runs near peak, the symmetric factorization is close to twice as fast as the general one for the same dimension.

One cost does not halve: the n square roots. On hardware where the square root is a multi-cycle instruction with limited throughput, a long sequence of them can become a bottleneck in the panel, which is one reason some implementations compute an L D L^T factorization with a unit diagonal and fold the square roots into a diagonal scaling applied later.

A failed square root as a verdict

The test s > 0 is the definiteness test. In exact arithmetic, s is the j-th pivot of symmetric elimination, and it is positive for every j if and only if A is positive definite. When the test fails, the matrix is not positive definite, and the algorithm stops at the first column where the property breaks.

The verdict is exact only in exact arithmetic. In floating point, s is computed with a relative error of order n * u times the magnitude of the terms being summed, and those terms can be much larger than s itself. For a matrix whose smallest eigenvalue is a small multiple of u * ||A||, the computed s can come out negative for a matrix that is positive definite, or positive for one that is not. The failure is confined to the boundary: for matrices with a comfortable margin of definiteness, the sign of the computed pivot is the sign of the true one.

Two consequences follow. First, a Cholesky factorization that completes is strong evidence of definiteness, and a breakdown is strong evidence against it, but neither is a proof at the edge of the representable range. Second, the factorization is often used as a cheap definiteness test in its own right, for example to check whether a Hessian is positive definite or whether a covariance matrix is valid, and in that role a breakdown at column j carries useful information about which leading submatrix is indefinite.

There is a common variant that avoids the square root entirely and makes the failure more informative: factor A = L D L^T with L unit lower triangular and D diagonal. The pivots d_j appear directly, without a square root, and the test becomes d_j > 0. The flop count is the same, the square roots are gone, and the pivots can be inspected individually. For a semidefinite matrix, the factorization proceeds until it reaches the rank of the matrix, at which point a pivot becomes zero and the corresponding column of L is undefined.

Indefinite matrices and the two-by-two pivot

Drop the definiteness assumption and keep symmetry, and the picture changes completely. A symmetric indefinite matrix can have a zero pivot while being perfectly invertible, so a pivot search becomes necessary. But a search restricted to single entries cannot bound the multipliers, because the symmetric structure forces the same permutation on rows and columns: swapping a large entry into the pivot position may require a column swap that moves an equally large entry out of the way.

Bunch and Kaufman resolved this by allowing 2 x 2 blocks on the diagonal of D. When the largest off-diagonal entry in the remaining column is comparable to the largest diagonal entry, the algorithm takes a 2 x 2 pivot involving both indices, which is always invertible for a symmetric matrix whose corresponding 2 x 2 submatrix is nonsingular, and which keeps the entries of L bounded by a fixed constant. The cost is bookkeeping: the factor is block diagonal rather than diagonal, and the triangular solves must handle 2 x 2 blocks.

The resulting factorization is backward stable with a growth bound independent of n, which is a stronger statement than partial pivoting provides for the general case. It is also the reason symmetric indefinite solvers expose an inertia computation: the numbers of positive, negative, and zero eigenvalues of D equal those of A by Sylvester’s law, so the factorization yields the inertia for free.

What the factorization does not report

Cholesky answers one question and no others. It says nothing about the conditioning of A, and the natural instinct to read the size of the pivots as a conditioning measure is unreliable: a matrix can have pivots of order one and a condition number of 1e12, because the small eigenvalue can be spread across many columns rather than concentrated in one. The reciprocal condition number must be estimated separately, from the factors, with a few triangular solves.

It also says nothing about structure that the factorization was not told about. A routine that reads only the lower triangle will happily accept an upper triangle that disagrees with it, which is a common source of silent errors when a matrix is assembled from blocks and only one triangle is filled correctly. And for a sparse matrix, the fill produced by the factorization can dwarf the original matrix; that cost is governed by the ordering, not by the arithmetic, and is the subject of Fill-In, Elimination Ordering, and Symbolic Factorization. Inside a sparse solver, the dense Cholesky kernel described here runs on the small dense fronts produced by the elimination tree, which is how Multifrontal Factorization and the Dense Front Matrix turns a sparse problem into a sequence of calls to this routine.