Подтвердить что ты не робот

Как оптимизировать чтение и запись в подразделы матрицы в R (возможно, используя data.table)

TL; DR

Какой самый быстрый метод в R для чтения и записи подмножества столбцы из очень большой матрицы. Я пытаюсь найти решение с data.table но нужен быстрый способ извлечения последовательности столбцов?

Короткий ответ: дорогой частью операции является присвоение. Таким образом, решение состоит в том, чтобы придерживаться матрицы и использовать Rcpp и С++ для изменения матрицы на месте. Ниже приведены примеры двух отличных ответов. [Для тех, кто обращается к другим проблемам, обязательно прочитайте отказы в решениях!]. Перейдите к нижней части вопроса для получения еще нескольких уроков.

Это мой первый вопрос о переполнении стека - я очень ценю ваше время, чтобы взглянуть, и я приношу свои извинения, если ничего не оставил. Я работаю над пакетом R, где у меня есть узкое место производительности из подмножества и записи в части матрицы (NB для статистиков, приложение обновляет достаточную статистику после обработки каждой точки данных). Индивидуальные операции невероятно быстры, но для их максимального количества требуется как можно быстрее. Простейшая версия идеи представляет собой матрицу размерности K на V, где K обычно находится между 5 и 1000, а V может быть между 1000 и 1,000,000.

set.seed(94253)
K <- 100
V <- 100000
mat <-  matrix(runif(K*V),nrow=K,ncol=V)

мы затем закончим выполнение вычисления по подмножеству столбцов и добавим его в полную матрицу. таким образом наивно он выглядит как

Vsub <- sample(1:V, 20)
toinsert <- matrix(runif(K*length(Vsub)), nrow=K, ncol=length(Vsub))
mat[,Vsub] <- mat[,Vsub] + toinsert
library(microbenchmark)
microbenchmark(mat[,Vsub] <- mat[,Vsub] + toinsert)

потому что это делается так много раз, что это может быть довольно медленным в результате семантики R copy-on-change (но см. уроки, полученные ниже, модификация может фактически произойти на месте в некоторых условиях).

Для моей задачи объект не обязательно должен быть матрицей (и я чувствителен к различию, описанному здесь Назначить матрицу подмножеству data.table). Мне всегда нужен полный столбец, поэтому структура списка кадра данных прекрасна. Мое решение состояло в том, чтобы использовать замечательный пакет данных. Запись может быть выполнена чрезвычайно быстро, используя set(). К сожалению, получение значения несколько сложнее. Мы должны вызвать переменные, заданные с помощью = FALSE, что резко замедляет работу.

library(data.table)
DT <- as.data.table(mat)  
set(DT, i=NULL, j=Vsub,DT[,Vsub,with=FALSE] + as.numeric(toinsert))

В функции set() с использованием я = NULL для ссылки на все строки невероятно быстро, но (по-видимому, из-за того, что вещи хранятся под капотом) нет сопоставимого варианта для j. @Roland отмечает в комментариях, что одним из вариантов было бы преобразование в тройное представление (номер строки, номер столбца, значение) и использование бинарного поиска data.tables для ускорения поиска. Я тестировал это вручную, и, хотя он быстрый, он примерно в три раза превышает требования к памяти для матрицы. Я хотел бы избежать этого, если это возможно.

Следуя этому вопросу: Время получения отдельных элементов из объектов data.table и data.frame. Хэдли Уикхэм дал невероятно быстрое решение для одного индекса

Vone <- Vsub[1]
toinsert.one <- toinsert[,1]
set(DT, i=NULL, j=Vone,(.subset2(DT, Vone) + toinsert.one))

однако, поскольку .subset2 (DT, i) является просто DT [[i]] без отправки методов, я, как мне известно, не имеет возможности захватить сразу несколько столбцов, хотя, похоже, это должно быть возможно. Как и в предыдущем вопросе, кажется, что, поскольку мы можем быстро перезаписать значения, мы должны быстро их прочитать.

Любые предложения? Также, пожалуйста, дайте мне знать, есть ли лучшее решение, чем data.table для этой проблемы. Я понял, что на самом деле это не предполагаемый прецедент во многих отношениях, но я стараюсь не переносить всю последовательность операций на C.

Ниже приведена последовательность таймингов обсуждаемых элементов: первые две - все столбцы, а вторая - только один столбец.

microbenchmark(mat[,Vsub] <- mat[,Vsub] + toinsert,
              set(DT, i=NULL, j=Vsub,DT[,Vsub,with=FALSE] + as.numeric(toinsert)),
              mat[,Vone] <- mat[,Vone] + toinsert.one,
              set(DT, i=NULL, j=Vone,(.subset2(DT, Vone) + toinsert.one)),
              times=1000L)

Unit: microseconds
                  expr      min       lq   median       uq       max neval
                Matrix   51.970   53.895   61.754   77.313   135.698  1000
            Data.Table 4751.982 4962.426 5087.376 5256.597 23710.826  1000
     Matrix Single Col    8.021    9.304   10.427   19.570 55303.659  1000
 Data.Table Single Col    6.737    7.700    9.304   11.549    89.824  1000

Ответ и извлеченные уроки:

Комментарии идентифицировали самую дорогую часть операции как процесс назначения. Оба решения дают ответы на основе кода C, которые изменяют матрицу на месте, нарушая правило R, не изменяя аргумент функции, но обеспечивая гораздо более быстрый результат.

Хэдли Уикхем остановился в комментариях, чтобы отметить, что матричная модификация фактически выполняется на месте, пока объектный мат не упоминается нигде (см. http://adv-r.had.co.nz/memory.html#modification-in-place). Это указывает на интересный и тонкий момент. Я выполнял свои оценки в RStudio. RStudio, как отмечает Хэдли в своей книге, создает дополнительную ссылку для каждого объекта, не входящего в функцию. Таким образом, хотя в случае производительности функции модификация будет происходить на месте, в командной строке он производит эффект копирования на смену. Hadley package pryr имеет несколько полезных функций для отслеживания ссылок и адресов памяти.

4b9b3361

ОТВЕТЫ

Ответ 1

Развлечения с Rcpp:

Вы можете использовать класс Eigen Map для изменения объекта R.

library(RcppEigen)
library(inline)

incl <- '
using  Eigen::Map;
using  Eigen::MatrixXd;
using  Eigen::VectorXi;
typedef  Map<MatrixXd>  MapMatd;
typedef  Map<VectorXi>  MapVeci;
'

body <- '
MapMatd              A(as<MapMatd>(AA));
const MapMatd        B(as<MapMatd>(BB));
const MapVeci        ix(as<MapVeci>(ind));
const int            mB(B.cols());
for (int i = 0; i < mB; ++i) 
{
A.col(ix.coeff(i)-1) += B.col(i);
}
'

funRcpp <- cxxfunction(signature(AA = "matrix", BB ="matrix", ind = "integer"), 
                       body, "RcppEigen", incl)

set.seed(94253)
K <- 100
V <- 100000
mat2 <-  mat <-  matrix(runif(K*V),nrow=K,ncol=V)

Vsub <- sample(1:V, 20)
toinsert <- matrix(runif(K*length(Vsub)), nrow=K, ncol=length(Vsub))
mat[,Vsub] <- mat[,Vsub] + toinsert

invisible(funRcpp(mat2, toinsert, Vsub))
all.equal(mat, mat2)
#[1] TRUE

library(microbenchmark)
microbenchmark(mat[,Vsub] <- mat[,Vsub] + toinsert,
               funRcpp(mat2, toinsert, Vsub))
# Unit: microseconds
#                                  expr    min     lq  median      uq       max neval
# mat[, Vsub] <- mat[, Vsub] + toinsert 49.273 49.628 50.3250 50.8075 20020.400   100
#         funRcpp(mat2, toinsert, Vsub)  6.450  6.805  7.6605  7.9215    25.914   100

Я думаю, что это в основном то, что предложил Джошуа Ульрих. Его предупреждения относительно нарушения функциональной парадигмы R применяются.

Я делаю добавление в С++, но тривиально изменять функцию только для назначения.

Очевидно, что если вы можете реализовать весь цикл в Rcpp, вы избежите повторных вызовов функций на уровне R и получите производительность.

Ответ 2

Вот что я имел в виду. Вероятно, это может быть намного сексуальнее с Rcpp и друзьями, но я не так хорошо знаком с этими инструментами.

#include <R.h>
#include <Rinternals.h>
#include <Rdefines.h>
SEXP addCol(SEXP mat, SEXP loc, SEXP matAdd)
{
  int i, nr = nrows(mat), nc = ncols(matAdd), ll = length(loc);
  if(ll != nc)
    error("length(loc) must equal ncol(matAdd)");
  if(TYPEOF(mat) != TYPEOF(matAdd))
    error("mat and matAdd must be the same type");
  if(nr != nrows(matAdd))
    error("mat and matAdd must have the same number of rows");
  if(TYPEOF(loc) != INTSXP)
    error("loc must be integer");
  int *iloc = INTEGER(loc);

  switch(TYPEOF(mat)) {
    case REALSXP:
      for(i=0; i < ll; i++)
        memcpy(&(REAL(mat)[(iloc[i]-1)*nr]),
               &(REAL(matAdd)[i*nr]), nr*sizeof(double));
      break;
    case INTSXP:
      for(i=0; i < ll; i++)
        memcpy(&(INTEGER(mat)[(iloc[i]-1)*nr]),
               &(INTEGER(matAdd)[i*nr]), nr*sizeof(int));
      break;
    default:
      error("unsupported type");
  }
  return R_NilValue;
}

Поместите вышеуказанную функцию в addCol.c, затем запустите R CMD SHLIB addCol.c. Тогда в R:

addColC <- dyn.load("addCol.so")$addCol
.Call(addColC, mat, Vsub, mat[,Vsub]+toinsert)

Небольшое преимущество этого подхода к Roland заключается в том, что это только назначение. Его функция делает дополнение для вас, что происходит быстрее, но также означает, что вам нужна отдельная функция C/С++ для каждой операции, которую вам нужно выполнить.