ускорить работу матрицы rowMeans

Question

Рассмотрим следующую матрицу:

nc <- 5000
nr <- 1024
m <- matrix(rnorm(nc*nr), ncol=nc)

Я хочу взять разницу между rowMeans двух групп одинакового размера, взятых случайным образом в этой матрице.

n <- 1000 # group size

system.time(replicate(100, {
   ind1 <- sample(seq.int(nc), n) 
   ind2 <- sample(seq.int(nc), n)
   rowMeans(m[, ind1]) - rowMeans(m[, ind2])
}))

Это довольно медленно, к сожалению, я не понял вывод Rprof (казалось, большая часть времени была потрачена на is.data.frame ??)

Предложения для чего-то более эффективного?

Я обдумал следующее:

• Rcpp: из моих онлайн-чтений я считаю, что rowMeans R довольно эффективен, поэтому не ясно, поможет ли он на этом этапе,Я хотел бы убедиться в том, где на самом деле узкое место, возможно, весь мой дизайн неоптимален.Если большая часть времени будет потрачена на изготовление копий для каждой из меньших матриц, будет ли Rcpp работать лучше?

• при обновлении до R-devel, похоже, появилась новая функция .rowMeansеще эффективнее.Кто-нибудь пробовал?

Спасибо.

Answer 1

Каждый вызов rowSums() для подмножества столбцов из m можно рассматривать как умножение матрицы между m и вектором 0 или 1, указывающим выбранные столбцы. Если вы сопоставите все эти векторы, вы получите умножение между двумя матрицами (что гораздо эффективнее):

ind1 <- replicate(100, seq.int(nc) %in% sample(seq.int(nc), n)) 
ind2 <- replicate(100, seq.int(nc) %in% sample(seq.int(nc), n))
output <- m %*% (ind1 - ind2)

Answer 2

Вам не нужно 2 звонка на rowMeans. Вы можете сначала выполнить вычитание и вызвать результат rowMeans.

x1 <- rowMeans(m[,ind1])-rowMeans(m[,ind2])
x2 <- rowMeans(m[,ind1]-m[,ind2])
all.equal(x1,x2)
# [1] TRUE

is.data.frame является частью проверок, выполненных в rowMeans.

ОБНОВЛЕНИЕ: относительно .rowMeans в R-devel, похоже, что это просто прямой вызов внутреннего кода (при условии, что do_colsum не изменился). Это определяется как:

.rowMeans <- function(X, m, n, na.rm = FALSE)
    .Internal(rowMeans(X, m, n, na.rm))

В вашем случае m=1024 и n=1000.