Can someone suggest a way to get Cholesky factorization of a singular covariance matrix? I need it to match Cholesky on full-rank matrices, ie coordinate order should be preserved. My attempt below was to use ldl routine in scipy, but that gives me factorization on a different ordering, any ideas?

import numpy as np
from scipy import linalg
def modified_cholesky(arr):
    """Use Cholesky to produce LDL' factorization of arr."""
    chol = linalg.cholesky(arr, lower=True)
    d1 = np.diag(np.diag(chol))
    L = chol@linalg.inv(d1)
    return L, d1

def modified_ldl(arr):
    lu, d, perm=linalg.ldl(arr, lower=True)
    lu2 = lu[perm,:]
    return lu2, np.sqrt(d), perm

# sanity checks
arr=np.array([[ 3.,  5.], [ 5., 11.]])
mchol, d = modified_cholesky(arr)
np.testing.assert_allclose(mchol @ d @ d @ mchol.T, arr, rtol=1e-6, atol=1e-7)
mchol2, d2, perm = modified_ldl(arr)
np.testing.assert_allclose(mchol2 @ d2 @ d2 @ mchol2.T, arr[perm,:][:,perm], rtol=1e-6, atol=1e-7)

np.testing.assert_allclose(mchol, mchol2)  # fails because linalg.ldl is permuted
modified_cholesky(np.array([[1,1],[1,1]]))  # fails with 2-th leading minor of the array is not positive definite


If your covariance matrix is singular, then you really should consider why the matrix is singular and come up with a higher-level approach that avoids the singularity.

However, if you insist on finding a Cholesky factorization somehow, you should look at modified Cholesky factorization algorithms that perturb the covariance as little as possible to make it positive definite and produce a Cholesky factorization for the perturbed matrix.

See the answer to this mathematics SE question for a pointer to recent research on this topic:


| cite | improve this answer | |
  • $\begingroup$ Specifically, consider doing a truncated SVD or QR decomposition that stops when you reach the zero eigenvalues. $\endgroup$ – Wolfgang Bangerth Sep 2 '19 at 22:23
  • $\begingroup$ One word of caution for modified Cholesky schemes - they typically run much slower than regular Cholesky since as far as I know the diagonal perturbation matrix makes it difficult to use Level 3 BLAS operations. This may be a consideration depending how large your matrix is. $\endgroup$ – vibe Sep 3 '19 at 4:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.