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

【问题标题】：R: Writing double summations in matrix formR：以矩阵形式编写双重求和
【发布时间】：2023-04-07 04:41:01
【问题描述】：

我想摆脱对 $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]))
  }

【问题讨论】：

你的问题是什么？
@deschen 如何摆脱 $t$ 的求和以避免使用 for 循环？

标签： r matrix vectorization linear-algebra

【解决方案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

【讨论】：