42 votos

Cómo optimizar las operaciones de Lectura y Escritura a las subsecciones de una matriz en R (posiblemente el uso de datos.de la tabla)

TL;DR

¿Cuál es el método más rápido en R para la lectura y la escritura de un subconjunto de las columnas de una matriz. Yo intento de una solución con los datos.tabla pero necesita una manera rápida de extraer una secuencia de columnas?

Respuesta corta: Las caras de la parte de la operación de asignación. Por lo tanto la solución es seguir con una matriz y el uso de Rcpp y C++ para modificar la matriz en su lugar. Hay dos excelentes respuestas con ejemplos.[para aquellos que se aplican a otros problemas, asegúrese de leer los avisos en las soluciones!]. Desplácese hasta la parte inferior de la pregunta para algunas de las lecciones aprendidas.


Este es mi primer Stack Overflow pregunta - aprecio mucho tu tiempo en echar un vistazo y me disculpo si he omitido algo. Estoy trabajando en un paquete de R, donde tengo un cuello de botella de rendimiento de subconjunto y de la escritura a las partes de una matriz (NOTA para los estadísticos de la aplicación es la actualización de las estadísticas suficientes después de procesar cada punto de datos). Las operaciones individuales son increíblemente rápido, pero el gran número de ellos requiere para ser tan rápido como sea posible. La versión más simple de la idea es una matriz de dimensión K de V donde K es por lo general entre 5 y 1000 V y puede ser entre 1000 y 1.000.000.

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

nos acaban de realizar un cálculo en un subconjunto de las columnas y la adición de este en la matriz completa. por lo tanto 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)

debido a esto se hace muchas veces puede ser bastante lento como resultado de R de la copia en cambio semántica (pero ver las lecciones aprendidas a continuación, de modificación de la realidad puede suceder en algunos cricumstances).

Para mi problema, el objeto no necesita ser una matriz (y soy sensible a la diferencia como se indica aquí Asignar una matriz a un subconjunto de datos.de la tabla). Siempre quiero que la columna completa y así la lista de la estructura de una trama de datos es buena. Mi solución fue usar Mateo Dowle impresionante de datos.tabla paquete. La escritura puede ser hecho extraordinariamente rápidamente a través de la set(). Por desgracia, obteniendo el valor es algo más complicado. Tenemos que llamar a las variables de configuración con=FALSE, el cual reduce 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() con i=NULL para hacer referencia a todas las filas es increíblemente rápido, pero (presumiblemente debido a la forma en que las cosas son almacenados bajo el capó), no es comparable opción para j. @Roland notas en los comentarios que una opción sería convertir a una triple representación (número de fila, col número, valor) y el uso de datos.tablas de búsqueda binaria para acelerar la recuperación. He probado esto de forma manual y si bien es rápido, hace aproximadamente el triple de los requisitos de memoria para la matriz. Me gustaría evitar esto si es posible.

Tras la pregunta aquí: Tiempo en llegar solo elemets de datos.tabla y los datos.marco objetos. Hadley Wickham dio una increíble solución rápida para un único índice

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

sin embargo, desde el .subset2(DT,i) es simplemente DT[[i]] sin los métodos de envío no hay ninguna manera (a mi conocimiento) para captar varias columnas a la vez aunque ciertamente parece que debería ser posible. Como en la pregunta anterior, parece que ya podemos sobrescribir los valores de la rapidez con la que debe ser capaz de leer con ellos rápidamente.

Alguna sugerencia? También por favor, hágamelo saber si hay una solución mejor que la de datos.tabla para este problema. Me di cuenta de que no es realmente la intención de caso de uso, en muchos aspectos, pero estoy tratando de evitar trasladar toda la serie de operaciones a. C.

Aquí hay una secuencia de tiempos de los elementos de discusión de los dos primeros son todas las columnas, las dos segundas son sólo 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 las Lecciones Aprendidas:

Comentarios identificado la parte más costosa de la operación, como el proceso de asignación. Ambas soluciones de dar respuestas basadas en el código C que modificar la matriz en lugar de romper la I convención de no modificar el argumento a una función, pero proporcionando un mucho más rápido resultado.

Hadley Wickham se detuvo en los comentarios a la nota que la matriz de la modificación es que realmente se hace en el lugar, siempre que el objeto de la estera no se hace referencia en otros lugares (ver http://adv-r.had.co.nz/memory.html#modification-in-place). Esto apunta a una interesante y sutil punto. Yo estaba realizando mi evaluaciones en RStudio. RStudio como Hadley notas en su libro crea una referencia adicional para cada objeto no dentro de una función. Así, mientras en el rendimiento de caso de una función de la modificación que iba a suceder en su lugar, en la línea de comando que se estaba produciendo una copia por efecto del cambio. Hadley del paquete pryr tiene algunas funciones interesantes para el seguimiento de referencias y direcciones de memoria.

19voto

Roland Puntos 37641

Diversión con Rcpp:

Usted puede utilizar el Nombre de la clase de Mapa para modificar un R objeto 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 @Josué Ulrich propuesto. Sus advertencias con respecto a romper R del paradigma funcional aplicar.

Tengo que hacer la suma en C++, pero es trivial para cambiar la función a sólo hacer la asignación.

Obviamente, si usted puede implementar su bucle en la Rcpp, evitar las repetidas llamadas a la función en el nivel I y aumentará su rendimiento.

11voto

Joshua Ulrich Puntos 68776

Aquí lo que yo tenía en mente. Esto probablemente podría ser mucho más sexy con Rcpp y amigos, pero yo no estoy tan familiarizado con estas 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;
}

Poner la función anterior en addCol.c, a continuación, ejecute R CMD SHLIB addCol.c. A continuación, 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 sólo se hace la asignación. Su función realiza la adición para usted, que es más rápido, pero también significa que usted necesita separar la C/C++ de la función para cada operación que usted necesita hacer.

Iteramos.com

Iteramos es una comunidad de desarrolladores que busca expandir el conocimiento de la programación mas allá del inglés.
Tenemos una gran cantidad de contenido, y también puedes hacer tus propias preguntas o resolver las de los demás.

Powered by:

X