比较 R 中的矩阵求逆 - Cholesky 方法有什么问题？答案

【问题标题】：Comparing matrix inversions in R - what is wrong with the Cholesky method?比较 R 中的矩阵求逆 - Cholesky 方法有什么问题？
【发布时间】：2014-12-30 08:25:11
【问题描述】：

我比较了计算对称矩阵逆的各种方法：

解决（来自 LAPCK 包）
求解（但使用更高的机器精度）
qr.solve（据说更快）
ginv（MASS 包，Moore-Penrose 算法的实现）
chol2inv（使用 Cholesky 分解）

通过它们的特征值比较逆矩阵：

R
library(MASS)

## Create the matrix
m = replicate(10, runif(n=10)) 
m[lower.tri(m)] = t(m)[lower.tri(m)]

## Inverse the matrix
inv1 = solve(m)
inv2 = solve(m, tol = .Machine$double.eps)
inv3 = qr.solve(m)
inv4 = ginv(m)
inv5 = chol2inv(m)

## Eigenvalues of the inverse
em1=eigen(inv1)
em2=eigen(inv2)
em3=eigen(inv3)
em4=eigen(inv4)
em5=eigen(inv5)

## Plot the abs of the eigenvalues (may be complex)
myPch=c( 20, 15, 17, 25, 3 )
plot(abs(em1$values),pch=myPch[1],cex=1.5)
points(abs(em2$values),pch=myPch[2], cex=1.5)
points(abs(em3$values),pch=myPch[3], cex=1.5)
points(abs(em4$values),pch=myPch[4], cex=1.5)
points(abs(em5$values),pch=myPch[5], cex=1.5)
legend( "topright", c("solve","solve-double","solve-fast","Moore-Penrose","Cholesky"), pch=myPch )

如您所见，Cholesky 方法给出的逆显然与其他方法不同。

根据这篇文章，如果矩阵是对称的（在我们的例子中是对称的），则首选 Cholesky 方法： Matrix inversion or Cholesky?

但是solve()是“官方普遍”的R方法来反转方法，我可能会误解一些东西......

有什么好的建议吗？

提前致谢，

【问题讨论】：

您的示例无法重现，因为您没有设置随机种子。但是，请检查矩阵是否为正定矩阵。此外，一般建议是避免反转矩阵。
另外，当然应该是inv5 = chol2inv(chol(m))。
@Roland Nitpicking;一般建议是避免显式反转矩阵。使用分解通常是一种避免对原始矩阵进行显式求逆的方法。
@Roland & 如果有人按照您的建议执行chol2inv(chol(m))，对于我刚刚运行的随机示例，R 声明：Error in chol.default(m) : the leading minor of order 3 is not positive definite，从而突出了 OP 面临的问题。 +1
改进我的代码以使矩阵为正定且对称的穷人方法：library(caper); t = rtree(n=10); m = vcv(t)

标签： r matrix matrix-inverse

【解决方案1】：

您需要将 Cholesky 分解传递给chol2inv：

inv5 = chol2inv(chol(m))

如果m 是正定的（它可能不适用于您不可重现的输入），这应该与其他方法给出相同的结果。

【讨论】：