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)]

Subset of Cholesky

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