Быстрое подмножество матрицы в R

Question

Я сталкиваюсь со следующей проблемой: мне нужно много подмножеств большой матрицы. На самом деле мне просто нужны представления в качестве входных данных для другой функции f(), поэтому мне не нужно менять значения. Однако кажется, что R ужасно медленный для этой задачи, или я делаю что-то не так (что кажется более вероятным). В игрушечном примере показано, сколько времени требуется, чтобы выбрать столбцы, а затем использовать их в другой функции (в этом игрушечном примере - примитивная функция sum()). В качестве "эталона" я также сравниваю время расчета с суммированием всей матрицы, что на удивление быстрее. Я также экспериментировал с пакетом ref, однако не смог добиться увеличения производительности. Таким образом, ключевой вопрос заключается в том, как поднастроить матрицу, не копируя ее? Я ценю любую помощь, спасибо!

library(microbenchmark)
library(ref)

m0 <- matrix(rnorm(10^6), 10^3, 10^3)
r0 <- refdata(m0)
microbenchmark(m0[, 1:900], sum(m0[, 1:900]), sum(r0[,1:900]), sum(m0))

Unit: milliseconds
             expr       min        lq      mean    median        uq
      m0[, 1:900] 10.087403 12.350751 16.697078 18.307475 19.054157
 sum(m0[, 1:900]) 11.067583 13.341860 17.286514 19.123748 19.990661
 sum(r0[, 1:900]) 11.066164 13.194244 16.869551 19.204434 20.004034
          sum(m0)  1.015247  1.040574  1.059872  1.049513  1.067142
       max neval
 58.238217   100
 25.664729   100
 23.505308   100
  1.233617   100

Контрольная задача суммирования всей матрицы занимает 1,059872 миллисекунды и примерно в 16 раз быстрее, чем другие функции.

Answer 1

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

У вас есть два решения:

Если время взято с sum в целом матрица в порядке с вами, вы могли бы использовать colSums на всю матрицу и подмножество результата:

sum(colSums(m0)[1:900])

Или вы можете использовать Rcpp для вычисления sum с поднабором без копирования матрицы.

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double sumSub(const NumericMatrix& x,
              const IntegerVector& colInd) {

  double sum = 0;

  for (IntegerVector::const_iterator it = colInd.begin(); it != colInd.end(); ++it) {
    int j = *it - 1;
    for (int i = 0; i < x.nrow(); i++) {
      sum += x(i, j);
    }
  }

  return sum;
}

    microbenchmark(m0[, 1:900], sum(m0[, 1:900]), sum(r0[,1:900]), sum(m0),
                   sum(colSums(m0)[1:900]),
                   sumSub(m0, 1:900))
Unit: milliseconds
                    expr      min       lq     mean   median       uq       max neval
             m0[, 1:900] 4.831616 5.447749 5.641096 5.675774 5.861052  6.418266   100
        sum(m0[, 1:900]) 6.103985 6.475921 7.052001 6.723035 6.999226 37.085345   100
        sum(r0[, 1:900]) 6.224850 6.449210 6.728681 6.705366 6.943689  7.565842   100
                 sum(m0) 1.110073 1.145906 1.175224 1.168696 1.197889  1.269589   100
 sum(colSums(m0)[1:900]) 1.113834 1.141411 1.178913 1.168312 1.201827  1.408785   100
       sumSub(m0, 1:900) 1.337188 1.368383 1.404744 1.390846 1.415434  2.459361   100

Вы можете использовать оптимизацию развертывания для дальнейшей оптимизации версии Rcpp.

Answer 2

С помощью compiler Я написал функцию, которая получает результат примерно в 2 раза быстрее, чем другие методы (в 8 раз sum(m0) вместо 16х):

require(compiler)

compiler_sum <- cmpfun({function(x) {
     tmp <- 0
     for (i in 1:900)
         tmp <- tmp+sum(x[,i])
     tmp
}})

microbenchmark( 
               sum(m0),
               compiler_sum(m0)
               )

Unit: milliseconds
             expr      min       lq     mean   median      uq       max
          sum(m0) 1.016532 1.056030 1.107263 1.084503 1.11173  1.634391
 compiler_sum(m0) 7.655251 7.854135 8.000521 8.021107 8.29850 16.760058
 neval
   100
   100