在 R 中高效使用 Choleski 分解
这个问题与这个相关< /a> 和这个一个
我有两个满秩矩阵 A1, A2 尺寸为 pxp 和 p 向量 y。
这些矩阵在某种意义上是密切相关的 矩阵A2是矩阵A1的一级更新。
我对矢量感兴趣
β2 | (β1, y, A1, A2, A1-1})
其中
β2 = (A2' A2)-1(A2'y)
和
β1 = (A1' A1)-1(A1' y)
现在，在之前的 问题在这里我被告知 自 Choleski 方法以来，通过 Choleski 方法估计 β2 使用 R 函数（例如 chud()）可以轻松更新分解 在SamplerCompare包中。
下面是 R 中求解线性系统的两个函数，第一个使用 solve() 函数和第二个 Choleski 方法 （我可以有效更新第二个）。
fx01 <- function(ll,A,y) chol2inv(chol(crossprod(A))) %*% crossprod(A,y)
fx03 <- function(ll,A,y) solve(A,y)
p <- 5
A <- matrix(rnorm(p^2),p,p)
y <- rnorm(p)
system.time(lapply(1:1000,fx01,A=A,y=y))
system.time(lapply(1:1000,fx03,A=A,y=y))
我的问题是：对于 p 小，两个函数似乎具有可比性 （实际上 fx01 甚至更快）。但当我增加 p 时， fx01 变得越来越慢，因此对于 p = 100， fx03 的速度是 fx01 的三倍。
是什么原因导致 fx01 性能下降？ 得到改进/解决（也许我对 Choleski 的实现太天真了？我不应该使用 Choleski 星座的函数，例如 backsolve，如果是的话，如何实现？
A %*% B是 A 与 B 矩阵乘法的 R 术语。crossprod(A,B)是 A' 的 R 术语B（即A矩阵的转置 乘以矩阵/向量 B)。solve(A,b)求解 x 线性方程组 A x=b。chol(A)是 PSD 矩阵 A 的 Choleski 分解。chol2inv根据 X' X)-1 计算 QR 分解的（R 部分） >X。
This question is related to this one and this one
I have two full rank matrices A1, A2 each
of dimension p x p and a p-vector y.
These matrices are closely related in the sense that
matrix A2 is a rank one update of matrix A1.
I'm interested in the vector
β2 | (β1, y, A1, A2, A1-1})
where
β2 = (A2' A2)-1(A2'y)
and
β1 = (A1' A1)-1(A1' y)
Now, in a previous question here I have been advised
to estimate β2 by the Choleski approach since the Choleski
decomposition is easy to update using R functions such as chud()
in package SamplerCompare.
Below are two functions to solve linear systems in R, the first one uses
the solve() function and the second one the Choleski approach
(the second one I can efficiently update).
fx01 <- function(ll,A,y) chol2inv(chol(crossprod(A))) %*% crossprod(A,y)
fx03 <- function(ll,A,y) solve(A,y)
p <- 5
A <- matrix(rnorm(p^2),p,p)
y <- rnorm(p)
system.time(lapply(1:1000,fx01,A=A,y=y))
system.time(lapply(1:1000,fx03,A=A,y=y))
My question is: for p small, both functions seems to be comparable
(actually fx01 is even faster). But as I increase p,fx01 becomes increasingly slower so that for p = 100,fx03 is three times as fast as fx01.
What is causing the performance deterioration of fx01 and can it
be improved/solved (maybe my implementation of the Choleski is too naive? Shouldn't I be using functions of the Choleski constellation such as backsolve, and if yes, how?
A %*% Bis the R lingo for matrix multiplication of A by B.crossprod(A,B)is the R lingo for A' B (ie transpose of A matrix
multiplying the matrix/vector B).solve(A,b)solves for x the linear system A x=b.chol(A)is the Choleski decomposition of a PSD matrix A.chol2invcomputes (X' X)-1 from the (R part) of the QR decomposition of X.
如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱，如果其他用户或者作者回复了您的评论，将不能在第一时间通知您！
发布评论
评论（1）
正如您所提到的，您的“fx01”实现有些幼稚，并且比“fx03”方法执行的工作要多得多。在线性代数中（我对主 StackOverflow 不支持 LaTeX 表示歉意！），“fx01”
因此，成本看起来与 2n^3 + 4n^2 非常相似，而您的“fx03”方法使用默认的“求解”例程，该例程可能会执行带有部分旋转（2/3 n^3 flops）和两个的 LU 分解三角形在 2n^2 次失败中求解（加上旋转）。因此，您的“fx01”方法渐进地执行了三倍的工作，这与您的实验结果惊人地一致。请注意，如果 A 是实对称或复厄米特矩阵，则 LDL^T 或 LDL' 因式分解和求解只需要一半的工作量。
话虽如此，我认为您应该将 A' A 的 Cholesky 更新替换为 A 的更稳定的 QR 更新，正如我刚刚回答的 在你之前的问题中。 QR 分解大约需要 4/3 n^3 次浮点运算，而 QR 分解的秩一更新仅为 O(n^2)，因此，当存在多个相关解时，这种方法仅对一般 A 有意义。仅仅是一级修改。
Your 'fx01' implementation is, as you mentioned, somewhat naive and is performing far more work than the 'fx03' approach. In linear algebra (my apologies for the main StackOverflow not supporting LaTeX!), 'fx01' performs:
Thus, the cost looks very similar to 2n^3 + 4n^2, whereas your 'fx03' approach uses the default 'solve' routine, which likely performs an LU decomposition with partial pivoting (2/3 n^3 flops) and two triangle solves (plus pivoting) in 2n^2 flops. Your 'fx01' approach therefore performs three times as much work asymptotically, and this amazingly agrees with your experimental results. Note that if A was real symmetric or complex Hermitian, that an LDL^T or LDL' factorization and solve would only require half as much work.
With that said, I think that you should replace your Cholesky update of A' A with a more stable QR update of A, as I just answered in your previous question. A QR decomposition costs roughly 4/3 n^3 flops and a rank-one update to a QR decomposition is only O(n^2), so this approach only makes sense for general A when there is more than just one related solve that is simply a rank-one modification.