6.3. Cholesky decomposition

Cholesky decomposition is an efficient matrix decomposition method that can be used when a square matrix is positive definite.

Definition 6.2 (Positive definite matrix)

A square matrix \(A\) is said to be positive definite if

\[ \mathbf{x}^\mathsf{T} A \mathbf{x} > 0, \]

for all \(\mathbf{x} \in \mathbb{R}^n \backslash \{\mathbf{0}\}\).

Theorem 6.2 (Sylvestor’s criterion)

A square matrix \(A\) is positive definite if it is symmetric and the determinants of all leading principal submatrices are positive, e.g.,

\[\begin{split} (a_{11}), \qquad \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix}, \qquad \ldots \end{split}\]

Example 6.5

Use the determinant test to show that the following matrix is positive definite

\[\begin{split} \begin{align*} A = \begin{pmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{pmatrix}. \end{align*} \end{split}\]

---

Solution

Matrix \(A\) is symmetric since \(A=A^\mathsf{T}\). Checking the determinants of the upper-left sub-matrices

\[\begin{split} \begin{align*} \det(2) &= 2 > 0, \\ \det\begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix} &= 4 - 1 = 3 > 0, \\ \det\begin{pmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{pmatrix} &= 6 - 2 + 0 = 4 > 0 . \end{align*} \end{split}\]

Since \(A\) is a symmetric matrix and all the determinants of the upper-left sub-matrices are positive then \(A\) is a positive definite matrix.

Given a positive definite matrix \(A\) the Cholesky decomposition decomposes \(A\) into the product of a lower triangular matrix \(L\) and its transpose \(L^\mathsf{T}\), i.e.,

\[ \begin{align*} A = LL^\mathsf{T}. \end{align*} \]

Since \(LL^\mathsf{T}\) is completely determined by LLL, only one triangular matrix needs to be stored and computed. Consequently, Cholesky decomposition requires roughly half the computational work and storage of a general LU decomposition.

Consider the Cholesky decomposition of a \(3\times 3\) matrix

\[\begin{split} \begin{align*} \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{pmatrix} &= \begin{pmatrix} \ell_{11} & 0 & 0 \\ \ell_{21} & \ell_{22} & 0 \\ \ell_{31} & \ell_{32} & \ell_{33} \end{pmatrix} \begin{pmatrix} \ell_{11} & \ell_{21} & \ell_{31} \\ 0 & \ell_{22} & \ell_{32} \\ 0 & 0 & \ell_{33} \end{pmatrix}\\ &= \begin{pmatrix} \ell_{11}^2 & \ell_{11} \ell_{21} & \ell_{11} \ell_{31} \\ \ell_{11} \ell_{21} & \ell_{21}^2 +\ell_{22}^2 & \ell_{21} \ell_{31} +\ell_{22} \ell_{33} \\ \ell_{11} \ell_{31} & \ell_{21} \ell_{31} +\ell_{22} \ell_{33} & \ell_{31}^2 +\ell_{32}^2 +\ell_{33}^2 \end{pmatrix}. \end{align*} \end{split}\]

The elements on the main diagonal are

\[ \begin{align*} a_{jj} =\ell_{jj}^2 +\sum_{k=1}^{j-1} \ell_{jk}^2 , \end{align*} \]

and the other elements are

\[ \begin{align*} a_{ij} =\sum_{k=1}^j \ell_{ik} \ell_{jk} = \ell_{jj} \ell_{ij} +\sum_{k=1}^{j-1} \ell_{ik} \ell_{jk}. \end{align*} \]

Rearranging these expressions gives the following definition

Theorem 6.3 (Cholesky decomposition)

The Cholesky decomposition of an \(n \times n\) positive definite matrix \(A\) results in an \(n \times n\) matrix \(L\) such that \(A = LL^\mathsf{T}\). The elements of \(L\) are calculated using

(6.3)\[\begin{split} \begin{align} \ell_{jj} &= \sqrt{a_{jj} - \sum_{k=1}^{j-1} \ell_{jk}^2 }, & j &= 1, \ldots, n,\\ \ell_{ij} &= \dfrac{1}{\ell_{jj} }\left(a_{ij} -\displaystyle \sum_{k=1}^{j-1} \ell_{ik} \ell_{jk} \right), & i &= j + 1,\ldots ,n. \end{align} \end{split}\]

For an \(n\times n\) matrix, Cholesky decomposition requires approximately \(\frac13 n^3\) operations, compared with approximately \(\frac23 n^3\) operations for LU decomposition. For symmetric positive definite matrices, Cholesky decomposition is therefore usually preferred.

Algorithm 6.4 (Cholesky decomposition)

Inputs: An \(n \times n\) positive definite matrix \(A\).

Outputs: An \(n \times n\) lower triangular matrix \(L\) such that \(A = LL^\mathsf{T}\).

• \(L \gets \mathbf{0}_{n \times n}\)

• For \(j = 1, \ldots n\) do

  • \(\ell_{jj} \gets \sqrt{a_{jj} - \displaystyle\sum_{k=1}^{j-1} \ell_{jk}^2}\)

  • \(\ell_{ij} \gets \dfrac{1}{\ell_{jj} }\left(a_{ij} -\displaystyle \sum_{k=1}^{j-1} \ell_{ik} \ell_{jk} \right), \qquad i = j + 1, \ldots, n\)

• Return \(L\)

Example 6.6

Calculate the Cholesky decomposition of the following matrix

\[\begin{split} \begin{align*} A = \begin{pmatrix} 4 & -2 & -4\\ -2 & 10 & 5\\ -4 & 5 & 14 \end{pmatrix}. \end{align*} \end{split}\]

---

Solution

Check that \(A\) is positive definite

\[\begin{split} \begin{align*} \det(4) &= 4, \\ \det \begin{pmatrix} 4 & -2 \\ -2 & 10 \end{pmatrix} &= 44, \\ \det \begin{pmatrix} 4 & -2 & -4 \\ -2 & 10 & 5 \\ -4 & 5 & 14 \end{pmatrix} &= 4 \det \begin{pmatrix} 10 & 5 \\ 5 & 14 \end{pmatrix} + 2\det \begin{pmatrix} -2 & 5 \\ -4 & 14 \end{pmatrix} -4 \det \begin{pmatrix} -2 & 10 \\ -4 & 5 \end{pmatrix} \\ &= 4(115) + 2(-8) - 4(30) = 324. \end{align*} \end{split}\]

Since all leading principal minors are positive, \(A\) is a positive definite matrix.

Computing the Cholesky decomposition of \(A\)

\[\begin{split} \begin{align*} j &= 1: & \ell_{11} &= \sqrt{a_{11}} = \sqrt{4} = 2, \\ && \ell_{21} &= \frac{1}{\ell_{11}}(a_{11}) = \frac{1}{2}(-2) = -1, \\ && \ell_{31} &= \frac{1}{\ell_{11}}(a_{31}) = \frac{1}{2}(-4) = -2, \\ j &= 2: & \ell_{22} &= \sqrt{a_{22} - \ell_{21}^2} = \sqrt{10 - (-1)^2}= 3, \\ && \ell_{32} &= \frac{1}{\ell_{22}}(a_{32} - \ell_{31}\ell_{21}) = \frac{1}{3}(5 - (-2)(-1)) = 1, \\ j &= 3: & \ell_{33} &= \sqrt{a_{33} - \ell_{31}^2 - \ell_{32}^2} = \sqrt{14 - (-2)^2 - 1^2} = 3, \end{align*} \end{split}\]

therefore

\[\begin{split}L=\begin{pmatrix} 2 & 0 & 0 \\ -1 & 3 & 0 \\ -2 & 1 & 3 \end{pmatrix}. \end{split}\]

Checking that \(A = LL^\mathsf{T}\)

\[\begin{split} \begin{align*} \begin{pmatrix} 2 & 0 & 0\\ -1 & 3 & 0\\ -2 & 1 & 3 \end{pmatrix} \begin{pmatrix} 2 & -1 & -2\\ 0 & 3 & 1\\ 0 & 0 & 3 \end{pmatrix} = \begin{pmatrix} 4 & -2 & -4\\ -2 & 10 & 5\\ -4 & 5 & 14 \end{pmatrix}. \end{align*} \end{split}\]

---

6.3.1. Solving systems of linear equations using Cholesky decomposition

Given a system of linear equations of the form \(A \mathbf{x} = \mathbf{b}\) where Cholesky decomposition has been applied to the coefficient matrix then since \(A = LL^\mathsf{T}\) we have

\[ LL^\mathsf{T} \mathbf{x} = \mathbf{b}.\]

Similar to the approach of solving systems using LU decomposition, we introduce an intermediate vector \(\mathbf{y}\) such that

\[ L \mathbf{y} = \mathbf{b}, \]

which is computed using forward substitution. Once \(\mathbf{y}\) has been found, the solution vector \(\mathbf{x}\) is obtained from

\[ L^\mathsf{T} \mathbf{x} = \mathbf{y}, \]

using back substitution.

Algorithm 6.5 (Solving systems of linear equations using Cholesky decomposition)

Inputs: A system of linear equations expressed using the coefficient matrix \(A\) and variable vector \(\mathbf{b}\) such that \(A \mathbf{x} = \mathbf{b}\).

Outputs: The variable vector \(\mathbf{x}\).

• Calculate the Cholesky decomposition of \(A\)

• \(y_i \gets b_i - \displaystyle\sum_{j=1}^i \ell_{ij} y_j, \qquad i = 1, \ldots, n\)

• \(x_n \gets \dfrac{y_n}{\ell_{nn}}\)

• \(x_i \gets \dfrac{1}{\ell_{ii}} \left( y_i - \displaystyle \sum_{j=i+1}^{n} \ell_{ji} x_j \right), \qquad i = n - 1, \ldots, 1\)

• Return \(\mathbf{x}\)

Example 6.7

Solve the following system of linear equations using the Cholesky decomposition.

\[\begin{split} \begin{align*} \begin{pmatrix} 4 & -2 & -4\\ -2 & 10 & 5\\ -4 & 5 & 14 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} -2 \\ 49 \\ 27 \end{pmatrix}. \end{align*} \end{split}\]

---

Solution

We saw in Example 6.6 that the Cholesky decomposition of the matrix \(A\) is

\[\begin{split} L=\begin{pmatrix} 2 & 0 & 0 \\ -1 & 3 & 0 \\ -2 & 1 & 3 \end{pmatrix}. \end{split}\]

Solving \(L\mathbf{y}=\mathbf{b}\) using forward substitution

\[\begin{split} \begin{align*} \begin{pmatrix} 2 & 0 & 0\\ -1 & 3 & 0\\ -2 & 1 & 3 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} = \begin{pmatrix} -2 \\ 49 \\ 27 \end{pmatrix}, \end{align*} \end{split}\]

gives

\[\begin{split} \begin{align*} y_1 &=\frac{-2}{2}=-1,\\ y_2 &=\frac{1}{3}(49+y_1 )=\frac{1}{3}(49-1)=16,\\ y_3 &=\frac{1}{3}(27+2y_1 -y_2 )=\frac{1}{3}(27+2(-1)-16)=3. \end{align*} \end{split}\]

Solving \(L^\mathsf{T} \mathbf{x}=\mathbf{y}\) using back substitution

\[\begin{split} \begin{align*} \begin{pmatrix} 2 & -1 & -2\\ 0 & 3 & 1\\ 0 & 0 & 3 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} -1 \\ 16 \\ 3 \end{pmatrix}, \end{align*} \end{split}\]

gives

\[\begin{split} \begin{align*} x_3 &=\frac{3}{3}=1,\\ x_2 &=\frac{1}{3}(16-x_3 )=\frac{1}{3}(16-1)=5,\\ x_1 &=\frac{1}{2}(-1+x_2 +2x_3 )=\frac{1}{2}(-1+5+2(1))=3. \end{align*} \end{split}\]

So the solution is \(\mathbf{x}=(3,5,1)\).