递归矩阵的快速交叉积

Question

我有一个矩阵A的递归计算(这将是一个帽子矩阵)，例如：

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

对于每一步i，我都需要A(i)的跟踪。在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是对称的和幂等的(因为它们是线性模型的hat矩阵)，并且可以是非常大的。我想并行计算在这里会失败，因为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)了。