3.7.2 The Cholesky factorisation

A square root of a matrix

Every positive number $a$ has a square root: $a = \sqrt{a}\cdot\sqrt{a}$ . A positive definite matrix has a matrix version of this, built from a lower triangular matrix: one with only zeros above the diagonal.

Take

\begin{gathered} A = \begin{pmatrix} 4 & 2 \\ 2 & 5 \end{pmatrix}, \quad L = \begin{pmatrix} 2 & 0 \\ 1 & 2 \end{pmatrix} \\[6pt] L^\top = \begin{pmatrix} 2 & 1 \\ 0 & 2 \end{pmatrix} \end{gathered}

Multiply $L$ by its transpose, row by column (lesson 3.2.1):

entry of $LL^\top$	row of $L$ · column of $L^\top$	value
$(1, 1)$	$(2, 0) \cdot (2, 0)$	4
$(1, 2)$ and $(2, 1)$	$(2, 0) \cdot (1, 2)$	2
$(2, 2)$	$(1, 2) \cdot (1, 2)$	5

So $LL^\top = A$ . The two off-diagonal entries are equal because $LL^\top$ is always symmetric: $(LL^\top)^\top = LL^\top$ by the rule of lesson 3.2.3.

Rule: The Cholesky factorisation

\begin{gathered} A = LL^\top \\[4pt] L \text{ lower triangular, } l_{ii} > 0 \end{gathered}

A symmetric matrix $A$ has such an $L$ exactly when it is positive definite, and then $L$ is unique. $L$ is the Cholesky factor of $A$ , named after André-Louis Cholesky.

It is the completed square. Lesson 3.7.1 wrote a quadratic form as a sum of squares. With $A = LL^\top$ ,

\begin{aligned} \mathbf{x}^\top A\,\mathbf{x} &= (L^\top\mathbf{x})^\top(L^\top\mathbf{x}) \\ &= \lVert L^\top\mathbf{x} \rVert^2 \end{aligned}

Here $L^\top\mathbf{x} = (2x_1 + x_2,\ 2x_2)$ , so $4x_1^2 + 4x_1x_2 + 5x_2^2 = (2x_1 + x_2)^2 + (2x_2)^2$ . That is lesson 3.7.1's $a\big(x_1 + \frac{b}{a}x_2\big)^2 + \frac{ac - b^2}{a}x_2^2 = 4\big(x_1 + \frac{1}{2}x_2\big)^2 + 4x_2^2$ , written with a matrix. At $(1, -1)$ both sides give 5.

The order matters: $L^\top L = \begin{pmatrix} 5 & 2 \\ 2 & 4 \end{pmatrix}$ , which is not $A$ .

LL^\top

import numpy as np A = np.array([[4, 2], [2, 5]]) L = np.array([[2, 0], [1, 2]]) print(L @ L.T) print(L.T @ L) x = np.array([1, -1]) print(x @ A @ x, L.T @ x, (L.T @ x) @ (L.T @ x))

⌘/Ctrl + Enter runs the cell

The 2 × 2 formulas

Write $L = \begin{pmatrix} l_{11} & 0 \\ l_{21} & l_{22} \end{pmatrix}$ and multiply out:

LL^\top = \begin{pmatrix} l_{11}^2 & l_{11}l_{21} \\ l_{11}l_{21} & l_{21}^2 + l_{22}^2 \end{pmatrix}

Match this with $A = \begin{pmatrix} a & b \\ b & c \end{pmatrix}$ , one entry at a time: $l_{11}^2 = a$ first, then $l_{11}l_{21} = b$ , then $l_{21}^2 + l_{22}^2 = c$ .

Rule: Cholesky for 2 × 2

\begin{gathered} l_{11} = \sqrt{a}, \quad l_{21} = \frac{b}{l_{11}} \\[4pt] l_{22} = \sqrt{c - l_{21}^2} \end{gathered}

Each step uses only what came before. The number under the last square root is $c - \frac{b^2}{a} = \frac{ac - b^2}{a}$ , so both roots are real and positive exactly when $a > 0$ and $ac - b^2 > 0$ : lesson 3.7.1's leading-minor test.

Three examples:

$A$	$l_{11}$	$l_{21}$	$l_{22}$
$\begin{pmatrix} 4 & 2 \\ 2 & 5 \end{pmatrix}$	$\sqrt{4} = 2$	$2 \div 2 = 1$	$\sqrt{5 - 1} = 2$
$\begin{pmatrix} 1 & 0.6 \\ 0.6 & 1 \end{pmatrix}$	1	0.6	$\sqrt{1 - 0.36} = 0.8$
$\begin{pmatrix} 3 & 1 \\ 1 & 2 \end{pmatrix}$	$\sqrt{3} \approx 1.7321$	$\frac{1}{\sqrt{3}} \approx 0.5774$	$\sqrt{5/3} \approx 1.291$

The second row has 1 on its diagonal and $\rho = 0.6$ off it. For any such matrix, $L = \begin{pmatrix} 1 & 0 \\ \rho & \sqrt{1 - \rho^2} \end{pmatrix}$ .

When it fails. For lesson 3.7.1's indefinite $N = \begin{pmatrix} 1 & 2 \\ 2 & 1 \end{pmatrix}$ , $l_{11} = 1$ and $l_{21} = 2$ , and $l_{22}$ would be $\sqrt{1 - 4} = \sqrt{-3}$ : no real number. For the semidefinite $J = \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}$ , $l_{22} = \sqrt{1 - 1} = 0$ , which is not positive, and $L$ has no inverse.

L

import math import numpy as np def chol2(a, b, c): l11 = math.sqrt(a) l21 = b / l11 l22 = math.sqrt(c - l21**2) # fails if c - l21**2 < 0 return np.array([[l11, 0], [l21, l22]]) for a, b, c in [(4, 2, 5), (1, 0.6, 1), (3, 1, 2)]: L = chol2(a, b, c) print(np.round(L, 4), np.allclose(L @ L.T, [[a, b], [b, c]])) print(1 - (2 / math.sqrt(1))**2) # N: what l22 would be the square root of

⌘/Ctrl + Enter runs the cell

A 3 × 3 example by hand

For larger matrices, match $LL^\top = A$ one column of $L$ at a time, from left to right. Each diagonal entry is a square root, and each entry below it a division.

Rule: Cholesky, column j

\begin{gathered} l_{jj} = \sqrt{a_{jj} - \textstyle\sum_{k<j} l_{jk}^2} \\[6pt] l_{ij} = \dfrac{a_{ij} - \sum_{k<j} l_{ik}\,l_{jk}}{l_{jj}} \\[4pt] \text{for each } i > j \end{gathered}

The sums run over the columns already found, to the left of column $j$ ; for $j = 1$ they are empty. If a number under a square root is 0 or negative, $A$ is not positive definite and the method stops.

Factor

A = \begin{pmatrix} 4 & 2 & -2 \\ 2 & 10 & 2 \\ -2 & 2 & 6 \end{pmatrix}

step	entry	working	value
1	$l_{11}$	$\sqrt{4}$	2
2	$l_{21}$	$2 \div 2$	1
3	$l_{31}$	$-2 \div 2$	$-1$
4	$l_{22}$	$\sqrt{10 - 1^2}$	3
5	$l_{32}$	$(2 - (-1)(1)) \div 3$	1
6	$l_{33}$	$\sqrt{6 - (-1)^2 - 1^2}$	2

L = \begin{pmatrix} 2 & 0 & 0 \\ 1 & 3 & 0 \\ -1 & 1 & 2 \end{pmatrix}

Check one entry: row 3 of $L$ dotted with row 2 of $L$ is $(-1)(1) + (1)(3) + (2)(0) = 2 = a_{32}$ .

The link with elimination. Eliminate on $A$ as in lesson 3.3.1, with no row swaps: $R_2 - \frac{1}{2}R_1$ , $R_3 + \frac{1}{2}R_1$ , then $R_3 - \frac{1}{3}R_2$ . The pivots are 4, 9 and 4, the squares of $L$ 's diagonal.

A symmetric matrix is positive definite exactly when elimination without swaps gives all positive pivots. Since $\det A$ is the product of the pivots (lesson 3.4.1), $\det A = (2 \cdot 3 \cdot 2)^2 = 144$ .

Cholesky does about half the arithmetic of elimination, because $A$ 's symmetry means only one triangle needs computing. For a positive definite matrix it also needs no row swaps.

j

import math import numpy as np def cholesky(A): A = np.array(A, dtype=float) n = len(A) L = np.zeros((n, n)) for j in range(n): L[j, j] = math.sqrt(A[j, j] - L[j, :j] @ L[j, :j]) for i in range(j + 1, n): L[i, j] = (A[i, j] - L[i, :j] @ L[j, :j]) / L[j, j] return L A = np.array([[4, 2, -2], [2, 10, 2], [-2, 2, 6]]) L = cholesky(A) print(L) print(np.allclose(L @ L.T, A), np.allclose(L, np.linalg.cholesky(A))) print(round(np.linalg.det(A), 6), np.prod(np.diag(L))**2)

⌘/Ctrl + Enter runs the cell

L

import numpy as np A = np.array([[4, 2, -2], [2, 10, 2], [-2, 2, 6]], dtype=float) U = A.copy() U[1] -= (U[1, 0] / U[0, 0]) * U[0] U[2] -= (U[2, 0] / U[0, 0]) * U[0] U[2] -= (U[2, 1] / U[1, 1]) * U[1] print(U) print(np.diag(U), np.sqrt(np.diag(U)))

⌘/Ctrl + Enter runs the cell

In NumPy: a test for positive definiteness

np.linalg.cholesky(A) returns $L$ . When $A$ is not positive definite it stops with LinAlgError: Matrix is not positive definite, as inv stopped on a singular matrix in lesson 3.2.3.

That makes it the quickest test on a computer: try the factorisation and catch the error. It costs less than finding the eigenvalues. It fails for the indefinite $N$ and for the semidefinite $J$ .

Rule: Cholesky as a test

\begin{gathered} \text{the factorisation succeeds} \\ \iff A \text{ is positive definite} \end{gathered}

The equivalence holds for symmetric $A$ , so check symmetry first. A failure says only "not positive definite": the eigenvalues tell whether $A$ is indefinite or semidefinite.

Two traps in the code:

NumPy does not check symmetry. It reads only the lower triangle and the diagonal. $\begin{pmatrix} 4 & 100 \\ 2 & 5 \end{pmatrix}$ returns the same $L$ as $A$ , with no warning. Test np.allclose(A, A.T) before factoring.
SciPy returns the other triangle. scipy.linalg.cholesky(A) returns the upper factor $L^\top$ unless you pass lower=True.

A semidefinite matrix sits on the edge. For $J$ the last step is exactly $\sqrt{0}$ , so the test fails. When rounding is involved, a zero can come out as a tiny positive number, and the factorisation may succeed with a useless, near-zero diagonal entry.

L

import numpy as np from scipy.linalg import cholesky as scipy_cholesky print(np.linalg.cholesky(np.array([[4, 2], [2, 5]]))) for M in ([[1, 2], [2, 1]], [[1, 1], [1, 1]]): try: np.linalg.cholesky(np.array(M)) except np.linalg.LinAlgError as err: print(M, "LinAlgError:", err) typo = np.array([[4, 100], [2, 5]]) print(np.allclose(typo, typo.T), np.linalg.cholesky(typo)) print(scipy_cholesky(np.array([[4, 2], [2, 5]]))) def is_positive_definite(M): M = np.asarray(M, dtype=float) if M.shape[0] != M.shape[1] or not np.allclose(M, M.T): return False try: np.linalg.cholesky(M) return True except np.linalg.LinAlgError: return False for M in ([[3, 1], [1, 2]], [[1, 2], [2, 1]], [[1, 1], [1, 1]], [[4, 100], [2, 5]], [[4, -1, 0], [-1, 3, 2], [0, 2, 5]]): print(is_positive_definite(M))

⌘/Ctrl + Enter runs the cell

Solving Ax = b with two triangular solves

With $A = LL^\top$ , the system $A\mathbf{x} = \mathbf{b}$ is $L(L^\top\mathbf{x}) = \mathbf{b}$ . Call the bracket $\mathbf{y}$ and solve in two steps, each with a triangular matrix.

Rule: Solving with Cholesky

\begin{gathered} L\mathbf{y} = \mathbf{b} \quad \text{(forward)} \\[4pt] L^\top\mathbf{x} = \mathbf{y} \quad \text{(back)} \end{gathered}

Forward substitution solves $L\mathbf{y} = \mathbf{b}$ from the top: row 1 of $L$ has one unknown, and each row below adds one more. Back substitution (lesson 3.3.1) then solves $L^\top\mathbf{x} = \mathbf{y}$ from the bottom.

Solve $A\mathbf{x} = (-2, -4, 8)$ with the $3 \times 3$ example. Forward, with the rows of $L$ , $(2, 0, 0)$ , $(1, 3, 0)$ and $(-1, 1, 2)$ :

\begin{aligned} y_1 &= -2 \div 2 = -1 \\ y_2 &= (-4 - y_1) \div 3 = -1 \\ y_3 &= (8 + y_1 - y_2) \div 2 = 4 \end{aligned}

Back, with $L^\top = \begin{pmatrix} 2 & 1 & -1 \\ 0 & 3 & 1 \\ 0 & 0 & 2 \end{pmatrix}$ :

\begin{aligned} x_3 &= 4 \div 2 = 2 \\ x_2 &= (-1 - x_3) \div 3 = -1 \\ x_1 &= (-1 - x_2 + x_3) \div 2 = 1 \end{aligned}

So $\mathbf{x} = (1, -1, 2)$ . Check it in row 3 of $A$ : $-2(1) + 2(-1) + 6(2) = 8$ .

Factoring is the expensive part. Once $L$ is known, each new right-hand side $\mathbf{b}$ costs only the two triangular solves, far less than a new elimination.

Try it solve_triangular from SciPy does each step; lower=True says the matrix is lower triangular. The answer agrees with np.linalg.solve . The last lines are the forward substitution written out, as in the text.

import numpy as np from scipy.linalg import solve_triangular A = np.array([[4, 2, -2], [2, 10, 2], [-2, 2, 6]]) L = np.linalg.cholesky(A) b = np.array([-2, -4, 8]) y = solve_triangular(L, b, lower=True) # forward: L y = b x = solve_triangular(L.T, y, lower=False) # back: L^T x = y print(y, x) print(np.allclose(A @ x, b), np.allclose(x, np.linalg.solve(A, b))) y1 = b[0] / L[0, 0] y2 = (b[1] - L[1, 0] * y1) / L[1, 1] y3 = (b[2] - L[2, 0] * y1 - L[2, 1] * y2) / L[2, 2] print(y1, y2, y3)

⌘/Ctrl + Enter runs the cell

L as a square root

$A = LL^\top$ is the matrix version of $a = \sqrt{a}\cdot\sqrt{a}$ . For a $1 \times 1$ matrix $(a)$ , $L = (\sqrt{a})$ . For a diagonal matrix, $L$ takes the square root of each diagonal entry: $\begin{pmatrix} 0.04 & 0 \\ 0 & 0.09 \end{pmatrix}$ gives $L = \begin{pmatrix} 0.2 & 0 \\ 0 & 0.3 \end{pmatrix}$ .

And $\det A = (\det L)^2$ , with $\det L$ the product of $L$ 's diagonal.

As a map, $L$ sends the unit circle to an ellipse, as lesson 3.6.3 found for symmetric matrices. Its first column $L\mathbf{e}_1 = (2, 1)$ and its second $L\mathbf{e}_2 = (0, 2)$ both land on that ellipse. $L\mathbf{e}_2$ points straight up, because the first row of $L$ is $(2, 0)$ .

$L$ is not the only square root. Lesson 3.6.3 writes $A = Q\Lambda Q^\top$ . Taking $\Lambda^{1/2}$ , the diagonal matrix of $\sqrt{\lambda_i}$ , gives the symmetric square root $B = Q\Lambda^{1/2}Q^\top$ , with $BB^\top = A$ too. For $A$ above, $B \approx \begin{pmatrix} 1.9403 & 0.4851 \\ 0.4851 & 2.1828 \end{pmatrix}$ .

$L$ and $B$ send the unit circle to the same ellipse, because $LL^\top = BB^\top = A$ , though each point lands somewhere different. Its semi-axes are $\sqrt{\lambda}$ for the eigenvalues of $A$ : about 2.56 and 1.56.

Any $C$ with $CC^\top = A$ gives that ellipse. Cholesky picks the one that is lower triangular with a positive diagonal, which is the cheapest to compute and to solve with. Case study 3.7.4 uses $L$ to turn independent random numbers into correlated ones, with a covariance matrix as $A$ .

B

import numpy as np A = np.array([[4, 2], [2, 5]]) L = np.linalg.cholesky(A) lam, Q = np.linalg.eigh(A) B = Q @ np.diag(np.sqrt(lam)) @ Q.T # the symmetric square root print(np.round(B, 4)) print(np.allclose(B @ B.T, A), np.allclose(L @ L.T, A)) print(np.round(lam, 4), np.round(np.sqrt(lam), 4)) print(np.linalg.cholesky(np.diag([0.04, 0.09]))) print(round(np.linalg.det(L), 6), round(np.linalg.det(A), 6))

⌘/Ctrl + Enter runs the cell

Common mistakes

$LL^\top$ , not $L^\top L$ . For $L = \begin{pmatrix} 2 & 0 \\ 1 & 2 \end{pmatrix}$ , $L^\top L = \begin{pmatrix} 5 & 2 \\ 2 & 4 \end{pmatrix}$ , not $A$ .
Take the square root on the diagonal. $l_{11} = \sqrt{a_{11}}$ , not $a_{11}$ .
Below the diagonal, divide by the diagonal entry of the same column, $l_{jj}$ , after subtracting the products already found.
Check symmetry first. np.linalg.cholesky reads only the lower triangle and never complains about the upper one.
Know which triangle you get. NumPy returns $L$ ; scipy.linalg.cholesky returns $L^\top$ unless lower=True.
Forward first, then back. Solve $L\mathbf{y} = \mathbf{b}$ from the top, then $L^\top\mathbf{x} = \mathbf{y}$ from the bottom.
A failure means "not positive definite", nothing more. The matrix may be indefinite, or semidefinite with a zero eigenvalue; ask eigvalsh which.

The Cholesky factorisation

0 of 12 correct

1.
Find $l_{11}$ in the Cholesky factor of $\begin{pmatrix} 9 & 3 \\ 3 & 5 \end{pmatrix}$
2.
For the same matrix, find $l_{21}$
3.
For the same matrix, find $l_{22}$
4.
$L = \begin{pmatrix} 2 & 0 \\ 3 & 1 \end{pmatrix}$ . Find the entry in row 2, column 2 of $LL^\top$
5.
Find $l_{22}$ for $\begin{pmatrix} 1 & 0.8 \\ 0.8 & 1 \end{pmatrix}$
6.
$\begin{pmatrix} 4 & 6 \\ 6 & c \end{pmatrix}$ is positive definite only when $c$ is above which number?
7.
Find $l_{32}$ for $\begin{pmatrix} 1 & 2 & 1 \\ 2 & 8 & 4 \\ 1 & 4 & 6 \end{pmatrix}$
8.
For that matrix, find $l_{33}$
9.
A Cholesky factor has diagonal 1, 2 and 2. Find $\det A$
10.
Solve $L\mathbf{y} = (4, 11)$ forward for $L = \begin{pmatrix} 2 & 0 \\ 1 & 3 \end{pmatrix}$ . What is $y_2$ ?
11.
Then solve $L^\top\mathbf{x} = (2, 3)$ back, with $L^\top = \begin{pmatrix} 2 & 1 \\ 0 & 3 \end{pmatrix}$ . What is $x_1$ ?
12.
For $L = \begin{pmatrix} 2 & 0 \\ 1 & 2 \end{pmatrix}$ , find the entry in row 1, column 1 of $L^\top L$ , the wrong order

How to type answers

Powers with ^: x^7, 2^(1/3), a^-2.
Multiply with * or by writing side by side: 3x^2, 2(x+1).
ln(x) is the natural log, log(x, 2) is log base 2, e is Euler's number, sqrt(x) is a square root.
sin(x), cos(x) and tan(x) take radians; write pi for π and sin(x)^2 for sin²x.
Your first answer to each question is saved: it counts towards your score and brings the question back for review.
Exact answers like ln(20)/ln(3) are best. Decimals need about 4 significant figures.

A

The Cholesky factorisation

A square root of a matrix

Rule: The Cholesky factorisation

The 2 × 2 formulas

Rule: Cholesky for 2 × 2

A 3 × 3 example by hand

Rule: Cholesky, column j

In NumPy: a test for positive definiteness

Rule: Cholesky as a test

Solving Ax = b with two triangular solves

Rule: Solving with Cholesky

L as a square root

Common mistakes

The Cholesky factorisation

Summary card

Free further reading