Benchmark de funções

Introdução

Estamos interessados em analisar o desempenho de funções que façam multiplicação matricial, para isso iremos comparar diversas funções implementadas em R e em C. A seguir definiremos as funções. Estamos interessados na que seja mais rápida, exigindo menos tempo computacional.

Funções em R

Utilizaremos as funções implementadas em R como vemos a seguir

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

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

m_2 <- 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)
}


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

Funções em C

Utilizaremos as funções em C implementadas a seguir

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

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

// [[Rcpp::export]]
Rcpp::NumericVector c_1(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 c_2(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 c_2inmp(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;
  #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 c_2outmp(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);
}

Conclusões

Para analisar o desempenho de cada função utilizaremos o pacote “microbenchmark”, que estima o tempo que cada método levou para ser executado. Multiplicaremos uma matriz de tamanho nxn por um vetor de tamanho n. Primeiramente trabalharemos com o n = 1000 com 10 repetições.

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


base <- microbenchmark(M0 = m_0(M, x),
                      M1 = m_1(M, x),
                      M2 = m_2(M, x),
                      M3 = m_3(M, x), 
                      C0 = c_0(M, x),
                      C1 = c_1(M, x),
                      C2 = c_2(M, x),
                      C2in = c_2inmp(M, x),
                      C2out = c_2outmp(M, x),times = 10)

Podemos observar os resultados a partir do gráfico de boxplot a seguir.

base %>% 
  ggplot() + 
  geom_boxplot(aes(x = reorder(expr, time, median), y = time, fill = expr, alpha = 0.5))

Como podemos ver, o método M2 está muito distante dos outros, dificultando a análise. Iremos portanto eliminá-lo do gráfico para facilitar a visualização

base %>% filter(expr != "M2") %>%
  ggplot() +
  geom_boxplot(aes(x = reorder(expr, time, median), y = time, fill = expr, alpha = 0.5))

Notamos que as funções implementadas em C tendem a ser mais rápidas do que as implementadas em R, com exceção da função M0 (função do R que utiliza %*% para multiplicação matricial). Investigaremos se esse comportamento permanece ao aumentarmos o n para n = 10000.

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


base2 <- microbenchmark(M0 = m_0(M, x),
                      M1 = m_1(M, x),
                      M2 = m_2(M, x),
                      M3 = m_3(M, x), 
                      C0 = c_0(M, x),
                      C1 = c_1(M, x),
                      C2 = c_2(M, x),
                      C2in = c_2inmp(M, x),
                      C2out = c_2outmp(M, x),times = 10)

base2 %>% filter(expr != "M2") %>%
  ggplot() +
  geom_boxplot(aes(x = reorder(expr, time, median), y = time, fill = expr, alpha = 0.5))

Ao aumentar o n para 10.000 o método M0 se consolida como sendo realmente o mais rápido, se destacando mais ainda em relação aos outros visto que a diferença de tempo de execução para os outros métodos aumentou ainda mais. Portanto podemos concluir que para os n testados, as funções definidas aqui tem esse comportamento. Sendo as funções para multiplicação matricial em C geralmente mais rápidas que as funções em R, com exceção do método M0 (que utiliza %*%) que se destaca como o método mais rápido.