首先，让我展示我的笔记本电脑的时序统计数据。我使用一个 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接口，并使用
R CMD SHLIB进行编译。
colSums和rowSums之间是否存在实质性差异？
如果我们有一个比CPU缓存容量大得多的n x n矩阵，colSums会从RAM加载n x n数据到缓存中，但rowSums会加载两倍的数据，即2 x n x n。
想一想保存总和的向量：这个长度为n的向量从RAM加载到缓存中的次数是多少？对于colSums，它只被加载一次，但对于rowSums，它被加载了n次。每次将矩阵列添加到其中时，它都会被加载到缓存中，但然后被驱逐，因为它太大了。
对于大的n:
colSums会从RAM加载n x n + n数据到缓存中;rowSums会从RAM加载n x n + n x n数据到缓存中。
换句话说，rowSums在理论上不太内存高效，并且可能会更慢。
如何改善colSums的性能？
由于RAM和缓存之间的数据流已经趋于最优，唯一的改进是循环展开。将内部循环（求和循环）展开2个深度就足够了，我们将会看到2倍的提升。
循环展开的作用在于启用CPU的指令管道。如果我们每次迭代只进行一次加法，那么无法进行流水线处理；而进行两次加法时，这种指令级并行性就开始起作用了。我们也可以将循环展开深度增加到4，但我的经验是深度为2的展开已足以获得大部分的循环展开好处。
如何改善rowSums的性能？
优化数据流是第一步。我们需要首先进行缓存块处理，以将数据传输从2 x n x n降至n x n。
将这个n x n矩阵分成多个行块：每个块大小为2040 x n（最后一个块可能较小），然后逐个应用普通的rowSums块。对于每个块，累加器向量的长度为2040，大约是32KB CPU缓存的一半。另一半则用于矩阵列添加到这个累加器向量中。这样，累加器向量就可以在缓存中保持，直到此块中的所有矩阵列都被处理完毕。因此，累加器向量只加载到缓存中一次，因此总体内存性能与colSums相当。
现在我们可以进一步对每个块中的
rowSums应用循环展开。将外部循环和内部循环都展开2层，我们会看到一个提升。一旦外部循环展开，块大小应该减小到1360，因为现在我们需要在缓存中保持每个外部循环迭代的三个长度为1360的向量的空间。
C代码：让Armadillo垮掉
使用循环展开编写代码可能是一项糟糕的工作，因为现在我们需要为函数编写多个不同版本。
对于colSums，我们需要两个版本：
colSums_1x1：内部和外部循环都展开深度为1，即没有循环展开的版本；colSums_2x1：没有外部循环展开，而内部循环深度为2。
对于rowSums，我们最多可以有四个版本，rowSums_sxt，其中s = 1或2是内部循环的展开深度，t = 1或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;
if (nrc > 0) a0 = A[0];
for (i = nrc; i < nr; i += s) {
a0 += A[i];
if (s > 1) a1 += A[i + 1];
}
if (s > 1) a0 += a1;
*sum++ = a0;
}
}
#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)
该模板函数计算矩阵A（具有LDA行）的列和，表达式为
sum <- colSums(A[1:nr, 1:nc])。参数
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;
if (ncr > 0) {
if (nrr > 0) sum[0] += A[0];
for (i = nrr; i < nr; i += s) {
sum[i] += A[i];
if (s > 1) sum[i + 1] += A[i + 1];
}
A += LDA;
}
for (; A < A_end; A += t * LDA) {
if (t > 1) B = A + LDA;
if (nrr > 0) {
a0 = A[0]; if (t > 1) a0 += A[LDA];
sum[0] += a0;
}
for(i = nrr; i < nr; i += 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
在这种情况下，我们观察到了缓存阻塞所带来的收益。
最后的想法
基本上，这个答案已经解决了所有问题，除了以下几点：
- 为什么R的
rowSums比应该更高效。
- 为什么没有任何优化，
rowSums（“row1”）比colSums（“col1”）更快。
再次强调，我无法解释第一个问题，实际上我不在意，因为我们可以很容易地编写比R内置版本更快的版本。
第二个问题值得追求。 我在我们的讨论室中复制了我的评论以备记录。
这个问题归结为：“为什么单个向量的加法比逐元素添加两个向量要慢？”我经常看到类似的现象。第一次遇到这种奇怪的行为是几年前，当时我编写了自己的矩阵乘法。我发现DAXPY比DDOT快。DAXPY执行以下操作：
y += a * x，其中
x和
y是向量，
a是标量；DDOT执行以下操作：
a += x * y。鉴于DDOT是一种约简操作，我预计它比DAXPY快。但不，DAXPY更快。然而，一旦我展开矩阵乘法的三重循环嵌套中的循环，DDOT比DAXPY快得多。你的问题也会发生非常相似的情况。一个约简操作：
a = x [1]+x [2]+...+x [n]比逐元素添加：
y [i] += x [i]慢。但是一旦展开循环，后者的优势就会丧失。我不确定以下解释是否正确，因为我没有证据。约简操作具有依赖链，因此计算是严格串行的；另一方面，逐元素操作没有依赖链，因此CPU可能更好地处理它。一旦我们展开循环，每次迭代都需要更多的算术运算，CPU可以在两种情况下进行流水线处理。然后可以观察到约简操作的真正优势。
回复Jaap关于使用matrixStats中的rowSums2和colSums2
仍然使用上面的5000 x 5000示例。
A <- matrix(0, 5000, 5000)
library(microbenchmark)
source("matSums.R")
library(matrixStats)
microbenchmark("col0" = base::colSums(A),
"col*" = matrixStats::colSums2(A),
"col1" = colSums_zheyuan(A, 1),
"col2" = colSums_zheyuan(A, 2),
"row0" = base::rowSums(A),
"row*" = matrixStats::rowSums2(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 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在性能上并没有优势。
