R：以矩阵形式写双重求和

Question

我想摆脱 $t$ 的求和并在 R 中实现它。到目前为止，我有这个

Z_diff = mapply(FUN = function(i) {Z-Z[i]}, i=1:n, SIMPLIFY = TRUE)
rss = 0
for(t in 1:n){
    rss = rss + t(Y - X %*% B[,t]) %*% (diag(c(Z_diff[,t])) %*% (Y - X %*% B[,t]))
  }

Answer 1

下面我们定义一个测试输入，并使用问题中的原始公式计算 rss。然后我们展示给出相同结果的矩阵公式。如果你在最后一行的总和中写出矩阵分量的公式，结果的原因应该很清楚。

# input
set.seed(123)
n <- 4
Z <- rnorm(n)
Y <- rnorm(n)
B <- matrix(rnorm(n*n), n)
X <- matrix(rnorm(n*n), n)

# scalar sum
rss <- 0
for(i in 1:n) 
  for(t in 1:n) 
    rss <- rss + (Y[i] - c( crossprod(X[i, ], B[, t])))^2 * (Z[t] - Z[i])
rss
## [1] 43.39269

# matrix sum
ZZ <- t(outer(Z, Z, "-"))
sum(ZZ * (Y - X%*%B)^2)
## [1] 43.39269

性能

使用上面的数据，矩阵公式在我的机器上运行速度提高了大约 300 倍。

library(microbenchmark)

microbenchmark(
scalar = {
  rss <- 0
  for(i in 1:n) 
    for(t in 1:n) 
      rss <- rss + (Y[i] - c( crossprod(X[i, ], B[, t])))^2 * (Z[t] - Z[i])
  rss
},
matrix = {
  ZZ <- t(outer(Z, Z, "-"))
  sum(ZZ * (Y - X%*%B)^2)
})
+ })
## Unit: microseconds
##    expr     min       lq      mean   median      uq      max neval cld
##  scalar 34428.2 36966.35 44244.423 38257.80 40968.6 407679.6   100   b
##  matrix    80.8    82.45   148.764   156.45   177.8    721.7   100  a

R：以矩阵形式写双重求和

R: Writing double summations in matrix form

r

matrix

linear-algebra

vectorization

性能