Nous présentons des méthodes simples pour accélérer des gros calculs sous R en utilisant de la parallélisation. Afin d’évaluer différentes librairies de calcul, nous utilisons le script R-benchmark-25.
Un BLAS contient les routines de bas niveau permettant de réaliser les opérations basiques sur les vecteurs et les matrices (addition, produits scalaire et matriciel, etc.). Le BLAS de R natif est single-thread et des alternatives plus performantes existent telles que OpenBLAS qui est multi-thread, ATLAS qui s’adapte automatiquement à l’architecture locale de manière optimisée ou Intel MKL pour les processeurs Intel (aussi multi-thread).
Voici un tutoriel pour installer de nouvelles librairies BLAS/LAPACK optimisées sur les principaux système d’exploitations: [https://csantill.github.io/RPerformanceWBLAS/]
Benchmark sur Intel® Core™ i5-8250U CPU @ 1.60GHz × 8 avec 8GO de RAM sous Ubuntu.
BLAS/LAPLACK
Creation, transp., deformation of a 2500x2500 matrix (sec): 0.735666666666667
2400x2400 normal distributed random matrix ^1000____ (sec): 0.54
Sorting of 7,000,000 random values__________________ (sec): 0.772333333333333
2800x2800 cross-product matrix (b = a' * a)_________ (sec): 16.5476666666667
Linear regr. over a 3000x3000 matrix (c = a \ b')___ (sec): 8.70833333333333 ATLAS
Creation, transp., deformation of a 2500x2500 matrix (sec): 0.687
2400x2400 normal distributed random matrix ^1000____ (sec): 0.530333333333333
Sorting of 7,000,000 random values__________________ (sec): 0.762
2800x2800 cross-product matrix (b = a' * a)_________ (sec): 1.93433333333333
Linear regr. over a 3000x3000 matrix (c = a \ b')___ (sec): 0.966666666666666 OpenBLAS (avec parallélisation)
Creation, transp., deformation of a 2500x2500 matrix (sec): 0.709333333333333
2400x2400 normal distributed random matrix ^1000____ (sec): 0.544666666666667
Sorting of 7,000,000 random values__________________ (sec): 0.780666666666667
2800x2800 cross-product matrix (b = a' * a)_________ (sec): 0.835333333333333
Linear regr. over a 3000x3000 matrix (c = a \ b')___ (sec): 0.274000000000001 Mais suivant la manière dont les différents package sont codés, les gains ne sont pas forcément évident. Voici un benchmark avec différents package d’analyse de réseaux: sbm, missSBM, GREMLINS:
BLAS/LAPLACK
Data frenchblog2007
==================== missSBM 1 core: ==================
Time difference of 28.6378 secs
==================== sbm 1 core: ==================
Time difference of 22.2913 secs
==================== missSBM 4 cores: ==================
Time difference of 21.01817 secs
==================== sbm 4 cores: ==================
Time difference of 17.59729 secs
Data MPEcoNetwork
==================== Multipartite 1 core: ==================
Time difference of 1.683621 mins
==================== Multipartite 4 cores: ==================
Time difference of 1.008362 minsAtlas
Data frenchblog2007
==================== missSBM 1 core: ==================
Time difference of 28.85101 secs
==================== sbm 1 core: ==================
Time difference of 16.92881 secs
==================== missSBM 4 cores: ==================
Time difference of 21.15375 secs
==================== sbm 4 cores: ==================
Time difference of 13.93179 secs
Data MPEcoNetwork
==================== Multipartite 1 core: ==================
Time difference of 1.722041 mins
==================== Multipartite 4 cores: ==================
Time difference of 56.11809 secsopenBLAS
Data frenchblog2007
==================== missSBM 1 core: ==================
Time difference of 1.282785 mins
==================== sbm 1 core: ==================
Time difference of 17.01022 secs
==================== missSBM 4 cores: ==================
Time difference of 1.094572 mins
==================== sbm 4 cores: ==================
Time difference of 15.65691 secs
Data MPEcoNetwork
==================== Multipartite 1 core: ==================
Time difference of 1.66076 mins
==================== Multipartite 4 cores: ==================
Time difference of 55.89039 secsMicrosoft R Open est une distribution améliorée de R conçue par Microsoft. Elle permet aussi d’utiliser facilement le BLAS d’Intel MKL avec R, ce qui accélère considérablement les calculs, notamment sous Windows (mais aussi sous Mac et Linux).
Voici les résultats du benchmark sur Intel(R) Core(TM) i3-7100U CPU @ 2.40GHz x 2 sous Windows 10 :
Sous R 4.1.1 et le BLAS natif :
Creation, transp., deformation of a 2500x2500 matrix (sec): 0.913333333333333
2400x2400 normal distributed random matrix ^1000____ (sec): 1.22333333333333
Sorting of 7,000,000 random values__________________ (sec): 1.05
2800x2800 cross-product matrix (b = a' * a)_________ (sec): 20.0033333333333
Linear regr. over a 3000x3000 matrix (c = a \ b')___ (sec): 9.58333333333333 Sous Microsoft R Open 4.0.2 et Intel MKL :
Creation, transp., deformation of a 2500x2500 matrix (sec): 0.923333333333333
2400x2400 normal distributed random matrix ^1000____ (sec): 1.16666666666667
Sorting of 7,000,000 random values__________________ (sec): 1.04
2800x2800 cross-product matrix (b = a' * a)_________ (sec): 0.513333333333332
Linear regr. over a 3000x3000 matrix (c = a \ b')___ (sec): 0.213333333333334 Un autre avantage de Microsoft R Open est la reproductibilité. En effet, chaque version de MROpen est associée à un snapshot du CRAN. Par exemple, la dernière version de MROpen est la 4.0.2 est fixée au 16/07/2020. Par conséquent, les packages qu’on essaie installer sous MROpen 4.0.2 sont récupérés dans leur état au 16/07/2020.
Sous Linux et MacOS, il est possible de paralléliser un code en faisant du multiprocess (en plus du multi-thread). Cela n’est pas possible sous Windows. Quelques librairies R pour le calcul parallèle : parallel, foreach, doParallel. Ici, on discute d’une curiosité rencontrée avec la fonction mclapply de la librairie parallel.
Les systèmes d’exploitation Linux et macOS ont un système de “shared memory” via POSIX qui évite de recopier certains objets. Toutefois il est difficile de complètement comprendre le comportement.
Just the tip what might have been going on R-devel Digest, Vol 149, Issue 22
Radford Neal’s answer from Jul 26, 2015:
“When mclapply forks to start a new process, the memory is initially shared with the parent process. However, a memory page has to be copied whenever either process writes to it. Unfortunately, R’s garbage collector writes to each object to mark and unmark it whenever a full garbage collection is done, so it’s quite possible that every R object will be duplicated in each process, even though many of them are not actually changed (from the point of view of the R programs).”
Voici un example:
library(lobstr)
library(parallel)
library(bettermc)##
## Attaching package: 'bettermc'## The following object is masked from 'package:parallel':
##
## mclapplyIci, on voit bien que la première case de la liste B pointe vers la case mémoire de la matrice A.
A <- rnorm(1e6)
B <- list(A)
print(lobstr::obj_addr(A))## [1] "0x5574ba1127d0"print(lobstr::obj_addrs(B))## [1] "0x5574ba1127d0"print(lobstr::obj_sizes(A, B))## * 8,000,048 B
## * 56 BOn a le même comportement avec un lapply:
list_A <- lapply(
seq(10),
function(x) {
B[[1]]
}
)
print(lobstr::obj_addrs(list_A))## [1] "0x5574ba1127d0" "0x5574ba1127d0" "0x5574ba1127d0" "0x5574ba1127d0"
## [5] "0x5574ba1127d0" "0x5574ba1127d0" "0x5574ba1127d0" "0x5574ba1127d0"
## [9] "0x5574ba1127d0" "0x5574ba1127d0"print(lobstr::obj_sizes(A, list_A, B))## * 8,000,048 B
## * 176 B
## * 56 BPar contre dès que l’on parallélise, on perd se comportement:
list_A <- parallel::mclapply(
seq(10),
function(x) {
return(B[[1]])
}, mc.cores = 2L
)
print(lobstr::obj_addrs(list_A))## [1] "0x5574b690cc90" "0x5574bc534780" "0x5574b70aded0" "0x5574bccd59c0"
## [5] "0x5574bae510c0" "0x5574bd476c00" "0x5574bb5f2300" "0x5574bdc17e40"
## [9] "0x5574bbd93540" "0x5574be3b9080"print(lobstr::obj_sizes(A, list_A, B))## * 8,000,048 B
## * 80,000,656 B
## * 56 BL’objet semble être recopier lors du return:
list_A <- bettermc::mclapply(
seq(10),
function (x) {
C <- B[[1]] + x
lobstr::obj_addr(C)
}, mc.cores = 2L
)
unique(unlist(list_A))## [1] "0x5574beb5a2c0" "0x5574bf2fb500" "0x5574bfa9c740" "0x5574c023d980"
## [5] "0x5574c09debc0"Et le comportement semble différent suivant que l’on preschedule ou non :
list_A <- bettermc::mclapply(
seq(10),
function (x) {
C <- B[[1]] + x
lobstr::obj_addr(C)
}, mc.cores = 2L, mc.preschedule = FALSE
)
unique(unlist(list_A))## [1] "0x5574beb5a2c0"Le package bettermc [https://github.com/gfkse/bettermc] donne plus de controle sur le comportement, en donnant des options pour controler les retours des processus enfants en utilisant la “shared memory” POSIX.
f <- function(i) A
microbenchmark::microbenchmark(
bettermc1 = bettermc::mclapply(1:2, f, mc.share.copy = FALSE),
bettermc2 = bettermc::mclapply(1:2, f),
bettermc3 = bettermc::mclapply(1:2, f, mc.share.vectors = FALSE),
bettermc4 = bettermc::mclapply(1:2, f, mc.share.vectors = FALSE, mc.shm.ipc = FALSE),
parallel = parallel::mclapply(1:2, f),
times = 10, setup = gc()
)## Unit: milliseconds
## expr min lq mean median uq max neval
## bettermc1 29.0073 30.0416 31.83080 31.57230 32.9582 35.2849 10
## bettermc2 48.5614 48.6860 51.37020 49.68985 53.5721 59.6321 10
## bettermc3 72.4206 72.8788 74.87865 73.97910 75.3169 81.5929 10
## bettermc4 59.5879 66.2118 69.01306 68.88600 70.4146 76.3639 10
## parallel 60.7242 61.1275 64.46019 62.80345 67.5300 74.1411 10Le package permet également de faire des parallélisation reproductible et donne beaucoup plus de controle sur la gestion des erreurs en permettant entre autre de relancer les processus qui ont renvoyé des erreurs.
L’extention bettermcExt [https://github.com/gfkse/bettermcExt] non autorisé sur le CRAN permet d’overloader la fonction mclapply de parallel d’un package donné, permettant ainsi un gain de temps et de rendre les calculs reproductibles.