R中的scipy.linalg.cho_solve对应物？

Question

我想知道在R中是否有一个与scipy.linalg.cho_solve相对应的函数，这个函数的作用是给定A (A = LL')和b的乔勒斯基因子L，它解决了原始问题Ax = b (而不是Lx = b)。

(所以它不同于backsolve/forwardsolve)

谢谢你，Joon

Answer 1

我意识到这个问题有点陈旧，但我看到了答案

forwardsolve(L, forwardsolve(L, b), transp=TRUE)

还没有给出。这使用了三角形结构，同时保持了原始问题。对于较大的矩阵，这应该更快、更准确。还值得注意的是，由于chol返回的是上三角矩阵，所以使用L <- t(chol(A))。

A <- matrix(c(1,1,1,1,5,5,1,5,14), nrow=3)
# Cholesky decomposition A = LL'
L <- t(chol(A))

# Make some b with known x
x <- c(1, 2, 3)
b <- A %*% x

# Solve
forwardsolve(L, forwardsolve(L, b), transp=TRUE)

给出答案：

> forwardsolve(L, forwardsolve(L, b), transp=TRUE)
     [,1]
[1,]    1
[2,]    2
[3,]    3

Answer 2

我想不出有什么函数能自动为你做到这一点，但假设你有乔勒斯基因子L，很容易通过重构分解A=LL定义的A矩阵在一行中完成：

 A=matrix(c(1,1,1,1,5,5,1,5,14),nrow=3)
 # Cholesky decomposition A = LL'
 L <- chol(A)

 # Make some b with known x
 x <- c(1,2,3)
 b <- A%*%x

 # Solve
 solve( t(L) %*% L, b)

编辑:请注意，在R中，Cholesky因子的定义与A= L'L相关，这就是为什么您必须将转置放在求解的第一位。

edit2 :读了贝茨的文章后，我意识到它应该是：

> solve(crossprod(L),b)
     [,1]
[1,]    1
[2,]    2
[3,]    3

Answer 3

如果我没理解错的话，那么Doug Bates在2004年为R新闻写的一篇article中就谈到了这一点(参见linK的第18页)。

相关的位是：

ch <- chol(crossprod(X))
chol.sol <- backsolve(ch, forwardsolve(ch, crossprod(X, y),
                                       upper = TRUE, trans = TRUE))

其中X是预测变量的矩阵。

Doug的文章继续展示了如何使用Matrix包(随R提供)中的功能来快速解决相同的系统。