accélérer une opération matricielle rowMeans

Question

Considérons la matrice suivante,

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

Je souhaite faire la différence entre le rowMeans de deux groupes de taille identique pris au hasard dans cette 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])
}))

C'est assez lent, malheureusement je n'ai pas compris le résultat de Rprof (il semblait que la plupart du temps était consacré à is.data.frame??)

Des suggestions pour quelque chose de plus efficace ?

J'ai envisagé ce qui suit :

• Rcpp:d'après mes lectures en ligne, je pense que rowMeans de R est assez efficace, il n'est donc pas clair que cela aiderait à cette étape.J'aimerais d'abord être convaincu de l'endroit où se situe réellement le goulot d'étranglement, peut-être que toute ma conception n'est pas optimale.Si la plupart du temps était consacré à faire des copies pour chacune des matrices plus petites, Rcpp fonctionnerait-il mieux ?

• mise à jour vers R-devel, il semble y avoir un nouveau .rowMeans fonction encore plus efficace.Quelqu'un l'a-t-il essayé ?

Merci.

Answer 1

Chaque rowSums() faire appel à un sous-ensemble de colonnes de m peut être vu comme la multiplication matricielle entre m et un vecteur de 0 ou 1 indiquant les colonnes sélectionnées.Si vous juxtaposez tous ces vecteurs, vous vous retrouvez avec une multiplication entre deux matrices (ce qui est bien plus efficace) :

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

Vous n'avez pas besoin de 2 appels pour rowMeans.Vous pouvez d'abord faire la soustraction et appeler rowMeans sur le résultat.

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

is.data.frame fait partie des contrôles effectués dans rowMeans.

MISE À JOUR:concernant .rowMeans dans R-devel, il semble que ce soit juste un appel direct au code interne (en supposant que do_colsum n'a pas changé).Il est défini comme :

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

Dans ton cas, m=1024 et n=1000.