Paralelização – {RcppParallel}

Como fazer seu código Rcpp ser ainda mais rápido

Jose Storopoli https://scholar.google.com/citations?user=xGU7H1QAAAAJ&hl=en (UNINOVE)https://www.uninove.br
February 2, 2021

Vimos que o {Rcpp} faz o seu código R ficar muito mais rápido. Tudo o que mostramos até agora foi usando apenas um único core/processador (single thread) do computador1. Agora imaginem o quão rápido seu código R pode ficar se você conseguir rodar {Rcpp} em paralelo 🤯!

Código `{Rcpp}` rodando em paralelo: muito rápido!

Figure 1: Código {Rcpp} rodando em paralelo: muito rápido!

C++ em Paralelo – Biblioteca Intel TBB

O Pacote {RcppParallel} usa a biblioteca TBB da Intel. TBB (Threading Building Blocks) é uma biblioteca de C++ desenvolvida pela Intel para programação paralela em processadores multi-core. Usando TBB, um cálculo é dividido em tarefas que podem ser executadas em paralelo. A biblioteca gerencia e agenda threads para executar essas tarefas.

Como usar {Rcpp} em paralelo – {RcppParallel}

Primeiro, certifique-se que você possui a biblioteca TBB da Intel instalada:

Segundo, instale o pacote {RcppParallel} para R.

Terceiro, coloque em todo código que deseja paralelizar com {RcppParallel} a seguinte síntaxe:

#include <Rcpp.h>
#include <RcppParallel.h>
using namespace Rcpp;
using namespace RcppParallel;

//[[Rcpp::depends(RcppParallel)]]

Pronto! É isso.

Como usar {RcppParallel}

É possível implementar paralelização em diversas partes do seu código {Rcpp} com o {RcppParallel}. Aqui eu vou cobrir apenas os dois algoritmos paralelos do {RcppParallel}2:

Ambos os algortimos usam o struct Worker definido no código do {RcppParallel} que é uma interface para a biblioteca TBB.

Além disso {RcppParallel} usa duas classes, uma para vetores e outra para matrizes:

Quantos Threads?

Ao carregar o {RcppParallel} é importante você designar o número de threads/cores que deseja que o código {RcppParallel} use para paralelização. Caso queira usar todos os seus threads/cores disponíveis, coloque como argumento parallel::detectCores() que retorna um número inteiro com todos os threads/cores disponíveis no seu computador. Aqui estou usando todos os threads/cores disponíveis: 12 threads/cores.

Se não especificado, por padrão {RcppParallel} usa todos os threads/cores disponíveis.

library(Rcpp)
library(RcppParallel)
setThreadOptions(parallel::detectCores())
print(parallel::detectCores())

[1] 12

parallelFor – Paralelizando Loops for

Para usar o parallelFor você deve criar um objeto Worker e definir um operador operator() desse objeto que será invocado pelo {RcppParallel} e roda em paralelo. Isso cria uma função com com o nome do objeto Worker que você criou4. Essa função toma como argumento um intervalo [começo, fim) e lida com todas as questões de segurança e travas de threads que são um porre bem complicadas de maneira automática. Note que o elemento fim do intervalo não é incluído no intervalo (mesmo padrão de comportamento dos iteradores end da biblioteca padrão C++11 STL).

Para mais detalhes, consulte a documentação do parallelFor no site do {RcppParallel}.

Exemplo parallelFor – Raiz Quadrada de Elementos da Matriz.

Aqui vou usar o exemplo da documentação do parallelFor no site do {RcppParallel} de uma função paralela que cacula a raiz quadrada dos elementos de uma matriz. Adicionei alguns comentários para você entender o que está sendo feito. Além disso, há uma versão single-thread também que vamos testar desempenho.

#include <Rcpp.h>
#include <RcppParallel.h>
#include <algorithm>

using namespace Rcpp;
using namespace RcppParallel;

// [[Rcpp::depends(RcppParallel)]]

// Criando um objeto Worker chamado SquareRoot
struct SquareRoot : public Worker
{
   // Variáveis Membro públicas
   const RMatrix<double> input;
   RMatrix<double> output;

   // Construtor do Objeto Worker SquareRoot
   SquareRoot(const Rcpp::NumericMatrix input, Rcpp::NumericMatrix output)
      : input(input), output(output) {}

   // Overload do operador () -- functor
   void operator()(std::size_t begin, std::size_t end) {
      std::transform(input.begin() + begin,
                     input.begin() + end,
                     output.begin() + begin,
                     ::sqrt);
   }
};

// Função que chama o Objeto Worker SquareRoot
// [[Rcpp::export]]
NumericMatrix parallelMatrixSqrt(NumericMatrix x) {

  // Variável local output inicializada
  NumericMatrix output(x.nrow(), x.ncol());

  // Invocação do operador() do Objeto Worker SquareRoot
  SquareRoot squareRoot(x, output);

  // Paralelização do loop for
  parallelFor(0, x.length(), squareRoot);

  return output;
}

// Versão single-thread
// [[Rcpp::export]]
NumericMatrix matrixSqrt(NumericMatrix orig) {
  NumericMatrix mat(orig.nrow(), orig.ncol());
  std::transform(orig.begin(), orig.end(), mat.begin(), ::sqrt);
  return mat;
}
set.seed(123)
b1 <- bench::press(
  n = 10^c(2:3),
  {
    X = matrix(rnorm(n * n), nrow = n)
    bench::mark(
      Rcpp = matrixSqrt(X),
      RcppParallel = parallelMatrixSqrt(X),
      check = FALSE,
      relative = TRUE
    )
})
b1

# A tibble: 4 x 7
  expression       n   min median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>   <dbl> <dbl>  <dbl>     <dbl>     <dbl>    <dbl>
1 Rcpp           100  3.14   3.33      1            1     1   
2 RcppParallel   100  1      1         3.01         1     2.88
3 Rcpp          1000  2.36   2.90      1            1     1   
4 RcppParallel  1000  1      1         2.78         1     2.83
Benchmarks de Raiz Quadrada de Elementos de Matriz: `Rcpp` vs `RcppParallel`

Figure 2: Benchmarks de Raiz Quadrada de Elementos de Matriz: Rcpp vs RcppParallel

Um ganho de 3x com paralelização.

parallelReduce – Paralelizando operações Reduce

Reduce é um algoritmo bem conhecido em ciências da computação. Reduce aplica um operação binária (como adição) em uma sequência definida de elementos, resultando em um único valor. O exemplo sum_of_squares do tutorial 2. Como incorporar C++ no R - {Rcpp} é uma aplicação de um Reduce5. Toda vez que você tiver essa situação você pode paralelizar com parallelReduce.

A lógica do parallelReduce é similar ao parallelFor. Primeiro ambos usam objetos Worker, com algumas diferenças:

Para mais detalhes, consulte a documentação do parallelReduce no site do {RcppParallel}

Exemplo parallelReduce – Soma dos Quadrados

Vamos reutilizar o exemplo sum_of_squares do tutorial 2. Como incorporar C++ no R - {Rcpp}.

Soma dos quadrados é algo que ocorre bastante em computação científica, especialmente quando estamos falando de regressão, mínimos quadrados, ANOVA etc. Vamos paralelizar a implementação ingênua que fizemos no tutorial 2. Como incorporar C++ no R - {Rcpp} com dois loops for. Lembrando que esta implementação será uma função que aceita como parâmetro um vetor de números reais (C++ double / R numeric) e computa a soma de todos os elementos do vetor elevados ao quadrado.

Aqui vamos inserir um std::accumulate() do header <numeric>.

Novamente vou incluir comentários para o entendimento do que estamos fazendo no {RcppParallel}. Além disso, há uma versão single-thread também que vamos testar desempenho.

#include <Rcpp.h>
#include <RcppParallel.h>
#include <algorithm>
using namespace RcppParallel;
using namespace Rcpp;

// [[Rcpp::depends(RcppParallel)]]

// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::plugins("cpp2a")]]

// Criando um objeto Worker chamado sum_of_squares
struct sum_of_squares : public Worker
{
  // Variáveis Membro públicas
  const RVector<double> input;
  double value;

  // Construtor padrão do Objeto Worker
  sum_of_squares(const NumericVector input) : input(input), value(0) {}

  // Construtor "divisor"
  sum_of_squares(const sum_of_squares& sum, Split) : input(sum.input), value(0) {}

  // Overload do operador ()
  void operator()(std::size_t begin, std::size_t end) {
      value += std::accumulate(input.begin() + begin,
                               input.begin() + end,
                               0.0,
                               [] (auto i, auto j) {return i + (j * j);});
   }

  void join(const sum_of_squares& rhs) {
      value += rhs.value;
   }
};

// Função que chama o Objeto Worker sum_of_squares
// [[Rcpp::export]]
double parallel_sum_of_squares(NumericVector x) {
   // variável local inicializada
   sum_of_squares sum(x);

   // Paralelização do Reduce
   parallelReduce(0, x.length(), sum);

   return sum.value;
}

// Versão single-thread
// [[Rcpp::export]]
double sum_of_squares(NumericVector x) {
   return std::accumulate(x.begin(),
                          x.end(),
                          0.0,
                          [] (auto i, auto j) {return i + (j * j);});
}
b2 <- bench::press(
  n = 10^c(4:6),
  {
    v = rnorm(n)
    bench::mark(
      Rcpp = sum_of_squares(v),
      RcppParallel = parallel_sum_of_squares(v),
      check = FALSE,
      relative = TRUE
    )
  })
b2

# A tibble: 6 x 7
  expression         n   min median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>     <dbl> <dbl>  <dbl>     <dbl>     <dbl>    <dbl>
1 Rcpp           10000  1.27   1         1.13         1      NaN
2 RcppParallel   10000  1      1.15      1            1      Inf
3 Rcpp          100000  4.82   3.84      1            1      NaN
4 RcppParallel  100000  1      1         3.81         1      Inf
5 Rcpp         1000000  8.96   8.13      1            1      NaN
6 RcppParallel 1000000  1      1         8.34         1      NaN
Benchmarks de Soma dos Quadrados: `Rcpp` vs `RcppParallel`

Figure 3: Benchmarks de Soma dos Quadrados: Rcpp vs RcppParallel

Mais um sucesso! Ganho de 8x 🤯 com paralelização.

Usar {RcppParallel} no seu Pacote R

As instruções abaixo foram retiradas da documentação do {RcppParallel}.

Se você deseja usar {RcppParallel} de dentro de um pacote R, você precisa editar vários arquivos para criar os links de construção e tempo de execução necessários. As seguintes adições devem ser feitas:

Ambiente

sessionInfo()

R version 4.0.4 (2021-02-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.10

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
[1] RcppParallel_5.0.3 Rcpp_1.0.6        

loaded via a namespace (and not attached):
 [1] tidyselect_1.1.0  xfun_0.22         bslib_0.2.4      
 [4] purrr_0.3.4       colorspace_2.0-0  vctrs_0.3.6      
 [7] generics_0.1.0    htmltools_0.5.1.1 emo_0.0.0.9000   
[10] yaml_2.2.1        utf8_1.1.4        rlang_0.4.10     
[13] jquerylib_0.1.3   pillar_1.5.1      glue_1.4.2       
[16] DBI_1.1.1         lifecycle_1.0.0   stringr_1.4.0    
[19] munsell_0.5.0     gtable_0.3.0      ragg_1.1.1       
[22] bench_1.1.1       evaluate_0.14     knitr_1.31       
[25] parallel_4.0.4    fansi_0.4.2       profmem_0.6.0    
[28] highr_0.8         scales_1.1.1      jsonlite_1.7.2   
[31] debugme_1.1.0     farver_2.1.0      systemfonts_1.0.1
[34] textshaping_0.3.2 distill_1.2       ggplot2_3.3.3    
[37] digest_0.6.27     stringi_1.5.3     dplyr_1.0.5      
[40] grid_4.0.4        cli_2.3.1         tools_4.0.4      
[43] magrittr_2.0.1    sass_0.3.1        tibble_3.1.0     
[46] tidyr_1.1.3       crayon_1.4.1      pkgconfig_2.0.3  
[49] downlit_0.2.1     ellipsis_0.3.1    lubridate_1.7.10 
[52] assertthat_0.2.1  rmarkdown_2.7     rstudioapi_0.13  
[55] R6_2.5.0          compiler_4.0.4

  1. tecnicamente, Eigen e Armadillo podem, dependendo da configuração do sistema operacional, automaticamente se beneficiar de paralelizações usando o OpenMP.↩︎

  2. a biblioteca TBB tem muito mais algoritmos complexos caso necessario. Recomendo você olhar este link da documentação do {RcppParallel}.↩︎

  3. tecnicamente é um MapReduce.↩︎

  4. mais questões técnicas: quando você define um operador operator() de um objeto em C++ você dá um overload no operador parenthesis do objeto e o resultado é uma síntaxe similar à uma função com o nome do objeto.↩︎

  5. tecnicamente é um MapReduce.↩︎

Corrections

If you see mistakes or want to suggest changes, please create an issue on the source repository.

Reuse

Text and figures are licensed under Creative Commons Attribution CC BY-SA 4.0. Source code is available at https://github.com/storopoli/Rcpp, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

Citation

For attribution, please cite this work as

Storopoli (2021, Feb. 2). Rcpp - A interface entre R e C++: Paralelização -- `{RcppParallel}`. Retrieved from https://storopoli.github.io/Rcpp/4-RcppParallel.html

BibTeX citation

@misc{storopoli2021rcppparallel,
  author = {Storopoli, Jose},
  title = {Rcpp - A interface entre R e C++: Paralelização -- `{RcppParallel}`},
  url = {https://storopoli.github.io/Rcpp/4-RcppParallel.html},
  year = {2021}
}

© Copyright 2021 Jose Storopoli.

Content licensed under the CC BY-SA 4.0.