Concours estival

[Logez Maxime]

Bonjour à tous,

Pour "s'amuser" un peu je vous propose un petit concours. On part tous de l'objet suivant :

Code : Tout sélectionner

set.seed(198)
dta <- as.data.frame(matrix(rnorm(1e6), 1e3))

Et on essaye de calculer le plus rapidement possible la somme en ligne.
L'arbitre étant la fonction microbenchmark du package éponyme.

En plus du côté jeu, je trouve ça intéressant de voir les différentes possibilités, logiques de code, et ça pourrait donné une idée des possibilités pour les gens qui démarrent.

Bon été et bonnes vacances aux chanceux (oui j'en ai besoin aussi ...)
Maxime

[Michaël Delorme]

Le plus évident pour commencer (si j'ai bien compris) :-)

Code : Tout sélectionner

> microbenchmark::microbenchmark(rowSums(dta))
Unit: milliseconds
         expr      min       lq     mean   median       uq      max neval
 rowSums(dta) 7.936299 8.505044 9.234051 8.773643 9.123363 39.63762   100

[Logez Maxime]

re,

presque il s'agissait de la somme en lignes et non en colonne mais le principe y est !
Une autre proposition ?

Maxime

[Serge Rapenne]

allez, à mon tour,
deux pour le prix d'une mais je ne dois pas maitriser la parallélisation car là ce n'est ouvertement pas rentable :

Code : Tout sélectionner

library(microbenchmark)
set.seed(198)
dta <- as.data.frame(matrix(rnorm(1e6), 1e3))
library(parallel)

colSums_par<-function(dta){
  no_cores <- detectCores() - 1
  features <- 1:ncol(dta)
  parts <- split(features, cut(features, no_cores))
  cl <- makeCluster(no_cores, type = "FORK")
 
  tmp<-do.call("c",parLapply(cl = cl, X = 1:n, fun = function(z) colSums(dta[,parts[[z]]])))
  stopCluster(cl)
  tmp
}

microbenchmark(apply(dta,2,sum),
                         colSums(dta),
                         colSums_par(dta))
Unit: milliseconds
               expr       min        lq      mean    median        uq       max neval
 apply(dta, 2, sum)  19.01546  23.76075  34.89021  25.38630  29.12751 614.56589   100
       colSums(dta)  10.94021  12.33058  14.93149  13.10856  15.98018  43.99783   100
   colSums_par(dta) 215.78403 252.59157 305.25291 303.75229 355.96283 441.04322   100


Serge

[Serge Rapenne]

C'est bcp mieux en sortant la mise en oeuvre de la parallelisation de la fonction :

Code : Tout sélectionner

library(microbenchmark)
set.seed(198)
dta <- as.data.frame(matrix(rnorm(1e6), 1e3))
library(parallel)
no_cores <- detectCores() - 1
features <- 1:ncol(dta)
cl <- makeCluster(no_cores, type = "FORK")

colSums_par<-function(dta){
  parts <- split(features, cut(features, no_cores))
do.call("c",parLapply(cl = cl, X = 1:n, fun = function(z) colSums(dta[,parts[[z]]])))
}

microbenchmark(apply(dta,2,sum),colSums(dta),colSums_par(dta))
Unit: milliseconds
               expr      min       lq     mean   median       uq       max neval
 apply(dta, 2, sum) 18.62092 18.81474 23.27708 18.98158 29.53223  43.06510   100
       colSums(dta) 10.37398 10.51229 11.46364 10.62497 10.76803  28.01406   100
   colSums_par(dta) 12.77959 85.57806 85.29144 86.95222 87.96718 168.41382   100

stopCluster(cl)

Serge

[Logez Maxime]

Arff mon idée de départ était de faire la somme en ligne (et pas la somme des lignes), l'équivalent d'un rowSums et non en colonne.
Selon le sens les méthodes les plus optimales changent ...

Mais on peut envisager les deux.

Je trouve intéressant de montrer que dans un cas comme celui-ci ou le temps de calcul est court la parallélisation n'a pas de sens et que le temps de calcul lié au cluster ralentit le calcul initial.

Maxime

[Serge Rapenne]

Une petite précision, le code a été lancé sur une machine virtuelle sous centos 7 avec 4 cœurs. Je n'ai pas le matériel pour tester l'effet du nombre de coeur sur la vitesse de calcul.

Sauf erreur de ma part mon code ne fonctionnera pas sous Windows qui n'accepte pas le type "FORK" pour le cluster

Serge

[Logez Maxime]

Bonjour,

Effectivement sous windows ça ne fonctionne pas, mais on peut utiliser makePSOCKcluster à la place et lancer les calculs via parLapply et compagnie ou via foreach.
Çà reste de toute manière plus lent.
Après pour évaluer le temps de calcul d'une parallélisation c'est un peu compliqué parce que ça dépend de pas mal de paramètres, dont le type de clusters, le type de parallélisation, le nombre de cœurs utilisés et bien entendu du code lancé.
Mais dans un cas comme celui là, à mon avis elle sera toujours plus lente, même si je ne suis pas un expert.

Cordialement,
Maxime

[Serge Rapenne]

Deux nouvelles versions pour la forme. Je n'ai pas fait le benchmark car il bloquait R trop longtemps et j'en ai quand même un peu besoin pour travail ;-)

Uniquement des boucles :

Code : Tout sélectionner

rowSums_loop<-function(dta){
resu <- rep(0.0,nrow(dta))
for (i in seq_along(resu)){
tmp<-0
   for (j in 1:ncoldta)){
   tmp<-tmp+dta[i,j]
   }
mat[i]<-tmp
  }
mat
}

1 boucle

Code : Tout sélectionner

rowSums_loop2<-function(dta){
  resu <- rep(0.0,nrow(dta))
  for (i in seq_along(resu)){
    mat[i]<-sum(dta[i,])
  }
  mat
}

par contre la 1ere version appelle une question

Code : Tout sélectionner

A<-rowSums(dta)
B<-rowSums_loop2(dta)
C<-rowSums_loop(dta)

 identical(A,B)
[1] TRUE
> identical(A,C)
[1] FALSE

les écarts entre A et C sont de l'ordre de 10^(-13) mais d'ou viennent ils ?

Serge

[Bastien Gamboa]

Bonjour,

Ma contribution (sachant que j'ai le même problème avec identical, c'est pour ça que je l'ai changé par all.equal(). Les différences venant du stockage des integers si je ne m'abuse) :

Code : Tout sélectionner

require(microbenchmark)
require(Rcpp)
require(compiler)
set.seed(198)
dta <- matrix(rnorm(1e6), 1e3)
dta2 <- as.data.frame(dta)

# rowSums en C++
cppFunction("NumericVector F1(NumericMatrix w) {
              NumericVector ww=rowSums(w);
              return(ww);
            }")
F1c <- cmpfun(F1) # la même fonction compilée

# double boucle en C++
cppFunction("RObject F2(NumericMatrix x) {
            int nrow = x.nrow(), ncol = x.ncol();
            NumericVector out(nrow);
            for(int i = 0; i < nrow; i++) {
              double total = 0;
              for(int j = 0; j < ncol; j++) {
                total += x(i, j);
              }
              out[i] = total;
            }
            return out;
            }")
F2c <- cmpfun(F2)

F3 <- function(x) rowSums(x)
F3c <- cmpfun(F3)
F4 <- function(x) apply(x, 1, sum)
F4c <- cmpfun(F4)
F5 <- function(x) { SUM <- rep(0, times=nrow(x)) ; for(i in 1:nrow(x)) SUM[i] <- sum(x[i,]) ; SUM }
F5c <- cmpfun(F5)

all.equal(F1(dta), F1c(dta)) # TRUE
all.equal(F1(dta), F2(dta)) # TRUE
all.equal(F1(dta), F2c(dta)) # TRUE
all.equal(F1(dta), F3(dta)) # TRUE
all.equal(F1(dta), F3c(dta)) # TRUE
all.equal(F1(dta), F4(dta)) # TRUE
all.equal(F1(dta), F4c(dta)) # TRUE
all.equal(F1(dta), F5(dta)) # TRUE
all.equal(F1(dta), F5c(dta)) # TRUE

# Idem mais pour colSums
cppFunction("NumericVector F1COL(NumericMatrix w) {
              NumericVector ww=colSums(w);
            return(ww);
            }")
F1cCOL <- cmpfun(F1COL)

cppFunction("RObject F2COL(NumericMatrix x) {
              int nrow = x.nrow(), ncol = x.ncol();
              NumericVector out(ncol);
              for(int i = 0; i < ncol; i++) {
                double total = 0;
                for(int j = 0; j < nrow; j++) {
                  total += x(j, i);
                }
                out[i] = total;
              }
              return out;
            }")
F2cCOL <- cmpfun(F2COL)

F3COL <- function(x) colSums(x)
F3cCOL <- cmpfun(F3COL)
F4COL <- function(x) apply(x, 2, sum)
F4cCOL <- cmpfun(F4COL)
F5COL <- function(x) { SUM <- rep(0, times=ncol(x)) ; for(i in 1:ncol(x)) SUM[i] <- sum(x[,i]) ; SUM }
F5cCOL <- cmpfun(F5COL)

all.equal(F1COL(dta), F1cCOL(dta)) # TRUE
all.equal(F1COL(dta), F2COL(dta)) # TRUE
all.equal(F1COL(dta), F2cCOL(dta)) # TRUE
all.equal(F1COL(dta), F3COL(dta)) # TRUE
all.equal(F1COL(dta), F3cCOL(dta)) # TRUE
all.equal(F1COL(dta), F4COL(dta)) # TRUE
all.equal(F1COL(dta), F4cCOL(dta)) # TRUE
all.equal(F1COL(dta), F5COL(dta)) # TRUE
all.equal(F1COL(dta), F5cCOL(dta)) # TRUE

microbenchmark(times=10,
  # rowSum sur la matrix
  F1(dta),
  F1c(dta),
  F2(dta),
  F2c(dta),
  F3(dta),
  F3c(dta),
  F4(dta),
  F4c(dta),
  F5(dta),
  F5c(dta),
  # rowSum sur le data.frame
  F3(dta2),
  F3c(dta2),
  F4(dta2),
  F4c(dta2),
  F5(dta2),
  F5c(dta2),

  # colSum sur la matrix
  F1COL(dta),
  F1cCOL(dta),
  F2COL(dta),
  F2cCOL(dta),
  F3COL(dta),
  F3cCOL(dta),
  F4COL(dta),
  F4cCOL(dta),
  F5COL(dta),
  F5cCOL(dta),
  # colSum sur le data.frame
  F3COL(dta2),
  F3cCOL(dta2),
  F4COL(dta2),
  F4cCOL(dta2),
  F5COL(dta2),
  F5cCOL(dta2)
)
Unit: microseconds
         expr          min           lq         mean       median           uq          max neval cld
      F1(dta)      756.083      767.570 8.078558e+02 7.749540e+02      778.236     1014.947    10  a
     F1c(dta)      758.954      774.134 8.061326e+02 7.927995e+02      810.235      921.822    10  a
      F2(dta)     2114.403     2174.299 2.323629e+03 2.197274e+03     2348.653     3249.962    10  a
     F2c(dta)     2109.481     2175.940 2.255036e+03 2.189889e+03     2343.731     2541.058    10  a
      F3(dta)     2709.669     2710.900 3.066418e+03 2.757463e+03     3022.276     4689.512    10  a
     F3c(dta)     2705.157     2707.208 2.759760e+03 2.723207e+03     2789.667     2919.304    10  a
      F4(dta)     9532.043     9922.186 1.123185e+04 1.025264e+04    11023.694    17553.973    10  a
     F4c(dta)     9638.707     9905.776 1.078452e+04 1.010433e+04    10345.149    13775.617    10  a
      F5(dta)     8642.632     8915.855 9.432148e+03 9.158720e+03     9440.148    12239.250    10  a
     F5c(dta)     8625.401     8752.988 9.088240e+03 8.959546e+03     9338.407    10288.945    10  a
     F3(dta2)    10619.602    10839.083 1.186966e+04 1.096913e+04    11338.761    15889.609    10  a
    F3c(dta2)    10677.857    10681.959 1.109352e+04 1.114759e+04    11384.299    11645.625    10  a
     F4(dta2)    22271.379    22763.674 2.315345e+04 2.307833e+04    23516.473    23984.152    10  a
    F4c(dta2)    21641.244    22683.676 2.405332e+04 2.326130e+04    23497.602    30801.193    10  a
     F5(dta2) 11754632.897 11774772.644 1.184984e+07 1.182246e+07 11836917.782 12228561.010    10   b
    F5c(dta2) 11585246.034 11718607.646 1.187673e+07 1.174873e+07 11843983.841 12994136.370    10   b
   F1COL(dta)     1357.092     1359.963 1.403162e+03 1.382732e+03     1419.448     1549.907    10  a
  F1cCOL(dta)     1354.630     1355.861 1.370055e+03 1.360374e+03     1396.064     1401.398    10  a
   F2COL(dta)     1309.914     1316.888 1.471837e+03 1.319349e+03     1401.398     2561.571    10  a
  F2cCOL(dta)     1316.067     1317.298 1.370670e+03 1.357502e+03     1419.038     1481.806    10  a
   F3COL(dta)     1013.306     1018.229 1.041449e+03 1.022332e+03     1056.382     1135.560    10  a
  F3cCOL(dta)     1016.178     1021.101 1.038823e+03 1.038331e+03     1058.023     1061.305    10  a
   F4COL(dta)     7704.812     7759.374 8.474267e+03 8.189721e+03     8636.479    11579.165    10  a
  F4cCOL(dta)     7667.889     8165.927 9.053492e+03 8.182747e+03     9135.746    12186.738    10  a
   F5COL(dta)     3797.228     3972.813 8.890666e+03 4.103477e+03     7522.253    42051.341    10  a
  F5cCOL(dta)     3705.334     3949.840 5.076372e+03 4.200910e+03     7140.726     7804.911    10  a
  F3COL(dta2)     8666.836     8862.523 9.214062e+03 9.152361e+03     9337.587    10117.052    10  a
 F3cCOL(dta2)     8579.454     9118.926 9.348130e+03 9.258205e+03     9817.574    10078.899    10  a
  F4COL(dta2)    15607.361    16549.694 1.795015e+04 1.683420e+04    18982.036    22909.311    10  a
 F4cCOL(dta2)    15852.687    16052.067 1.788619e+04 1.695337e+04    19903.446    21275.715    10  a
  F5COL(dta2)     8867.036     9167.335 9.520679e+03 9.388662e+03     9612.041    11065.128    10  a
 F5cCOL(dta2)     8837.088     8947.444 9.423246e+03 9.537582e+03     9669.475    10089.566    10  a

Je n'ai fait que 10 itérations car ça prend déjà pas mal de minutes sur mon petit pc portable ! En faisant plus d'itérations, les séparations entre fonctions sont plus fiables.
Je trouve étonnant que pour le data.frame la boucle en R soit rapide en colonne et très lente en ligne ! Ca viendrait du fait que le data.frame est une liste et qu'il serait plus facile d'extraire une colonne qu'une ligne ?

HTH,
Bastien

[François Bonnot]

Bonjour,
Une solution:

Code : Tout sélectionner

microbenchmark(as.matrix(dta) %*% rep(1,ncol(dta)))

[Logez Maxime]

Bonjour,

mes propositions :

Code : Tout sélectionner

set.seed(198)
dta <- as.data.frame(matrix(rnorm(1e6), 1e3))
fun <- function(x) {
  res <- numeric(nrow(x))
  for (i in seq_len(ncol(x))) {
    res <- res + x[[i]]
    }
  res
}
fun <- compiler:::cmpfun(fun) # pas nécessaire dans les dernières version de R déjà fait

fun(dta)
fun2 <- function(x) {
  res <- numeric(nrow(x))
  for (i in x) {
    res <- res + i
    }
  res
}
fun2 <- compiler:::cmpfun(fun2)
fun2(dta)

Reduce("+", dta)

tapply(unlist(dta), rep(1:nrow(dta), ncol(dta)), sum)

microbenchmark(fun(dta), fun2(dta), Reduce("+", dta),
  tapply(unlist(dta), rep(1:nrow(dta), ncol(dta)), sum),
  rowSums(dta), map(dta, sum),
  times = 1e3)


Pour les sommes en colonne :

Code : Tout sélectionner

fun2 <- function(x) {
  res <- numeric(ncol(x))
  for (i in seq_along(res))
    res[i] <- sum(i)
  res
}
fun2 <- compiler:::cmpfun(fun2)

require(purrr)
map(dta, sum)

microbenchmark(fun2(dta), colSums(dta),
  map(dta, sum), lapply(dta, sum),
  sapply(dta, sum), vapply(dta, sum, numeric(1)),
  colSums(dta),
  times = 1e3)


Maxime

[François Bonnot]

Hum... il semble difficile de battre fun2 (qui arrive en tête) et Reduce, parce que ces deux fonctions travaillent directement sur les colonnes du data frame avec un minimum d'opérations.

[Logez Maxime]

re,

Pour la parallélisation en ligne j'ai essayé avec le package doSNOW et foreach (pour un exemple reproductible aussi sous Windows) et ma fois c'est une catastrophe dans un cas comme celui-là :

Code : Tout sélectionner

require(doSNOW)
cl <- makeSOCKcluster(5)
registerDoSNOW(cl)

microbenchmark(
  foreach(i=iter(dta, "row")) %dopar% {
    sum(i)
  }, times = 100) 

Unit: seconds
                                                 expr      min       lq     mean   median       uq      max neval
 foreach(i = iter(dta, "row")) %dopar% {     sum(i) } 6.683085 6.970803 7.156997 7.135713 7.309011 7.934822   100


Pour la parallélisation en colonne avec mclapply (sous linux donc), plus le nombre de cœurs augmente et plus le temps de calcul est long. Encore une fois surement parce que le calcul ici est très rapide et que c'est le temps de communication qui est limitant.

Cordialement,
Maxime

[Eric Wajnberg]

Je n'ai pas le temps de m'y coller en ce moment, mais si je devais résoudre ce problème, je coderais ça en C est l’intégrerait dans R avec la fonction .C(). Je mets ma main à couper que ce bâterait toutes les autres propositions utilisant du R natif. C'est autorisé dans le concours ?

Eric.