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

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

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

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

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

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

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

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

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

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

Ответы

Ответ 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/С++ для каждой операции, которую вам нужно выполнить.

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

TL; DR

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

Ответы

Ответ 1

Ответ 2