r - 为什么乔列斯基分解不能得到与简单矩阵求逆相同的结果？

Question

设置重现我的最小工作示例

我有以下矩阵

K <- matrix(c(1.250000e+00, 3.366892e-07, 4.641930e-10, 1.455736e-08, 1.049863e-06, 
              3.366892e-07, 1.250000e+00, 5.482775e-01, 8.606555e-01, 9.776887e-01,
              4.641930e-10, 5.482775e-01, 1.250000e+00, 8.603413e-01, 4.246732e-01,
              1.455736e-08, 8.606555e-01, 8.603413e-01, 1.250000e+00, 7.490100e-01,
              1.049863e-06, 9.776887e-01, 4.246732e-01, 7.490100e-01, 1.250000e+00), nrow=5)

以及以下向量

y <- matrix(c(39.13892, 12.34428, 12.38426, 14.71951, 10.52160), nrow=5)

问题

我想计算 K 的倒数与向量 y 之间的乘积。

简单的方法 - 它有效

天真的方法有效(我有一种检查方法，但在这里并不重要)

solve(K) %*% y

           [,1]
[1,] 31.3111308
[2,]  3.0620869
[3,]  3.7383357
[4,]  6.6257060
[5,]  0.7820081

Cholesky 分解 - 失败

然而“聪明”的方法失败了。我使用乔列斯基分解，它给出了上三角矩阵。然后，我通过后向替换求解系统 L z = y，并通过前向替换求解系统 L^T x = L^{-1} y。

L <- chol(K)  ## upper triangular
forwardsolve(t(L), backsolve(L, y))

          [,1]
[1,]  31.31112
[2,] -14.16259
[3,]   9.84534
[4,]  39.67900
[5,]  33.54842

发生了什么事？这个矩阵K和这个向量“y”只是一个例子。许多其他类似的向量和矩阵也会发生这种情况。为什么？

Answer 1

关键点是，当取乘积的逆时，必须反转逆的乘积:

solve(A %*% B) = solve(B) %*% solve(A)

在问题中，顺序没有颠倒。

如果R = chol(K)，其中我们使用R来强调它是右上三角形，那么:

solve(K, y)
= solve(t(R) %*% R, y)   since K = t(R) %*% R
= solve(t(R) %*% R) %*% y
= solve(R) %*% solve(t(R)) %*% y  note that we have reversed the order
= solve(R) %*% solve(t(R), y)
= backsolve(R, forwardsolve(t(R), y))

在最后一行中，我们使用了 R 的转置是左下三角矩阵的事实，并且 forwardsolve 适用于此类矩阵，而 backsolve 适用于右上三角矩阵。

我们可以检查这确实给出了与使用 solve direclty 相同的答案:

R = chol(K)
all.equal(backsolve(R, forwardsolve(t(R), y)), solve(K, y))
# [1] TRUE