对于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](最高速缓存适用 "最近最少未使用" 高速缓存行替换策略).当我们回来扫描第二行时,我们必须A[2, 1]再次从RAM中取出.因此,行式扫描给出了高速缓存未命中率100%.相比之下,逐列扫描仅具有高速缓存未命中率1 / 16 = 6.25%.在这个例子中,我们将观察到列式扫描要快得多.

总之,使用10000 x 3矩阵,无论我们是按行还是按列扫描,我们都具有相同的缓存性能.我没有看到这rowWise比报告columnWise的中位时间更快microbenchmark.他们的执行时间可能不完全相同,但差异太小,不足以引起我们的关注.

对于5000 x 5000矩阵,rowWise比慢得多columnWise.

谢谢你的验证.

备注

我们应该确保最内层循环中的顺序内存访问的"黄金法则"是效率的一般准则.但是从狭义上说不明白.

事实上,如果你把三列A作为三个向量x,y和z,并考虑元素加法(即行的和A):, z[i] = x[i] + y[i]我们没有对所有三个向量的顺序访问吗？这不属于"黄金法则"吗？10000 x 3逐行扫描矩阵与顺序交替读取三个矢量没有区别.这非常有效.