Cómo optimizar las secciones de Lectura y Escritura de una matriz en R(posiblemente usando data.table)

performance

Question

TL; DR

¿Cuál es el método más rápido en R para leer y escribir un subconjunto de columnas de una matriz muy grande? Intento una solución con data.table pero necesito una forma rápida de extraer una secuencia de columnas?

Respuesta corta: la parte costosa de la operación es asignación. Por lo tanto, la solución es quedarse con una matriz y usar Rcpp y C ++ para modificar la matriz en su lugar. Hay dos respuestas excelentes a continuación con ejemplos. [Para aquellos que se aplican a otros problemas, ¡asegúrese de leer los descargos de responsabilidad en las soluciones!]. Desplácese hasta el final de la pregunta para obtener más lecciones aprendidas.

Esta es mi primera pregunta sobre el desbordamiento de pila: agradezco mucho su tiempo para echar un vistazo y me disculpo si he omitido algo. Estoy trabajando en un paquete R donde tengo un cuello de botella de rendimiento de subconjunto y escritura en partes de una matriz (NB para estadísticos la aplicación está actualizando suficientes estadísticas después de procesar cada punto de datos). Las operaciones individuales son increíblemente rápidas, pero la gran cantidad de ellas requiere que sea lo más rápido posible. La versión más simple de la idea es una matriz de dimensión K por V, donde K es generalmente entre 5 y 1000 y V puede estar entre 1000 y 1,000,000.

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

luego terminamos realizando un cálculo en un subconjunto de las columnas y agregando esto a la matriz completa. así ingenuamente parece

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)

porque esto se hace tantas veces puede ser bastante lento como resultado de la semántica de copia en el cambio de R (pero vea las lecciones aprendidas a continuación, la modificación puede ocurrir en algunas circunstancias).

Para mi problema, el objeto no necesita ser una matriz (y soy sensible a la diferencia como se describe aquí Asignar una matriz a un subconjunto de una tabla de datos ). Siempre quiero la columna completa y la estructura de lista de un marco de datos está bien. Mi solución fue utilizar el increíble paquete de data.table de Matthew Dowle. La escritura se puede hacer extraordinariamente rápido usando set (). Lamentablemente, obtener el valor es algo más complicado. Tenemos que llamar a la configuración de variables con = FALSO que ralentiza drásticamente las cosas.

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

Dentro de la función set () que usa i = NULL para hacer referencia a todas las filas es increíblemente rápido pero (presumiblemente debido a la forma en que se almacenan las cosas debajo del capó) no hay una opción comparable para j. @Roland señala en los comentarios que una opción sería convertir a una representación triple (número de fila, número de columna, valor) y usar la búsqueda binaria de data.tables para acelerar la recuperación. Probé esto manualmente y, aunque es rápido, aproximadamente triplica los requisitos de memoria para la matriz. Me gustaría evitar esto si es posible.

Siguiendo la pregunta aquí: Tiempo en obtener elementos individuales de objetos data.table y data.frame . Hadley Wickham dio una solución increíblemente rápida para un solo índice

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

sin embargo, dado que el .subset2 (DT, i) es simplemente DT [[i]] sin el envío de métodos, no hay forma (que yo sepa) de tomar varias columnas a la vez, aunque parece que debería ser posible. Al igual que en la pregunta anterior, parece que como podemos sobreescribir los valores rápidamente, deberíamos poder leerlos rápidamente.

¿Alguna sugerencia? También, avíseme si hay una solución mejor que data.table para este problema. Me di cuenta de que no es realmente el caso de uso previsto en muchos aspectos, pero estoy tratando de evitar transferir toda la serie de operaciones a C.

Aquí hay una secuencia de tiempos de los elementos discutidos: los primeros dos son todas las columnas, los segundos son solo una columna.

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

Respuesta y Lecciones Aprendidas:

Los comentarios identificaron la parte más costosa de la operación como el proceso de asignación. Ambas soluciones dan respuestas basadas en el código C que modifican la matriz en su lugar, rompiendo la convención de R de no modificar el argumento para una función, pero proporcionando un resultado mucho más rápido.

Hadley Wickham se detuvo en los comentarios para observar que la modificación de la matriz se realiza realmente en el lugar siempre y cuando no se haga referencia a la estera del objeto en otro lugar (ver http://adv-r.had.co.nz/memory.html#modification-in-place ). Esto apunta a un punto interesante y sutil. Estaba realizando mis evaluaciones en RStudio. RStudio, como lo señala Hadley en su libro, crea una referencia adicional para cada objeto que no está dentro de una función. Por lo tanto, mientras que en el caso de rendimiento de una función la modificación ocurriría en su lugar, en la línea de comando producía un efecto de copia en el cambio. El paquete pryr de Hadley tiene algunas bonitas funciones para rastrear referencias y direcciones de memoria.

Answer 1

Esto es lo que tenía en mente. Esto podría ser mucho más sexy con Rcpp y sus amigos, pero no estoy tan familiarizado con esas herramientas.

#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;
}

Coloque la función anterior en addCol.c , luego ejecute R CMD SHLIB addCol.c. Luego en R:

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

La ligera ventaja de este enfoque sobre Roland es que esto solo cumple la tarea. Su función hace la adición para usted, que es más rápida, pero también significa que necesita una función C / C ++ por separado para cada operación que necesita hacer.

Answer 2

Diversión con Rcpp:

Puede usar la clase del Mapa de Eigen para modificar un objeto R en su lugar.

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

Creo que esto es básicamente lo que propuso Joshua Ulrich. Sus advertencias sobre romper el paradigma funcional de R se aplican.

Hago la adición en C ++, pero es trivial cambiar la función para hacer solo la asignación.

Obviamente, si puede implementar todo su ciclo en Rcpp, evitará llamadas a funciones repetidas en el nivel R y obtendrá rendimiento.