centrar eficientemente una gran matriz en I

Question

Tengo una matriz grande que me gustaría al centro:

X <- matrix(sample(1:10, 5e+08, replace=TRUE), ncol=10000) 

Encontrar el del medio es rápido y eficiente con colMeans:

means <- colMeans(X) 

Pero lo que es una buena (rápida y eficiente en la memoria) forma de restar la media respectiva de cada columna? Esto funciona, pero no se siente bien:

for (i in 1:length(means)){ 
    X[,i] <- X[,i]-means[i] 
} 

¿Hay una manera mejor?

/editar: He aquí una modificación de los puntos de referencia los diversos Dwin escribió, en una matriz más grande, incluyendo el otro publicado sugerencias:

require(rbenchmark) 
X <- matrix(sample(1:10, 5e+07, replace=TRUE), ncol=10000) 
frlp.c <- compiler:::cmpfun(function(mat){ 
    means <- colMeans(mat) 
    for (i in 1:length(means)){ 
    mat[,i] <- mat[,i]-means[i] 
    } 
    return(mat) 
}) 

mat.c <- compiler:::cmpfun(function(mat){ 
    t(t(X) - colMeans(X)) 
}) 

swp.c <- compiler:::cmpfun(function(mat){ 
    sweep(mat, 2, colMeans(mat), FUN='-') 
}) 

scl.c <- compiler:::cmpfun(function(mat){ 
    scale(mat, scale=FALSE) 
}) 

matmult.c <- compiler:::cmpfun(function(mat){ 
    mat-rep(1, nrow(mat)) %*% t(colMeans(mat)) 
}) 

benchmark( 
    frlp.c=frlp.c(X), 
    mat=mat.c(X), 
    swp=swp.c(X), 
    scl=scl.c(X), 
    matmult=matmult.c(X), 
    replications=10, 
    order=c('replications', 'elapsed')) 

La función matMult parece ser nuevo ganador! Realmente quiero probar esto en una matriz de elementos 5e + 08, pero me quedo sin memoria RAM.

 test replications elapsed relative user.self sys.self user.child sys.child 
5 matmult   10 11.98 1.000  7.47  4.47   NA  NA 
1 frlp.c   10 35.05 2.926  31.66  3.32   NA  NA 
2  mat   10 50.56 4.220  44.52  5.67   NA  NA 
4  scl   10 58.86 4.913  50.26  8.42   NA  NA 
3  swp   10 61.25 5.113  51.98  8.64   NA  NA 

Answer 1

Ésta parece ser aproximadamente el doble de rápido que sweep().

X - rep(1, nrow(X)) %*% t(colMeans(X)) 

X <- matrix(sample(1:10, 5e+06, replace=TRUE), ncol=10000) 
system.time(sweep(X, 2, colMeans(X))) 
    user system elapsed 
    0.33 0.00 0.33 
system.time(X - rep(1, nrow(X)) %*% t(colMeans(X))) 
    user system elapsed 
    0.15 0.03 0.19 

Dwin edición: Cuando hice esto con una matriz más pequeña que la OP utilizada (sólo 5e + 07) consigo estos tiempos, donde Josh es mat2 (El más grande se desbordó en la memoria virtual en mi Mac w/32 GB y necesario para ser terminados):

test replications elapsed relative user.self sys.self user.child sys.child 
2 mat2   1 0.546 1.000000  0.287 0.262   0   0 
3 mat   1 2.372 4.344322  1.569 0.812   0   0 
1 frlp   1 2.520 4.615385  1.720 0.809   0   0 
4 swp   1 2.990 5.476190  1.959 1.043   0   0 
5 scl   1 3.019 5.529304  1.984 1.046   0   0 

Answer 2

¿Será útil para usted?

sweep(X, 2, colMeans(X)) # this substracts the colMean to each col 
scale(X, center=TRUE, scale=FALSE) # the same 

sweep(X, 2, colMeans(X), FUN='/') # this makes division 

Si desea acelerar su código basado en el bucle for puede utilizar cmpfun de compiler paquete. Ejemplo

X <- matrix(sample(1:10, 500000, replace=TRUE), ncol=100) # some data 
means <- colMeans(X) # col means 

library(compiler) 

# One of your functions to be compiled and tested 
Mean <- function(x) { 
    for (i in 1:length(means)){ 
     X[,i] <- X[,i]-means[i] 
    } 
    return(X) 
} 



CMean <- cmpfun(Mean) # compiling the Mean function 

system.time(Mean(X)) 
    user system elapsed 
    0.028 0.016 0.101 
system.time(CMean(X)) 
    user system elapsed 
    0.028 0.012 0.066 

Quizás esta sugerencia pueda ayudarlo.

Answer 3

Entiendo por qué Jilber no estaba segura de lo que quería, ya que en un momento se solicita la división, pero en el código se usa la resta. La operación de barrido que sugiere es superflua aquí. Sólo mediante la escala lo haría:

cX <- scale(X, scale=FALSE) # does the centering with subtraction of col-means 
sX <- scale(X, center=FALSE) # does the scaling operation 
csX <- scale(X) # does both 

(Es difícil de creer que scale es más lento Mira que es código utiliza sweep en columnas..)

scale.default # since it's visible. 

Un enfoque de la matriz:

t(t(X)/colMeans(X)) 

Editar: Algunos tiempos (me equivoqué acerca de que scale es equivalente a sweep-colMeans):

require(rbenchmark) 
benchmark(
    mat={sX <- t(t(X)/colMeans(X)) }, 
    swp ={swX <- sweep(X, 2, colMeans(X), FUN='/')}, 
    scl={sX <- scale(X, center=FALSE)}, 
    replications=10^2, 
    order=c('replications', 'elapsed')) 
#----------- 
    test replications elapsed relative user.self sys.self user.child sys.child 
1 mat   100 0.015 1.000000  0.015  0   0   0 
2 swp   100 0.015 1.000000  0.015  0   0   0 
3 scl   100 0.025 1.666667  0.025  0   0   0 

Algunas cosas divertidas suceden cuando se amplía esto. Los timigns anteriores estaban enojados con una pequeña matriz-X. A continuación se muestra con algo más cercano a lo que estaba utilizando:

 benchmark( 
     frlp ={means <- colMeans(X) 
         for (i in 1:length(means)){ 
           X[,i] <- X[,i]-means[i] 
           } 
         }, 
     mat={sX <- t(t(X) - colMeans(X)) }, 
     swp ={swX <- sweep(X, 2, colMeans(X), FUN='-')}, 
     scl={sX <- scale(X, scale=FALSE)}, 
    replications=10^2, 
    order=c('replications', 'elapsed')) 
#  
    test replications elapsed relative user.self sys.self user.child sys.child 
2 mat   100 2.075 1.000000  1.262 0.820   0   0 
3 swp   100 2.964 1.428434  1.917 1.058   0   0 
4 scl   100 2.981 1.436627  1.935 1.059   0   0 
1 frlp   100 3.651 1.759518  2.540 1.128   0   0 

Answer 4

Quizás compilar su función frlp() aceleraría las cosas un poco?

frlp.c <- compiler:::cmpfun(function(mat){ 
       means <- colMeans(mat) 
       for (i in 1:length(means)){ 
       mat[,i] <- mat[,i]-means[i] 
       } 
       mat 
      } 
     ) 

[Editar]: para mí no acelerar las cosas, pero tuve que reducir en gran medida X a trabajar en mi equipo.Se puede escalar bien, no sé

También es posible que desee comparar con JIT:

frlp.JIT <- function(mat){ 
       means <- colMeans(mat) 
       compiler::enableJIT(2) 
       for (i in 1:length(means)){ 
       mat[,i] <- mat[,i]-means[i] 
       } 
       mat 
      } 

Answer 5

Éstos son algunos más, ninguno tan rápido como Josh:

X <- matrix(runif(1e6), ncol = 1000) 
matmult <- function(mat) mat - rep(1, nrow(mat)) %*% t(colMeans(mat)) 
contender1 <- function(mat) mat - colMeans(mat)[col(mat)] 
contender2 <- function(mat) t(apply(mat, 1, `-`, colMeans(mat))) 
contender3 <- function(mat) mat - rep(colMeans(mat), each = nrow(mat)) 
contender4 <- function(mat) mat - matrix(colMeans(mat), nrow(mat), ncol(mat), 
             byrow = TRUE) 
benchmark(matmult(X), 
      contender1(X), 
      contender2(X), 
      contender3(X), 
      contender4(X), 
      replications = 100, 
      order=c('replications', 'elapsed')) 
#  test replications elapsed relative user.self sys.self 
# 1 matmult(X)   100 1.41 1.000000  1.39  0.00 
# 5 contender4(X)   100 1.90 1.347518  1.90  0.00 
# 4 contender3(X)   100 2.69 1.907801  2.69  0.00 
# 2 contender1(X)   100 2.74 1.943262  2.73  0.00 
# 3 contender2(X)   100 6.30 4.468085  6.26  0.03 

Tenga en cuenta que estoy probando en una matriz de valores numéricos, no enteros; Creo que a mucha gente le resultará útil (si hace alguna diferencia)