Эффективное применение или сопоставление для нескольких аргументов матрицы по строке

У меня есть две матрицы, к которым я хочу применить функцию, по строкам:

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 для матриц, Множественная матричная строка -поле 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)) )
}

Ответы

Ответ 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

UPDATE
Здесь "более чистая" версия с использованием функции 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) )

Ответ 2

Это решение позволяет избежать расщепления и списков, поэтому, возможно, это будет быстрее, чем ваша версия:

## original data:
tmp1 <- matrix(sample(1:100, 20), nrow = 5)
tmp2 <- matrix(sample(1:100, 20), nrow = 5)

## combine them together
tmp3 <- cbind(tmp1, tmp2)

## calculate t.stats:
t.stats <- apply(tmp3, 1, function(x) t.test(x[1:ncol(tmp1)], 
  x[(1 + ncol(tmp1)):ncol(tmp3)])$statistic)

Изменить: просто протестировал его на двух матрицах из 9000 строк и по 5 столбцов каждый, и он завершился менее чем за 6 секунд:

tmp1 <- matrix(rnorm(5 * 9000), nrow = 9000)
tmp2 <- matrix(rnorm(5 * 9000), nrow = 9000)
tmp3 <- cbind(tmp1, tmp2)
system.time(t.st <- apply(tmp3, 1, function(x) t.test(x[1:5], x[6:10])$statistic))

- > Истекшая пользовательская система

- > 5,640 0,012 5,705