Cholesky for ill-conditioned/singular covariance matrices

Question

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

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:

https://math.stackexchange.com/questions/2946069/what-can-be-added-to-a-non-positive-definite-matrix-to-make-it-positive-definite