Lapack：Cholesky矩阵分解问题

0 投票

2 回答

1387 浏览

提问于 2025-04-18 03:00

问题 1

有没有人能推荐一种更简单的方法来在Python中进行Cholesky分解？特别是最后一行让我觉得很别扭。

SigmaSqrt = matrix(Sigma)
cvxopt.lapack.potrf(SigmaSqrt)
SigmaSqrt = matrix(np.tril(SigmaSqrt))

问题 2

我遇到一个问题，就是有一整行和一整列（比如第一行的所有元素和第一列的所有元素）都是零，这样的话，lapack会报错说矩阵不是正定的。处理这个问题的最佳方法是什么？

目前我这样做：（感觉超级别扭……）

try:
    SigmaSqrt = matrix(Sigma)
    cvxopt.lapack.potrf(SigmaSqrt)
    SigmaSqrt = matrix(np.tril(SigmaSqrt))
except ArithmeticError:
    SigmaSqrt = matrix(Sigma.ix[1:,1:])
    cvxopt.lapack.potrf(SigmaSqrt)
    SigmaSqrt = matrix(np.tril(SigmaSqrt))
    SigmaSqrt = sparse([[v0],[v0[1:].T, SigmaSqrt]])

error handling numerical methods lapack cholesky decomposition positive definite matrix matrix factorization

2 个回答

另一个选择是 chompack 模块：chompack 的主页，里面有 chompack.cholesky 这个功能。我自己也在用这个模块，通常是和 cvxopt 模块一起使用。它和 cvxopt 里的（稀疏）矩阵配合得非常好。

回答于 2025-04-18 由 Python大师

分享举报

你可以直接使用 numpy.linalg.cholesky。另外，如果某一列或某一行全是零，那么这个矩阵就会是奇异的，至少会有一个特征值为零，因此它就不是正定的。因为Cholesky分解只适用于“厄米（如果是实数就是对称的）且正定”的矩阵，所以在这种情况下是无法使用的。

补充：要“解决”你的问题，取决于你想要什么。你所做的任何让它工作的操作，得到的Cholesky分解都不会是原始矩阵的Cholesky分解。如果你在做一个迭代过程，并且可以稍微调整一下，如果矩阵已经是对称的，可以使用 numpy.linalg.eigvalsh 来找到矩阵的特征值。设d为最小的特征值。然后设置 A += (abs(d) + 1e-4) * np.indentity(len(A))。这样可以让它变成正定的。

补充：这是在Levenberg–Marquardt算法中使用的一个技巧。这里有一篇关于牛顿法的维基百科文章提到过这个内容，因为关于Levenberg–Marquardt的文章没有详细说明。此外，这里还有一篇相关的论文 here。基本上，这个方法会将所有特征值都向 (abs(d) + 1e-4) 的方向移动，从而使它们都变为正值，这样就满足了矩阵是正定的条件。