Ответ 1
Две возможности spring на первый взгляд, как в функции R MatrixSubset в строке 265.
Это может быть ни один из них. Просто гадать.
1. Кажется, что петля в кэше неэффективна.
for (i = 0; i < nrs; i++) { // rows
...
for (j = 0; j < ncs; j++) { // columns
...
В вашем примере много столбцов (8000). Каждый раз, когда внутренний цикл извлекает новый столбец, ему нужно получить страницу ОЗУ, содержащую это значение из ОЗУ в кеш (скорее всего, L2). Следующий выбор - это другой столбец, и поэтому он менее вероятно, сможет повторно использовать страницу уже в L2. A matrix является внутренним одним огромным смежным вектором: весь столбец 1, за которым следует весь столбец 2 и т.д. Извлечение страницы относительно дорого. Переход в "неправильное" направление приводит к слишком большому набору страниц. Подробнее о кэше CPU здесь.
Хороший компилятор должен выполнить Loop interchange автоматически, например gcc -floop-interchange, который включен по умолчанию. Подробнее здесь. В этом случае эта оптимизация может не произойти из-за сложности того, что внутри циклов for; возможно, в этом случае операторы switch. Или, возможно, версия R, которую вы используете на вашей ОС, не была скомпилирована с компилятором с этой опцией или не была включена.
2. Переключатель() слишком глубокий
Тип включения происходит в каждом элементе matrix. Даже если matrix - это единственный тип! Так что это расточительно. Даже если переключатель оптимизирован с помощью таблицы переходов, таблица перехода, вероятно, все еще происходит для каждого элемента в матрице (возможно, потому, что CPU может предсказать переключатель). Поскольку ваш пример matrix крошечный на 61 МБ, я склоняюсь больше к тому, чтобы быть виновником, а не идти в неправильном направлении.
Предлагаемое исправление для обоих выше (непроверенных)
// 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. Это уже продемонстрировало, потому что я ленился в том, что я не использовал приложение LGLSXP специально для LOGICAL. Я знаю, что LOGICAL точно так же, как INTEGER в настоящее время в базе R. Но это может измениться в будущем, поэтому моя лень (из-за раздувания кода) может привести к ошибке в R в будущем, если LOGICAL действительно изменится (скажем char, а не int для эффективности ОЗУ).
Один из возможных вариантов решения проблемы раздувания кода, обратите внимание, что все, что действительно происходит, перемещает память. Таким образом, все типы (кроме STRSXP, VECSXP и EXPRSXP) могут быть выполнены с помощью одного двойного цикла, используя memcpy с размером шрифта. SET_STRING_ELT и SET_VECTOR_ELT по-прежнему должны использоваться для поддержания ссылок на эти объекты. Таким образом, это должно быть всего лишь 3 повторения двойных циклов for для поддержки. Альтернативно, эта идиома может быть завернута в #define, которая выполняется в других частях R.
Наконец, есть ли какие-либо NA в строке или столбцах, переданных (очень обычный случай, чтобы не запрашивать строку NA'th или столбец NA'th!), можно обнаружить в первом цикле проверки границ. Если нет NA, то самый глубокий тройной ((ii == NA_INTEGER || jj == NA_INTEGER) ? :) (2000 * 8000 звонков на эту ветку) можно сохранить, подняв эту ветку за пределы. Но со стоимостью более сложного повторяющегося кода. Однако, возможно, предсказание ветвей будет надежно защищаться от всех архитектур, и мы не должны беспокоиться об этом.
data.table выполняет как трюк memcpy, так и глубину сохранения в некоторых, но не во всех местах. Он также начал подмножество параллельно, по столбцу. Но не в этом случае еще только потому, что он новый и все еще выкатившийся (setkey очень похож и уже параллелен). Главный поток обрабатывает столбцы character и list один за другим, хотя (не параллельно), поскольку SET_STRING_ELT и SET_VECTOR_ELT не являются потокобезопасными в R. Другие потоки обрабатывают все целочисленные, реальные, сложные и необработанные столбцы параллельно. Затем он идет так же быстро, как память io может идти.
Я действительно не вижу разницы, которую вы видите на 61 МБ, но масштабируясь до (все еще малого) 610 МБ, увеличивая количество столбцов с 10х до 80 000. Я вижу разницу.
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
У меня 128 Мбайт кэша 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