代码之家  ›  专栏  ›  技术社区  ›  Mikko Marttila

R vs Rcpp vs犰狳中矩阵行和与列和的效率

  •  7
  • Mikko Marttila  · 技术社区  · 6 年前

    背景

    从R编程开始,我正在以C/C++的形式扩展到编译代码的过程中。 Rcpp公司 . 作为一个手练习循环交换的效果(一般是C/C++),我实现了R的等价物。 rowSums() colSums() 矩阵的函数 Rcpp公司 (我知道它们存在于RCPP糖和犰狳中,这只是一个练习)。

    问题

    我的C++实现 行和() 列和() 随着 Rcpp sugar arma::sum() 版本 this matsums.cpp file . 我的只是这样的简单循环:

    NumericVector Cpp_colSums(const NumericMatrix& x) {
      int nr = x.nrow(), nc = x.ncol();
      NumericVector ans(nc);
      for (int j = 0; j < nc; j++) {
        double sum = 0.0;
        for (int i = 0; i < nr; i++) {
          sum += x(i, j);
        }
        ans[j] = sum;
      }
      return ans;
    }
    
    NumericVector Cpp_rowSums(const NumericMatrix& x) {
      int nr = x.nrow(), nc = x.ncol();
      NumericVector ans(nr);
      for (int j = 0; j < nc; j++) {
        for (int i = 0; i < nr; i++) {
          ans[i] += x(i, j);
        }
      }
      return ans;
    }
    

    ( R矩阵存储列major,因此外部循环中的列应该是更有效的方法。这就是我最初测试的。 )

    在运行这些基准测试时,我遇到了一些出乎意料的事情:行和列和之间存在明显的性能差异(参见下面的基准测试):

    1. 使用内置的R函数, 列和() 大约是 行和() .
    2. 有了我自己的Rcpp和sugar版本,这是相反的:它是 行和() 速度大约是 列和() .
    3. 最后,加上犰狳的实现,操作大致相等(col sum可能会快一点,正如我在R中所期望的那样)。

    所以我的首要问题是: 为什么 Cpp_rowSums() 明显快于 Cpp_colSums() ?

    作为第二个兴趣,我也很好奇为什么相同的差异在R实现中是相反的。(我浏览了一下 the C source (第三,犰狳为什么会有同样的表现?)

    基准

    我在 10,000 x 10,000 对称矩阵:

    Rcpp::sourceCpp("matsums.cpp")
    
    set.seed(92136)
    n <- 1e4 # build n x n test matrix
    x <- matrix(rnorm(n), 1, n)
    x <- crossprod(x, x) # symmetric
    
    bench::mark(
      rowSums(x),
      colSums(x),
      Cpp_rowSums(x),
      Cpp_colSums(x),
      Sugar_rowSums(x),
      Sugar_colSums(x),
      Arma_rowSums(x),
      Arma_colSums(x)
    )[, 1:7]
    
    #> # A tibble: 8 x 7
    #>   expression            min     mean   median      max `itr/sec` mem_alloc
    #>   <chr>            <bch:tm> <bch:tm> <bch:tm> <bch:tm>     <dbl> <bch:byt>
    #> 1 rowSums(x)        192.2ms  207.9ms  194.6ms  236.9ms      4.81    78.2KB
    #> 2 colSums(x)         93.4ms   97.2ms   96.5ms  101.3ms     10.3     78.2KB
    #> 3 Cpp_rowSums(x)     73.5ms   76.3ms     76ms   80.4ms     13.1     80.7KB
    #> 4 Cpp_colSums(x)    126.5ms  127.6ms  126.8ms  130.3ms      7.84    80.7KB
    #> 5 Sugar_rowSums(x)   73.9ms   75.6ms   74.3ms   79.4ms     13.2     80.7KB
    #> 6 Sugar_colSums(x)  124.2ms  125.8ms  125.6ms  127.9ms      7.95    80.7KB
    #> 7 Arma_rowSums(x)    73.2ms   74.7ms   73.9ms   79.3ms     13.4     80.7KB
    #> 8 Arma_colSums(x)    62.8ms   64.4ms   63.7ms   69.6ms     15.5     80.7KB
    

    (再次,您可以找到C++源文件。 垫子.cpp here .)

    平台:

    > sessioninfo::platform_info()
     setting  value                       
     version  R version 3.5.1 (2018-07-02)
     os       Windows >= 8 x64            
     system   x86_64, mingw32             
     ui       RStudio                     
     language (EN)                        
     collate  English_United States.1252  
     tz       Europe/Helsinki             
     date     2018-08-09  
    

    更新

    进一步研究,我还使用R的传统C接口编写了相同的函数:源代码是 here . 我 compiled the functions 具有 R CMD SHLIB ,并再次测试:行和比列和快的相同现象仍然存在( benchmarks ). 我又看了看 the disassembly with objdump ,但在我看来(对asm的理解非常有限),编译器并没有真正优化主循环体( rows , cols )离C码还有远吗?

    2 回复  |  直到 6 年前
        1
  •  11
  •   Zheyuan Li    6 年前

    首先,让我在笔记本电脑上显示计时统计数据。我使用的5000 x 5000矩阵足以进行基准测试,并且 microbenchmark 包用于100次评估。

    Unit: milliseconds
                 expr       min        lq      mean    median        uq       max
           colSums(x)  71.40671  71.64510  71.80394  71.72543  71.80773  75.07696
       Cpp_colSums(x)  71.29413  71.42409  71.65525  71.48933  71.56241  77.53056
     Sugar_colSums(x)  73.05281  73.19658  73.38979  73.25619  73.31406  76.93369
      Arma_colSums(x)  39.08791  39.34789  39.57979  39.43080  39.60657  41.70158
           rowSums(x) 177.33477 187.37805 187.57976 187.49469 187.73155 194.32120
       Cpp_rowSums(x)  54.00498  54.37984  54.70358  54.49165  54.73224  64.16104
     Sugar_rowSums(x)  54.17001  54.38420  54.73654  54.56275  54.75695  61.80466
      Arma_rowSums(x)  49.54407  49.77677  50.13739  49.90375  50.06791  58.29755
    

    R内核中的C代码并不总是比我们自己编写的代码好。那个 Cpp_rowSums rowSums 显示这个。我觉得自己没有能力解释为什么R的版本比应该的慢。我将专注于: 如何进一步优化自身 colSums 行和 打败犰狳 . 请注意,我编写C,使用R的旧C接口并使用 R CMD SHLIB .


    有什么实质性的区别吗 科尔苏姆 行和 ?

    如果我们有一个 n x n 比CPU缓存容量大得多的矩阵, 科尔苏姆 荷载 n×n 从RAM到缓存的数据,但是 行和 装载量是原来的两倍。, 2 x n x n .

    想想得到的向量,它包含了和:这个长度的多少倍- n 向量从RAM加载到缓存中?为了 科尔苏姆 ,它只加载一次,但是 行和 ,它已加载 n个 时代。每次向它添加矩阵列时,它都会加载到缓存中,但由于太大而被逐出。

    对于一个大的 n个 :

    • 科尔苏姆 原因 n x n + n 从RAM到cache的数据加载;
    • 行和 原因 n x n + n x n 数据从RAM加载到缓存。

    换句话说, 行和 从理论上讲,它的记忆效率较低,而且速度可能较慢。


    如何提高 科尔苏姆 ?

    由于RAM和cache之间的数据流很容易达到最优,因此惟一的改进是循环展开。将内循环(总和循环)展开2就足够了,我们将看到2倍的提升。

    循环展开工作,因为它启用CPU的指令管道。如果每次迭代只进行一次加法运算,那么就不可能进行流水线操作;通过两次加法运算,指令级并行开始工作。我们还可以将循环展开4的深度,但我的经验是,深度2展开足以获得循环展开的大部分好处。

    如何提高 行和 ?

    优化数据流是第一步。我们需要首先执行缓存阻塞以减少 2 x n x n 下降到 n×n .

    切掉这个 n×n 矩阵成若干行块:每个 2040 x n (最后一块可能更小),然后应用普通 行和 一块一块。对于每个块,累加器向量的长度为2040,大约是32KB CPU缓存容量的一半。对于添加到此累加器向量的矩阵列,另一半是相反的。这样,累加器向量就可以保存在缓存中,直到处理该块中的所有矩阵列。因此,累加器向量只加载到缓存中一次,因此总体内存性能与 科尔苏姆 .

    现在我们可以进一步应用循环展开 行和 在每一块中。展开外环和内环的深度为2,我们将看到一个提升。一旦外部循环被展开,块大小应该减少到1360,因为现在我们需要缓存中的空间来保存三个长度-每个外部循环迭代1360个向量。


    C代码:让我们打败犰狳

    使用循环展开来编写代码可能是一件麻烦的工作,因为我们现在需要为一个函数编写几个不同的版本。

    为了 科尔苏姆 ,我们需要两个版本:

    • colSums_1x1 :内部和外部循环都以深度1展开,即这是一个没有循环展开的版本;
    • colSums_2x1 :不展开外部循环,而展开内部循环时深度为2。

    为了 行和 我们最多可以有四个版本, rowSums_sxt ,其中 s = 1 or 2 是内环和 t = 1 or 2 是外环的展开深度。

    如果我们一个接一个地编写每个版本,代码编写可能会非常乏味。在经历了多年的挫折之后,我开发了一个使用内联模板函数和宏的“自动代码/版本生成”技巧。

    #include <stdlib.h>
    #include <Rinternals.h>
    
    static inline void colSums_template_sx1 (size_t s,
                                             double *A, size_t LDA,
                                             size_t nr, size_t nc,
                                             double *sum) {
    
      size_t nrc = nr % s, i;
      double *A_end = A + LDA * nc, a0, a1;
    
      for (; A < A_end; A += LDA) {
        a0 = 0.0; a1 = 0.0;  // accumulator register variables
        if (nrc > 0) a0 = A[0];  // is there a "fractional loop"?
        for (i = nrc; i < nr; i += s) {  // main loop of depth-s
          a0 += A[i];  // 1st iteration
          if (s > 1) a1 += A[i + 1];  // 2nd iteration
          }
        if (s > 1) a0 += a1;  // combine two accumulators
        *sum++ = a0;  // write-back
        }
    
      }
    
    #define macro_define_colSums(s, colSums_sx1) \
    SEXP colSums_sx1 (SEXP matA) { \
      double *A = REAL(matA); \
      size_t nrow_A = (size_t)nrows(matA); \
      size_t ncol_A = (size_t)ncols(matA); \
      SEXP result = PROTECT(allocVector(REALSXP, ncols(matA))); \
      double *sum = REAL(result); \
      colSums_template_sx1(s, A, nrow_A, nrow_A, ncol_A, sum); \
      UNPROTECT(1); \
      return result; \
      }
    
    macro_define_colSums(1, colSums_1x1)
    macro_define_colSums(2, colSums_2x1)
    

    template函数计算(在R语法中) sum <- colSums(A[1:nr, 1:nc]) 对于矩阵 A 具有 LDA (A)行的前导维度。参数 s 是内部循环展开的版本控制。模板函数包含许多 if . 但是,它被宣布 static inline . 如果通过将已知常数1或2传递给 s公司 ,优化编译器能够计算 如果 在编译时,消除不可访问的代码,删除“set but not used”变量(注册初始化、修改但不写回RAM的变量)。

    宏用于函数声明。接受常数 s公司 以及函数名,它生成具有所需循环展开版本的函数。

    以下为 行和 .

    static inline void rowSums_template_sxt (size_t s, size_t t,
                                             double *A, size_t LDA,
                                             size_t nr, size_t nc,
                                             double *sum) {
    
      size_t ncr = nc % t, nrr = nr % s, i;
      double *A_end = A + LDA * nc, *B;
      double a0, a1;
    
      for (i = 0; i < nr; i++) sum[i] = 0.0;  // necessary initialization
    
      if (ncr > 0) {  // is there a "fractional loop" for the outer loop?
        if (nrr > 0) sum[0] += A[0];  // is there a "fractional loop" for the inner loop?
        for (i = nrr; i < nr; i += s) {  // main inner loop with depth-s
          sum[i] += A[i];
          if (s > 1) sum[i + 1] += A[i + 1];
          }
        A += LDA;
        }
    
      for (; A < A_end; A += t * LDA) {  // main outer loop with depth-t
        if (t > 1) B = A + LDA;
        if (nrr > 0) {  // is there a "fractional loop" for the inner loop?
          a0 = A[0]; if (t > 1) a0 += A[LDA];
          sum[0] += a0;
          }
        for(i = nrr; i < nr; i += s) {  // main inner loop with depth-s
          a0 = A[i]; if (t > 1) a0 += B[i];
          sum[i] += a0;
          if (s > 1) {
            a1 = A[i + 1]; if (t > 1) a1 += B[i + 1];
            sum[i + 1] += a1;
            }
          }
        }
    
      }
    
    #define macro_define_rowSums(s, t, rowSums_sxt) \
    SEXP rowSums_sxt (SEXP matA, SEXP chunk_size) { \
      double *A = REAL(matA); \
      size_t nrow_A = (size_t)nrows(matA); \
      size_t ncol_A = (size_t)ncols(matA); \
      SEXP result = PROTECT(allocVector(REALSXP, nrows(matA))); \
      double *sum = REAL(result); \
      size_t block_size = (size_t)asInteger(chunk_size); \
      size_t i, block_size_i; \
      if (block_size > nrow_A) block_size = nrow_A; \
      for (i = 0; i < nrow_A; i += block_size_i) { \
        block_size_i = nrow_A - i; if (block_size_i > block_size) block_size_i = block_size; \
        rowSums_template_sxt(s, t, A, nrow_A, block_size_i, ncol_A, sum); \
        A += block_size_i; sum += block_size_i; \
        } \
      UNPROTECT(1); \
      return result; \
      }
    
    macro_define_rowSums(1, 1, rowSums_1x1)
    macro_define_rowSums(1, 2, rowSums_1x2)
    macro_define_rowSums(2, 1, rowSums_2x1)
    macro_define_rowSums(2, 2, rowSums_2x2)
    

    注意,模板函数现在接受 s公司 t ,宏定义的函数已应用行分块。

    即使我在代码中留下了一些注释,代码可能仍然不容易理解,但是我不能花更多的时间来解释更详细的内容。

    要使用它们,请将它们复制并粘贴到一个名为“matSums.C”的C文件中,并使用 R CMD SHLIB -c matSums.c .

    对于R端,在“matSums.R”中定义以下函数。

    colSums_zheyuan <- function (A, s) {
      dyn.load("matSums.so")
      if (s == 1) result <- .Call("colSums_1x1", A)
      if (s == 2) result <- .Call("colSums_2x1", A)
      dyn.unload("matSums.so")
      result
      }
    
    rowSums_zheyuan <- function (A, chunk.size, s, t) {
      dyn.load("matSums.so")
      if (s == 1 && t == 1) result <- .Call("rowSums_1x1", A, as.integer(chunk.size))
      if (s == 2 && t == 1) result <- .Call("rowSums_2x1", A, as.integer(chunk.size))
      if (s == 1 && t == 2) result <- .Call("rowSums_1x2", A, as.integer(chunk.size))
      if (s == 2 && t == 2) result <- .Call("rowSums_2x2", A, as.integer(chunk.size))
      dyn.unload("matSums.so")
      result
      }
    

    现在让我们有一个基准,同样有一个5000 x 5000矩阵。

    A <- matrix(0, 5000, 5000)
    
    library(microbenchmark)
    source("matSums.R")
    
    microbenchmark("col0" = colSums(A),
                   "col1" = colSums_zheyuan(A, 1),
                   "col2" = colSums_zheyuan(A, 2),
                   "row0" = rowSums(A),
                   "row1" = rowSums_zheyuan(A, nrow(A), 1, 1),
                   "row2" = rowSums_zheyuan(A, 2040, 1, 1),
                   "row3" = rowSums_zheyuan(A, 1360, 1, 2),
                   "row4" = rowSums_zheyuan(A, 1360, 2, 2))
    

    在我的笔记本电脑上我得到:

    Unit: milliseconds
     expr       min        lq      mean    median        uq       max neval
     col0  65.33908  71.67229  71.87273  71.80829  71.89444 111.84177   100
     col1  67.16655  71.84840  72.01871  71.94065  72.05975  77.84291   100
     col2  35.05374  38.98260  39.33618  39.09121  39.17615  53.52847   100
     row0 159.48096 187.44225 185.53748 187.53091 187.67592 202.84827   100
     row1  49.65853  54.78769  54.78313  54.92278  55.08600  60.27789   100
     row2  49.42403  54.56469  55.00518  54.74746  55.06866  60.31065   100
     row3  37.43314  41.57365  41.58784  41.68814  41.81774  47.12690   100
     row4  34.73295  37.20092  38.51019  37.30809  37.44097  99.28327   100
    

    注意循环展开如何加快 科尔苏姆 行和 . 通过完全优化(“col2”和“row4”),我们击败了犰狳(见本答案开头的时间表)。

    在这种情况下,行分块策略并不能明显地产生效益。让我们试试有数百万行的矩阵。

    A <- matrix(0, 1e+7, 20)
    microbenchmark("row1" = rowSums_zheyuan(A, nrow(A), 1, 1),
                   "row2" = rowSums_zheyuan(A, 2040, 1, 1),
                   "row3" = rowSums_zheyuan(A, 1360, 1, 2),
                   "row4" = rowSums_zheyuan(A, 1360, 2, 2))
    

    我明白了

    Unit: milliseconds
     expr      min       lq     mean   median       uq      max neval
     row1 604.7202 607.0256 617.1687 607.8580 609.1728 720.1790   100
     row2 514.7488 515.9874 528.9795 516.5193 521.4870 636.0051   100
     row3 412.1884 413.8688 421.0790 414.8640 419.0537 525.7852   100
     row4 377.7918 379.1052 390.4230 379.9344 386.4379 476.9614   100
    

    在本例中,我们观察缓存阻塞带来的好处。


    最后的想法

    基本上,这个答案解决了所有问题,除了以下问题:

    • 为什么是R 行和 没有它应有的效率。
    • 为什么没有任何优化, 行和 (“row1”)比 科尔苏姆 (“col1”)。

    再说一遍,我不能解释第一个,实际上我不在乎,因为我们可以很容易地编写一个比R的内置版本快的版本。

    第二个绝对值得追求。我在我们的讨论室里把我的评论抄下来存档。

    这个问题可以归结为:“为什么添加一个向量要比添加两个向量慢?”

    我时常看到类似的现象。我第一次遇到这种奇怪的行为是,几年前,我编码我自己的矩阵乘法。我发现DAXPY比DDOT快。

    DAXPY这样做: y += a * x ,其中 x y 是向量和 a 是标量;DDOT执行以下操作: a += x * y .

    假设DDOT是一个还原操作,我希望它比DAXPY快。但不,达斯比更快。

    然而,只要我在矩阵乘法的三环嵌套中展开循环,DDOT就比DAXPY快得多。

    你的问题也发生了类似的事情。还原操作: a = x[1] + x[2] + ... + x[n] 比按元素添加慢: y[i] += x[i] . 但一旦完成循环展开,后者的优势就丧失了。

    我不确定下面的解释是否正确,因为我没有证据。

    归约操作有一个依赖链,因此计算是严格串行的;另一方面,元素级操作没有依赖链,因此CPU可以更好地使用它。

    一旦我们展开循环,每次迭代都有更多的算法要做,CPU在这两种情况下都可以做流水线。然后可以观察到还原操作的真正优势。


    答复 Jaap 关于使用 rowSums2 colSums2 matrixStats

    仍然使用上面的5000 x 5000示例。

    A <- matrix(0, 5000, 5000)
    
    library(microbenchmark)
    source("matSums.R")
    library(matrixStats)  ## NEW
    
    microbenchmark("col0" = base::colSums(A),
                   "col*" = matrixStats::colSums2(A),  ## NEW
                   "col1" = colSums_zheyuan(A, 1),
                   "col2" = colSums_zheyuan(A, 2),
                   "row0" = base::rowSums(A),
                   "row*" = matrixStats::rowSums2(A),  ## NEW
                   "row1" = rowSums_zheyuan(A, nrow(A), 1, 1),
                   "row2" = rowSums_zheyuan(A, 2040, 1, 1),
                   "row3" = rowSums_zheyuan(A, 1360, 1, 2),
                   "row4" = rowSums_zheyuan(A, 1360, 2, 2))
    
    Unit: milliseconds
     expr       min        lq      mean    median        uq       max neval
     col0  71.53841  71.72628  72.13527  71.81793  71.90575  78.39645   100
     col*  75.60527  75.87255  76.30752  75.98990  76.18090  87.07599   100
     col1  71.67098  71.86180  72.06846  71.93872  72.03739  77.87816   100
     col2  38.88565  39.03980  39.57232  39.08045  39.16790  51.39561   100
     row0 187.44744 187.58121 188.98930 187.67168 187.86314 206.37662   100
     row* 158.08639 158.26528 159.01561 158.34864 158.62187 174.05457   100
     row1  54.62389  54.81724  54.97211  54.92394  55.04690  56.33462   100
     row2  54.15409  54.44208  54.78769  54.59162  54.76073  60.92176   100
     row3  41.43393  41.63886  42.57511  41.73538  41.81844 111.94846   100
     row4  37.07175  37.25258  37.45033  37.34456  37.47387  43.14157   100
    

    我不认为 游泳池2 阴道炎2 .

        2
  •  1
  •   Jesper Juhl    6 年前

    “为什么Cpp_rowSums()比Cpp_colSums()快得多?”-当获取“row major”时,CPU预取器可以预测您正在做什么,并在需要之前将下一批您需要的数据从主存提取到CPU缓存。这将加快您对数据的访问。

    当你访问“column major”时,预取器很难预测你下一步需要什么,因此它不会提前高效地(如果有的话)将东西填充到缓存中,这会减慢你的速度。

    中央处理器 线性访问数据。如果你不做他们喜欢的事情,你就要付出缓存丢失和主内存访问延迟的代价。