efficiente centrare una grande matrice R

Question

Ho un grande matrice che vorrei centro:

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

Trovare il il mezzo è rapido ed efficiente con colMeans:

means <- colMeans(X) 

Ma che cosa è un buon modo (veloce ed efficiente in termini di memoria) per sottrarre la media corrispondente da ciascuna colonna? Funziona, ma non sembra giusto:

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

C'è un modo migliore?

/EDIT: Ecco una modifica dei vari parametri di riferimento DWin ha scritto, su una matrice più grande, tra cui l'Altro Pubblicato suggerimenti:

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 funzione matmult sembra essere nuovo vincitore! Voglio davvero provarli su una matrice di elementi 5e + 08, ma continuo a rimanere senza 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

Questo sembra essere circa due volte più veloce 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 edit: Quando ho fatto questo con una matrice più piccola della OP utilizzato (solo 5e + 07) ottengo questi tempi, in cui Josh è mat2 (Il più grande straripato nella memoria virtuale sul mio Mac w/32GB e aveva bisogno di essere terminati):

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

Potrebbe essere utile per voi?

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 

Se si vuole accelerare il vostro codice basato sul for ciclo è possibile utilizzare cmpfun da compiler pacchetto. Esempio

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 

Forse questo suggerimento potrebbe aiutarti.

Answer 3

Posso capire perché Jilber era incerto riguardo a ciò che volevi, dal momento che a un certo punto chiedi la divisione ma nel tuo codice usi la sottrazione. L'operazione di spazzata che suggerisce è superflua qui. Proprio utilizzando la scala avrebbe fatto:

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 

(E 'difficile credere che scale è più lento Guarda è codice utilizza sweep su colonne..)

scale.default # since it's visible. 

Un approccio a matrice:

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

Modifica: alcuni tempi (mi sono sbagliato su scale equivale 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 

Alcune cose divertenti accadono quando si scala questo. I tempi di cui sopra erano pazzi di samall matrix-X. Qui di seguito è con qualcosa di più vicino a quello che si sta utilizzando:

 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

compilazione Forse la vostra funzione frlp() sarebbe accelerare le cose un po '?

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

[Edit]: per me non accelerare le cose, ma ho dovuto ridimensionare notevolmente X a lavorare sul mio computer.Può scala bene, non so

Puoi anche confrontare 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

Qui ci sono un paio di più, nessuno veloce come 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 

Nota che sto testando su una matrice di valori numerici, non interi; Penso che più persone lo troveranno utile (se fa alcuna differenza)