为什么在计算矩阵的所有 rows/columns 的范数时使用 for 循环比应用更快？

Question

考虑以下因素

n <- 10^4
p <- 2
foo <- matrix(runif(p*n), n, p)

我想计算矩阵每一行的范数，即计算 sqrt(crossprod(a_i))，其中 a_i 是 foo 的第 i 行。我可以用 apply 或 for 循环来做到这一点：

for_loop <- function(x){
  range <- seq_along(x[,1])
  foo <- range
  for (i in range){
    foo[i] <- sqrt(crossprod(x[i,]))
    }
  foo
}

use_apply <- function(x){
  apply(x, 1, function(r) sqrt(crossprod(r)))
}

我认为更简单的 apply 代码会更快，但是：

> microbenchmark(for_loop(foo), use_apply(foo), times = 1000)
Unit: milliseconds
           expr      min       lq     mean   median       uq      max neval
  for_loop(foo) 16.07111 18.87690 24.25369 20.78997 27.66441 179.8374  1000
 use_apply(foo) 24.77948 29.05891 35.98689 31.89625 40.30085 205.1632  1000

请注意，times = 1000 可能需要相当长的时间，如果您的机器速度不快，您可能需要使用 microbenchmark 默认值。为什么 apply 比 for 循环代码慢？ purrr 中是否有一些函数会更快？

EDIT 我无法相信 crossprod(x) 会比 sum(x*x) 慢这么多，所以我想检查 Emmanuel-Lin 的结果。我得到非常不同的时间：

my_loop <- function(x){
  range <- seq_along(x[,1])
  foo <- range
  for (i in range){
    foo[i] <- sqrt(sum((x[i,] *x[i,])))
    }
  foo
}

my_apply <- function(x){
  apply(x, 1, function(r) sqrt(sum(r*r)))
}

for_loop <- function(x){
  range <- seq_along(x[,1])
  foo <- range
  for (i in range){
    foo[i] <- sqrt(crossprod(x[i,]))
  }
  foo
}

use_apply <- function(x){
  apply(x, 1, function(r) sqrt(crossprod(r)))
}
> microbenchmark(for_loop(foo), my_loop(foo), use_apply(foo), my_apply(foo))
Unit: milliseconds
           expr       min       lq     mean   median       uq       max neval
  for_loop(foo) 16.299758 17.77176 21.59988 19.04428 22.44558 131.33819   100
   my_loop(foo)  9.950813 12.02106 14.43540 12.66142 15.26865  45.42030   100
 use_apply(foo) 25.480019 27.95396 31.98351 29.85244 36.41599  60.88678   100
  my_apply(foo) 13.277354 14.98329 17.60356 15.98103 19.70325  34.07097   100

好的，my_apply 和 my_loop 更快（我仍然不敢相信！什么，crossprod 是针对慢速优化的吗？:-/）但是 不像 Emmanuel-lin 发现的那样快。这可能与 crossprod 执行的一些维度一致性检查有关。

Answer 1

如果您检查代码，

apply 实际上是一个 R for 循环：

   #only the for-loop code shown here
   if (length(d.call) < 2L) {
        if (length(dn.call)) 
            dimnames(newX) <- c(dn.call, list(NULL))
        for (i in 1L:d2) {
            tmp <- forceAndCall(1, FUN, newX[, i], ...)
            if (!is.null(tmp)) 
                ans[[i]] <- tmp
        }
    }
    else for (i in 1L:d2) {
        tmp <- forceAndCall(1, FUN, array(newX[, i], d.call, 
            dn.call), ...)
        if (!is.null(tmp)) 
            ans[[i]] <- tmp
    }

除上述之外，apply 还将运行进行一系列检查，以确保您提供的参数正确无误。是上面的让它慢了一点。

但是，lapply、sapply 和 vapply 是基于 C 的 for 循环，因此比基于 R 的 for 循环快得多。

Answer 2

要在 RAM 上完成@LyzandeR 的回答。

您可以通过自己编写乘法来更快地执行计算：

将 crossprod 替换为 sum(r * r)

my_loop <- function(x){
  range <- seq_along(x[,1])
  foo <- range
  for (i in range){
    foo[i] <- sqrt(sum(x[i,] * x[i,]))
  }
  foo
}


my_sapply <- function(x){
  apply(x, 1, function(r) sqrt(sum(r * r)))
}

microbenchmark(for_loop(X), 
               use_apply(X),
               my_loop(X),
               my_sapply(X),
               times = 100)

结果：

Unit: milliseconds
         expr       min        lq      mean    median        uq       max neval
  for_loop(X) 122.45210 145.67150 179.84469 177.63446 199.10468 460.73182   100
 use_apply(X) 141.99250 169.11596 198.82019 198.11953 223.50906 296.94566   100
   my_loop(X)  10.38776  11.61263  16.47609  14.24066  19.07957  58.50008   100
 my_sapply(X)  13.21431  15.32081  23.23124  18.39573  26.08099 222.57685   100

所以快了10倍多！

你也可以注意到你的机器比我的快得多:/

为什么在计算矩阵的所有 rows/columns 的范数时使用 for 循环比应用更快？

Why is using a for loop faster than apply, in computing the norm of all rows/columns of a matrix?

for-loop

r

apply

purrr