Accelerare un'operazione di matrice di remi

Question

Considera la seguente matrice,

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

Desidero prendere la differenza tra il rowMeans di due gruppi di dimensioni identiche prese a caso in questa matrice.

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

È abbastanza lento, purtroppo non ho capito l'output di RProf (sembrava che la maggior parte del tempo fosse trascorso is.data.frame??)

Suggerimenti per qualcosa di più efficiente?

Ho contemplato quanto segue:

• Rcpp: Dalle mie letture online credo che i remi di R siano abbastanza efficienti, quindi non è chiaro che aiuterebbe a questo passaggio. Mi piacerebbe essere convinto di dove il collo di bottiglia è davvero il primo, forse il mio intero design è non ottimale. Se la maggior parte del tempo viene trascorso per creare copie per ciascuna delle matrici più piccole, RCPP funzionerebbe meglio?

• Aggiornando a R-devel, sembra esserci un nuovo .rowMeans funzionare ancora più efficiente. Qualcuno l'ha provato?

Grazie.

Answer 1

A testa rowSums() chiamare un sottoinsieme di colonne da m può essere visto come la moltiplicazione della matrice tra m e un vettore di 0 o 1 indicando le colonne selezionate. Se giustappondi tutti quei vettori, finisci con una moltiplicazione tra due matrici (che è molto più efficiente):

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

Non hai bisogno di 2 chiamate a rowMeans. Puoi prima fare la sottrazione e chiamare rowMeans sul risultato.

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

is.data.frame fa parte dei controlli effettuati in rowMeans.

Aggiornamento: per quanto riguarda .rowMeans A R-devel, sembra che sia solo una chiamata diretta al codice interno (supponendo do_colsum non è cambiato). È definito come:

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

Nel tuo caso, m=1024 e n=1000.