Introduction
Theory
HOWTO
Error Analysis
Examples
Questions
Applications in Engineering
Matlab
Maple
Introduction
In this topic, we see that under certain circumstances, we may factor
a matrix M in the form M = LL^{T}. Such
a decomposition is called a Cholesky decomposition. It requires half
the memory, and half the number operations of an PLU decomposition, but
it may only be applied in restricted circumstances, namely when the
matrix M is real, symmetric, and positive definite.
Background
References
Interactive Maplet
A Cholesky Decomposition Tutorial
Theory
Solving Systems of Linear Equations
Suppose that the n × n matrix M given by a system of linear equations
is both real and symmetric. The symmetry suggests that we can store the
matrix in half the memory required by a full nonsymmetric matrix of the same
size. Consequentially, one may suspect that it may also be possible to
write M = LL^{T} instead of
write M = LU, as we saw in the topic on the LU decomposition.
First, if the matrix M is 1 × 1, then M contains a single
entry. In this case, we can write
This immediately suggests that if we want to keep L real, then
m_{1, 1} ≥ 0.
Let us take a more general lower triangular matrix, for example, the 4 × 4 matrix:
Calculating LL^{T}, we get:
By inspection, we see that this matrix is symmetric. Thus, if we wanted to
write a general symmetric matrix M as LL^{T}, from the
first column, we get that:
 l_{1, 1} = sqrt( m_{1, 1} )
 l_{2, 1} = m_{2, 1}/l_{1, 1}
 l_{3, 1} = m_{3, 1}/l_{1, 1}
 l_{4, 1} = m_{4, 1}/l_{1, 1}
Having calculated these values from the entries of the matrix M, we may
go to the second column, and we note that, because we have already solved for
the entries of the form l_{i, 1}, we may continue to solve:
 l_{2, 2} = sqrt( m_{2, 2}  l_{2, 1}^{2})
 l_{3, 2} = (m_{3, 2}  l_{2, 1}l_{3, 1})/l_{2, 2}
 l_{4, 2} = (m_{4, 2}  l_{2, 1}l_{4, 1})/l_{2, 2}
Having solved these three, we find that we can solve for
l_{3, 3} and l_{4, 3}:
 l_{3, 3} = sqrt( m_{3, 3}  l_{3, 1}^{2}  l_{3, 2}^{2})
 l_{4, 3} = (m_{4, 3}  l_{3, 1}l_{4, 1}  l_{3, 2}l_{4, 2})/l_{3, 3}
and thus, finally we may solve for the last entry:
 l_{4, 4} = sqrt( m_{4, 4}  l_{4, 1}^{2}  l_{4, 2}^{2}  l_{4, 3}^{2})
This may seem exceptionally complex, but by using dot products, we can
simplify this algorithm significantly, as is covered in the howto.
Three Observations
How can we ensure that all of the square roots are positive? Without proof,
we will state that the Cholesky decomposition is real if the matrix M is
positive definite.
Solving a problem Mx = b where M is real and positive
definite may be reduced to finding the Cholesky decomposition and then
setting y = L^{T}x, solving Ly = b and
then solving L^{T}y = b.
We have not discussed pivoting. If the matrix is diagonally dominant, then
pivoting is not required for the PLU decomposition, and consequentially, not required
for Cholesky decomposition, either.
HOWTO
Problem
Given the n × n real, symmetric, and diagonally dominant matrix M, find a decomposition
M = LL^{T} where
L is a real
lowertriangular matrix. Next, given
a column vector b, solve the system of linear equations defined
by M x = b for the vector of unknowns x.
Assumptions
We will assume that M is real, symmetric, and diagonally dominant, and
consequently, it must be invertible.
Cholesky Decomposition
Start with the candidate matrix L = 0_{n, n}
where 0_{m, n} is the
m × n zero matrix.
Then, for i = 1, 2, ..., n do:
Subtract from m_{i, i} the dot
product of the ith row of L (as constructed so far) with itself and
set l_{i, i} to be the square root of
this result.
For j = i + 1, ..., n, subtract
the dot product of the ith and jth rows of L (as constructed so far) and
set l_{j, i} to be the this result
divided by l_{i, i}.

This is so simple to program in Matlab that we should cover it here:
M = [5 1.2 0.3 0.6; 1.2 6 0.4 0.9; 0.3 0.4 8 1.7; 0.6 0.9 1.7 10];
n = length( M );
L = zeros( n, n );
for i=1:n
L(i, i) = sqrt( M(i, i)  L(i, :)*L(i, :)' );
for j=(i + 1):n
L(j, i) = ( M(j, i)  L(i, :)*L(j, :)' )/L(i, i);
end
end
to get the result:
>> format long
>> L
L =
2.236067977499790 0.000000000000000 0.000000000000000 0.000000000000000
0.536656314599949 2.389979079406345 0.000000000000000 0.000000000000000
0.134164078649987 0.197491268466351 2.818332343581848 0.000000000000000
0.268328157299975 0.436823907370487 0.646577012719190 3.052723872310221
Solving with Forward and Backward Substitution
Calculate use forward substitution to solve
for y in Ly = b.
Then, having found y, use backward substitution to solve for
x in L^{T}x = y.
Error Analysis
The error analysis for the Cholesky decomposition is similar to that for
the PLU decomposition, which we will look at
when we look at matrix and vector norms.
Examples
Example 1
Find the Cholesky decomposition of the matrix M = (m_{i, j}):
To begin, we note that M is real, symmetric, and diagonally dominant, and therefore positive definite, and thus a real Cholesky decomposition
exists.
To do this, we first set up the matrix L = 0_{3} (the 3 × 3 zero
matrix):
For the first column, we may make the following simplification:
the first entry l_{1,1} is the square root of m_{1,1}:
Next, for all other entries in the first column,
l_{i, 1} = m_{i, 1}/l_{1, 1}:
Next, we go to the 2nd column: we subtract the dot product of the 2nd row
of L with itself from the entry m_{2, 2} and set l_{2, 2}
to the square root of this result:
For the 3rd row of the 2nd column, we subtract the dot product
of the 2nd and 3rd rows of L from m_{3,2} and
set l_{3,2} to this result divided by l_{2, 2}.
Finally, to complete our Cholesky decomposition, we subtract the dot product of the 3rd row
of L with itself from the entry m_{3, 3} and set l_{3, 3}
to the square root of this result:
This is the Cholesky decomposition of M, and a quick test shows
that L⋅L^{T} = M.
Example 2
Use the Cholesky decomposition from Example 1 to solve Mx = b
for x when b = (55, 19, 114)^{T}.
We rewrite Mx = b as LL^{T}x = b and
let L^{T}x = y.
First we solve Ly = b using forward substitution to get y = (11, 2, 14)^{T}.
Next, using this, we solve L^{T}x = y using backward substitution to get
x = (1, 2, 2)^{T}.
Example 3
Find the Cholesky decomposition of the matrix M = (m_{i, j}):
To begin, we note that M is real, symmetric, and diagonally dominant, and therefore positive definite, and thus a real Cholesky decomposition
exists.
Set L to be the 4 × 4 zero matrix. Starting with the first column, set l_{1,1} to be
the square root of m_{1, 1} = 9, and then set
l_{i, 1} = m_{i, 1}/l_{1,1}:
Next, for the 2nd column, we subtract off the dot product of the 2nd row of L
with itself from m_{2, 2} and set l_{2, 2} to be the square
root of this result:
For the entry l_{3, 2}, we subtract off the dot product of
rows 3 and 2 of L from m_{3, 2} and divide this by l_{2,2}.
Similarly, for the entry l_{4, 2}, we subtract off the dot product of
rows 4 and 2 of L from m_{4, 2} and divide this by l_{2,2}:
Next, for the 3rd column, we subtract off the dot product of the 3rd row of L
with itself from m_{3, 3} and set l_{3, 3} to be the square
root of this result:
For the entry l_{4, 3}, we subtract off the dot product of
rows 4 and 3 of L from m_{4, 3} and divide this by l_{3,3}:
Finally, for the 4th column, we subtract off the dot product of the 4th row of L
with itself from m_{4, 4} and set l_{4, 4} to be the square
root of this result:
In Matlab, we see that LL^{T} = M:
>> L = [3 0 0 0; 0.2 4 0 0; 0.1 0.3 2 0; 0.5 0.4 0.2 5]
L =
3.00000 0.00000 0.00000 0.00000
0.20000 4.00000 0.00000 0.00000
0.10000 0.30000 2.00000 0.00000
0.50000 0.40000 0.20000 5.00000
>> L * L'
ans =
9.0000 0.6000 0.3000 1.5000
0.6000 16.0400 1.1800 1.5000
0.3000 1.1800 4.1000 0.5700
1.5000 1.5000 0.5700 25.4500
Example 4
Use the Cholesky decomposition from Example 3 to solve Mx = b
for x when b = (2.49, 0.566, 0.787, 2.209)^{T}.
We rewrite Mx = b as LL^{T}x = b and
let L^{T}x = y.
First we solve Ly = b using forward substitution to get y = (0.83, 0.1, 0.42, 0.5)^{T}.
Next, using this, we solve L^{T}x = y using backward substitution to get
x = (0.3, 0, 0.2, 0.1)^{T}.
Questions
Question 1
Find the Cholesky decomposition of the matrix M:
Answer:
L = [7 0 0 0; 2 5 0 0; 1 2 6 0; 1 0 3 5]
Question 2
Using the Cholesky decomposition of M from Question 1, use
backward and forward substitution
to solve Mx = b
where b = (14, 46, 36, 32)^{T}.
Answer: x = (1, 2, 0, 1)^{T}.
Question 3
Find the Cholesky decomposition of the matrix M:
Answer:
L = [2 0 0 0;0.2 1 0 0;0.4 0.2 3 0;0.1 0.3 0.5 2]
Question 4
Using the Cholesky decomposition of M from Question 3, use
backward and forward substitution
to solve Mx = b
where b = (0.2, 0.32, 13.52, 14.17)^{T}.
Answer: x = (0, 1, 1 3)^{T}.
Applications to Engineering
The conductance matrix formed by a circuit is positive definite, as
are the matrices required to solve a leastsquares linear regression.
Matlab
The follow Matlab code finds the Cholesky decomposition of the
matrix M:
M = [0.9 0.06 0.39 0.24; 0.06 1.604 0.134 0.464; 0.39 0.134 2.685 0.802; 0.24 0.464 0.802 1.977];
L = chol( M )' % chol returns an upper triangular matrix
To do this manually, consider:
M = [0.9 0.06 0.39 0.24; 0.06 1.604 0.134 0.464; 0.39 0.134 2.685 0.802; 0.24 0.464 0.802 1.977]
n = length( M );
L = zeros( n, n );
for i=1:n
L(i, i) = sqrt(M(i, i)  L(i, :)*L(i, :)')
for j=(i + 1):n
L(j, i) = (M(j, i)  L(i,:)*L(j ,:)')/L(i, i)
end
end
To solve the system M x = b where
b = [0.063, 0.6358, 0.5937, 0.1907]'
let Matlab use forward and backward substitution as follows:
y = L \ b
x = L' \ y
Maple
The following commands in Maple finds the Cholesky decomposition
of the given matrix M:
M := Matrix([[9, 0.6, 0.3, 1.5], [0.6, 16.04, 1.1800, 1.5], [0.3, 1.18, 4.1, 0.57], [1.5, 1.5, 0.57, 25.45]] );
LinearAlgebra:LUDecomposition( M, method = Cholesky );
For more help on either of these routines or on the LinearAlgebra
package, enter:
?LinearAlgebra
?LinearAlgebra,LUDecomposition
?LinearAlgebra,LinearSolve