Effizientes eine große Matrix in R Zentrum

Question

Ich habe eine große Matrix, die ich Mitte möchte: die

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

Finden der Mittel ist schnell und effizient mit colMeans:

means <- colMeans(X) 

Aber was ist ein gute (schnelle und speichereffiziente) Möglichkeit, den jeweiligen Mittelwert aus jeder Spalte zu subtrahieren? Das funktioniert, aber es fühlt sich nicht richtig an:

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

Gibt es einen besseren Weg?

/edit: Hier ist eine Modifikation die die verschiedenen Benchmarks Dwin schrieb, auf eine größere Matrix, einschließlich der anderen Vorschläge gepostet:

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')) 

Die matmult Funktion scheint neue Gewinner zu sein! Ich möchte das wirklich auf einer 5e + 08-Element-Matrix ausprobieren, aber mir geht immer noch der Arbeitsspeicher aus.

 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

Hier geht es um zu sein scheint doppelt so schnell wie 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 bearbeiten: Wenn ich mit einer kleineren Matrix tat dies als die OP verwendet (nur 5e + 07) ich diese Zeiten kommen, wo Joshs ist mat2 (Je größer ein auf meinem Mac in den virtuellen Speicher überschwemmt w/32GB und musste beendet werden):

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

Könnte das für Sie nützlich sein?

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 

Wenn Sie Ihren Code beschleunigen basierend auf dem for Schleifen Sie cmpfun von compiler Paket verwenden können. Beispiel

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 

Vielleicht könnte dieser Vorschlag Ihnen helfen.

Answer 3

ich sehen kann, warum Jilber unsicher war in Bezug auf was Sie wollten, an einem Punkt, da Sie für die Teilung fragen, aber in Ihrem Code Sie Subtraktion verwenden. Die von ihm vorgeschlagene Sweep-Operation ist hier überflüssig. Nur mit Skala es tun würde:

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 ist schwer zu glauben, dass scale langsamer ist Betrachtet es Code ist Verwendet sweep auf Spalten..)

scale.default # since it's visible. 

Ein Matrix-Ansatz:

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

Edit: Einige Timings (ich war falsch über scale gleichwertig fegen-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 

Einige lustige Dinge passieren, wenn Sie diese Scal auf. Die obigen Timigns waren mit Samall Matrix-X verrückt. Im Folgenden finden Sie mit etwas näher an, was Sie verwendet haben: beschleunigen würde die Dinge ein wenig

Vielleicht

 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

Ihre frlp() Funktion kompilieren?

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

[Bearbeiten]: für mich ist es nicht die Dinge beschleunigen, aber ich hatte X maßstäblich stark nach unten auf meinem Computer zu arbeiten.Es kann gut skalieren, weiß nicht

Sie können auch mit JIT vergleichen möchten:

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

Hier sind ein paar mehr, keine so schnell wie Joshs:

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 

Bitte beachte, dass ich auf einer Matrix aus Numerik mich entschieden, nicht ganze Zahlen sind; Ich denke, mehr Leute werden das nützlich finden (wenn es einen Unterschied macht.)