Cholesky decomposition

Definition

For a real symmetric positive-definite (SPD) matrix $A\in\mathbb{R}^{n\times n}$, the Cholesky decomposition, or Cholesky factorization, is

where $L$ is lower triangular with strictly positive diagonal entries. This factor exists and is unique. Equivalently, $A=R^TR$ with $R=L^T$ upper triangular. SPD means

See van de Geijn’s definition and existence theorem. For complex Hermitian positive-definite matrices, the form is $A=LL^{\ast}$, where $L^{\ast}$ is the conjugate transpose. LAPACK’s ZPOTRF documentation gives the complex convention.

This is a structured matrix decomposition. Unlike the unit lower triangular factor in LU decomposition, the Cholesky factor generally has diagonal entries other than $1$.

Intuition: positive energy becomes ordinary squared length

The question behind Cholesky is: can we rewrite a positive quadratic expression as a sum of squares, computing one square at a time? We will do that for a small example, then collect the coefficients into a matrix. That construction will give us the factor $L$.

A quadratic form is the scalar $q(x)=x^TAx$. Think of it as an energy or penalty assigned to a displacement $x$. For the real symmetric matrix \(A=\begin{pmatrix}4&2\\2&3\end{pmatrix}\),

What a cross term means

A cross term multiplies two different variables together. Here it is $4x_1x_2$. In contrast, $4x_1^2$ involves only $x_1$, and $3x_2^2$ involves only $x_2$.

The cross term describes an interaction: its contribution depends on both variables and their relative signs. Compare two inputs:

The squared terms have the same values in both cases. Only the cross term changes sign. This is what it means to say that the expression couples the coordinates: their joint contribution cannot be described using only separate squares of $x_1$ and $x_2$.

You can also see where the coefficient comes from in the matrix product. The two off-diagonal entries contribute $2x_1x_2$ and $2x_2x_1$, which add to $4x_1x_2$. The quadratic form page explains this correspondence between matrix entries and terms.

What completing the square means

Completing the square means rewriting a quadratic expression using a squared quantity, while keeping its value unchanged. Start with

For example, to rewrite $x^2+6x$, choose $t=3$ so that $2tx=6x$. The square supplies an extra $9$, so subtract it:

Adding and subtracting the same amount preserves the expression. See OpenStax’s explanation of completing the square.

For our two-variable example, notice that

This square accounts for the first two terms and one of the three copies of $x_2^2$. Two copies remain, giving

The expression is now a sum of squares of combinations of the original variables.

Why completing squares in order creates a triangular matrix

The important choice is to finish with one variable before moving to the next. We first collect every term containing $x_1$ into the square $(2x_1+x_2)^2$. Subtracting that square from the original expression leaves

There is no $x_1$ left in this remainder. Consequently, the next square, $(\sqrt2x_2)^2$, has an $x_1$ coefficient of zero. The coefficient rows of the two squares therefore have the pattern

With three variables, the same procedure first absorbs all terms containing $x_1$, leaving a quadratic expression in $x_2,x_3$. Completing the next square absorbs all terms containing $x_2$ from that remainder, leaving only a multiple of $x_3^2$. The resulting combinations have the form

Here $r_{ij}$ denotes the coefficient of $x_j$ in the $i$th combination. Each later combination omits the variables already absorbed into earlier squares. Those missing coefficients are precisely the zeros below the diagonal of an upper triangular coefficient matrix. Transposing it gives the lower triangular factor $L$.

Thus triangularity comes from completing squares in variable order. An arbitrary sum-of-squares representation need not have this pattern. Positive definiteness ensures that each successive remainder supplies a positive coefficient for the next square; the proof below explains why this process can continue.

How Cholesky records these combinations

Define two new quantities:

Then $q(x)=y_1^2+y_2^2$, the ordinary squared length of the vector $y$. For instance, $x=(1,1)^T$ gives $y=(3,\sqrt2)^T$, so the same value $11$ becomes $3^2+(\sqrt2)^2$.

Collect the coefficients of these two combinations into a matrix, which we will call $R$:

The entries are the coefficients in the two equations for $y$. The zero below the diagonal records the step we just took: after absorbing all $x_1$ terms into the first square, the second square uses only $x_2$. This is the upper triangular pattern described above.

Now multiply its transpose by itself:

We have performed the factorization by completing the square. To use the lower triangular convention in the definition, name the transpose $L=R^T$:

This is where $L$ comes from: its columns contain the coefficients of the squares we just constructed. Cholesky carries out this idea systematically for an SPD matrix, choosing positive diagonal coefficients. The general algorithm below makes that sequential construction precise; see the FLAME definition and algorithm.

This gives an order for recovering the original variables. First recover $x_2$ from the second equation, then substitute it into the first:

The pair $(y_1,y_2)$ is another way to describe the same input, so we call it a new set of coordinates. Calculating $y$ from $x$ is a coordinate transformation. The word triangular describes the matrix and its dependency pattern; it does not describe a special kind of coordinate or a triangular shape.

In matrix notation, the sum of squares is

Here $L^Tx=y$, and the squared-norm notation on the right means $y_1^2+y_2^2$. Since we can recover $x$ from $y$, a nonzero $x$ cannot give $y=(0,0)^T$. At least one square is therefore positive. This explains the SPD condition in this example; see Higham’s discussion of positive-definite quadratic forms.

This also connects the factorization to the contour plots below. In the original $x$ coordinates, $q(x)=4$ is an ellipse. After applying $y=L^Tx$, the same condition becomes $y_1^2+y_2^2=4$, a circle of radius $2$. The factor supplies the change of coordinates that turns the elliptical energy levels into ordinary circles. The contour shape motivates this goal; completing the square tells us how to achieve it.

Positivity concerns directions rather than individual entries

Positive definite does not mean that every entry of the matrix is positive. For example, \(\begin{pmatrix}2&-1\\-1&2\end{pmatrix}\) gives

Conversely, \(\begin{pmatrix}1&2\\2&1\end{pmatrix}\) has positive entries but gives $x^TAx=-2$ for $x=(1,-1)^T$. The condition concerns every direction, not the signs of individual entries.

The semidefinite example has a direction that costs no energy. The indefinite example has a direction with negative energy, so it cannot equal a squared Euclidean length for all $x$. These geometric differences correspond to zero or invalid square-root pivots in the usual algorithm.

Another way to read the factor we constructed

We already found $L$ from the completed squares. We can now read its entries in another way to understand the computation. Matrix multiplication takes a row of $L$ and a column of $L^T$, which is a row of $L$ turned sideways. Thus each entry of $LL^T$ is a dot product of two rows of $L$ (multiply corresponding coordinates and add). Calling those rows $r_1,\ldots,r_n$, we have $a_{ij}=r_i\cdot r_j$. In our example,

Their squared lengths are $4$ and $3$, and their dot product is $2$. These are precisely the two diagonal entries and the off-diagonal entry of $A$. Triangularity lets us construct these vectors one coordinate at a time: give the first row its length, match the next row’s dot product with it, then choose the remaining coordinate to supply the remaining squared length. This is the Gram-matrix interpretation of $LL^T$; see the FLAME definition and algorithm.

The connection to orthogonalization

If $B=QR$ is a reduced QR factorization with independent columns and positive diagonal $R$, then

Thus $R^T$ is the lower Cholesky factor of the Gram matrix $B^TB$. Gram–Schmidt constructs $Q$ by removing overlapping components and records their coefficients in $R$. This identity links the factorizations in exact arithmetic; the QR accuracy discussion explains why forming the Gram matrix first can lose information through rounding.

A triangular square root has a specific meaning

$A=LL^T$ resembles $a=\sqrt a\sqrt a$, but generally $L^2\ne A$. The transpose matters. In our example, the upper-right entry of $L^2$ is zero, while $a_{12}=2$. A symmetric matrix square root satisfying $S^2=A$ is a different factorization. Cholesky chooses a triangular factor because it supports sequential computation and solves.

Deriving the algorithm

Equating entries in $A=LL^T$ gives

Compute columns in order, for $j=1,\ldots,n$:

then, for $i=j+1,\ldots,n$,

Empty sums are zero. Entries above the diagonal remain zero. These recurrences are the lower triangular counterpart of the elementwise Cholesky algorithm presented by Higham.

Each diagonal entry supplies what is still missing

For a two-by-two matrix, write

Matching $LL^T$ to $A$ yields three requirements:

First choose $\ell_{11}=\sqrt a$. Next choose $\ell_{21}=b/\sqrt a$ to match the shared component. Finally choose $\ell_{22}=\sqrt{c-b^2/a}$ to supply the squared length that the first coordinate has not supplied.

For $a=4$, $b=2$, $c=3$, the second row needs total squared length $3$. Its first coordinate is $1$, already contributing $1$. The remaining coordinate must contribute $2$, so it is $\sqrt2$. This makes the subtraction inside the recurrence concrete: subtract the contribution already accounted for, then take the square root of what remains.

The same remainder appears when completing the square:

For $a>0$, a positive remainder means the second direction still has positive energy. A zero remainder leaves a flat direction. A negative remainder gives a negative-energy direction. In higher dimensions, the Schur complement is the matrix version of this remainder.

Why the square roots are positive

Partition an SPD matrix as

Here $\alpha>0$. For nonzero $z$, set $x=(-c^Tz/\alpha,z)^T$. Then

Thus the Schur complement $B-cc^T/\alpha$ is SPD. Inductively factor it as $KK^T$, and obtain

The positive square root and each successive column are uniquely determined. This gives an inductive existence and uniqueness proof; see the FLAME notes for a full treatment.

Worked example

Let

It is SPD because, for $x\ne0$,

The recurrences give

Therefore

To solve $Ax=b$ for $b=(6,5)^T$, first solve $Ly=b$:

Then solve $L^Tx=y$ to obtain $x_2=1$ and $x_1=1$. Direct substitution confirms $A(1,1)^T=b$.

A second example: shared and independent variation

Let $z_1,z_2$ be independent random variables with mean zero and variance one, and set

Then $\operatorname{Var}(u_1)=4$, $\operatorname{Var}(u_2)=1+2=3$, and $\operatorname{Cov}(u_1,u_2)=2$. Both outputs share $z_1$, which produces covariance. The independent $z_2$ supplies the remaining variance of the second output. In matrix notation, $u=Lz$ and $\operatorname{Cov}(u)=LL^T=A$. See Stata’s covariance-transformation explanation.

Reading the same factor in different ways gives different intuitions: $L^T$ removes the coupling from a quadratic energy, while $L$ introduces a desired covariance into independent inputs. The algebraic identity is the same, but the direction of the transformation and the problem being solved differ.

What a nearly zero remainder tells us

For

positive $\varepsilon$ gives an SPD matrix. At $\varepsilon=0$, the two rows represent the same underlying direction and the matrix is singular. For small positive $\varepsilon$, the second independent contribution is tiny: the system is close to losing a direction. A triangular solve divides by $\sqrt\varepsilon$, explaining why small changes can matter greatly. This example illustrates the distinction between existence and numerical sensitivity discussed below.

Applications and computational cost

For an SPD linear system, solve $Ly=b$ by forward substitution, then $L^Tx=y$ by backward substitution. Reuse $L$ for additional right-hand sides. LAPACK’s DPOTRS routine implements this factor-and-solve workflow.

The determinant follows from the product:

See Higham’s determinant identity; the logarithmic formula follows by taking logs of positive diagonal entries.

Conventional dense Cholesky costs $\frac13 n^3+O(n^2)$ floating-point operations, counting multiplication and addition separately, and two triangular solves cost $O(n^2)$ per right-hand side. The cubic leading term is half that of conventional dense LU. Only one triangular factor needs storage. See the cost discussion in the FLAME notes.

Numerical behavior and limitations

For SPD input, standard Cholesky needs no pivoting and has small backward error: the computed factor $\widehat L$ satisfies $A+\Delta A=\widehat L\widehat L^T$ for a small perturbation $\Delta A$. An ill-conditioned system can still have a large solution error. In floating-point arithmetic, a nearly singular SPD matrix may cause a nonpositive computed pivot. Higham discusses stability and breakdown.

Symmetry alone is insufficient. An indefinite matrix does not have the positive-diagonal Cholesky factor defined above. Positive-semidefinite matrices can have triangular square-root factors with zero diagonal entries, but uniqueness and the standard algorithm’s divisions no longer follow. See Higham’s semidefinite discussion. For symmetric indefinite systems, LAPACK uses pivoted factorizations with block diagonal $D$ instead.

History

The original problem was adjusting a measured network

André-Louis Cholesky was a French artillery officer working in the army’s geographic service. When the French triangulation network was being revised, officers needed ways to adjust the measured network with greater speed, convenience, and precision. Benoît’s account of Cholesky’s work identifies this practical problem as the setting for his least-squares calculation procedure. See Benoît’s account, translated by Richard W. Cottle.

What does adjusting a network mean? Consider a simplified surveying example in a flat plane. Three measured triangle angles might add up to $180^\circ+\varepsilon$, although the geometric condition requires $180^\circ$. Corrections $v_1,v_2,v_3$ must satisfy

There are infinitely many possible corrections. If the observations have equal weight, minimizing $v_1^2+v_2^2+v_3^2$ gives $v_1=v_2=v_3=-\varepsilon/3$. A network has many conditions sharing the same observations, so the corrections must be calculated together. This example illustrates the adjustment problem; an actual geodetic model also accounts for the Earth’s curvature and the geometry of the network.

The original question was therefore: how can we calculate small corrections that make a whole measured network consistent, with a procedure people can carry out reliably? Least squares specifies the desired corrections. Cholesky’s contribution concerns the numerical work needed to obtain them.

Why this produced normal equations

The title of Benoît’s 1924 paper refers to fewer equations than unknowns. This describes the condition equations for observation corrections, rather than the overdetermined regression problem commonly used to introduce least squares today. De Falguerolles, §5.2, explains this formulation.

In modern notation, let $v\in\mathbb R^m$ contain corrections and let $C\in\mathbb R^{p\times m}$ contain $p<m$ independent linearized conditions. For equal observation weights, solve

Introduce Lagrange multipliers $\lambda\in\mathbb R^p$:

Stationarity and the constraints give

Thus an underdetermined correction problem becomes a square system for the multipliers. Its coefficient matrix $N=CC^T$ is symmetric. Full row rank makes it positive definite, since

For positive observation weights collected in a diagonal matrix $W$, minimizing $\frac12v^TWv$ instead gives $N=CW^{-1}C^T$ and $v=W^{-1}C^T\lambda$. This follows from the same stationarity calculation.

In regression, the familiar route is instead

Both routes produce a Gram matrix, positive definite under the appropriate rank assumption. The historical distinction is what the unknowns represent: observation corrections and their multipliers in network adjustment, or fitted parameters in regression. See de Falguerolles, §§5.1–5.2.

The 1910 method organized the solve into two stages

Cholesky’s surviving manuscript is dated 2 December 1910, establishing that the method existed by then. The exact invention date is unknown. See Brezinski (2005), p. 217.

The following uses modern notation for the construction described in Brezinski (2006), pp. 283–286. Instead of solving $N\lambda=d$ directly, find a lower triangular $L$ such that

Introduce an intermediate vector $y=L^T\lambda$. The calculation becomes

The first equation gives $y_1$, then $y_2$, and so on. The second gives $\lambda_p$, then $\lambda_{p-1}$, and so on. Matching entries of $N$ and $LL^T$ supplies the square-root recurrences given earlier in this article. The factorization organizes a coupled system into successive scalar calculations.

Cholesky also described the calculation layout, use of a Dactyle mechanical calculator, and checks for mistakes. These details mattered because the work was performed by people, including soldiers with limited mathematical training. His discussion already included finite-precision errors; concern for numerical reliability was present before electronic computers. See Brezinski (2006), pp. 285–286.

Publication in 1924 and wider adoption after 1944

Cholesky died in 1918. Commandant Benoît published the method in 1924, six years later. The publication date differs from the manuscript date. See Cottle’s translation and publication note.

Henry Jensen’s 1944 comparison of methods for normal equations reached John Todd, who used it in his 1946 course at King’s College London. Leslie Fox brought the method to the National Physical Laboratory. Work by Fox, Huskey, Wilkinson, and Turing in 1948 developed its numerical analysis, followed by dissemination through textbooks. See Brezinski (2005), pp. 216–217. This connects with the postwar accuracy questions discussed in the LU history.

The scope broadened and the accuracy question became more precise

Cholesky’s manuscript already considers general square systems. Its construction, written in modern notation for nonsingular $B$, is

Solve the symmetric system for $\lambda$, then recover $x$. See Cholesky’s manuscript, systems I–III, pp. 221–223.

Modern numerical analysis distinguishes algebraic equivalence from computational accuracy. For nonsingular square $B$,

Forming normal equations can therefore make the computation much more sensitive to rounding. A stable Cholesky solve of the resulting matrix does not undo this effect. QR or SVD can avoid forming the Gram matrix when that loss of accuracy matters. See ALAFF, §4.2.5 and LAPACK’s least-squares methods.

The accuracy question also became more explicit: can a computed solution be interpreted as the exact solution of a nearby problem, and how does conditioning turn that perturbation into solution error? LAPACK’s backward-error discussion expresses these as separate quantities. This refines the original aim of reliable calculation rather than replacing it.

Electronic computing added another implementation problem: moving data through memory efficiently. LAPACK reorganizes algorithms into blocks so that much of the work uses matrix-matrix operations. The identity $A=LL^T$ remains the same, while the arrangement of the computation changes with the machine.

The factor became useful beyond solving equations

A modern Cholesky factor can also represent a covariance structure. Suppose $\Sigma=LL^T$ is a positive-definite covariance matrix and $z$ has independent standard normal components. Then

Here the factor generates correlated observations rather than solving for unknowns. Stata’s explanation of multivariate-normal simulation gives this construction. It demonstrates how the same triangular factor acquired additional uses without changing its mathematical definition.

The historical thread is thus a continuing problem of usable computation: first, obtain least-squares corrections for measured networks; then, establish precise accuracy guarantees and efficient implementations; and, in additional applications, use the factor itself as a representation of a positive-definite matrix.

Interactive example: vary the remaining energy

The experiment uses symmetric \(A=\begin{pmatrix}a&b\\b&c\end{pmatrix}\). Predict the sign of $c-b^2/a$ before changing the shared coefficient $b$. When the factorization fails, find a direction with zero or negative energy. A failure here is a consequence of the chosen matrix, rather than a reason to take the square root of a negative number. The calculations follow the Cholesky recurrences in the FLAME notes.

Practice problems

Try each problem before opening its hint or solution. The first two build hand-calculation fluency, the next two ask you to interpret or choose an approach, the next two develop logical reasoning, and the final two connect the concept to practical work. These are original exercises; the references support the mathematical concepts they use.

1. Hand calculation: factor and solve

Find $L$ with positive diagonal, solve $Ax=b$, and compute $\det A$ for

2. Hand calculation: construct three rows

Compute the Cholesky factor and determinant of

3. Thinking: diagnose the failed square root

Compare \(A_1=\begin{pmatrix}1&2\\2&5\end{pmatrix}\) and \(A_2=\begin{pmatrix}1&2\\2&3\end{pmatrix}\). Which admits positive-diagonal Cholesky? For the other, exhibit a negative-energy direction.

4. Thinking: interpret correlation near one

For the covariance matrix \(A_\rho=\begin{pmatrix}1&\rho\\\rho&1\end{pmatrix}\) with $0\le\rho<1$, find $L$. Interpret its second diagonal when $\rho=0.9$ and as $\rho\to1$.

5. Logic: prove positivity and examine uniqueness

If $L$ is an invertible real lower triangular matrix, prove $LL^T$ is SPD. Why does the uniqueness theorem require a positive diagonal?

6. Logic: identify the factor of a leading block

If $A=LL^T$ is an SPD Cholesky factorization, show that the upper-left $k\times k$ block of $L$ is the Cholesky factor of the upper-left $k\times k$ block of $A$.

7. Applied: generate correlated sensor errors

Two sensors need zero-mean errors with covariance \(\Sigma=\begin{pmatrix}4&1\\1&1\end{pmatrix}\). Construct the errors from independent zero-mean, unit-variance variables $z_1,z_2$. Verify the two variances, covariance, and correlation.

8. Applied: minimize a quadratic penalty

An optimization step minimizes

Derive the system for its minimizer and solve it with Cholesky. Why is the stationary point the unique minimum?

Further learning

• FLAME: Notes on Cholesky Factorization: study the full proof, algorithm, operation count, and accompanying exercises.
• Higham: What Is a Symmetric Positive Definite Matrix?: connect energy, eigenvalues, and positivity conditions.
• ALAFF: why normal equations can lose accuracy: understand why a stable factorization is only one part of a reliable least-squares computation.
• Stata: multivariate-normal simulation: connect exercise 7 to a complete covariance-generation example.

Related pages

• QR factorization

• Gram–Schmidt orthogonalization

• Quadratic form

• Symmetric positive definite matrix

• Matrix factorization and decomposition
• LU decomposition

References

1. van de Geijn, R. A. (2011). Notes on Cholesky Factorization. University of Texas at Austin. Definition, proof, algorithm, and operation count.
2. Higham, N. J. (2020, August 11). What Is a Cholesky Factorization?. Determinant identity, numerical stability, and semidefinite extensions.
3. Anderson, E., et al. (1999). LAPACK Users’ Guide, third edition. SIAM. Linear Equations. Real and complex factorizations and triangular solves.
4. LAPACK (2011). ZPOTRF and DPOTRS, version 3.3.1. Complex Cholesky factorization and solves using an existing factor.
5. Cholesky, A.-L. (1910, December 2). Sur la résolution numérique des systèmes d’équations linéaires. Manuscript reproduced in Brezinski (2005), §4, pp. 221–229.
6. Benoît, Commandant (1924). Note sur une méthode de résolution des équations normales provenant de l’application de la méthode des moindres carrés à un système d’équations linéaires en nombre inférieur à celui des inconnues (Procédé du Commandant Cholesky). Bulletin géodésique, 2, 67–77. DOI. Publication history described in Brezinski (2005).
7. Benoît, Commandant. Account of Cholesky’s life and geodetic work, translated by R. W. Cottle in 1975 and circulated in NA Digest, February 18, 1990.
8. Brezinski, C. (2005). La méthode de Cholesky. Revue d’histoire des mathématiques, 11(2), 205–238. Full text. Historical investigation and reproduction of the original manuscript.
9. Brezinski, C. (2006). The life and work of André Cholesky. Numerical Algorithms, 43, 279–288. Author-uploaded full text.
10. de Falguerolles, A. (2019). Cholesky and the Cholesky decomposition: a commemoration by an applied statistician. Journal de la Société Française de Statistique, 160(2), 83–96. Sections 5.1–5.2 distinguish fitting and compensation.
11. van de Geijn, R. A., and Myers, M. Advanced Linear Algebra: Foundations to Frontiers. §4.2.5, Why using the Method of Normal Equations could be bad.
12. Anderson, E., et al. (1999). LAPACK Users’ Guide. Linear Least Squares Problems and Further Details: Error Bounds for Linear Equation Solving. See also the LAPACK project description for blocked computation.

13. StataCorp. Stata 6: Simulating multivariate normal observations. Cholesky factors as covariance transformations.

14. Higham, N. J. (2020, July 21). What Is a Symmetric Positive Definite Matrix?. Quadratic-form interpretation and positivity conditions.
15. OpenStax. College Algebra 2e, §2.5, Quadratic Equations, subsection Completing the Square. Algebraic background for rewriting a quadratic using a squared expression.