r - 加速二次形式的评估

Question

我的问题是另一个“矢量化!”。类似的问题出现在其他地方( Efficient way of calculating quadratic forms: avoid for loops? )，但不知何故我似乎无法让它适用于我的情况。

我想为大小为 x'Sx 的样本中的每个 p 维观察 x 计算二次形式 n 。我想不出一个漂亮的矢量化代码，所以我最后的办法是通过 for loop 来完成。以下玩具示例适用于 p=2 、 n=100 。

set.seed(123)
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
x <- cbind(x1,x2)
Sigma <- matrix(c(1, 2, 3, 4), ncol = 2)
z  <- rep(0, n)
for (i in 1:n) {
   z[i]  <- x[i, ] %*% solve(Sigma, x[i, ]) #quadratic form of x'S^{-1}x
}

像许多其他崇拜矢量化代码的 R 用户一样，for 循环的使用引起了情绪上的痛苦。为了减轻痛苦，我使用几种常见的矢量化技术修改了我的代码。

ap <- function(Sigma, x) apply(x, 1, function(x) x %*% solve(Sigma, x))
lap <- function(Sigma, x) unlist(lapply(1:n, function(i) x[i, ] %*% solve(Sigma, x[i, ])))
loop <- function(Sigma, x){
  z  <- rep(0, n)
  for (i in 1:n) {
    z[i]  <- x[i, ] %*% solve(Sigma, x[i, ])
  }
  z
}

但是速度比较显示没有什么收获。

library(microbenchmark)
microbenchmark(lap(Sigma, x), ap(Sigma, x), loop(Sigma, x))

# Unit: milliseconds
#           expr      min       lq     mean   median       uq       max neval
#  lap(Sigma, x) 4.207434 4.444895 5.092389 4.616912 5.283504  8.440802   100
#   ap(Sigma, x) 4.360204 4.523306 5.317304 4.685396 5.412771 10.168674   100
# loop(Sigma, x) 4.518645 4.679317 6.204626 4.827831 5.438908 94.115144   100

是否有任何改进的余地，或者我应该去 Rcpp 以摆脱使用 for loops 的罪恶？

Answer 1

如果将 x 的行存储在列表中并使用 vapply 而不是 lapply ，则可以按如下方式加快速度

# First, make a list of the rows of x
xl <- vector("list",nrow(x))
for (i in seq_along(xl)) xl[[i]] <- x[i, ] 

# Apply solve
solve.mat <- vapply(xl, solve, numeric(2), a = Sigma)
# Take the dot product of each pair of elements
result <- colSums(solve.mat * t(x))
all(result == lap(Sigma, x))
# [1] TRUE

一步写出来并比较

library(microbenchmark)
microbenchmark(lap = lap(Sigma, x),
    csums = colSums(vapply(xl, solve, numeric(2), a = Sigma) * t(x)))
# Unit: milliseconds
#   expr      min       lq     mean   median       uq      max neval
#    lap 3.013343 3.050855 3.164558 3.097901 3.136355 4.206923   100
#  csums 2.224350 2.263772 2.354349 2.289751 2.317672 3.660294   100