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

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？.我现在的代码已经按行将矩阵拆分成列表,但是必须拆分它会使它相当慢并且不比for循环快得多,因为我在每个矩阵中有近9000行:

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