La famille «* apply» n'est-elle vraiment pas vectorisée?

138

Nous avons donc l'habitude de dire à chaque nouvel utilisateur R que " applyn'est pas vectorisé, regardez le Patrick Burns R Inferno Circle 4 " qui dit (je cite):

Un réflexe courant consiste à utiliser une fonction de la famille apply. Ce n'est pas de la vectorisation, c'est un masquage de boucle . La fonction apply a une boucle for dans sa définition. La fonction lapply enterre la boucle, mais les temps d'exécution ont tendance à être à peu près égaux à une boucle for explicite.

En effet, un rapide coup d'œil sur le applycode source révèle la boucle:

grep("for", capture.output(getAnywhere("apply")), value = TRUE)
## [1] "        for (i in 1L:d2) {"  "    else for (i in 1L:d2) {"

Ok jusqu'à présent, mais un regard lapplyou vapplyrévèle en fait une image complètement différente:

lapply
## function (X, FUN, ...) 
## {
##     FUN <- match.fun(FUN)
##     if (!is.vector(X) || is.object(X)) 
##        X <- as.list(X)
##     .Internal(lapply(X, FUN))
## }
## <bytecode: 0x000000000284b618>
## <environment: namespace:base>

Donc, apparemment, il n'y a pas de forboucle R cachée là-bas, ils appellent plutôt une fonction écrite interne en C.

Un rapide coup d' oeil dans le lapin trou révèle à peu près la même image

De plus, prenons la colMeansfonction par exemple, qui n'a jamais été accusée de ne pas être vectorisée

colMeans
# function (x, na.rm = FALSE, dims = 1L) 
# {
#   if (is.data.frame(x)) 
#     x <- as.matrix(x)
#   if (!is.array(x) || length(dn <- dim(x)) < 2L) 
#     stop("'x' must be an array of at least two dimensions")
#   if (dims < 1L || dims > length(dn) - 1L) 
#     stop("invalid 'dims'")
#   n <- prod(dn[1L:dims])
#   dn <- dn[-(1L:dims)]
#   z <- if (is.complex(x)) 
#     .Internal(colMeans(Re(x), n, prod(dn), na.rm)) + (0+1i) * 
#     .Internal(colMeans(Im(x), n, prod(dn), na.rm))
#   else .Internal(colMeans(x, n, prod(dn), na.rm))
#   if (length(dn) > 1L) {
#     dim(z) <- dn
#     dimnames(z) <- dimnames(x)[-(1L:dims)]
#   }
#   else names(z) <- dimnames(x)[[dims + 1]]
#   z
# }
# <bytecode: 0x0000000008f89d20>
#   <environment: namespace:base>

Hein? Il appelle aussi juste .Internal(colMeans(...que nous pouvons également trouver dans le terrier du lapin . Alors, en quoi est-ce différent .Internal(lapply(..?

En fait, un benchmark rapide révèle que les sapplyperformances ne sont pas pires colMeanset bien meilleures qu'une forboucle pour un ensemble de données volumineuses

m <- as.data.frame(matrix(1:1e7, ncol = 1e5))
system.time(colMeans(m))
# user  system elapsed 
# 1.69    0.03    1.73 
system.time(sapply(m, mean))
# user  system elapsed 
# 1.50    0.03    1.60 
system.time(apply(m, 2, mean))
# user  system elapsed 
# 3.84    0.03    3.90 
system.time(for(i in 1:ncol(m)) mean(m[, i]))
# user  system elapsed 
# 13.78    0.01   13.93 

En d'autres termes, est-il correct de dire cela lapplyet vapply sont en fait vectorisés (par rapport à applylaquelle est une forboucle qui appelle également lapply) et que voulait vraiment dire Patrick Burns?

David Arenburg
la source
8
Tout cela est dans la sémantique, mais je ne les considérerais pas comme vectorisés. Je considère une approche vectorisée si une fonction R est appelée une seule fois et peut se voir passer un vecteur de valeurs. *applyles fonctions appellent à plusieurs reprises les fonctions R, ce qui les rend boucles. En ce qui concerne les bonnes performances de sapply(m, mean): Peut-être le code C de la lapplyméthode envoie-t-il une seule fois et appelle-t-il ensuite la méthode à plusieurs reprises? mean.defaultest assez optimisé.
Roland
4
Excellente question, et merci de vérifier le code sous-jacent. Je cherchais s'il avait été récemment modifié, mais rien à ce sujet dans les notes de version R à partir de la version 2.13.0.
ilir
1
Dans quelle mesure les performances dépendent-elles à la fois de la plate-forme et des indicateurs du compilateur C et de l'éditeur de liens utilisés?
smci
3
@DavidArenburg En fait, je ne pense pas que ce soit bien défini. Au moins, je ne connais pas de référence canonique. La définition du langage mentionne les opérations "vectorisées", mais ne définit pas la vectorisation.
Roland du
3
Très lié: la famille Apply de R est-elle plus que le sucre syntaxique? (Et, comme ces réponses, aussi une bonne lecture.)
Gregor Thomas

Réponses:

73

Tout d'abord, dans votre exemple vous faites des tests sur un "data.frame" qui n'est pas juste pour colMeans, applyet "[.data.frame"puisqu'ils ont une surcharge:

system.time(as.matrix(m))  #called by `colMeans` and `apply`
#   user  system elapsed 
#   1.03    0.00    1.05
system.time(for(i in 1:ncol(m)) m[, i])  #in the `for` loop
#   user  system elapsed 
#  12.93    0.01   13.07

Sur une matrice, l'image est un peu différente:

mm = as.matrix(m)
system.time(colMeans(mm))
#   user  system elapsed 
#   0.01    0.00    0.01 
system.time(apply(mm, 2, mean))
#   user  system elapsed 
#   1.48    0.03    1.53 
system.time(for(i in 1:ncol(mm)) mean(mm[, i]))
#   user  system elapsed 
#   1.22    0.00    1.21

En reprenant la partie principale de la question, la principale différence entre lapply/ mapply/ etc et les boucles R simples est l'endroit où la boucle est effectuée. Comme le note Roland, les boucles C et R doivent évaluer une fonction R à chaque itération, ce qui est le plus coûteux. Les fonctions C vraiment rapides sont celles qui font tout en C, donc, je suppose, cela devrait être ce que signifie "vectorisé"?

Un exemple où nous trouvons la moyenne dans chacun des éléments d'une "liste":

( EDIT 11 mai 16 : Je crois que l'exemple de recherche de la "moyenne" n'est pas une bonne configuration pour les différences entre l'évaluation d'une fonction R de manière itérative et le code compilé, (1) en raison de la particularité de l'algorithme de moyenne de R sur "numérique" s sur un simple sum(x) / length(x)et (2), il devrait avoir plus de sens de tester sur les «listes» avec length(x) >> lengths(x). Ainsi, l'exemple «moyen» est déplacé à la fin et remplacé par un autre.)

A titre d'exemple simple, nous pourrions considérer la découverte de l'opposé de chaque length == 1élément d'une "liste":

Dans un tmp.cfichier:

#include <R.h>
#define USE_RINTERNALS 
#include <Rinternals.h>
#include <Rdefines.h>

/* call a C function inside another */
double oppC(double x) { return(ISNAN(x) ? NA_REAL : -x); }
SEXP sapply_oppC(SEXP x)
{
    SEXP ans = PROTECT(allocVector(REALSXP, LENGTH(x)));
    for(int i = 0; i < LENGTH(x); i++) 
        REAL(ans)[i] = oppC(REAL(VECTOR_ELT(x, i))[0]);

    UNPROTECT(1);
    return(ans);
}

/* call an R function inside a C function;
 * will be used with 'f' as a closure and as a builtin */    
SEXP sapply_oppR(SEXP x, SEXP f)
{
    SEXP call = PROTECT(allocVector(LANGSXP, 2));
    SETCAR(call, install(CHAR(STRING_ELT(f, 0))));

    SEXP ans = PROTECT(allocVector(REALSXP, LENGTH(x)));     
    for(int i = 0; i < LENGTH(x); i++) { 
        SETCADR(call, VECTOR_ELT(x, i));
        REAL(ans)[i] = REAL(eval(call, R_GlobalEnv))[0];
    }

    UNPROTECT(2);
    return(ans);
}

Et en face R:

system("R CMD SHLIB /home/~/tmp.c")
dyn.load("/home/~/tmp.so")

avec des données:

set.seed(007)
myls = rep_len(as.list(c(NA, runif(3))), 1e7)

#a closure wrapper of `-`
oppR = function(x) -x

for_oppR = compiler::cmpfun(function(x, f)
{
    f = match.fun(f)  
    ans = numeric(length(x))
    for(i in seq_along(x)) ans[[i]] = f(x[[i]])
    return(ans)
})

Analyse comparative:

#call a C function iteratively
system.time({ sapplyC =  .Call("sapply_oppC", myls) }) 
#   user  system elapsed 
#  0.048   0.000   0.047 

#evaluate an R closure iteratively
system.time({ sapplyRC =  .Call("sapply_oppR", myls, "oppR") }) 
#   user  system elapsed 
#  3.348   0.000   3.358 

#evaluate an R builtin iteratively
system.time({ sapplyRCprim =  .Call("sapply_oppR", myls, "-") }) 
#   user  system elapsed 
#  0.652   0.000   0.653 

#loop with a R closure
system.time({ forR = for_oppR(myls, "oppR") })
#   user  system elapsed 
#  4.396   0.000   4.409 

#loop with an R builtin
system.time({ forRprim = for_oppR(myls, "-") })
#   user  system elapsed 
#  1.908   0.000   1.913 

#for reference and testing 
system.time({ sapplyR = unlist(lapply(myls, oppR)) })
#   user  system elapsed 
#  7.080   0.068   7.170 
system.time({ sapplyRprim = unlist(lapply(myls, `-`)) }) 
#   user  system elapsed 
#  3.524   0.064   3.598 

all.equal(sapplyR, sapplyRprim)
#[1] TRUE 
all.equal(sapplyR, sapplyC)
#[1] TRUE
all.equal(sapplyR, sapplyRC)
#[1] TRUE
all.equal(sapplyR, sapplyRCprim)
#[1] TRUE
all.equal(sapplyR, forR)
#[1] TRUE
all.equal(sapplyR, forRprim)
#[1] TRUE

(Suit l'exemple original de la recherche moyenne):

#all computations in C
all_C = inline::cfunction(sig = c(R_ls = "list"), body = '
    SEXP tmp, ans;
    PROTECT(ans = allocVector(REALSXP, LENGTH(R_ls)));

    double *ptmp, *pans = REAL(ans);

    for(int i = 0; i < LENGTH(R_ls); i++) {
        pans[i] = 0.0;

        PROTECT(tmp = coerceVector(VECTOR_ELT(R_ls, i), REALSXP));
        ptmp = REAL(tmp);

        for(int j = 0; j < LENGTH(tmp); j++) pans[i] += ptmp[j];

        pans[i] /= LENGTH(tmp);

        UNPROTECT(1);
    }

    UNPROTECT(1);
    return(ans);
')

#a very simple `lapply(x, mean)`
C_and_R = inline::cfunction(sig = c(R_ls = "list"), body = '
    SEXP call, ans, ret;

    PROTECT(call = allocList(2));
    SET_TYPEOF(call, LANGSXP);
    SETCAR(call, install("mean"));

    PROTECT(ans = allocVector(VECSXP, LENGTH(R_ls)));
    PROTECT(ret = allocVector(REALSXP, LENGTH(ans)));

    for(int i = 0; i < LENGTH(R_ls); i++) {
        SETCADR(call, VECTOR_ELT(R_ls, i));
        SET_VECTOR_ELT(ans, i, eval(call, R_GlobalEnv));
    }

    double *pret = REAL(ret);
    for(int i = 0; i < LENGTH(ans); i++) pret[i] = REAL(VECTOR_ELT(ans, i))[0];

    UNPROTECT(3);
    return(ret);
')                    

R_lapply = function(x) unlist(lapply(x, mean))                       

R_loop = function(x) 
{
    ans = numeric(length(x))
    for(i in seq_along(x)) ans[i] = mean(x[[i]])
    return(ans)
} 

R_loopcmp = compiler::cmpfun(R_loop)


set.seed(007); myls = replicate(1e4, runif(1e3), simplify = FALSE)
all.equal(all_C(myls), C_and_R(myls))
#[1] TRUE
all.equal(all_C(myls), R_lapply(myls))
#[1] TRUE
all.equal(all_C(myls), R_loop(myls))
#[1] TRUE
all.equal(all_C(myls), R_loopcmp(myls))
#[1] TRUE

microbenchmark::microbenchmark(all_C(myls), 
                               C_and_R(myls), 
                               R_lapply(myls), 
                               R_loop(myls), 
                               R_loopcmp(myls), 
                               times = 15)
#Unit: milliseconds
#            expr       min        lq    median        uq      max neval
#     all_C(myls)  37.29183  38.19107  38.69359  39.58083  41.3861    15
#   C_and_R(myls) 117.21457 123.22044 124.58148 130.85513 169.6822    15
#  R_lapply(myls)  98.48009 103.80717 106.55519 109.54890 116.3150    15
#    R_loop(myls) 122.40367 130.85061 132.61378 138.53664 178.5128    15
# R_loopcmp(myls) 105.63228 111.38340 112.16781 115.68909 128.1976    15
alexis_laz
la source
10
Excellent point sur les coûts de conversion du data.frame en matrice, et merci de fournir des repères.
Joshua Ulrich
C'est une très bonne réponse, même si je n'ai pas pu compiler vos fonctions all_Cet C_and_R. J'ai également trouvé dans les documentations d' compiler::cmpfunune ancienne version R de lapply qui contient une véritable forboucle R , je commence à soupçonner que Burns faisait référence à cette ancienne version qui a été vectorisée depuis et c'est la vraie réponse à ma question. ..
David Arenburg
@DavidArenburg: l'analyse comparative la1de ?compiler::cmpfunsemble, toujours, donner la même efficacité avec tout sauf les all_Cfonctions. Je suppose que cela - en fait - est une question de définition; "vectorisé" signifie-t-il toute fonction qui accepte non seulement des scalaires, toute fonction qui a du code C, toute fonction qui utilise des calculs en C seulement?
alexis_laz
1
Je suppose que toutes les fonctions de R contiennent du code C, simplement parce que tout dans R est une fonction (qui devait être écrite dans une langue). Donc, fondamentalement, si je comprends bien, vous dites que ce lapplyn'est pas vectorisé simplement parce qu'il évalue une fonction R à chaque itération avec son code C?
David Arenburg
5
@DavidArenburg: Si je dois définir la "vectorisation" d'une certaine manière, je suppose que je choisirais une approche linguistique; c'est à dire une fonction qui accepte et sait gérer un "vecteur", qu'il soit rapide, lent, écrit en C, en R ou autre. Dans R, l'importance de la vectorisation réside dans le fait que de nombreuses fonctions sont écrites en C et gèrent des vecteurs alors que dans d'autres langues, les utilisateurs bouclent généralement sur l'entrée pour -eg- trouver la moyenne. Cela fait que la vectorisation est liée, indirectement, à la vitesse, l'efficacité, la sécurité et la robustesse.
alexis_laz
65

Pour moi, la vectorisation consiste avant tout à rendre votre code plus facile à écrire et à comprendre.

Le but d'une fonction vectorisée est d'éliminer la comptabilité associée à une boucle for. Par exemple, au lieu de:

means <- numeric(length(mtcars))
for (i in seq_along(mtcars)) {
  means[i] <- mean(mtcars[[i]])
}
sds <- numeric(length(mtcars))
for (i in seq_along(mtcars)) {
  sds[i] <- sd(mtcars[[i]])
}

Tu peux écrire:

means <- vapply(mtcars, mean, numeric(1))
sds   <- vapply(mtcars, sd, numeric(1))

Cela permet de voir plus facilement ce qui est identique (les données d'entrée) et ce qui est différent (la fonction que vous appliquez).

Un avantage secondaire de la vectorisation est que la boucle for est souvent écrite en C plutôt qu'en R. Cela présente des avantages substantiels en termes de performances, mais je ne pense pas que ce soit la propriété clé de la vectorisation. La vectorisation consiste fondamentalement à sauver votre cerveau, pas à sauver le travail de l'ordinateur.

hadley
la source
5
Je ne pense pas qu'il y ait une différence de performance significative entre les forboucles C et R. OK, une boucle C peut être optimisée par le compilateur, mais le point principal pour les performances est de savoir si le contenu de la boucle est efficace. Et évidemment, le code compilé est généralement plus rapide que le code interprété. Mais c'est probablement ce que vous vouliez dire.
Roland
3
@Roland ouais, ce n'est pas la boucle for en elle-même, c'est tout ce qui l'entoure (le coût d'un appel de fonction, la possibilité de modifier sur place, ...).
hadley
10
@DavidArenburg "La cohérence inutile est le hobgoblin des petits esprits";)
hadley
6
Non, je ne pense pas que la performance soit le point principal de la vectorisation de votre code. Réécrire une boucle dans un lapply est bénéfique même si ce n'est pas plus rapide. Le point principal de dplyr est qu'il facilite l'expression de la manipulation des données (et c'est juste très bien que ce soit aussi rapide).
hadley
12
@DavidArenburg c'est parce que vous êtes un utilisateur R expérimenté. La plupart des nouveaux utilisateurs trouvent les boucles beaucoup plus naturelles et doivent être encouragés à vectoriser. Pour moi, utiliser une fonction comme colMeans n'est pas nécessairement une question de vectorisation, il s'agit de réutiliser du code rapide que quelqu'un a déjà écrit
hadley
49

Je suis d'accord avec le point de vue de Patrick Burns selon lequel il s'agit plutôt de masquage de boucles et non de vectorisation de code . Voici pourquoi:

Considérez cet Cextrait de code:

for (int i=0; i<n; i++)
  c[i] = a[i] + b[i]

Ce que nous aimerions faire est assez clair. Mais comment la tâche est effectuée ou comment elle pourrait être exécutée n'est pas vraiment. Une boucle for par défaut est une construction série. Il n'indique pas si ou comment les choses peuvent être faites en parallèle.

Le moyen le plus évident est que le code est exécuté de manière séquentielle . Chargez a[i]et continuez b[i]dans les registres, ajoutez-les, stockez le résultat dans c[i]et faites-le pour chacun i.

Cependant, les processeurs modernes ont un jeu d' instructions vectorielles ou SIMD qui est capable de fonctionner sur un vecteur de données pendant la même instruction lors de l'exécution de la même opération (par exemple, en ajoutant deux vecteurs comme indiqué ci-dessus). Selon le processeur / l'architecture, il peut être possible d'ajouter, disons, quatre nombres à partir de aet bsous la même instruction, au lieu d'un à la fois.

Nous aimerions exploiter les données multiples à instruction unique et effectuer un parallélisme au niveau des données , c'est-à-dire charger 4 choses à la fois, ajouter 4 choses à la fois, stocker 4 choses à la fois, par exemple. Et c'est la vectorisation du code .

Notez que ceci est différent de la parallélisation de code - où plusieurs calculs sont effectués simultanément.

Ce serait formidable si le compilateur identifie ces blocs de code et les vectorise automatiquement , ce qui est une tâche difficile. La vectorisation automatique du code est un sujet de recherche difficile en informatique. Mais au fil du temps, les compilateurs se sont améliorés. Vous pouvez vérifier les capacités de vectorisation automatique d' GNU-gcc ici . De même pour LLVM-clang ici . Et vous pouvez également trouver quelques points de repère dans le dernier lien par rapport à gccet ICC(compilateur Intel C ++).

gcc(I'm on v4.9) par exemple, ne vectorise pas le code automatiquement au -O2niveau de l'optimisation. Donc, si nous devions exécuter le code ci-dessus, il serait exécuté séquentiellement. Voici le moment pour ajouter deux vecteurs entiers d'une longueur de 500 millions.

Nous devons soit ajouter le drapeau, -ftree-vectorizesoit changer l'optimisation au niveau -O3. (Notez que cela -O3effectue également d' autres optimisations supplémentaires ). Le drapeau -fopt-info-vecest utile car il informe quand une boucle a été vectorisée avec succès).

# compiling with -O2, -ftree-vectorize and  -fopt-info-vec
# test.c:32:5: note: loop vectorized
# test.c:32:5: note: loop versioned for vectorization because of possible aliasing
# test.c:32:5: note: loop peeled for vectorization to enhance alignment    

Cela nous indique que la fonction est vectorisée. Voici les timings comparant les versions non vectorisées et vectorisées sur des vecteurs entiers de longueur 500 millions:

x = sample(100L, 500e6L, TRUE)
y = sample(100L, 500e6L, TRUE)
z = vector("integer", 500e6L) # result vector

# non-vectorised, -O2
system.time(.Call("Csum", x, y, z))
#    user  system elapsed 
#   1.830   0.009   1.852

# vectorised using flags shown above at -O2
system.time(.Call("Csum", x, y, z))
#    user  system elapsed 
#   0.361   0.001   0.362

# both results are checked for identicalness, returns TRUE

Cette partie peut être ignorée en toute sécurité sans perdre la continuité.

Les compilateurs n'auront pas toujours suffisamment d'informations pour vectoriser. Nous pourrions utiliser la spécification OpenMP pour la programmation parallèle , qui fournit également une directive de compilateur simd pour demander aux compilateurs de vectoriser le code. Il est essentiel de s'assurer qu'il n'y a pas de chevauchements de mémoire, de conditions de course, etc. lors de la vectorisation manuelle du code, sinon cela entraînera des résultats incorrects.

#pragma omp simd
for (i=0; i<n; i++) 
  c[i] = a[i] + b[i]

En faisant cela, nous demandons spécifiquement au compilateur de le vectoriser quoi qu'il arrive. Nous devrons activer les extensions OpenMP en utilisant l'indicateur de temps de compilation -fopenmp. En faisant cela:

# timing with -O2 + OpenMP with simd
x = sample(100L, 500e6L, TRUE)
y = sample(100L, 500e6L, TRUE)
z = vector("integer", 500e6L) # result vector
system.time(.Call("Cvecsum", x, y, z))
#    user  system elapsed 
#   0.360   0.001   0.360

qui est genial! Cela a été testé avec gcc v6.2.0 et llvm clang v3.9.0 (tous deux installés via homebrew, MacOS 10.12.3), tous deux prenant en charge OpenMP 4.0.


En ce sens, même si la page Wikipedia sur Array Programming mentionne que les langages qui opèrent sur des tableaux entiers appellent généralement cela comme des opérations vectorisées , c'est vraiment une boucle qui cache l' OMI (à moins qu'elle ne soit réellement vectorisée).

Dans le cas de R, pair rowSums()ou colSums()code en C n'exploitent pas la vectorisation de code IIUC; c'est juste une boucle en C. Il en va de même lapply(). Dans le cas de apply(), c'est dans R. Tous ces éléments se cachent donc en boucle .

En bref, envelopper une fonction R par:

écrire simplement une boucle for dans C! = vectoriser votre code.
écrire simplement une boucle for dans R! = vectoriser votre code.

Intel Math Kernel Library (MKL) par exemple implémente des formes vectorisées de fonctions.

HTH


Références:

  1. Conférence de James Reinders, Intel (cette réponse est surtout une tentative de résumer cet excellent discours)
Arun
la source
35

Donc, pour résumer les bonnes réponses / commentaires en une réponse générale et fournir un arrière-plan: R a 4 types de boucles ( dans l'ordre non vectorisé à vectorisé )

  1. forBoucle R qui appelle à plusieurs reprises des fonctions R dans chaque itération ( non vectorisée )
  2. Boucle C qui appelle à plusieurs reprises des fonctions R dans chaque itération ( non vectorisée )
  3. Boucle C qui n'appelle la fonction R qu'une seule fois (un peu vectorisée )
  4. Une boucle C simple qui n'appelle aucune fonction R du tout et utilise ses propres fonctions compilées ( vectorisées )

La *applyfamille est donc le deuxième type. Sauf applyqui est plus du premier type

Vous pouvez comprendre cela à partir du commentaire dans son code source

/ * .Interne (lapply (X, FUN)) * /

/ * Ceci est un .Internal spécial, donc des arguments non évalués. Il est
appelé à partir d'un wrapper de fermeture, donc X et FUN sont des promesses. FUN doit être non évalué pour être utilisé par exemple dans bquote. * /

Cela signifie que le lapplycode C accepte une fonction non évaluée de R et l'évalue plus tard dans le code C lui-même. C'est fondamentalement la différence entre l' appel de lapplys.Internal

.Internal(lapply(X, FUN))

Qui a un FUNargument qui contient une fonction R

Et l' colMeans .Internalappel qui n'a pas d' FUNargument

.Internal(colMeans(Re(x), n, prod(dn), na.rm))

colMeans, contrairement à lapplysait exactement quelle fonction il doit utiliser, il calcule donc la moyenne en interne dans le code C.

Vous pouvez clairement voir le processus d'évaluation de la fonction R à chaque itération dans le lapplycode C

 for(R_xlen_t i = 0; i < n; i++) {
      if (realIndx) REAL(ind)[0] = (double)(i + 1);
      else INTEGER(ind)[0] = (int)(i + 1);
      tmp = eval(R_fcall, rho);   // <----------------------------- here it is
      if (MAYBE_REFERENCED(tmp)) tmp = lazy_duplicate(tmp);
      SET_VECTOR_ELT(ans, i, tmp);
   }

Pour résumer, lapplyn'est pas vectorisé , bien qu'il présente deux avantages possibles par rapport à la forboucle R simple

  1. L'accès et l'affectation dans une boucle semblent être plus rapides en C (c'est-à-dire dans lapplyune fonction). Bien que la différence semble grande, nous restons toujours au niveau de la microseconde et le plus coûteux est la valorisation d'une fonction R à chaque itération. Un exemple simple:

    ffR = function(x)  {
        ans = vector("list", length(x))
        for(i in seq_along(x)) ans[[i]] = x[[i]]
        ans 
    }
    
    ffC = inline::cfunction(sig = c(R_x = "data.frame"), body = '
        SEXP ans;
        PROTECT(ans = allocVector(VECSXP, LENGTH(R_x)));
        for(int i = 0; i < LENGTH(R_x); i++) 
               SET_VECTOR_ELT(ans, i, VECTOR_ELT(R_x, i));
        UNPROTECT(1);
        return(ans); 
    ')
    
    set.seed(007) 
    myls = replicate(1e3, runif(1e3), simplify = FALSE)     
    mydf = as.data.frame(myls)
    
    all.equal(ffR(myls), ffC(myls))
    #[1] TRUE 
    all.equal(ffR(mydf), ffC(mydf))
    #[1] TRUE
    
    microbenchmark::microbenchmark(ffR(myls), ffC(myls), 
                                   ffR(mydf), ffC(mydf),
                                   times = 30)
    #Unit: microseconds
    #      expr       min        lq    median        uq       max neval
    # ffR(myls)  3933.764  3975.076  4073.540  5121.045 32956.580    30
    # ffC(myls)    12.553    12.934    16.695    18.210    19.481    30
    # ffR(mydf) 14799.340 15095.677 15661.889 16129.689 18439.908    30
    # ffC(mydf)    12.599    13.068    15.835    18.402    20.509    30
  2. Comme mentionné par @Roland, il exécute une boucle C compilée plutôt qu'une boucle R interprétée


Cependant, lors de la vectorisation de votre code, vous devez prendre en compte certains éléments.

  1. Si votre ensemble de données (appel Let It df) est de classe data.frame, certaines fonctions vectorisés ( par exemple colMeans, colSums, rowSums, etc.) devront convertir en une matrice d' abord, simplement parce que c'est la façon dont ils ont été conçus. Cela signifie que pour un gros, dfcela peut créer une surcharge énorme. Bien que lapplycela ne soit pas nécessaire car il extrait les vecteurs réels de df(comme data.framec'est juste une liste de vecteurs) et donc, si vous n'avez pas autant de colonnes mais de nombreuses lignes, lapply(df, mean)peut parfois être une meilleure option que colMeans(df).
  2. Une autre chose à retenir est que R a une grande variété de types de fonctions différents, tels que .Primitive, et générique ( S3, S4) voir ici pour quelques informations supplémentaires. La fonction générique doit faire un envoi de méthode qui est parfois une opération coûteuse. Par exemple, meanest une S3fonction générique tandis que sumest Primitive. Ainsi, certains moments lapply(df, sum)pourraient être très efficaces par colSumsrapport aux raisons énumérées ci-dessus
David Arenburg
la source
1
Résumé très cohérent. Juste quelques notes: (1) C sait comment gérer les "data.frame", car ce sont des "list" avec des attributs; c'est ça colMeansetc. qui sont construits pour ne gérer que des matrices. (2) Je suis un peu confus par votre troisième catégorie; Je ne peux pas dire à quoi -exaclty- vous parlez. (3) Puisque vous faites allusion spécifiquement à lapply, je pense que cela ne fait pas de différence entre "[<-"R et C; ils pré-allouent tous les deux une "liste" (un SEXP) et la remplissent à chaque itération ( SET_VECTOR_ELTen C), sauf si je manque votre point.
alexis_laz
2
Je comprends votre point de vue do.callen ce qu'il construit un appel de fonction dans l'environnement C et l'évalue simplement; même si j'ai du mal à le comparer à la boucle ou à la vectorisation car il fait une chose différente. Vous avez, en fait, raison d'accéder et d'attribuer des différences entre C et R, bien que les deux restent au niveau de la microseconde et n'affectent pas énormément le résultat, le résultat étant énormément coûteux, l'appel itératif de la fonction R (comparer R_loopet R_lapplydans ma réponse ). (Je vais éditer votre message avec un repère; j'espère que cela ne vous dérangera toujours pas)
alexis_laz
2
Je n'essaye pas d'être en désaccord - et je suis confus, honnêtement, sur ce avec quoi vous n'êtes pas d'accord. Mon commentaire précédent aurait pu être mieux formulé. J'essaie d'affiner la terminologie utilisée parce que le terme «vectorisation» a deux définitions qui sont souvent confondues. Je ne pense pas que ce soit discutable. Burns et vous semblez vouloir l'utiliser uniquement dans le sens de l'implémentation, mais Hadley et de nombreux membres de R-Core (en prenant Vectorize()comme exemple) l'utilisent également dans le sens de l'interface utilisateur. Je pense qu'une grande partie du désaccord dans ce fil est causée par l'utilisation d'un terme pour deux concepts distincts mais liés.
Gregor Thomas
3
@DavidArenburg et n'est-ce pas une vectorisation au sens de l'interface utilisateur, qu'il y ait une boucle for dans R ou C en dessous?
Gregor Thomas
2
@DavidArenburg, Gregor, je pense que la confusion est entre "vectorisation de code" et "fonctions vectorisées". En R, l'usage semble incliné vers ce dernier. La "vectorisation de code" décrit le fonctionnement sur un vecteur de longueur 'k' dans la même instruction. Emballage d'un fn. autour du code en boucle aboutit à des «fonctions vectorisées» (oui, cela n'a pas de sens et est déroutant, je suis d'accord, mieux serait le masquage de boucle ou les fonctions vectorielles i / p ) et n'a rien à voir avec la vectorisation du code . Dans R, apply serait une fonction vectorisée , mais elle ne vectorise pas votre code, mais opère plutôt sur des vecteurs.
Arun