r - 对于R中的矩阵，为什么逐列运算不比逐行运算快(应该如此)

Question

考虑以下函数，它们按行和列存储值。

 #include <Rcpp.h>
using namespace Rcpp;

const int m = 10000;
const int n = 3;

// [[Rcpp::export]]
SEXP rowWise() {
    SEXP A = Rf_allocMatrix(INTSXP, m, n);
    int* p = INTEGER(A);
    int i, j;
    for (i = 0; i < m; i++){
        for(j = 0; j < n; j++) {
            p[m * j + i] = j;
        }
    }
    return A;
}

// [[Rcpp::export]]
SEXP columnWise() {
  SEXP A = Rf_allocMatrix(INTSXP, n, m);
  int* p = INTEGER(A);
  int i, j;
  for(j = 0; j < m; j++) {
    for (i = 0; i < n; i++){
      p[n * j + i] = i;
    }
  }
  return A;
}


/*** R
library(microbenchmark)
gc()
microbenchmark(
  rowWise(),
  columnWise(),
  times = 1000
)
*/

上面的代码产生

Unit: microseconds
         expr    min     lq     mean  median      uq       max neval
    rowWise() 12.524 18.631 64.24991 20.4540 24.8385 10894.353  1000
 columnWise() 11.803 19.434 40.08047 20.9005 24.1585  8590.663  1000

按行分配值比按列分配值更快(如果不是更慢)，这与我所相信的是相反的。

但是，它确实神奇地取决于 m和 n的值。所以我想我的问题是:为什么 columnWise的 速度不比rowWise快很多？

Answer 1

矩阵的尺寸(形状)会产生影响。

当我们对10000 x 3整数矩阵A进行逐行扫描时，我们仍然可以有效地进行缓存。为了简化说明，我假设A的每一列都与高速缓存行对齐。

--------------------------------------
A[1, 1] A[1, 2] A[1, 3]        M  M  M
A[2, 1] A[2, 2] A[2, 3]        H  H  H
   .        .       .          .  .  .
   .        .       .          .  .  .
A[16,1] A[16,2] A[16,3]        H  H  H
--------------------------------------
A[17,1] A[17,2] A[17,3]        M  M  M
A[18,1] A[18,2] A[18,3]        H  H  H
   .        .       .          .  .  .
   .        .       .          .  .  .
A[32,1] A[32,2] A[32,3]        H  H  H
--------------------------------------
A[33,1] A[33,2] A[33,3]        M  M  M
A[34,1] A[34,2] A[34,3]        H  H  H
   .        .       .          .  .  .
   .        .       .          .  .  .

64位高速缓存行可以容纳16个整数。当我们访问 A[1, 1]时，将填充完整的缓存行，即 A[1, 1]到 A[16, 1]均已加载到缓存中。当我们扫描一行 A[1, 1], A[1, 2], A[1, 3]时，一个 16 x 3矩阵被加载到缓存中，它比缓存容量(32 KB)小得多。虽然第一行中的每个元素都有一个缓存未命中(M)，但是当我们开始扫描第二行时，每个元素都有一个缓存命中(H)。因此，我们有一个周期性的模式:

[3 Misses] -> [45 Hits] -> [3 Misses] -> [45 Hits] -> ...

也就是说，我们的缓存未命中率平均为 3 / 48 = 1 / 16 = 6.25%。实际上，如果我们按列扫描 A，则它等于高速缓存未命中率，在这种情况下，我们具有以下周期性模式:

[1 Miss] -> [15 Hits] -> [1 Miss] -> [15 Hits] -> ...

尝试使用 5000 x 5000矩阵。在那种情况下，在读取第一行之后，将 16 x 5000元素提取到缓存中，但它比缓存容量大得多，因此缓存逐出将 A[1, 1]踢到了 A[16, 1](大多数缓存应用 "least recently unused"缓存行替换策略)。回到扫描第二行时，我们必须再次从RAM中获取 A[2, 1]。因此，按行扫描给出的缓存未命中率为 100%。相反，逐列扫描仅具有 1 / 16 = 6.25%的缓存未命中率。在此示例中，我们将观察到逐列扫描要快得多。

总而言之，使用 10000 x 3矩阵，无论按行还是按列进行扫描，我们都具有相同的缓存性能。从 rowWise报告的 中位数时间中，我看不出 columnWise比 microbenchmark快。它们的执行时间可能并不完全相等，但是差异太小而无法引起我们的关注。

For a 5000 x 5000 matrix, rowWise is much slower than columnWise.

感谢您的验证。

评论

我们应该确保在最内层循环中顺序访问内存的“黄金法则”是提高效率的一般准则。但不要狭义地理解它。

实际上，如果将 A的三列视为三个向量 x， y， z，并考虑逐元素加法(即 A的按行总和): z[i] = x[i] + y[i]，我们是否无法对所有对象进行顺序访问三个向量？这不属于“黄金法则”吗？逐行扫描 10000 x 3矩阵与依次交替读取三个向量没有区别。这是非常有效的。