Accélérez le fonctionnement de la boucle dans R

193

J'ai un gros problème de performances dans R. J'ai écrit une fonction qui itère sur un data.frameobjet. Il ajoute simplement une nouvelle colonne à a data.frameet accumule quelque chose. (opération simple). Le data.framea environ 850K lignes. Mon PC fonctionne toujours (environ 10h maintenant) et je n'ai aucune idée du runtime.

dayloop2 <- function(temp){
    for (i in 1:nrow(temp)){    
        temp[i,10] <- i
        if (i > 1) {             
            if ((temp[i,6] == temp[i-1,6]) & (temp[i,3] == temp[i-1,3])) { 
                temp[i,10] <- temp[i,9] + temp[i-1,10]                    
            } else {
                temp[i,10] <- temp[i,9]                                    
            }
        } else {
            temp[i,10] <- temp[i,9]
        }
    }
    names(temp)[names(temp) == "V10"] <- "Kumm."
    return(temp)
}

Des idées pour accélérer cette opération?

Kay
la source

Réponses:

435

Le plus gros problème et la racine de l'inefficacité est l'indexation de data.frame, je veux dire toutes ces lignes où vous utilisez temp[,].
Essayez d'éviter cela autant que possible. J'ai pris ta fonction, change d'indexation et ici version_A

dayloop2_A <- function(temp){
    res <- numeric(nrow(temp))
    for (i in 1:nrow(temp)){    
        res[i] <- i
        if (i > 1) {             
            if ((temp[i,6] == temp[i-1,6]) & (temp[i,3] == temp[i-1,3])) { 
                res[i] <- temp[i,9] + res[i-1]                   
            } else {
                res[i] <- temp[i,9]                                    
            }
        } else {
            res[i] <- temp[i,9]
        }
    }
    temp$`Kumm.` <- res
    return(temp)
}

Comme vous pouvez le voir, je crée des vecteurs resqui rassemblent les résultats. À la fin, je l'ajoute data.frameet je n'ai pas besoin de jouer avec les noms. Alors, comment est-ce mieux?

Je lance chaque fonction pour data.frameavec nrowde 1.000 à 10.000 par 1000 et la mesure du temps avecsystem.time

X <- as.data.frame(matrix(sample(1:10, n*9, TRUE), n, 9))
system.time(dayloop2(X))

Le résultat est

performance

Vous pouvez voir que votre version dépend de façon exponentielle de nrow(X). La version modifiée a une relation linéaire, et le lmmodèle simple prévoit que pour 850 000 lignes, le calcul prend 6 minutes et 10 secondes.

Puissance de vectorisation

Comme Shane et Calimo le déclarent dans leurs réponses, la vectorisation est la clé d'une meilleure performance. À partir de votre code, vous pouvez sortir de la boucle:

  • conditionnement
  • initialisation des résultats (qui sont temp[i,9])

Cela conduit à ce code

dayloop2_B <- function(temp){
    cond <- c(FALSE, (temp[-nrow(temp),6] == temp[-1,6]) & (temp[-nrow(temp),3] == temp[-1,3]))
    res <- temp[,9]
    for (i in 1:nrow(temp)) {
        if (cond[i]) res[i] <- temp[i,9] + res[i-1]
    }
    temp$`Kumm.` <- res
    return(temp)
}

Comparez le résultat pour ces fonctions, cette fois nrowde 10 000 à 100 000 par 10 000.

performance

Accorder l'écoute

Un autre ajustement consiste à changer dans une boucle l'indexation temp[i,9]vers res[i](qui sont exactement les mêmes dans l'itération de boucle i-ème). C'est encore une différence entre l'indexation d'un vecteur et l'indexation d'un data.frame.
Deuxième chose: lorsque vous regardez la boucle, vous pouvez voir qu'il n'est pas nécessaire de boucler sur tout i, mais seulement pour celles qui correspondent à la condition.
Alors on y va

dayloop2_D <- function(temp){
    cond <- c(FALSE, (temp[-nrow(temp),6] == temp[-1,6]) & (temp[-nrow(temp),3] == temp[-1,3]))
    res <- temp[,9]
    for (i in (1:nrow(temp))[cond]) {
        res[i] <- res[i] + res[i-1]
    }
    temp$`Kumm.` <- res
    return(temp)
}

Les performances que vous gagnez dépendent fortement d'une structure de données. Précisément - sur le pourcentage des TRUEvaleurs de la condition. Pour mes données simulées, il faut du temps de calcul pour 850 000 lignes en dessous d'une seconde.

performance

Je veux que tu puisses aller plus loin, je vois au moins deux choses qui peuvent être faites:

  • écrire un Ccode pour faire du sperme conditionnel
  • si vous savez que dans votre séquence de données max n'est pas grande, vous pouvez changer la boucle en while vectorisé, quelque chose comme

    while (any(cond)) {
        indx <- c(FALSE, cond[-1] & !cond[-n])
        res[indx] <- res[indx] + res[which(indx)-1]
        cond[indx] <- FALSE
    }

Le code utilisé pour les simulations et les figures est disponible sur GitHub .

Marek
la source
2
Comme je ne trouve pas le moyen de demander à Marek en privé, comment ces graphiques ont-ils été générés?
carbontwelve
@carbontwelve Vous posez des questions sur les données ou les graphiques? Les parcelles ont été faites avec un paquet de treillis. Si j'ai le temps, je mets le code quelque part sur le Web et je vous en informe.
Marek
@carbontwelve Ooops, je me suis trompé :) Ce sont des graphiques standard (à partir de la base R).
Marek
@Gregor Malheureusement non. Il est cumulatif donc vous ne pouvez pas le vectoriser. Exemple simple: res = c(1,2,3,4)et condest tout TRUE, alors résultat final devrait être: 1, 3(cause 1+2), 6(deuxième cause est maintenant 3, et le troisième est 3aussi), 10( 6+4). Faire simple somme que vous avez 1, 3, 5, 7.
Marek
Ah, j'aurais dû y réfléchir plus attentivement. Merci de m'avoir montré l'erreur.
Gregor Thomas
132

Stratégies générales pour accélérer le code R

Tout d'abord, déterminez se trouve réellement la partie lente. Il n'est pas nécessaire d'optimiser le code qui ne s'exécute pas lentement. Pour de petites quantités de code, le simple fait d'y réfléchir peut fonctionner. Si cela échoue, RProf et des outils de profilage similaires peuvent être utiles.

Une fois que vous avez déterminé le goulot d'étranglement, pensez à des algorithmes plus efficaces pour faire ce que vous voulez. Les calculs ne doivent être exécutés qu'une seule fois si possible, donc:

L' utilisation de fonctions plus efficaces peut produire des gains de vitesse modérés ou importants. Par exemple, paste0produit un petit gain d'efficacité mais .colSums()et ses parents produisent des gains un peu plus prononcés. meanest particulièrement lent .

Ensuite, vous pouvez éviter certains problèmes particulièrement courants :

  • cbind vous ralentira très rapidement.
  • Initialisez vos structures de données, puis remplissez-les plutôt que de les étendre à chaque fois .
  • Même avec la pré-allocation, vous pouvez passer à une approche passe par référence plutôt qu'à une approche passe par valeur, mais cela ne vaut peut-être pas la peine.
  • Jetez un œil au R Inferno pour d'autres pièges à éviter.

Essayez une meilleure vectorisation , ce qui peut souvent mais pas toujours aider. À cet égard, des commandes intrinsèquement vectorisées comme ifelse,diff et autres apporteront plus d'améliorations que la applyfamille de commandes (qui ne fournissent que peu ou pas d'augmentation de vitesse sur une boucle bien écrite).

Vous pouvez également essayer de fournir plus d' informations aux fonctions de R . Par exemple, utilisez vapplyplutôt quesapply et spécifiez colClasseslors de la lecture de données textuelles . Les gains de vitesse seront variables en fonction de la quantité de devinettes que vous éliminez.

Ensuite, considérez les packages optimisés : le data.tablepackage peut produire des gains de vitesse massifs là où son utilisation est possible, dans la manipulation de données et dans la lecture de grandes quantités de données (fread ).

Ensuite, essayez de gagner de la vitesse grâce à des moyens plus efficaces d'appeler R :

  • Compilez votre script R. Ou utilisez les packages Raet jiten concert pour la compilation juste à temps (Dirk a un exemple dans cette présentation ).
  • Assurez-vous d'utiliser un BLAS optimisé. Ceux-ci fournissent des gains de vitesse à tous les niveaux. Honnêtement, c'est dommage que R n'utilise pas automatiquement la bibliothèque la plus efficace lors de l'installation. Espérons que Revolution R contribuera le travail qu'ils ont fait ici à l'ensemble de la communauté.
  • Radford Neal a fait un tas d'optimisations, dont certaines ont été adoptées dans R Core, et beaucoup d'autres qui ont été transformées en pqR .

Et enfin, si tout ce qui précède ne vous permet toujours pas de vous rendre aussi rapidement que nécessaire, vous devrez peut-être passer à un langage plus rapide pour l'extrait de code lent . La combinaison de Rcppet inlineici permet de ne remplacer que la partie la plus lente de l'algorithme par du code C ++ particulièrement facile. Ici, par exemple, est ma première tentative de le faire , et cela épate même les solutions R hautement optimisées.

Si vous avez encore des problèmes après tout cela, vous avez juste besoin de plus de puissance de calcul. Examinez la parallélisation ( http://cran.r-project.org/web/views/HighPerformanceComputing.html ) ou même les solutions basées sur GPU ( gpu-tools).

Liens vers d'autres conseils

Ari B. Friedman
la source
36

Si vous utilisez des forboucles, vous codez probablement R comme s'il s'agissait de C ou Java ou autre chose. Un code R correctement vectorisé est extrêmement rapide.

Prenons par exemple ces deux simples bits de code pour générer une liste de 10000 entiers en séquence:

Le premier exemple de code est de savoir comment coder une boucle en utilisant un paradigme de codage traditionnel. Cela prend 28 secondes pour terminer

system.time({
    a <- NULL
    for(i in 1:1e5)a[i] <- i
})
   user  system elapsed 
  28.36    0.07   28.61 

Vous pouvez obtenir une amélioration presque 100 fois par la simple action de pré-allocation de mémoire:

system.time({
    a <- rep(1, 1e5)
    for(i in 1:1e5)a[i] <- i
})

   user  system elapsed 
   0.30    0.00    0.29 

Mais en utilisant l'opération de vecteur de base R utilisant l'opérateur deux-points, :cette opération est pratiquement instantanée:

system.time(a <- 1:1e5)

   user  system elapsed 
      0       0       0 
Andrie
la source
+1 bien que je considère votre deuxième exemple comme peu convaincant car a[i]ne change pas. Mais system.time({a <- NULL; for(i in 1:1e5){a[i] <- 2*i} }); system.time({a <- 1:1e5; for(i in 1:1e5){a[i] <- 2*i} }); system.time({a <- NULL; a <- 2*(1:1e5)})a un résultat similaire.
Henry
@Henry, juste commentaire, mais comme vous le faites remarquer, les résultats sont les mêmes. J'ai modifié l'exemple pour initialiser un to rep(1, 1e5)- les horaires sont identiques.
Andrie
17

Cela pourrait être rendu beaucoup plus rapide en sautant les boucles en utilisant des index ou des ifelse()instructions imbriquées .

idx <- 1:nrow(temp)
temp[,10] <- idx
idx1 <- c(FALSE, (temp[-nrow(temp),6] == temp[-1,6]) & (temp[-nrow(temp),3] == temp[-1,3]))
temp[idx1,10] <- temp[idx1,9] + temp[which(idx1)-1,10] 
temp[!idx1,10] <- temp[!idx1,9]    
temp[1,10] <- temp[1,9]
names(temp)[names(temp) == "V10"] <- "Kumm."
Shane
la source
Merci d'avoir répondu. J'essaye de comprendre vos déclarations. La ligne 4: "temp [idx1,10] <- temp [idx1,9] + temp [which (idx1) -1,10]" a provoqué une erreur car la longueur de l'objet plus long n'est pas un multiple de la longueur du objet plus court. "temp [idx1,9] = num [1: 11496]" et "temp [which (idx1) -1,10] = int [1: 11494]" donc 2 lignes sont manquantes.
Kay
Si vous fournissez un échantillon de données (utilisez dput () avec quelques lignes), je le réparerai pour vous. En raison de quel () - 1 bit, les index sont inégaux. Mais vous devriez voir comment cela fonctionne à partir d'ici: il n'y a pas besoin de boucle ou d'application; utilisez simplement des fonctions vectorisées.
Shane
1
Hou la la! Je viens de changer un bloc de fonction imbriqué if..else et mapply, en une fonction ifelse imbriquée et j'ai un speedup 200x!
James
Votre conseil général est correct, mais dans le code, vous avez manqué un fait, cette ivaleur -th dépend de i-1-th donc elle ne peut pas être définie comme vous le faites (en utilisant which()-1).
Marek
8

Je n'aime pas la réécriture du code ... Bien sûr, ifelse et lapply sont de meilleures options, mais parfois il est difficile de faire cela.

J'utilise fréquemment data.frames comme on utiliserait des listes telles que df$var[i]

Voici un exemple inventé:

nrow=function(x){ ##required as I use nrow at times.
  if(class(x)=='list') {
    length(x[[names(x)[1]]])
  }else{
    base::nrow(x)
  }
}

system.time({
  d=data.frame(seq=1:10000,r=rnorm(10000))
  d$foo=d$r
  d$seq=1:5
  mark=NA
  for(i in 1:nrow(d)){
    if(d$seq[i]==1) mark=d$r[i]
    d$foo[i]=mark
  }
})

system.time({
  d=data.frame(seq=1:10000,r=rnorm(10000))
  d$foo=d$r
  d$seq=1:5
  d=as.list(d) #become a list
  mark=NA
  for(i in 1:nrow(d)){
    if(d$seq[i]==1) mark=d$r[i]
    d$foo[i]=mark
  }
  d=as.data.frame(d) #revert back to data.frame
})

Version data.frame:

   user  system elapsed 
   0.53    0.00    0.53

version de la liste:

   user  system elapsed 
   0.04    0.00    0.03 

17 fois plus rapide d'utiliser une liste de vecteurs qu'un data.frame.

Des commentaires sur les raisons pour lesquelles les data.frames internes sont si lents à cet égard? On pourrait penser qu'ils fonctionnent comme des listes ...

Pour un code encore plus rapide, faites ceci class(d)='list'au lieu de d=as.list(d)etclass(d)='data.frame'

system.time({
  d=data.frame(seq=1:10000,r=rnorm(10000))
  d$foo=d$r
  d$seq=1:5
  class(d)='list'
  mark=NA
  for(i in 1:nrow(d)){
    if(d$seq[i]==1) mark=d$r[i]
    d$foo[i]=mark
  }
  class(d)='data.frame'
})
head(d)
Chris
la source
1
C'est probablement grâce à la surcharge de [<-.data.frame, qui est en quelque sorte appelée lorsque vous le faites d$foo[i] = market qui peut finir par faire une nouvelle copie du vecteur de éventuellement tout le data.frame à chaque <-modification. Cela ferait une question intéressante sur SO.
Frank
2
@Frank Il (i) doit s'assurer que l'objet modifié est toujours un data.frame valide et (ii) afaik en fait au moins une copie, éventuellement plus d'une. La sous-attribution de Dataframe est connue pour être lente et si vous regardez le long code source, ce n'est pas vraiment surprenant.
Roland
@Frank, @Roland: La df$var[i]notation passe- t -elle par la même [<-.data.framefonction? J'ai remarqué que c'était assez long en effet. Sinon, quelle fonction utilise-t-il?
Chris
@Chris, je crois, d$foo[i]=markest grossièrement traduit en d <- `$<-`(d, 'foo', `[<-`(d$foo, i, mark)), mais avec une certaine utilisation de variables temporaires.
Tim Goodman
7

Comme Ari l'a mentionné à la fin de sa réponse, les packages Rcppet inlinefacilitent incroyablement la rapidité des choses. À titre d'exemple, essayez ce inlinecode (avertissement: non testé):

body <- 'Rcpp::NumericMatrix nm(temp);
         int nrtemp = Rccp::as<int>(nrt);
         for (int i = 0; i < nrtemp; ++i) {
             temp(i, 9) = i
             if (i > 1) {
                 if ((temp(i, 5) == temp(i - 1, 5) && temp(i, 2) == temp(i - 1, 2) {
                     temp(i, 9) = temp(i, 8) + temp(i - 1, 9)
                 } else {
                     temp(i, 9) = temp(i, 8)
                 }
             } else {
                 temp(i, 9) = temp(i, 8)
             }
         return Rcpp::wrap(nm);
        '

settings <- getPlugin("Rcpp")
# settings$env$PKG_CXXFLAGS <- paste("-I", getwd(), sep="") if you want to inc files in wd
dayloop <- cxxfunction(signature(nrt="numeric", temp="numeric"), body-body,
    plugin="Rcpp", settings=settings, cppargs="-I/usr/include")

dayloop2 <- function(temp) {
    # extract a numeric matrix from temp, put it in tmp
    nc <- ncol(temp)
    nm <- dayloop(nc, temp)
    names(temp)[names(temp) == "V10"] <- "Kumm."
    return(temp)
}

Il existe une procédure similaire pour #includeles choses, où vous passez simplement un paramètre

inc <- '#include <header.h>

à la fonction cxx, comme include=inc . Ce qui est vraiment cool à ce sujet, c'est qu'il fait tout le lien et la compilation pour vous, donc le prototypage est vraiment rapide.

Avertissement: je ne suis pas totalement sûr que la classe de tmp doive être numérique et non numérique ou autre. Mais j'en suis surtout sûr.

Edit: si vous avez encore besoin de plus de vitesse après cela, OpenMP est une fonction de parallélisation idéale C++. Je n'ai pas essayé de l'utiliser à partir de inline, mais cela devrait fonctionner. L'idée serait, dans le cas des ncœurs, de faire effectuer une itération de boucle kpar k % n. Une introduction convenable en Matloff est l'art de la programmation R , disponible ici , au chapitre 16, Recourir à C .

jclancy
la source
3

Les réponses ici sont excellentes. Un aspect mineur non couvert est que la question déclare " Mon PC fonctionne toujours (environ 10h maintenant) et je n'ai aucune idée de l'exécution ". Je mets toujours le code suivant dans des boucles lors du développement pour avoir une idée de la façon dont les changements semblent affecter la vitesse et aussi pour surveiller le temps qu'il faudra pour terminer.

dayloop2 <- function(temp){
  for (i in 1:nrow(temp)){
    cat(round(i/nrow(temp)*100,2),"%    \r") # prints the percentage complete in realtime.
    # do stuff
  }
  return(blah)
}

Fonctionne également avec lapply.

dayloop2 <- function(temp){
  temp <- lapply(1:nrow(temp), function(i) {
    cat(round(i/nrow(temp)*100,2),"%    \r")
    #do stuff
  })
  return(temp)
}

Si la fonction dans la boucle est assez rapide mais que le nombre de boucles est important, envisagez d'imprimer de temps en temps, car l'impression sur la console elle-même a une surcharge. par exemple

dayloop2 <- function(temp){
  for (i in 1:nrow(temp)){
    if(i %% 100 == 0) cat(round(i/nrow(temp)*100,2),"%    \r") # prints every 100 times through the loop
    # do stuff
  }
  return(temp)
}
débutant
la source
Une option similaire, imprimez la fraction i / n. J'ai toujours quelque chose comme cat(sprintf("\nNow running... %40s, %s/%s \n", nm[i], i, n))puisque je suis généralement en boucle sur des choses nommées (avec des noms nm).
Frank
2

Dans R, vous pouvez souvent accélérer le traitement en boucle en utilisant les applyfonctions de famille (dans votre cas, ce serait probablement le cas replicate). Jetez un œil au plyrpackage qui fournit des barres de progression.

Une autre option consiste à éviter complètement les boucles et à les remplacer par de l'arithmétique vectorisée. Je ne sais pas exactement ce que vous faites, mais vous pouvez probablement appliquer votre fonction à toutes les lignes à la fois:

temp[1:nrow(temp), 10] <- temp[1:nrow(temp), 9] + temp[0:(nrow(temp)-1), 10]

Ce sera beaucoup plus rapide et vous pourrez ensuite filtrer les lignes avec votre condition:

cond.i <- (temp[i, 6] == temp[i-1, 6]) & (temp[i, 3] == temp[i-1, 3])
temp[cond.i, 10] <- temp[cond.i, 9]

L'arithmétique vectorisée nécessite plus de temps et de réflexion sur le problème, mais vous pouvez parfois économiser plusieurs ordres de grandeur en temps d'exécution.

Calimo
la source
14
vous êtes sur le fait que les fonctions vectorielles seront plus rapides que les boucles ou apply () mais ce n'est pas vrai que apply () soit plus rapide que les boucles. Dans de nombreux cas, apply () extrait simplement la boucle de l'utilisateur mais continue de boucler. Voir cette question précédente: stackoverflow.com/questions/2275896/…
JD Long
0

Le traitement avec data.tableest une option viable:

n <- 1000000
df <- as.data.frame(matrix(sample(1:10, n*9, TRUE), n, 9))
colnames(df) <- paste("col", 1:9, sep = "")

library(data.table)

dayloop2.dt <- function(df) {
  dt <- data.table(df)
  dt[, Kumm. := {
    res <- .I;
    ifelse (res > 1,             
      ifelse ((col6 == shift(col6, fill = 0)) & (col3 == shift(col3, fill = 0)) , 
        res <- col9 + shift(res)                   
      , # else
        res <- col9                                 
      )
     , # else
      res <- col9
    )
  }
  ,]
  res <- data.frame(dt)
  return (res)
}

res <- dayloop2.dt(df)

m <- microbenchmark(dayloop2.dt(df), times = 100)
#Unit: milliseconds
#       expr      min        lq     mean   median       uq      max neval
#dayloop2.dt(df) 436.4467 441.02076 578.7126 503.9874 575.9534 966.1042    10

Si vous ignorez les gains possibles du filtrage des conditions, c'est très rapide. Évidemment, si vous pouvez faire le calcul sur le sous-ensemble de données, cela aide.

Bulat
la source
2
Pourquoi répétez-vous la suggestion d'utiliser data.table? Cela a déjà été fait plusieurs fois dans les réponses précédentes.
IRTFM