Cholesky decomposition

1 Matrix prerequisites and definitions

Cholesky decomposition applies to a restricted class of matrices. The input matrix is typically real, symmetric, and at least positive semidefinite; with these properties, it can be represented as a product of a triangular matrix and its transpose (or an equivalent “square-root” form).

1.1 Symmetry and transpose structure

Let \(A\) be a real matrix. The classic Cholesky factorization requires that \(A\) equals its transpose, \(A=A^\top\). Under this condition, the triangular factor inherits a consistent structure: the factor multiplies with its transpose to reproduce \(A\), preserving both diagonal contributions and cross terms in a way that matches symmetry.

1.2 Positive definiteness vs. positive semidefiniteness

A symmetric matrix \(A\) is positive definite if \(x^\top A x>0\) for all nonzero vectors \(x\). It is positive semidefinite if \(x^\top A x\ge 0\) for all \(x\). Standard Cholesky, without modifications, corresponds to positive definiteness. When \(A\) is only semidefinite, some diagonal pivots can become zero, motivating pivoted or generalized variants.

1.3 Uniqueness of the Cholesky factor

For a symmetric positive definite matrix \(A\), the Cholesky factor \(L\) (taken lower triangular) is unique when its diagonal entries are constrained to be positive. If the diagonal sign convention is altered, multiple equivalent triangular factors can exist, but they differ only by predictable sign changes on diagonal elements.

1.4 Notation conventions (\(LL^\top\) vs. \(R^\top R\))

Two equivalent conventions are common. One writes \[ A=LL^\top \] with \(L\) lower triangular. Another writes \[ A=R^\top R \] with \(R\) upper triangular. The choice mostly affects how algorithms and back-substitution steps are presented; computationally, they correspond to the same underlying factor.

2 Statement of the decomposition

Cholesky decomposition is a structured factorization that converts a symmetric (semi)definite matrix into triangular factors, enabling stable computations for solving and related tasks.

2.1 The standard Cholesky factorization

For a symmetric positive definite matrix \(A\in \mathbb{R}^{n\times n}\), Cholesky decomposition finds a lower triangular matrix \(L\) with strictly positive diagonal entries such that \[ A = LL^\top. \] The entries of \(L\) can be constructed recursively, using previously computed elements to determine each new diagonal term and the corresponding subdiagonal row/column entries.

2.2 Strict vs. pivoted variants

When \(A\) is positive definite, the strict (non-pivoted) algorithm proceeds without breakdown because each computed diagonal pivot is positive. For semidefinite or nearly indefinite matrices, pivoted or modified methods may reorder rows/columns, or adjust the factorization approach, to avoid failure or reduce the impact of small pivots.

2.3 Conditions under which the factor exists

A Cholesky factorization \(A=LL^\top\) with real \(L\) and lower triangular form exists precisely under definiteness conditions aligned with the algorithm’s assumptions. For positive definite matrices it exists without pivoting. For positive semidefinite matrices, an exact factor may exist but may require generalized formulations (e.g., allowing zero diagonals) and, in practice, careful numerical handling.

3 Computational algorithm (factor-by-factor construction)

Cholesky decomposition can be viewed as building the factor \(L\) one column (or row) at a time, using algebraic identities that mirror the structure of the product \(LL^\top\).

3.1 Deriving the update formulas

Write \(A=LL^\top\) and consider the \((i,j)\) entry: \[ A_{ij} = \sum_{k=1}^{\min(i,j)} L_{ik}L_{jk}. \] When constructing \(L\) in a forward manner, previously computed terms allow one to solve for an unknown diagonal element \(L_{ii}\) and subdiagonal elements \(L_{ij}\) for \(j<i\). Diagonal entries satisfy \[ L_{ii} = \sqrt{A_{ii}-\sum_{k=1}^{i-1} L_{ik}^2}, \] and for \(j<i\), \[ L_{ij} = \frac{1}{L_{jj}}\left(A_{ij}-\sum_{k=1}^{j-1} L_{ik}L_{jk}\right). \] These formulas are the backbone of standard implementations.

3.2 Pseudocode for the basic algorithm

A common schematic for the non-pivoted factorization (lower triangular \(L\)) is:

For \(i=1\) to \(n\):

Compute \(L_{ii}\) from \(A_{ii}\) minus the sum of squares of previously computed \(L_{ik}\) for \(k<i\).
For \(j=i+1\) to \(n\):
Compute \(L_{j i}\) using \(A_{j i}\) minus the dot product of previously computed parts, divided by \(L_{ii}\).

In practice, the loops are arranged to improve memory access patterns.

3.3 Complexity and operation counts

The cost of Cholesky for an \(n\times n\) dense matrix is on the order of \(O(n^3)\) floating-point operations, with a leading term comparable to roughly one-third of the work required by general LU factorization for similar dimensions. Because only the lower triangular part of \(L\) is stored and computed, constants can be favorable relative to less structured factorizations.

3.4 Numerical stability considerations

Cholesky is often more numerically stable than elimination methods for symmetric positive definite problems because it leverages positivity of pivots and avoids certain forms of cancellation. Still, numerical errors can be influenced by conditioning: if \(A\) has eigenvalues that differ greatly in magnitude, rounding errors may be amplified. Ensuring appropriate scaling and using stable floating-point arithmetic are key.

4 Pivoting and generalized Cholesky

When the matrix fails to be strictly positive definite or when floating-point effects lead to small or negative pivots, generalized approaches broaden applicability.

4.1 Need for pivoting (rank deficiency and semidefiniteness)

If \(A\) is positive semidefinite, some diagonal terms in the Cholesky construction can become zero, corresponding to directions in which \(x^\top A x=0\). If the matrix is close to semidefinite or contaminated by numerical noise, naive factorization may encounter negative values under square roots. Pivoting and modifications attempt to find a stable decomposition by reordering or relaxing strictness.

4.2 Incomplete and modified Cholesky approaches

For large sparse problems, incomplete Cholesky computes only an approximation by dropping certain fill-in terms, producing factors used as preconditioners. Modified variants adjust the factorization to handle indefiniteness or improve robustness, sometimes by shifting diagonals or applying heuristics to maintain nonnegative pivots.

4.3 Relation to LDLᵀ factorization

A closely related family is \(A=LDL^\top\), where \(L\) is unit lower triangular and \(D\) is block diagonal (or diagonal in the simplest case). This representation can naturally accommodate semidefinite structure through zero entries in \(D\). Cholesky can be recovered as a special case where \(D\) is positive and diagonal, allowing square roots to be absorbed into the triangular factor.

4.4 Rank-revealing variants

Rank-revealing methods aim to identify numerical rank and expose nearly singular behavior. In practice, these algorithms may pivot to maximize diagonal growth or track pivot magnitudes, producing a factorization that is stable even when the effective rank is less than \(n\).

5 Solving linear systems using Cholesky

Once \(A\) is factored as \(LL^\top\), solving systems becomes a structured pair of triangular solves, typically faster and more reliable than direct elimination for appropriate matrices.

5.1 Forward substitution with \(L\)

To solve \(Ax=b\) with \(A=LL^\top\), first set \[ LL^\top x=b. \] Compute an intermediate vector \(y\) from \[ Ly=b \] using forward substitution. Because \(L\) is lower triangular, each component of \(y\) depends only on earlier computed entries.

5.2 Back substitution with \(L^\top\)

Then solve \[ L^\top x=y \] by back substitution. Since \(L^\top\) is upper triangular, each component of \(x\) depends on already computed later components, proceeding from the last index to the first.

5.3 Multiple right-hand sides

If multiple vectors \(b^{(1)},\dots,b^{(m)}\) share the same coefficient matrix \(A\), the factorization \(A=LL^\top\) can be reused. The additional cost is then primarily the sequence of triangular solves for each right-hand side, which can be batched efficiently in optimized implementations.

5.4 Error propagation in computed solutions

Rounding errors originate both in the factorization and in the triangular solves. For well-conditioned positive definite matrices, the overall backward error is often modest. For ill-conditioned problems, small perturbations in \(A\) (as reflected in the computed factor) can lead to larger deviations in the solution, aligning with standard conditioning theory.

6 Least squares and Gaussian models

Cholesky frequently appears when rewriting least-squares systems into forms that are symmetric and suitable for triangular solvers, as well as in probabilistic models involving covariance matrices.

6.1 Cholesky for normal-equation solving

A common route from least squares to Cholesky is via the normal equations. For an overdetermined linear model minimizing \(\\|Ax-b\\|_2\), the normal equations produce

\[ (A^\top A)x = A^\top b, \] where \(A^\top A\) is symmetric positive semidefinite and often positive definite when columns of \(A\) are independent. If positive definite, Cholesky can factor \(A^\top A\) and enable efficient solution.

6.2 Connection to covariance matrices

In statistics, covariance matrices are symmetric positive semidefinite. When a covariance matrix \(\Sigma\) admits a Cholesky factor \(\Sigma=LL^\top\), it directly supports computation of Mahalanobis-type quantities and enables numerically stable transformations between standard normal variables and correlated Gaussian variables.

6.3 Solving constrained or structured least squares

Cholesky can be integrated into broader least-squares workflows when the problem yields a symmetric (often block) normal system or when constraints can be encoded so that the reduced system remains compatible with triangular factorization. Structured problems, such as those with banded or sparse patterns, benefit from storage and ordering tailored to sparsity.

6.4 Conditioning and practical solver workflow

While using normal equations can be efficient, it can also worsen conditioning because forming \(A^\top A\) squares the singular values. In practice, many workflows prefer QR-based or SVD-based least squares for improved numerical reliability, using Cholesky where structure or conditioning is favorable, or when the problem size and constraints make it a good trade-off.

Triangular factors provide convenient pathways to compute determinants, solve related linear systems, and evaluate functions that depend on \(A\).

7.1 Computing \(\det(A)\) from the Cholesky factor

If \(A=LL^\top\), then \[ \det(A)=\det(L)\det(L^\top) = (\prod_{i=1}^n L_{ii})^2. \] Thus the determinant can be obtained from the diagonal of \(L\) without forming \(A^{-1}\), though care is needed for overflow or underflow when products become large or small.

7.2 Matrix inversion strategies via \(L\)

Although explicitly forming \(A^{-1}\) is often unnecessary and can be numerically less efficient, it can be computed by solving \(A X = I\) using the Cholesky factor: for each column of \(X\), solve two triangular systems with \(L\) and \(L^\top\). In many applications, more economical alternatives exist, such as computing \(A^{-1}b\) for specific vectors.

Quantities such as \(\mathrm{tr}(A^{-1}B)\) arise in uncertainty quantification and sensitivity analysis. Using triangular solves, one can compute products involving \(A^{-1}\) without explicit inversion, sometimes combined with stochastic trace estimation to avoid costly exact computations.

7.4 Log-determinants and numerical scaling

The log-determinant is widely used in optimization and probabilistic scoring. From \(A=LL^\top\), \[ \log\det(A) = 2\sum_{i=1}^n \log L_{ii}. \] Using sums of logarithms improves numerical robustness compared with directly multiplying diagonal entries, especially for high-dimensional matrices.

8 Applications in numerical methods

Beyond basic linear solves, Cholesky underpins algorithms where symmetric positive definite structure is exploited for speed or stability.

8.1 Iterative methods preconditioning (conceptual overview)

In iterative solvers, preconditioners approximate the inverse of the operator to accelerate convergence. Cholesky-based preconditioners, including incomplete variants for sparse systems, are popular because triangular solves are relatively cheap and the factorization captures the dominant curvature of the matrix.

8.2 Optimization and quadratic forms

Many optimization problems involve quadratic objectives or constraints leading to linear systems with symmetric positive definite Hessians (or approximate Hessians). Cholesky factorization provides a direct method to solve Newton-like systems or evaluate quadratic forms such as \(x^\top A x\) in a numerically controlled manner.

8.3 Sampling from multivariate normal distributions

To sample \(z\sim \mathcal{N}(0,\Sigma)\), a standard technique is to generate a vector \(u\sim \mathcal{N}(0,I)\) and set \(z=Lu\) when \(\Sigma=LL^\top\). This yields samples with the intended covariance while avoiding repeated expensive matrix operations.

8.4 Scientific computing workflows

Scientific workflows often require repeated solution of the same structural system across timesteps, parameter values, or multiple experimental conditions. When the coefficient matrices are symmetric positive definite (or can be regularized to be so), Cholesky offers a practical backbone for assembling pipelines that involve factorization once and many solves thereafter.

Cholesky sits within a broader ecosystem of matrix factorizations that share goals—stability, efficiency, and structural exploitation.

9.1 LDLᵀ decomposition

The \(LDL^\top\) factorization generalizes Cholesky by separating sign and scaling information into the diagonal (or block diagonal) matrix \(D\). This family is suited to symmetric indefinite matrices and supports semidefinite cases more directly, particularly in algorithms that avoid square roots.

9.2 Square-root and QR comparisons

A “square-root” viewpoint emphasizes that Cholesky computes a matrix square root in the form \(LL^\top\). QR factorization, by contrast, is often preferred for least squares because it avoids squaring condition numbers. The comparison is therefore application dependent: Cholesky is excellent for SPD systems, while QR or SVD may outperform it for general least squares.

9.3 Spectral and congruence viewpoints

Symmetric positive definite matrices admit eigenvalue decompositions and congruence transformations. Cholesky can be interpreted as producing a particular congruence factorization consistent with triangular structure, effectively providing a computationally efficient substitute for full spectral methods in many engineering contexts.

9.4 Cholesky of block matrices

When \(A\) has a block structure, block Cholesky can exploit conditional independence between parts of the system. Such approaches are common in sparse graphical models, domain decomposition, and systems where fill-in is limited by ordering strategies, allowing sub-block factorization and Schur complement updates.

10 Practical considerations and implementation details

Successful use depends on numerical and software-level details: data formats, ordering, and careful treatment of borderline cases.

10.1 Choice of data types and scaling

Floating-point precision (e.g., single vs. double) affects robustness, especially for ill-conditioned matrices. Scaling the problem—normalizing magnitudes of rows/columns or applying diagonal preconditioning—can reduce overflow risks and improve pivot quality, indirectly enhancing stability of the triangular factors.

10.2 Handling ill-conditioned matrices

When \(A\) is near singular, small pivots can lead to large relative errors in the computed \(L\). Practical responses include adding a regularization term (diagonal loading), using pivoted/generalized factorizations, or switching to alternative methods better suited to near-rank-deficiency.

10.3 Cache-friendly ordering and storage formats

For dense matrices, implementations aim to maximize cache reuse by arranging loops and data layout. For sparse matrices, storage formats such as compressed sparse row (CSR) or compressed sparse column (CSC) interact with fill-in patterns. Ordering strategies (e.g., permutations) help reduce fill-in and thus reduce both runtime and memory usage.

10.4 Common implementation pitfalls

Typical pitfalls include treating the matrix as exactly symmetric when it is only approximately so, neglecting to verify positive definiteness (or semidefiniteness) before factorization, and failing to account for cases where a zero or negative pivot occurs under floating-point arithmetic. Another frequent issue is using a factorization result in ways inconsistent with the chosen convention (lower vs. upper triangular).

11 Worked examples

Concrete calculations illustrate how the factorization is constructed and how it supports solving systems.

11.1 Decomposing a small 2×2 matrix

Consider a symmetric positive definite matrix \[ A=\begin{pmatrix} a & b \\ b & c \end{pmatrix}. \] Seek a lower triangular \[ L=\begin{pmatrix} \ell_{11} & 0 \\ \ell_{21} & \ell_{22} \end{pmatrix} \] such that \(A=LL^\top\). Matching entries gives: \[ \ell_{11}=\sqrt{a},\quad \ell_{21}=\frac{b}{\ell_{11}},\quad \ell_{22}=\sqrt{c-\ell_{21}^2}. \] The final square root requires \(c-b^2/a>0\), aligning with positive definiteness.

11.2 Decomposing a symmetric positive definite system

For a numeric SPD matrix, the procedure applies the same recursion: compute the first diagonal pivot from the corresponding diagonal entry, then fill in subdiagonal elements using division by that pivot, and continue inward. Each step uses only already computed entries, so the method is directly implementable with a nested loop structure.

11.3 Using the factor to solve \(Ax=b\)

After computing \(A=LL^\top\), solving \(Ax=b\) proceeds in two stages. First compute \(y\) from \(Ly=b\) via forward substitution. Then compute \(x\) from \(L^\top x=y\) via back substitution. This approach avoids forming \(A^{-1}\) and turns the dense system solve into a pair of triangular solves.

11.4 Example with semidefinite/pivoted behavior

If \(A\) is positive semidefinite, some diagonal pivot candidates can become zero, reflecting an exact or numerical rank deficiency. In such situations, a strict algorithm may fail due to a zero pivot appearing in denominators or due to roundoff producing a tiny negative value under a square root. Pivoted or generalized methods address this by reordering the system or using an \(LDL^\top\)-style framework that separates the problematic scaling into \(D\).