Lapack：Cholesky 矩阵分解问题答案

【问题标题】：Lapack: Cholesky matrix factorization problemsLapack：Cholesky 矩阵分解问题
【发布时间】：2014-04-16 16:32:03
【问题描述】：

问题 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]])

【问题讨论】：

标签： python matrix lapack cvxopt

【解决方案1】：

您可以使用numpy.linalg.cholesky。此外，如果一列或一行全部为零，则矩阵将是奇异的，至少有一个特征值为零，因此不是正定的。由于 Cholesky 仅针对“Hermitian（实值对称）和正定矩阵定义，因此它不适用于它。

编辑：“处理”您的问题取决于您想要什么。您为使其工作所做的任何事情都会产生一个不会是原始矩阵的 Cholesky 的 Cholesky。如果你正在做一个迭代过程并且它可能会有点糊涂，如果它已经是对称的，那么使用numpy.linalg.eigvalsh 来查找矩阵的特征值。令 d 为最负的特征值。然后设置A += (abs(d) + 1e-4) * np.indentity(len(A))。这将使它成为肯定的。

编辑：这是 Levenberg-Marquardt 算法中使用的一个技巧。这个链接到Newton's Method 上的一篇维基百科文章，其中提到它，因为关于 Levenberg–Marquard 的文章没有涉及到这个。此外，here 也是一篇关于它的论文。基本上，这会将所有特征值移动(abs(d) + 1e-4)，从而使它们全部为正，这是矩阵为正定的充分条件。

【讨论】：

感谢您的帮助。您在编辑中描述的过程的“名称”是什么，以便我可以阅读我的线性代数？
@ArikRaffaelFunke 添加了更多信息。这是 Levenberg-Marquardt 算法中使用的一个技巧。
另外，还有scipy.linalg.cho_factor 和scipy.linalg.cho_solve 可以串联使用，使用Cholesky 解决Ax=b。

【解决方案2】：

另一个选项是 chompack 模块：chompack homepage chompack.cholesky 我自己将它与 cvxopt 模块结合使用。它与来自 cvxopt 的（稀疏）矩阵完美配合。

【讨论】：

您是否知道chompack.cholesky 和numpy.linalg.cholesky 之间的密集矩阵是否存在速度差异？还是您使用chompack 仅仅是因为它使用cvxopt 矩阵很方便？