Discussion :

[Christophe Genolini]

Bonjour à tous,

Il est clair que R n'est pas très rapide, mais étant donnée notre objectif, il faudrait le montrer de manière un peu plus formelle. Pour cela, je pense qu'il faudrait comparer les performances relatives des différents langages. Un moyen simple serait de programmer le mini pseudo code suivant :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
### Lecture des données
A <- read.file("matrix5.csv") 

### Moyenne de chaque colonne, retourne un vecteur
meanCol(A) 

### Ecart type de chaque colonne, retourne un vecteur
sdCol(A) 

### Régression linéaire, la première colonne en fonction de toutes les autres
lm(A[,1]~A[,2]+A[,3]+A[,4]+A[,5]) 

### Régression linéaire, manuelle (détail dans le post suivant)
inv( t(B)*B ) * ( t(B)*a )

Voila pour l'idée générale. En pratique, ce qui nous intéresse est le temps d'exécution. De plus, on va faire REROLL=100 essais a chaque fois. Donc le code est :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
### Lecture des données
A <- read.file("matrix5.csv") 

### REROLL essais pour chaque opération.
tableauTemps <- matrice vide(4 colonnes et REROLL lignes)
boucle i variant de 1 à REROLL
   tableauTemps[i,1] <- temps de l'instruction( meanCol(A) )
   tableauTemps[i,2] <- temps de l'instruction( sdCol(A) )
   tableauTemps[i,3] <- temps de l'instruction( lm(A[,1]~A[,2]+A[,3]+A[,4]+A[,5]) )
   tableauTemps[i,4] <- temps de l'instruction( inv( t(B)*B ) * ( t(B)*a ) )
finBoucle

Et comme il est possible, pour les petites matrices, que cela soit trop rapide, on ajoute une boucle qui fait faire ITER fois chaque opération :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
### Lecture des données
A <- read.file("matrix5.csv") 

### 100 essais pour chaque opération.
tableauTemps <- matrice vide (4 colonnes et REROLL lignes)
boucle i variant de 1 à REROLL
   tableauTemps[i,1] <- temps de l'instruction(for(j variant de 1 à ITER) faire meanCol(A))
   tableauTemps[i,2] <- temps de l'instruction(for(j variant de 1 à ITER) faire sdCol(A))
   tableauTemps[i,3] <- temps de l'instruction(for(j variant de 1 à ITER) faire lm(A[,1]~A[,2]+A[,3]+A[,4]+A[,5]) )
   tableauTemps[i,4] <- temps de l'instruction(for(j variant de 1 à ITER) faire inv( t(B)*B ) * ( t(B)*a ) )
finBoucle

tableauTemps <- tableauTemps / ITER

Le but est de faire un code optimisé. Idéalement, il faudrait un exemple pour :

• R (Christophe)
• SAS(Laurent Orsy)
• SPSS
• STATA (Virginie)
• R, avec pqr (Mamoun, Christophe)
• R, avec compiler (Mamoun, Christophe)
• R, avec parallel (Mamoun, Christophe)
• C
• C++
• Calm
• Java (Malek)
• Phyton

A priori, il vous faut juste écrire le code, le tester sur 3 ou 4 matrices pour être sur qu'il marche, puis le poster ici. Les matrices sont là, soit en format wide (t colonnes, t² lignes), soit long (t² colonnes, t lignes) :

• Wide : matrix2.csv ; Long : matrixL2.csv (taille 1 Ko)
  
• Wide : matrix3.csv ; Long : matrixL3.csv (taille 1 Ko)
  
• Wide : matrix4.csv ; Long : matrixL4.csv (taille 1 Ko)
  
• Wide : matrix5.csv ; Long : matrixL5.csv (taille 1 Ko)
  
• Wide : matrix7.csv ; Long : matrixL7.csv (taille 2 Ko)
  
• Wide : matrix10.csv ; Long : matrixL10.csv (taille 6 Ko)
  
• Wide : matrix15.csv ; Long : matrixL15.csv (taille 19 Ko)
  
• Wide : matrix22.csv ; Long : matrixL22.csv (taille 57 Ko)
  
• Wide : matrix32.csv ; Long : matrixL32.csv (taille 174 Ko)
  
• Wide : matrix47.csv ; Long : matrixL47.csv (taille 548 Ko)
  
• Wide : matrix69.csv ; Long : matrixL69.csv (taille 1 732 Ko)
  
• Wide : matrix100.csv ; Long : matrixL100.csv (taille 5 266 Ko)
  
• Wide : matrix147.csv ; Long : matrixL147.csv (taille 16 717 Ko)
  
• Wide : matrix216.csv ; Long : matrixL216.csv (taille 53 015 Ko)
  
• Wide : matrix317.csv ; Long : matrixL317.csv (taille 167 527 Ko)
  
• Wide : matrix465.csv ; Long : matrixL465.csv (taille 528 673 Ko)
  
• Wide : matrix682.csv (taille 1 667 699 Ko)
  
• Wide : matrix1000.csv (taille 5 256 817 Ko)

Dans tous les cas, CA N'EST PAS LA PEINE de faire tous les tests puisque je referais tout moi-même sur ma machine (histoire de faire les comparaisons a machine constante).

Des volontaires ? Qui fait quoi ?.

Christophe

[Christophe Genolini]

Quelque précision sur la régression linéaire :

• Soit A une matrice
• Soit a la première colonne de A
• Soit B la matrice A dans laquelle la première colonne est remplacée par des 1
• soit inv, la fonction inverse de matrice
• soit t, la fonction transposée de matrice
• soit *, la multiplication matricielle

Alors, la régression linéaire de la première colonne de A par toutes les autres colonnes de A est :

lm(A) = inv( t(B)*B ) * ( t(B)*a )

Exemple (avec le résultat des calculs pour vérification) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
         (1 2 5) 
Soit A = (4 3 6) 
         (7 8 9) 

alors 

  (1) 
a=(4) 
  (7) 

    (1 2 5) 
B = (1 3 6) 
    (1 8 9) 

       (1 1 1) 
t(B) = (2 3 8) 
       (5 6 9) 

         ( 3   13   20) 
t(B)*B = (13   77  100) 
         (20  100  142) 

              (233.5  38.5 -60.0) 
inv(t(B)*B) = ( 38.5   6.5 -10.0) 
              (-60.0 -10.0  15.5) 

         (12) 
t(B)*a = (78) 
         (90) 

                       (-23) 
inv(t(B)*B)*(t(B)*a) = ( -3) 
                       (  6)

Donc le résultat de la fonction lm(A) doit être (-23,-3,6)

[Christophe Genolini]

Pour R, voila le résultat :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
test <- function(A,reroll,iter){
    tabTemps <- matrix(NA,reroll,4)
    colnames(tabTemps) <- c("colMeans","sdMean","lm")

    a = A[,1]
    B = A
    B[,1] = 1

    for(j in 1:reroll){
       tabTemps[j,1] <- system.time(for(i in 1:iter)apply(A,1,mean))["elapsed"]
       tabTemps[j,2] <- system.time(for(i in 1:iter)apply(A,1,sd))["elapsed"]
       tabTemps[j,3] <- system.time(for(i in 1:iter)lm(A[,1]~A[,-1]))["elapsed"]
       tabTemps[j,4] <- system.time(for(i in 1:iter)solve(t(B)%*%B)%*%(t(B)%*%a))["elapsed"]
    }
    return(tabTemps/iter)
}

### Pour matrix0
A0 <- as.matrix(read.csv2("matrix0.csv"))
temps0 <- test(A0,nbReroll=1000,nbIter=1000)

### Pour matrix2
A2 <- as.matrix(read.csv2("matrix2.csv"))
temps2 <- test(A2,nbReroll=1000,nbIter=1000)

### Pour matrix3
A3 <- as.matrix(read.csv2("matrix3.csv"))
temps3 <- test(A3,nbReroll=1000,nbIter=1000)

...

[Fintuiron]

Salut

Ton code me semble problématique. Pour moi, un for en R est très coûteux en termes de temps. Donc lorsque tu écris :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

for(i in 1:iter)apply(A,1,mean)

J'imagine qu'une grande partie du temps sera "dépensé" par le for et non apply. Donc le calcul sera faussé, ce n'est pas une mesure moyenne de apply mais de for+apply
Idem pour les autres

Il n'y a pas de fonctions de benchmark en R ? Si ce n'est pas le cas, il faut en créer une directement en bas niveau dans un module C

[Christophe Genolini]

Bonjour cher visiteur !

C'est une bonne remarque. Ceci étant dit, dans R, certaines boucles sont pénalisantes mais pas toutes (j'avais appris ça sur le GuR, le "Groupe des Utilisateurs de R"). Ici, si je remplace le for par un sapply ou un lapply, j'obtiens la même chose :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
> system.time(for(i in 1:10000)apply(A,1,mean))["elapsed"]
elapsed 
   0.76 
> system.time(sapply(1:10000,function(x)apply(A,1,mean)))["elapsed"]
elapsed 
   0.77 
> system.time(lapply(1:10000,function(x)apply(A,1,mean)))["elapsed"]
elapsed 
   0.78 
> f <- function(x)apply(A,1,mean)
> system.time(sapply(1:10000,f))["elapsed"]
elapsed 
   0.78

Qu'entends tu par "fonction de benchmark" ?

Christophe

[Fintuiron]

Pour ça que je disais "semble", sans benchmark sur le benchmark (), difficile de dire quelle est l'importance du coût du for.
Mais comme on voit souvent des benchmarks faussés parce que les évaluations posent problèmes, je suis très méfiant.

Pour fonctions de benchmark dans R, je pensais à des choses comme http://cran.r-project.org/web/packag...ark/index.html ou http://cran.r-project.org/web/packag...ark/index.html