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 = [email protected](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


1 Answer 1


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:


  • $\begingroup$ Specifically, consider doing a truncated SVD or QR decomposition that stops when you reach the zero eigenvalues. $\endgroup$ Sep 2, 2019 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, 2019 at 4:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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