用于病态/奇异协方差矩阵的 Cholesky

计算科学 线性代数 矩阵分解
2021-11-28 22:47:52

有人可以建议一种方法来获得奇异协方差矩阵的 Cholesky 分解吗?我需要它来匹配满秩矩阵上的 Cholesky,即应该保留坐标顺序。我在下面的尝试是在 scipy 中使用 ldl 例程,但这让我对不同的排序进行分解,有什么想法吗?

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

1个回答

如果您的协方差矩阵是奇异的,那么您真的应该考虑为什么矩阵是奇异的,并提出一种避免奇异性的更高级别的方法。

但是,如果您坚持以某种方式找到 Cholesky 分解,您应该查看修改后的 Cholesky 分解算法,该算法尽可能少地扰动协方差以使其为正定并为扰动矩阵生成 Cholesky 分解。

请参阅此数学 SE 问题的答案,以获取有关该主题的最新研究的指针:

https://math.stackexchange.com/questions/2946069/what-c​​an-be-added-to-a-non-positive-defined-matrix-to-make-it-positive-defined