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 diagonal matrix $\mat L$ such that

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

Computing the Cholesky decompostion

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.

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


Remco Bloemen
Math & Engineering