r - 按行有效地应用或映射多个矩阵参数

Question

我有两个矩阵要按行应用函数:

matrixA
           GSM83009  GSM83037  GSM83002  GSM83029  GSM83041
100001_at  5.873321  5.416164  3.512227  6.064150  3.713696
100005_at  5.807870  6.810829  6.105804  6.644000  6.142413
100006_at  2.757023  4.144046  1.622930  1.831877  3.694880

matrixB
          GSM82939 GSM82940 GSM82974 GSM82975
100001_at 3.673556 2.372952 3.228049 3.555816
100005_at 6.916954 6.909533 6.928252 7.003377
100006_at 4.277985 4.856986 3.670161 4.075533

我发现了几个类似的问题，但没有找到很多答案: mapply for matrices , Multi matrix row-wise mapply? .我现在的代码将矩阵逐行拆分为列表，但考虑到我在每个矩阵中有近 9000 行，因此必须将其拆分使其相当慢且比 for 循环快不了多少:

scores <- mapply(t.test.stat, split(matrixA, row(matrixA)), split(matrixB, row(matrixB)))

函数本身很简单，只要找到t值:

t.test.stat <- function(x, y)
{
    return( (mean(x) - mean(y)) / sqrt(var(x)/length(x) + var(y)/length(y)) )
}

Answer 1

拆分矩阵不是评估时间的最大贡献者。

set.seed(21)
matrixA <- matrix(rnorm(5 * 9000), nrow = 9000)
matrixB <- matrix(rnorm(4 * 9000), nrow = 9000)

system.time( scores <- mapply(t.test.stat,
    split(matrixA, row(matrixA)), split(matrixB, row(matrixB))) )
#    user  system elapsed 
#    1.57    0.00    1.58 
smA <- split(matrixA, row(matrixA))
smB <- split(matrixB, row(matrixB))
system.time( scores <- mapply(t.test.stat, smA, smB) )
#    user  system elapsed 
#    1.14    0.00    1.14 

查看 Rprof 的输出看到大部分时间——毫不奇怪——花在评估上 t.test.stat ( mean 、 var 等)。基本上，函数调用有相当多的开销。

Rprof()
scores <- mapply(t.test.stat, smA, smB)
Rprof(NULL)
summaryRprof()

您也许能够找到更快的通用解决方案，但没有一个能达到下面矢量化解决方案的速度。

由于您的函数很简单，您可以利用矢量化 rowMeans几乎立即执行此操作的函数(尽管有点困惑):

system.time({
ncA <- NCOL(matrixA)
ncB <- NCOL(matrixB)
ans <- (rowMeans(matrixA)-rowMeans(matrixB)) /
  sqrt( rowMeans((matrixA-rowMeans(matrixA))^2)*(ncA/(ncA-1))/ncA +
        rowMeans((matrixB-rowMeans(matrixB))^2)*(ncB/(ncB-1))/ncB )
})
#    user  system elapsed 
#      0       0       0 
head(ans)
# [1]  0.8272511 -1.0965269  0.9862844 -0.6026452 -0.2477661  1.1896181

更新
这是使用 rowVars 的“更清洁”版本功能:

rowVars <- function(x, na.rm=FALSE, dims=1L) {
  rowMeans((x-rowMeans(x, na.rm, dims))^2, na.rm, dims)*(NCOL(x)/(NCOL(x)-1))
}
ans <- (rowMeans(matrixA)-rowMeans(matrixB)) /
  sqrt( rowVars(matrixA)/NCOL(matrixA) + rowVars(matrixB)/NCOL(matrixB) )