Simulación

Cómputo paralelo en R

Cualquier operación donde no todos los pasos dependen de los resultados de los pasos anteriores es un candidato para ser paralelizado.

Un cluster es un grupo de dos o más núcleos que llevan a cabo tareas de forma conjunta en paralelo.

Aplicación de funciones a vectores de forma paralela

> sapply(1:10, function(x) 2^x)
[1] 2 4 8 16 32 64 128 256 512 1024

El paquete parallel [2] facilita la ejecución paralela en R. No requiere instalación manual ya que forma parte de la instalación base de R. Se carga para su uso con la instrucción library(parallel) que no produce ninguna salida en el caso de que haya sido exitosa su carga al sistema.

Se puede determinar el número de núcleos disponibles en un sistema con la llamada

> detectCores()
[1] 8

y se puede proceder a paralelizar cuando haya por lo menos dos. Si se quiere excluir los núcleos virtuales (en el caso que estén presentes), se llama con detectCores(logical=FALSE):

> detectCores(logical=FALSE)
[1] 4

No conviene utilizar todos los núcleos en las tareas de R, ya que el sistema operativo y su interfaz de usuario suelen requerir mantenerse responsivos durante los cálculos realizados, por lo cual es recomendable crear un cluster con un número menor de núcleos, por ejemplo cluster <- makeCluster(detectCores() - 1). Teniendo el cluster, se puede realizar llamadas que lo aprevechen, como por ejemplo

> parSapply(cluster, 1:10, function(x) 2^x)
[1] 2 4 8 16 32 64 128 256 512 1024

Una vez que ya no se requiera el cluster, se le libera con la instrucción stopCluster(cluster).

Cuando se ocupa utilizar una variable o una librería dentro de lo que se evalúa en el cluster, es necesario exportarlo de manera explícita al cluster:

> cluster <- makeCluster(detectCores() - 1)
> desde <- 3
> hasta <- 8
> parSapply(cluster, desde:hasta, function(x) 2^x)
[1] 8 16 32 64 128 256
> base <- 3
> parSapply(cluster, desde:hasta, function(x) base^x)
Error in checkForRemoteErrors(val) :
6 nodes produced errors; first error: object 'base' not found
> clusterExport(cluster, "base")
> parSapply(cluster, desde:hasta, function(x) base^x)
[1] 27 81 243 729 2187 6561
> stopCluster(cluster)

Procesamiento de sequencias en paralelo

Otra forma de realizar cálculos paralelizados es con el paquete foreach que no viene instalado por omisión sino requiere instalación manual:

> library(foreach)
Error in library(foreach) : there is no package called ‘foreach’
> install.packages("foreach")
also installing the dependency ‘iterators’

trying URL 'https://cran.itam.mx/bin/macosx/mavericks/contrib/3.3/iterators_1.0.8.tgz'
Content type 'application/x-gzip' length 310135 bytes (302 KB)
==================================================
downloaded 302 KB

trying URL 'https://cran.itam.mx/bin/macosx/mavericks/contrib/3.3/foreach_1.4.3.tgz'
Content type 'application/x-gzip' length 382270 bytes (373 KB)
==================================================
downloaded 373 KB


The downloaded binary packages are in
/var/folders/qj/1gw209fx1qb31lz61njz69k00000gn/T//RtmphnuSlq/downloaded_packages

Es importante tomar en cuenta que los números de los servidores de descarga no siempre son los mismos cuando uno instala poquetes en computoradoras diferentes o en momentos diferentes: siempre busco el más cercano a mi y tecleo la selección. Al utilizar R en un ambiente de interfaz gráfica, típicamente se abre una ventana para realizar la selección del servidor.

También se requiere un segundo paquete:

> install.packages("doParallel")
trying URL 'https://cran.itam.mx/bin/macosx/mavericks/contrib/3.3/doParallel_1.0.10.tgz'
Content type 'application/x-gzip' length 181113 bytes (176 KB)
==================================================
downloaded 176 KB

The downloaded binary packages are in
/var/folders/qj/1gw209fx1qb31lz61njz69k00000gn/T//RtmphnuSlq/downloaded_packages

Con el paquete foreach, la creación y liberación de un cluster es ligeramente diferente:

> library(doParallel)
Loading required package: foreach
foreach: simple, scalable parallel programming from Revolution Analytics
Use Revolution R for scalability, fault tolerance and more.
http://www.revolutionanalytics.com
Loading required package: iterators
> registerDoParallel(makeCluster(detectCores() - 1))
> foreach(potencia = desde:hasta, .combine=c) %dopar% base^potencia
[1] 27 81 243 729 2187 6561
> stopImplicitCluster()

En las llamadas a foreach, se puede proprocionar opciones .export y .packages para proporcionarle variables externas y librerías adicionales.

Cómputo con GPU

Tarjetas gráficas contemporaneas contienen procesdaores gráficos con bastante núcleos. El paquete gpuR permite aprovechar esos núcleos para cálculos. Los tipos de dato para realizar cálculos con un GPU son especiales y provienen de ese paquete. Finalizando la instalación del paquete, igual como al cargarlo al uso, se detectan GPUs presentes en la computadora en cuestión, por ejemplo, en una computadora de escritorio iMac detecta

Number of platforms: 1
- platform: Apple: OpenCL 1.2 (Apr  4 2017 19:07:42)
- gpu index: 0
    - AMD Radeon R9 M395X Compute Engine

o en una laptop Macbook detecta

Number of platforms: 1
- platform: Apple: OpenCL 1.2 (Apr  4 2017 19:07:42)
  - gpu index: 0
     - Intel(R) HD Graphics 5300

Por ejemplo, el producto de matrices que usa de prueba Charles Determan (el autor del paquete) en su blog post hace lo siguiente con el GPU de Radeon:

> suppressMessages(library(gpuR))
> n <- 2000
> M1 <- matrix(runif(n * n), nrow=n, ncol=n)
> M2 <- matrix(rnorm(n * n), nrow=n, ncol=n)
> system.time(M1 %*% M2)
   user  system elapsed
  4.587   0.002   4.591
> G1 <- gpuMatrix(M1, type="double")
> G2 <- gpuMatrix(M2, type="double")
> system.time(G1 %*% G2)
   user  system elapsed
  0.082   0.027   0.224

y esto con el GPU de Intel (hay que cambiar el double por float ya que es un GPU más modesto:

> suppressMessages(library(gpuR))
> n <- 2000
> M1 <- matrix(runif(n * n), nrow=n, ncol=n)
> M2 <- matrix(rnorm(n * n), nrow=n, ncol=n)
> system.time(M1 %*% M2)
   user  system elapsed
  7.016   0.033   7.092
> G1 <- gpuMatrix(M1, type="float")
> G2 <- gpuMatrix(M2, type="float")
> system.time(G1 %*% G2)
   user  system elapsed
  0.079   0.042   1.970 

En otras computadoras con GPUs de NVIDIA se obtiene lo siguiente:

> library(gpuR)
Number of platforms: 1
- platform: NVIDIA Corporation: OpenCL 1.2 CUDA 8.0.0
  - gpu index: 0
    - GeForce GT 720
checked all devices
completed initialization
gpuR 1.2.1
> n <- 2000
> M1 <- matrix(runif(n * n), nrow=n, ncol=n)
> M2 <- matrix(rnorm(n * n), nrow=n, ncol=n)
> system.time(M1 %*% M2)
   user system elapsed
   4.09   0.02    4.11
> G1 <- gpuMatrix(M1, type="double")
> G2 <- gpuMatrix(M2, type="double")
> system.time(G1 %*% G2)
   user system elapsed
   0.41   1.31    1.72

> library(gpuR)
Number of platforms: 1
- platform: NVIDIA Corporation: OpenCL 1.2 CUDA 8.0.0
  - gpu index: 0
    - Tesla K80
checked all devices
completed initialization
gpuR 1.2.1
> n <- 2000
> M1 <- matrix(runif(n * n), nrow=n, ncol=n)
> M2 <- matrix(rnorm(n * n), nrow=n, ncol=n)
> system.time(M1 %*% M2)
  user  system elapsed
5.424   0.000   5.424
> G1 <- gpuMatrix(M1, type="double")
> G2 <- gpuMatrix(M2, type="double")
> system.time(G1 %*% G2)
  user  system elapsed
 1.348   0.052   1.672

> library(gpuR)
Number of platforms: 1
- platform: NVIDIA Corporation: OpenCL 1.2 CUDA 11.1.102
  - context device index: 0
    - GeForce GTX 1650
checked all devices
completed initialization
gpuR 2.0.3
> n <- 2000
> M1 <- matrix(runif(n * n), nrow=n, ncol=n)
> M2 <- matrix(rnorm(n * n), nrow=n, ncol=n)
> system.time(M1 %*% M2)
  user  system elapsed
  4.577   0.000   4.582
> G1 <- gpuMatrix(M1, type="double")
> G2 <- gpuMatrix(M2, type="double")
> system.time(G1 %*% G2)
  user  system elapsed
  0.920   0.029   1.060

> library(gpuR)
library(gpuR)
Number of platforms: 1
- platform: NVIDIA Corporation: OpenCL 1.2 CUDA 11.1.102
  - context device index: 0
    - GeForce RTX 3090
checked all devices
completed initialization
gpuR 2.0.3
> n <- 2000
> M1 <- matrix(runif(n * n), nrow=n, ncol=n)
> M2 <- matrix(rnorm(n * n), nrow=n, ncol=n)
> system.time(M1 %*% M2)
  user  system elapsed
  4.770   0.000   4.771
> G1 <- gpuMatrix(M1, type="double")
> G2 <- gpuMatrix(M2, type="double")
> system.time(G1 %*% G2)
  user  system elapsed
  0.121   0.041   0.191

Participantes que cuentan con acceso a un GPU pueden consultar la documentación del paquete para mayor detalle en cuáles operaciones con matríces y vectores están disponibles para realizarse en el GPU. En todas las tareas, se aprecia el uso adecuado de cálculos GPU y son una alternativa válida a paralelismo en el CPU — también conviene incluir comparasiones del desempeño cuando ambas opciones están disponibles.

Actualizado el 16 de noviembre del 2021.
https://satuelisa.github.io/simulation/paralelo/R.html