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

Question

这个问题与this one有关而这个 one

我有两个满秩矩阵 A1、A2 每个
维度为 p x p 和 p 向量 y。

这些矩阵在某种意义上是密切相关的
矩阵 A2 是矩阵 的一阶更新A1 .

我对向量感兴趣

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

在哪里

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

和

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

现在，在之前的 question 中在这里我被告知
用 Choleski 方法估计 β2，因为 Choleski
使用 R 函数(如 chud())很容易更新分解。
包装内 采样器比较 .

下面是在 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)是 的 R 行话一个 ' 乙 (即 A 矩阵的转置
乘以矩阵/向量 B)。

solve(A,b)求解 x 线性系统 一个 x=b。

chol(A)是 PSD 矩阵的 Choleski 分解 一个 .

chol2inv从 的 QR 分解的(R 部分)计算 ( X ' X )-1 X .

Answer 1

正如您所提到的，您的“fx01”实现有些幼稚，并且比“fx03”方法执行的工作要多得多。在线性代数中(我为不支持 LaTeX 的主要 StackOverflow 道歉!)，'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' 例程，它可能执行带有部分旋转(2/3 n^3 触发器)和两个的 LU 分解三角形在 2n^2 次触发器中求解(加上旋转)。因此，您的“fx01”方法渐近执行了三倍的工作，这与您的实验结果惊人地一致。请注意，如果 A 是实对称或复 Hermitian，那么 LDL^T 或 LDL' 分解和求解只需要一半的工作。

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