在 R 中有效地居中一个大矩阵答案

【问题标题】：Efficiently center a large matrix in R在 R 中有效地居中一个大矩阵
【发布时间】：2012-09-02 04:13:09
【问题描述】：

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

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] 
}

有没有更好的办法？

/edit：这是 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

【问题讨论】：

也许scale 功能可以帮助你。见?scale。另一个有用的函数可能是sweep。
@Jiber：缩放功能比上面的 for 循环慢得多。扫描应该可以，谢谢！
谁是'wuber'？ benchmark 函数由 Wacek Kusnierczyk 编写。
@DWin：抱歉，我引用了您的帖子，但名称错误。我最近在交叉验证上阅读了很多 wuber 写的东西。

标签： r center data-manipulation

【解决方案1】：

这对你有用吗？

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

也许这个建议对你有帮助。

【讨论】：

【解决方案2】：

我明白为什么 Jilber 不确定您想要什么，因为有一次您要求除法，但在您的代码中您使用减法。他建议的扫描操作在这里是多余的。只需使用 scale 就可以了：

 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 等同于 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

当你扩大规模时会发生一些有趣的事情。上面的 timigns 对小矩阵-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

【讨论】：

其实sweep和scale好像都比我的for循环慢一倍左右。
我编辑了我的原始帖子。感谢基准代码。但是，似乎 for 循环实际上在较大的矩阵（5,000 行、10,000 列或 50,000 行和 10,000 列）上最快。

【解决方案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
            }

【讨论】：

【解决方案4】：

这似乎是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（较大的一个溢出到我的 Mac 上的虚拟内存中 w/32GB 并且需要被终止）：

  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

【讨论】：

我很着急，否则我会做一些更好的时间安排。如果您执行它们，任何人都可以随时将它们添加到我的答案中。
非常感谢@Dwin。看看简单的矩阵运算有多快真的很有趣。

【解决方案5】：

这里还有一些，没有一个比乔希的快：

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

请注意，我正在测试数字矩阵，而不是整数；我认为更多的人会发现这很有用（如果它有什么不同的话。）

【讨论】：