Benchmark de Funções em R e C

Funções em R

library(tidyverse)

R1 <- function(M, x){
  return(M %*% x)
}

R2 <- function(M, x){
  Mx <- matrix(NA, nrow = n, ncol = 1)
  for(i in 1:nrow(M))
    Mx[i,1] = M[i,] %*% x
  return(Mx)
}

R3 <- function(M,x){
  Mx <- matrix(NA, nrow = n, ncol = 1)
  for(i in 1:nrow(M)){
    val = 0
    for (j in 1:ncol(M)){
      val = val + M[i,j]*x[j,1]
    }
    Mx[i,1] = val
  }
  return(Mx)
}

R4 <- function(M, x)
{
  Mx = apply(M, 1, function(m){m %*% x})
  return(matrix(Mx, ncol = 1))
}

Como pode ser visto, foram criadas 4 funções do R, são elas: R1, R2, R3 e R4.

Funções em C

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
#include <omp.h>

// [[Rcpp::export]]
arma::vec C1(const arma::mat& M, const arma::vec& x){
  
  return(M*x);
}

// [[Rcpp::export]]
Rcpp::NumericVector C2(Rcpp::NumericMatrix M, Rcpp::NumericVector x)
{
  int r = M.nrow();
  int c = M.ncol();
  Rcpp::NumericVector Mx(r);
  for (int i = 0; i < r; i++)
  {
    double val = 0;
    for (int j = 0; j < c; j++)
    {
      val += M(i,j)*x(j);
    }
    Mx(i) = val;
  }
  return(Mx);
}

// [[Rcpp::export]]
arma::vec C3(const arma::mat M, const arma::vec x)
{
  arma::vec Mx(x.size(), arma::fill::zeros);
  int r = M.n_rows;
  int c = M.n_cols;
  for (int i = 0; i < r; i++)
  {
    double val = 0;
    for (int j = 0; j < c; j++)
    {
      val += M(i,j)*x(j);
    }
    Mx(i) = val;
  }
  return(Mx);
}

// [[Rcpp::export]]
arma::vec C3InMP(const arma::mat& M, const arma::vec& x, int cores = 4)
{
  arma::vec Mx(x.size(), arma::fill::zeros);
  omp_set_num_threads(cores);
  int r = M.n_rows;
  int c = M.n_cols;
#pragma omp parallel for schedule(dynamic)
  for (int i = 0; i < r; i++)
  {
    double val = 0;
    for (int j = 0; j < c; j++)
    {
      val += M(i,j)*x(j);
    }
    Mx(i) = val;
  }
  return(Mx);
}

// [[Rcpp::export]]
arma::vec C3OutMP(const arma::mat& M, const arma::vec& x, int cores = 4)
{
  arma::vec Mx(x.size(), arma::fill::zeros);
  omp_set_num_threads(cores);
  int r = M.n_rows;
  int c = M.n_cols;
#pragma omp parallel for schedule(static)
  for (int i = 0; i < r; i++)
  {
    double val = 0;
    for (int j = 0; j < c; j++)
    {
      val += M(i,j)*x(j);
    }
    Mx(i) = val;
  }
  return(Mx);
}

Por fim, foram criadas 5 funções em C, nomeadas de C1, C2, C3, C3InMP e C3OutMP, em que as duas últimas são adaptações da função C3.

Então, temos ao todo 9 funções a serem comparadas.

Benchmark

Realizando esta comparação por meio da função microbenchmark, obtemos o seguinte resultado:

library(microbenchmark)

n <- 1000
M <- matrix(runif(n*n), nrow = n)
x <- matrix(runif(n), nrow = n)

base <- microbenchmark(R1 = R1(M,x),
                       R2 = R2(M,x),
                       R3 = R3(M,x),
                       R4 = R4(M,x),
                       C1 = C1(M,x),
                       C2 = C2(M,x),
                       C3 = C3(M,x),
                       C3InMP = C3InMP(M,x),
                       C3OutMP = C3OutMP(M,x),
                       times = 10)

base %>% ggplot() + geom_boxplot(aes(x = reorder(expr, time, median), y = time, fill = expr)) + labs(fill = "Função", x = "Função", y = "Tempo de Execução") + ggtitle("Figura 1: Boxplots do Microbenchmark")

É bastante visível que as funções em C são, de modo geral, bem mais rápidas que as em R. Entretanto, a função R1 possui uma performance excepcional, sendo uma das melhores, pois utiliza uma multiplicação matricial direta.

Nos boxplots, notamos que a função R3 destoa bastante do restante, fazendo com que não possamos ver tantos detalhes nas comparações das demais funções. Retirando, então, esta função, e aumentando o tamanho das matrizes e das repetições feitas, obtemos os boxplots abaixo:

library(microbenchmark)

n <- 2000
M <- matrix(runif(n*n), nrow = n)
x <- matrix(runif(n), nrow = n)

base <- microbenchmark(R1 = R1(M,x),
                       R2 = R2(M,x),
                       R3 = R3(M,x),
                       R4 = R4(M,x),
                       C1 = C1(M,x),
                       C2 = C2(M,x),
                       C3 = C3(M,x),
                       C3InMP = C3InMP(M,x),
                       C3OutMP = C3OutMP(M,x),
                       times = 20)

base %>% filter(!(expr %in% c("R3"))) %>% ggplot() + geom_boxplot(aes(x = reorder(expr, time, median), y = time, fill = expr)) + labs(fill = "Função", x = "Função", y = "Tempo de Execução") + ggtitle("Figura 2: Boxplots do Microbenchmark com matrizes maiores")

Notamos que, claramente, as funções C1 e R1 tem uma perfomance muito melhor que as demais. Isto pode ser explicado pois utilizam diretamente a multiplicação matricial, sem loops em seus códigos.