递归矩阵的r快速叉积

Question

我有一个矩阵a的递归计算（这将是一个帽子矩阵），例如：

A(i) = A(i-1) + crossprod(B,A(i-1))

对于每个步骤CCD_ 2，我需要CCD_。有没有比以下实现更快的方法在R中实现这一点：

# define random matrices
set.seed(123)
n <- 7^2*10^4
steps <- 10
A <- matrix(rnorm(n), ncol=sqrt(n))
B <- matrix(rnorm(n), ncol=sqrt(n))
# preallocation
Amat <- traceA <- vector("list", steps)
Amat[[1]] <- A
# recursive computation for matrix A(i)
ptm <- proc.time()
for(i in 2:steps){
  Amat[[i]] <- Amat[[i-1]] + crossprod(B,Amat[[i-1]])
  traceA[[i]] <- sum(diag(Amat[[i]]))
}
proc.time() - ptm

我想提到的是，矩阵A（I）和矩阵B是对称的和幂等的（因为它们是线性模型的帽矩阵），并且可以非常大。我想并行计算在这里会失败，因为for循环需要之前步骤的矩阵A（I-1）。

这背后的想法是一种基于可能性的提升算法，其中我需要可以如上所述计算的帽子矩阵的每个提升迭代的轨迹。

Answer 1

看起来你的Amat_i可以写成Amat_i = (1+t(B))^(i-1) * A，既然你提到了B*B = B或t(B)*t(B) = t(B)，那么

 (1+B)^n = 1 + choose(n,1)*B + choose(n,2)*B^2 + ...
         = 1 + B * (choose(n,1) + choose(n,2) + ... + choose(n,n))
         = 1 + B * (2^n - 1)

把所有这些放在一起：

 tr(Amat_i) = tr(A) + (2^(i-1) - 1) * tr(t(B)*A)

所以只需要计算两个轨迹，就不需要再做任何矩阵乘法就可以得到所有的tr(Amat_i)。