3

我有一个矩阵 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 是对称且幂等的（因为它们是线性模型的帽子矩阵）并且可以非常大。我猜这里并行计算会失败，因为 for 循环需要前一步的矩阵 A(i-1)。

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

4

1 回答 1

3

看起来你Amat_i的 's 可以写成Amat_i = (1+t(B))^(i-1) * Aand 既然你提到了B*B = Bor 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)'s。

于 2013-07-30T17:05:50.887 回答