# Cholesky decomposition

$$\gdef\delim#1#2#3{\mathopen{}\mathclose{\left#1 #2 \right#3}} \gdef\p#1{\delim({#1})} \gdef\ps#1{\delim\{{#1}\}} \gdef\N{\mathcal N} \gdef\vec#1{\bm #1} \gdef\mat#1{\mathrm #1} \gdef\k{\mathrm k} \gdef\u{\mathrm u} \gdef\c{\mathrm c} \gdef\T{\mathrm T}$$

Positive definite matrices occur as covariance matrices. Given a positive definite matrix $\mat M$ it can be factored into a triangular matrix $\mat L$ such that

$$\mat M = \mat L ⋅ \mat L^\T$$

## Computing the Cholesky decompostion

https://en.wikipedia.org/wiki/Cholesky_decomposition#The_Cholesky_algorithm

https://en.wikipedia.org/wiki/Cholesky_decomposition#The_Cholesky%E2%80%93Banachiewicz_and_Cholesky%E2%80%93Crout_algorithms

A reasonably efficient way to compute it is:

def cholesky(M):
L = np.tril(M)
for i in range(L.shape[0]):
L[i:,i] -= L[i,:i] @ L[i:,:i].T
L[i:,i] /= np.sqrt(L[i,i])
return L


But it's more efficient to use the LAPACK dpotrf function, available as np.linalg.cholesky.

https://netlib.org/lapack/explore-html/d2/d09/group__potrf_ga84e90859b02139934b166e579dd211d4.html#ga84e90859b02139934b166e579dd211d4

def cholesky_update(L, x):
for k in range(L.shape[0]):
r = np.sqrt(L[k,k]**2 + x[k]**2)
c = r / L[k,k]
s = x[k] / L[k,k]
L[k,k] = r
L[k+1:, k] = (L[k+1:, k] + s * x[k+1:]) /c
x[k+1:] = c * x[k+1:] - s * L[k+1:, k]

def cholesky_subset(L, mask):
L = L.copy()
for i in range(L.shape[0]):
if not mask[i]:
cholesky_update(L[i+1:,i+1:], L[i+1:, i])
return L[np.ix_(mask, mask)]


## References

https://algowiki-project.org/en/Cholesky_decomposition

https://www.cs.utexas.edu/users/flame/Notes/NotesOnCholReal.pdf

https://www.cs.utexas.edu/users/flame/pubs/flawn41.pdf

https://christian-igel.github.io/paper/AMERCMAUfES.pdf

https://github.com/scipy/scipy/issues/8188

https://link.springer.com/article/10.1007/BF01933218

Remco Bloemen
Math & Engineering
https://2π.com