6

中央に配置したい大きなマトリックスがあります：

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

colMeans を使用すると、平均値をすばやく効率的に見つけることができます。

means <- colMeans(X)

しかし、各列からそれぞれの平均を減算する良い (高速でメモリ効率の良い) 方法は何ですか? これは機能しますが、正しくありません。

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

より良い方法はありますか？

/編集: これは、DWin が書いたさまざまなベンチマークを、他の投稿された提案を含む、より大きなマトリックスに変更したものです。

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

matmult 関数は新しい勝者のようです! これらを 5e+08 要素マトリックスで実際に試してみたいのですが、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
4

5 に答える 5

6

これはあなたにとって役に立ちますか？

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

forループに基づいてコードを高速化したい場合は、パッケージcmpfunから使用できます。compiler例

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

たぶん、この提案はあなたを助けることができます.

于 2012-09-08T16:16:26.997 に答える
3

ある時点で除算を要求しているのに、コードでは減算を使用しているため、Jilber があなたが何を望んでいるのかについて確信が持てなかった理由がわかります。彼が提案するスイープ操作は、ここでは不要です。スケールを使用するだけでそれができます：

 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

（それscaleが遅いとは信じがたいです。コードを見てください。sweep列での使用）

 scale.default # since it's visible.

マトリックスアプローチ:

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

編集：いくつかのタイミング（scalesweep-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

これをスケールアップすると、いくつかの面白いことが起こります。上記のティミンは、小さなマトリックス-Xで狂っていました。以下は、使用していたものに近いものです。

     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
于 2012-09-08T16:46:03.293 に答える
1

ここにいくつかありますが、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

整数ではなく数値のマトリックスでテストしていることに注意してください。より多くの人がそれが役立つと思うと思います（違いがあれば）。

于 2012-09-08T19:54:20.717 に答える