c++ - R vs Rcpp vs Armadillo 中矩阵 rowSums() 与 colSums() 的效率

Question

背景

来自 R 编程，我正在使用 扩展为 C/C++ 形式的编译代码。 Rcpp .作为关于循环交换效果的练习(一般来说只是 C/C++)，我实现了与 R 等效的 rowSums()和 colSums()矩阵函数 Rcpp (我知道这些存在于 Rcpp 糖和 Armadillo 中——这只是一个练习)。

问题

我有我的 C++ 实现 rowSums()和 colSums()连同 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 矩阵以列为主存储，因此外循环中的列应该是更有效的方法。这就是我最初测试的内容。)

在对这些进行基准测试时，我遇到了一些我没想到的事情:行总和和列总和之间存在明显的性能差异(请参阅下面的基准):

使用内置 R 函数，colSums()大约是 rowSums() 的两倍.

使用我自己的 Rcpp 和糖版本，这是相反的:它是 rowSums()这大约是 colSums() 的两倍.

最后，添加 Armadillo 实现，操作大致相同(col sum 可能会快一点，因为我也希望它们在 R 中)。

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

作为次要兴趣，我也很好奇为什么在 R 实现中颠倒了相同的差异。 (我浏览了 the C source ，但无法真正找出显着差异。)(第三， Armadillo 如何获得相同的性能？)

基准

我在 10,000 x 10,000 上测试了这两个函数的所有 4 个实现对称矩阵:

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++ 源文件 matsums.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 接口(interface)编写了相同的函数:来源是 here .我compiled the functions与 R CMD SHLIB ，并再次测试:行总和比列总和更快的相同现象持续存在( benchmarks )。然后我也看了the disassembly with objdump ，但在我看来(由于我对 asm 的理解非常有限)，编译器并没有真正优化主循环体( rows 、 cols )离 C 代码更远吗？

Answer 1

首先，让我在我的笔记本电脑上显示计时统计信息。我使用 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 和 rowSums 以击败 Armadillo 。请注意，我编写 C，使用 R 的旧 C 接口(interface)并使用 R CMD SHLIB 进行编译。
colSums 和 rowSums 之间有什么实质性的区别吗？
如果我们有一个远大于 CPU 缓存容量的 n x n 矩阵，则 colSums 将 n x n 数据从 RAM 加载到缓存，但 rowSums 加载的数据是 2 x n x n 的两倍。
想想保存总和的结果 vector :这个 length-n vector 从 RAM 加载到缓存中的次数是多少？对于 colSums ，它只加载一次，但对于 rowSums ，它被加载 n 次。每次向其中添加矩阵列时，它都会加载到缓存中，但由于太大而被逐出。
对于大型 n :

colSums 导致 n x n + n 数据从 RAM 加载到缓存；

rowSums 导致 n x n + n x n 数据从 RAM 加载到缓存。

换句话说， rowSums 理论上的内存效率较低，并且可能会更慢。

如何提高 colSums 的性能？
由于 RAM 和缓存之间的数据流很容易优化，唯一的改进是循环展开。将内部循环(求和循环)展开深度为 2 就足够了，我们将看到 2 倍的提升。
循环展开工作是因为它启用了 CPU 的指令管道。如果我们每次迭代只做一个加法，流水线是不可能的；通过两个添加，这个指令级并行开始工作。我们也可以将循环展开深度为 4，但我的经验是深度 2 展开足以从循环展开中获得大部分好处。
如何提高 rowSums 的性能？
数据流的优化是第一步。我们需要首先进行缓存阻塞以减少从 2 x n x n 到 n x n 的数据传输。
将这个n x n矩阵切成若干行块:每个块都是2040 x n(最后一个块可能更小)，然后逐块应用普通的rowSums块。对于每个块，累加器 vector 的长度为 2040，大约是 32KB CPU 缓存可以容纳的长度的一半。对于添加到此累加器 vector 的矩阵列，另一半被反转。这样，累加器 vector 可以保存在缓存中，直到处理完该块中的所有矩阵列。因此，累加器 vector 只加载到缓存中一次，因此整体内存性能与 colSums 一样好。
现在我们可以进一步对每个块中的 rowSums 应用循环展开。将外循环和内循环展开深度为 2，我们将看到提升。一旦外循环展开，块大小应该减少到 1360，因为现在我们需要缓存中的空间来保存每个外循环迭代的三个长度为 1360 的 vector 。

C 代码:让我们打败 Armadillo
编写带有循环展开的代码可能是一项令人讨厌的工作，因为我们现在需要为一个函数编写几个不同的版本。
对于 colSums ，我们需要两个版本:

colSums_1x1 : 内循环和外循环都以深度1展开，即这是一个没有循环展开的版本；

colSums_2x1 :没有展开外循环，而展开深度为 2 的内循环。

对于 rowSums 我们最多可以有四个版本， 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)

模板函数为具有 sum <- colSums(A[1:nr, 1:nc])(A 的前导维度)行的矩阵 A 计算(在 R 语法中)LDA。参数 s 是对内循环展开的版本控制。模板函数乍一看很可怕，因为它包含许多 if 。但是，它被声明为 static inline 。如果通过将已知常量 1 或 2 传递给 s 来调用它，则优化编译器能够在编译时评估这些 if，消除无法访问的代码并删除“设置但未使用”变量(已初始化的寄存器变量) ，已修改但未写回 RAM)。
宏用于函数声明。接受一个常量 s 和一个函数名，它生成一个具有所需循环展开版本的函数。
以下是 rowSums 。

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

请注意循环展开如何加速 colSums 和 rowSums 。通过全面优化(“col2”和“row4”)，我们击败了 Armadillo (请参阅本答案开头的时序表)。
在这种情况下，行分块策略并没有明显地产生好处。让我们尝试一个包含数百万行的矩阵。

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

在这种情况下，我们观察到缓存阻塞的 yield 。

最后的想法
基本上这个答案已经解决了所有问题，除了以下问题:

为什么 R 的 rowSums 效率低于应有的效率。

为什么没有任何优化，rowSums ("row1") 比 colSums ("col1") 快。

同样，我无法解释第一个，实际上我不在乎，因为我们可以轻松编写比 R 的内置版本更快的版本。
第二个绝对值得追求。我在我们的讨论室中复制了我的评论以作记录。

This issue is down to this: "why adding up a single vector is slower than adding two vectors element-wise?"

I see similar phenomenon from time to time. The first time I encountered this strange behavior was when I, a few years ago, coded my own matrix-matrix multiplication. I found that DAXPY is faster than DDOT.

DAXPY does this: y += a * x, where x and y are vectors and a is a scalar; DDOT does this: a += x * y.

Given than DDOT is a reduction operation I expect that it is faster than DAXPY. But no, DAXPY is faster.

However, as soon as I unroll the loop in the triple loop-nest of the matrix-multiplication, DDOT is much faster than DAXPY.

A very similar thing happens to your issue. A reduction operation: a = x[1] + x[2] + ... + x[n] is slower than element-wise add: y[i] += x[i]. But once loop unrolling is done, the advantage of the latter is lost.

I am not sure whether the following explanation is true as I have no evidence.

The reduction operation has a dependency chain so the computation is strictly serial; on the other hand, element-wise operation has no dependency chain, so that CPU may do better with it.

As soon as we unroll the loop, each iteration has more arithmetics to do and CPU can do pipelining in both cases. The true advantage of the reduction operation can then be observed.

回复 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

我没有看到 rowSums2 和 colSums2 的性能优势。