\(\newcommand{\re}{\mathbb{R}}\)

Cholesky factorization

Given a spd \(A \in \re^{m \times m}\), we want to write it as

\[ A = R^\top R \]

where \(R\) is upper triangular matrix.

%config InlineBackend.figure_format = 'svg'
import numpy as np
import matplotlib.pyplot as plt

The function returns the factor \(R\) which is an upper-triangular matrix.

# Algorithm 23.1
def cholesky(A):
    m = A.shape[0]
    R = np.triu(A)
    for k in range(m):
        for j in range(k+1,m):
            R[j,j:m] = R[j,j:m] - R[k,j:m] * (R[k,j] / R[k,k])
        R[k,k:m] = R[k,k:m] / np.sqrt(R[k,k])
    return R

Generate a random spd matrix as

\[ A = B^\top B \]

where \(B\) is a random \(m \times m\) matrix and hence very likely to be of full rank, so that \(A\) is positive definite.

m = 4
B = 2 * np.random.rand(m,m) - 1
A = B.T @ B

and compute its Cholesky decomposition.

R = cholesky(A)
print('max|A - R.T * R| = ', np.abs(A - R.T@R).max())

max|A - R.T * R| =  4.440892098500626e-16