logo标签列表

在R中高效地将一个大矩阵居中

rcenterdata-manipulation
7

我有一个很大的矩阵需要居中：

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

使用colMeans函数快速高效地找到均值：

使用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元素矩阵上尝试这些，但我总是内存不足。

     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
也许 scale 函数可以帮到你。请参考 ?scale。另一个有用的函数可能是 sweep。 - Jilber Urbina
@Jiber：比起上面的for循环，scale函数要慢得多。不过sweep应该可以用，谢谢！ - Zach
“wuber” 是谁？benchmark 函数是由 Wacek Kusnierczyk 编写的。 - IRTFM
@DWin：抱歉，我引用了你的帖子，但名字错了。最近我一直在阅读wuber在交叉验证上写的东西。 - Zach
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循环的代码，可以使用compiler包中的cmpfun。例如：

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

也许这个建议可以帮到您。
5

这看起来比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编辑：我使用比OP使用的矩阵更小的矩阵（仅为5e + 07）时，得到了以下时间记录，其中Josh的矩阵是mat2（更大的矩阵在我的32GB Mac上溢出到虚拟内存并需要终止）：

  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
我很匆忙，否则我会进行更好的计时。如果你能做到，请随意在我的答案中添加它们。 - Josh O'Brien
非常感谢，@Dwin。真的很有趣看到简单矩阵操作快了多少。 - Josh O'Brien
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) )

编辑：一些时间安排（关于scale等同于扫描-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
实际上，似乎无论是扫描还是缩放都比我的for循环慢两倍左右。 - Zach
我编辑了原始帖子。感谢提供基准代码。然而，似乎在更大的矩阵（5,000行，10,000列或50,000行和10,000列）上，for循环实际上是最快的。 - Zach
3
也许将您的frlp()函数编译一下会稍微加快速度？
frlp.c <- compiler:::cmpfun(function(mat){
              means <- colMeans(mat)
              for (i in 1:length(means)){
                mat[,i] <- mat[,i]-means[i] 
              }
              mat
            }
          )

[编辑]: 对我来说，它并没有加速事情的进行，但我不得不大大缩小X以在我的计算机上工作。它可能会很好地扩展，不知道。
您还可以与JIT进行比较：
frlp.JIT <- function(mat){
              means <- colMeans(mat)
              compiler::enableJIT(2)
              for (i in 1:length(means)){
                mat[,i] <- mat[,i]-means[i] 
              }
              mat
            }
1

这里还有几个，但没有乔希的快：

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

请注意，我正在对数字矩阵进行测试，而不是整数；如果有任何区别的话，我认为更多的人会发现这很有用。
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接