为什么采样矩阵行很慢？

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

乍一看有两种可能性,两者都在第265行的 R函数MatrixSubset中.

它们可能都不是.只是猜测.

它似乎在缓存中循环低效的方向.

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

你的例子有很多列(8,000).每次内部循环获取一个新列时,它需要获取RAM页面,该值将RAM中的值保存到缓存中(很可能是L2).下一次提取是一个不同的列,因此它不太可能重用已经在L2中的页面.A matrix在内部是一个巨大的连续向量:第1列全部后跟所有第2列,等等.页面提取相对昂贵.走向"错误"的方向会引发太多的页面提取.更多关于CPU缓存的信息.

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

2.开关()太深

开启类型发生在每个项目中matrix.即使a matrix是单一类型!所以这很浪费.即使使用跳转表优化开关,矩阵中的每个项目仍可能发生跳转表('可能'因为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完全相同.但是这可能会在将来发生变化,所以我的懒惰(由于代码膨胀)很可能在将来导致R中的错误,如果LOGICAL确实改变(比如说char而不是intRAM效率) ).

解决代码膨胀问题的一种可能选择,请注意所有真正发生的事情是移动内存.所以所有类型(STRSXP,VECSXP和EXPRSXP除外)都可以通过使用memcpy类型大小的单个双循环来完成.SET_STRING_ELT并且SET_VECTOR_ELT仍然必须用于维护这些对象的引用计数.所以这应该只需要重复3次双for循环来维护.或者,这个习语可以包裹成一个#define在R的其他部分完成的习语.

最后,在第一个边界检查循环中是否可以检测到传入的行或列中是否存在任何NA(非常见的情况是不请求第NA行或第NA列!).如果没有NA,则(ii == NA_INTEGER || jj == NA_INTEGER) ? :可以通过在外部提升该分支来保存最深的三元()(2000*8000调用该分支).但随着更复杂的重复代码的成本.但是,分支预测可能会在所有体系结构上可靠地发挥作用,我们不应该担心这一点.

data.table做到这些的memcpy伎俩和深支节省一些但不是所有的地方.它也开始逐列并行.但是在这种情况下还不是因为它是新的并且仍在推出(setkey非常相似且已经并行).主线程虽然(并非并行)逐个处理character和list列SET_STRING_ELT,SET_VECTOR_ELT因此在R中不是线程安全的.其他线程并行处理所有整数,实数,复数和原始列.然后它会像记忆io一样快.

我没有看到你在61MB上看到的差异,但通过增加10x到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