Эффективное применение или отображение нескольких матричных аргументов по строкам

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, учитывая, что в каждой матрице у меня почти 9000 строк:Сама функция очень проста, просто найти значение 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

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

Answer 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