为什么 cholesky 分解没有给我与简单地反转矩阵相同的结果?
Why is cholesky decomposition not giving me the same result as simply inverting the matrix?
设置重现我的最小工作示例
我有以下矩阵
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 分解 - 失败
但是 "clever" 方法失败了。我使用 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' 只是一个例子。许多其他类似的向量和矩阵都会发生这种情况。为什么?
这只是部分答案,但聊胜于无。基本上求K^{-1}y相当于求解下面的系统
使用我们的 cholesky 分解来写这个
基本上我们现在首先考虑Lz
作为我们的变量,称它为x
并且由于 L
是上三角,它的转置是下三角,所以我们可以使用 forwardsolve
找到 x
forwardsolve(t(L), y)
一旦我们找到 x
,我们就可以记住它的意思,即
这样我们就可以使用'backsolveto find
z`
backsolve(L, forwardsolve(t(L), y))
这给出了正确答案。不确定为什么反过来也行不通。
关键是,当取一个乘积的倒数时,必须反转倒数的乘积:
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
设置重现我的最小工作示例
我有以下矩阵
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 分解 - 失败
但是 "clever" 方法失败了。我使用 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' 只是一个例子。许多其他类似的向量和矩阵都会发生这种情况。为什么?
这只是部分答案,但聊胜于无。基本上求K^{-1}y相当于求解下面的系统
使用我们的 cholesky 分解来写这个
基本上我们现在首先考虑Lz
作为我们的变量,称它为x
并且由于 L
是上三角,它的转置是下三角,所以我们可以使用 forwardsolve
找到 x
forwardsolve(t(L), y)
一旦我们找到 x
,我们就可以记住它的意思,即
这样我们就可以使用'backsolveto find
z`
backsolve(L, forwardsolve(t(L), y))
这给出了正确答案。不确定为什么反过来也行不通。
关键是,当取一个乘积的倒数时,必须反转倒数的乘积:
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