为什么采样矩阵行很慢？答案

【问题标题】：Why sampling matrix row is very slow?为什么采样矩阵行很慢？
【发布时间】：2017-07-25 20:46:51
【问题描述】：

我尝试做一些引导和计算colMeans，自然我选择了矩阵来存储数据，但是采样速度很慢：

m[sample(n,replace=TRUE),]

原来data.table 是最快的。

require(microbenchmark)
require(data.table)
n = 2000
nc = 8000
m = matrix(1:(n*nc) ,nrow = n)
DF = as.data.frame(m)
DT = as.data.table(m)

s=sample(n, replace=TRUE)
microbenchmark(m[s,], DF[s,],DT[s,])

# Unit: milliseconds
    # expr      min       lq     mean   median       uq      max neval
  # m[s, ] 371.9271 402.3542 421.7907 420.8446 437.8251 506.1788   100
 # DF[s, ] 182.3189 199.0865 218.0746 213.9451 231.1518 409.8625   100
 # DT[s, ] 129.8225 139.1977 156.9506 150.4321 164.3104 254.2048   100

为什么采样矩阵比其他两个慢很多？

【问题讨论】：

@dww 这在我的系统上速度较慢（如我所料）。
请注意，此结果取决于m 的大小，尤其是列数。找出导致这种情况的原因需要分析internal C code。由于矩阵子集不应该更慢，你应该用 R 的开发版本确认这些时间，然后在 R-devel 邮件列表中提出这个问题。
请注意，在 m 的情况下，您正在用 length == nrow(m) * ncol(m) 置换 1 个向量（因为“矩阵”存储为具有“dim”属性的无量纲对象），而在DF/DT 您正在独立置换length == nrow(m) 的ncol(m) 向量（因为“data.frame”是向量的“列表”）。对于您的用例，我相信，一种有效的方法是将您的数据存储为tDF = as.data.frame(t(m))，因为这样可以避免大多数不必要的（深度）复制——microbenchmark(m[s, ], DF[s, ], tDF[, s], times = 50)（all.equal(colMeans(m[s, ]), rowMeans(tDF[, s]))）

标签： r matrix data.table

【解决方案1】：

乍一看有两种可能性，都在 R 的函数 MatrixSubset on line 265 中。

可能两者都不是。只是猜测。

1。它似乎朝着缓存低效的方向循环。

for (i = 0; i < nrs; i++) {    // rows
  ...
  for (j = 0; j < ncs; j++) {  // columns
    ...

您的示例有很多列 (8,000)。每次内部循环获取一个新列时，它都需要从 RAM 中获取保存该值的 RAM 页到缓存（很可能是 L2）。下一次提取是不同的列，因此不太可能重用已经在 L2 中的页面。 matrix 在内部是一个巨大的连续向量：第 1 列的所有内容，然后是第 2 列的所有内容，等等。页面获取相对昂贵。朝着“错误”的方向前进会导致过多的页面获取。更多关于 CPU 缓存here.

一个好的编译器应该自动执行Loop interchange，例如gcc -floop-interchange，它默认是开启的。更多here。由于 for 循环内部的复杂性，这种优化可能不会在这种情况下发生；也许在这种情况下是 switch 语句。或者，您在操作系统上使用的 R 版本可能没有使用带有该选项的编译器编译，或者没有打开。

2。 switch() 太深了

matrix 中的每个项目都在切换类型。即使 matrix 是单一类型！所以这是浪费。即使开关是optimized with a jump table，矩阵中的每个项目都可能仍在发生跳转表（“可能”，因为 CPU 可能会预测开关）。由于您的示例 matrix 很小，只有 61MB，我更倾向于这是罪魁祸首，而不是走错方向。

针对上述两者的建议修复（未经测试）

// Check the row numbers once up front rather than 8,000 times.
// This is a contiguous sweep and therefore almost instant
// Declare variables i and ii locally for safety and maximum compiler optimizations
for (int i = 0; i < nrs; i++) {
  int ii = INTEGER(sr)[i];
  if (ii != NA_INTEGER && (ii < 1 || ii > nr))
    errorcall(call, R_MSG_subs_o_b);
}

// Check the column numbers up front once rather than 2,000 times
for (int j = 0; j < ncs; j++) {
  int jj = INTEGER(sc)[j];
  if (jj != NA_INTEGER && (jj < 1 || jj > nc))
    errorcall(call, R_MSG_subs_o_b);
}

// Now switch once on type rather than 8,000 * 2,000 times
// Loop column-by-column not row-by-row

int resi=0;  // contiguous write to result (for page efficiency)
int ii, jj;  // the current row and column, bounds checked above
switch (TYPEOF(x)) {
  case LGLSXP:  // the INTSXP will work for LGLSXP too, currently
  case INTSXP:
    for (int j=0; j<ncs; j++) {  // column-by-column
      jj = INTEGER(sc)[j];
      for (int i=0; i<nrs; i++) {  // within-this-column
        ii = INTEGER(sr)[i];
        INTEGER(result)[resi++] = (ii == NA_INTEGER || jj == NA_INTEGER) ? NA_INTEGER : INTEGER(x)[ii + jj * nr];
      }
    }
    break;
  case REALSXP:
    for (int j=0; j<ncs; j++) {
      jj = INTEGER(sc)[j];
      for (int i=0; i<nrs; i++) {
        ii = INTEGER(sr)[i];
        REAL(result)[resi++] = (ii == NA_INTEGER || jj == NA_INTEGER) ? NA_REAL : REAL(x)[ii + jj * nr];
      }
    }
    break;
  case ...

如您所见，这种方式的代码更多，因为在 switch() 案例中必须一遍又一遍地重复相同的 for 循环。代码可读性和健壮性的原因可能是为什么原始代码是这样的：在 R 的实现中出现拼写错误的可能性较小。这已经证明了，因为我懒得没有专门为 LOGICAL 实现 LGLSXP 案例。我知道 LOGICAL 与当前在基础 R 中的 INTEGER 完全相同。但这可能会在未来发生变化，所以如果 LOGICAL 确实发生变化（比如char 而不是int RAM 效率）。

解决代码膨胀问题的一个可能选项，请注意真正发生的所有事情都是移动内存。因此，所有类型（STRSXP、VECSXP 和 EXPRSXP 除外）都可以通过使用 memcpy 和类型大小的单个双循环来完成。 SET_STRING_ELT 和 SET_VECTOR_ELT 仍然必须用于维护这些对象的引用计数。所以这应该只是要维护的双 for 循环的 3 次重复。或者，可以将该成语包装到 #define 中，这在 R 的其他部分完成。

最后，在第一个边界检查循环中可以检测到传入的行或列中是否有任何 NA（不请求第 NA 行或第 NA 列的非常常见的情况！）。如果没有 NA，则可以通过在外部提升该分支来保存最深的三元 ((ii == NA_INTEGER || jj == NA_INTEGER) ? :)（对该分支的 2000 * 8000 次调用）。但以更复杂的重复代码为代价。不过，branch prediction 可能会在所有架构中可靠地发挥作用，我们不应该担心这一点。

data.table 在某些但不是所有地方都执行memcpy 技巧和深度分支保存。它还开始逐列并行地进行子集化。但在这种情况下还不是因为它是新的并且仍在推出（setkey 非常相似并且已经是并行的）。主线程一一处理character 和list 列（不是并行处理），因为SET_STRING_ELT 和SET_VECTOR_ELT 在R 中不是线程安全的。其他线程处理所有整数、实数、复数和并行的原始列。然后它会尽可能快地运行内存 io。

我并没有真正看到您在 61MB 上看到的差异，但是通过将列数增加 10 倍到 80,000 来扩大到（仍然很小）610MB，我确实看到了差异。

n = 2000
nc = 8000    # same size as your example (61MB), on my laptop
microbenchmark(m[s,], DF[s,],DT[s,])
Unit: milliseconds
    expr       min        lq      mean    median        uq      max neval
  m[s, ] 108.75182 112.11678 118.60111 114.58090 120.07952 168.6079   100
 DF[s, ] 100.95019 105.88253 116.04507 110.84693 118.08092 163.9666   100
 DT[s, ]  63.78959  69.07341  80.72039  72.69873  96.51802 136.2016   100

n = 2000
nc = 80000     # 10x bigger (610MB)
microbenchmark(m[s,], DF[s,],DT[s,])
Unit: milliseconds
    expr       min        lq      mean    median        uq      max neval
  m[s, ] 1990.3343 2010.1759 2055.9847 2032.9506 2057.2498 2733.278   100
 DF[s, ] 1083.0373 1212.6633 1265.5346 1234.1558 1300.7502 2105.177   100
 DT[s, ]  698.1295  830.3428  865.5918  862.5773  907.7225 1053.393   100

不过，我有 128MB 的 L4 缓存。我猜你有更少的缓存。整个 61MB 都适合我的 L4 缓存，所以我并没有真正注意到这种大小的缓存效率低下。

$ lscpu
Architecture:          x86_64
CPU op-mode(s):        32-bit, 64-bit
Byte Order:            Little Endian
CPU(s):                8
On-line CPU(s) list:   0-7
Thread(s) per core:    2
Core(s) per socket:    4
Socket(s):             1
NUMA node(s):          1
Vendor ID:             GenuineIntel
CPU family:            6
Model:                 70
Model name:            Intel(R) Core(TM) i7-4980HQ CPU @ 2.80GHz
Stepping:              1
CPU MHz:               3345.343
CPU max MHz:           4000.0000
CPU min MHz:           800.0000
BogoMIPS:              5587.63
Virtualization:        VT-x
L1d cache:             32K
L1i cache:             32K
L2 cache:              256K
L3 cache:              6144K
L4 cache:              131072K
NUMA node0 CPU(s):     0-7

【讨论】：