r - 在 R 中高效使用 Choleski 分解

Question

此问题与 this one 相关还有这个one

我有两个满秩矩阵 A₁, A₂ 维度为 p x p 和 p 向量 y。

这些矩阵密切相关，因为 矩阵A₂是矩阵A₁的一级更新。

我对向量感兴趣

β₂ | (β₁, y, A₁, A₂, A₁^-1})

哪里

β₂ = (A₂' A₂)^-1(A₂'y)

和

β₁ = (A₁' A₁)^-1(A₁' y)

现在，在之前的 question 中在这里我被告知 自 Choleski 方法以来，通过 Choleski 方法估计 β₂ 使用 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，如果是，如何实现？

1. A %*% B 是 A 与 B 矩阵乘法的 R 语言。
2. crossprod(A,B) 是 A' B 的 R 术语(即 A 矩阵的转置 乘以矩阵/向量 B)。
3. solve(A,b) 求解 x 线性方程组 A x=b。
4. chol(A) 是 PSD 矩阵 A 的 Choleski 分解。
5. chol2inv 根据 QR 分解的(R 部分)计算 (X' X)^-1 X。

Answer 1

正如您所提到的，您的“fx01”实现有些幼稚，并且比“fx03”方法执行的工作要多得多。在线性代数中(我对主 StackOverflow 不支持 LaTeX 表示歉意!)，“fx01”执行:

• B := A' A 大约需要 n^3 次失败。
• L := chol(B) 大约需要 1/3 n^3 次失败。
• L := inv(L) 大约需要 1/3 n^3 次失败。
• B := L' L 大约需要 1/3 n^3 次失败。
• z := A y 大约需要 2n^2 次失败。
• x := B z 大约需要 2n^2 次失败。

因此，成本看起来与 2n^3 + 4n^2 非常相似，而您的 'fx03' 方法使用默认的 'solve' 例程，该例程可能会执行部分旋转的 LU 分解(2/3 n^3 flops) )和两个三角形在 2n^2 次失败中求解(加上旋转)。因此，您的“fx01”方法渐进地执行了三倍的工作，这与您的实验结果惊人地一致。请注意，如果 A 是实对称或复厄米特矩阵，则 LDL^T 或 LDL' 因式分解和求解只需要一半的工作量。

话虽如此，我认为您应该将 A' A 的 Cholesky 更新替换为 A 的更稳定的 QR 更新，正如我刚刚回答的 in your previous question 。 QR 分解大约需要 4/3 n^3 次浮点运算，而 QR 分解的秩一更新仅为 O(n^2)，因此，当存在多个相关解时，这种方法仅对一般 A 有意义。只是一个一级修改。