r - 为什么采样矩阵行很慢？

Question

我试着做一些引导并计算 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

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

Answer 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 ...

如您所见，这种方式有更多的代码，因为相同的 for循环必须在 switch() 中一遍又一遍地重复案件。代码可读性和健壮性的原因可能是原始代码为何如此:在 R 的实现中出现拼写错误的可能性较小。这已经证明了，因为我懒惰地没有专门为 LOGICAL 实现 LGLSXP 案例。我知道 LOGICAL 与当前在基础 R 中的 INTEGER 完全相同。但这可能会在 future 发生变化，所以如果 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