Integración de R y C++ – Acelerando los programas

En los últimos meses he estado tratando de aprender a optimizar mi código en R. En general, R funciona bastante bien, pero hay varias aplicaciones donde se puede poner muy, pero muy lento, en particular los loops. Un búsqueda rápida muestra que existen muchos artículos que discuten el tema con bastante profundidad.

Uno de los temas que encontré y llamaron mi atención fue la capacidad de integrar R con otros lenguajes de programación, en particular con algunos de más bajo nivel, como C, Fortran o C++. La integración con C y Fortran está ya solucionada directamente en R, a través de la función .Call, pero esta forma de llamar a código externo tiene la gran desventaja que todas las funciones deben ser declaradas como void, es decir, debemos modificar los valores entregados usando directamente punteros a los objetos en memoria. Esto es… engorroso, al menos lo fue en mis primeras experiencias programando.

La solución: usar C++. Es un lenguaje moderno, que corre eficientemente los loops, y que tiene muchas herramientas que facilitan la vida al programador. Para integrar R y C++ existe el excelente paquete Rcpp, el cual permite crear funciones de forma dinámica y sencilla en R, y que además permite usar objetos inherentes a R (vectores numéricos, vectores de caractéres, dataframes, etc.) en C++, pero aún este paquete tiene una falencia: no soporta operaciones matriciales eficientes. Para ello, un plugin del mismo paquete llamado RcppArmadillo ofrece una interfaz entre Rcpp y la librería Armadillo que entrega herramientas para operaciones matriciales en C++. Con esto tenemos todo lo que necesitamos.

En este artículo programaré un código sencillo para obtener un clúster de K-Medias y compararé las implementaciones en código en R, la implementación oficial (en C) y una implementación propia en C++. El algoritmo es bastante sencillo: inicializa vectores de forma aleatoria, para luego ir ajustando los clusters en cada iteración hasta alcanzar convergencia.

Implementación (ineficiente) en R.

kmeans.R <- function(x,k,eps = 10e-6){
  # k es la cantidad de clusters, x es la matriz de datos.
  # eps es la tolerancia máxima entre una iteración y otra.
  require(pdist)
  cluster <- rep(0,nrow(x)) # Cluster actual del caso.
  # Inicializa centroides.
  centroids <- matrix(runif(nrow(x) * k), nrow = k, ncol = ncol(x))
  # Guarda última iteración
  centroids.old <- matrix(0, nrow = k, ncol = ncol(x))
  eps.vec <- rep(eps,k)
  dists.centroids <- 1000
  while(any(dists.centroids > eps.vec)){
    centroids.old <- centroids
    # Matriz con las distancias entre casos y vectores.
    dists.x <- as.matrix(pdist(x, centroids)) 
    min.x <- apply(dists.x, 1, min)
    # Obtiene indice del menor valor, i.e., el cluster.
    cluster <- apply(cbind(dists.x,min.x), 1, 
                     function(x) which(x[1:k] == x[k+1])[1])     
    for(i in 1:k){
      # Calcula el centroide.
      centroids[i,] <- apply(x[cluster == i, ], 2, mean) 
    }
    dists.centroids <- diag(as.matrix(pdist(centroids, centroids.old)))
  }
  # Devuelve cluster y centroides.
  return(list(cluster = cluster, centroids = centroids)) 
}

Para probar el código, en R simularemos una matriz sencilla con dos distribuciones. La primera una normal de estándar y la segunda una normal de media 5 y desviación estándar 1. Posteriormente aplicamos la función y graficamos. Noten que estoy usando el paquete ‘pdist’, solo para facilitarme la vida cuando calculo las distancias entre matrices. Este paquete tira algunas warnings en la última iteración, simplemente para avisar que convergió, así que ignórenlas.

library(pdist)
x <- rbind(matrix(rnorm(1000 * 2, mean = 0), nrow = 1000, ncol = 2), 
           matrix(rnorm(1000 * 2, mean = 5), nrow = 1000, ncol = 2))
x.kmeans.R <- kmeans.R(x, 2)
x.kmeans.R$centroids
# Resultado:
#            [,1]       [,2]
# [1,] 0.04097412 0.04618219
# [2,] 5.02960813 4.97157822
plot(x, col=ifelse(x.kmeans.R$cluster == 1, "red", "green"))

El resultado es el siguiente:

Datos en dos segmentos.

Implementación en C++.
Implementemos ahora la misma función en C++. Para ello usaremos Rcpp y RcppArmadillo. Primero, es necesario instalar los paquetes:

install.packages(c('Rcpp', 'RcppArmadillo'))

.

Ahora, existen varias rutas para implementar la función. Es posible escribir directamente la función en un gran string de texto y compilarla en línea, o se puede guardar un archivo anexo que tenga la función programada. Seguiremos este camino. En un nuevo archivo (kmeansCpp.cpp) programamos la función de kmedias.

#include <RcppArmadillo.h>
using namespace Rcpp;

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

// [[Rcpp::export]]
List kmeansCpp(NumericMatrix xa, int k, double eps = 10e-3) {
  arma::mat x(xa.begin(), xa.nrow(), xa.ncol(), false);
  // Inicializa centroides de forma uniforme en [0-1]
  arma::mat centroids(k, xa.ncol(), arma::fill::randu);
  // Incializa con puros 0's.
  arma::mat centroids_old(k, xa.ncol(), arma::fill::zeros);
  // Crea vector fila tamaño 1xk.
  arma::rowvec eps_vec = eps * arma::ones<arma::rowvec>(k);
  arma::rowvec dists_centroids = 100 * arma::ones<arma::rowvec>(k);
  // Crea matriz tamaño casos x k.
  arma::mat dists_x(xa.nrow(), k);
  LogicalMatrix temp(xa.nrow(), k);
  arma::mat cluster(xa.nrow(), k);
  while(arma::any(dists_centroids > eps_vec)){
    centroids_old = centroids;
    // Calcula matriz de distancias
    for(int i = 0; i<k; i++){
      dists_x.col(i) = sqrt(sum((x - (arma::ones(xa.nrow(), 1) 
                                  * centroids.row(i))) % 
                      (x - (arma::ones(xa.nrow(), 1) * centroids.row(i))), 1));
    }
    // Obtiene indice del menor valor, i.e., el cluster.
    arma::colvec min_x = min(dists_x, 1);
    temp = wrap((min_x * arma::ones(1, k) == dists_x));
    cluster = as<arma::mat>(temp);
    // Calcula el centroide.
    for(int ia = 0; ia < k; ia++){
      centroids.row(ia) = mean(x.rows(find(cluster.col(ia) == 1)), 0);
    }
    dists_centroids = (sqrt(sum((centroids - centroids_old) % 
                      (centroids - centroids_old), 1))).t();
  }
  // Devuelve cluster y centroides.
  return List::create(Named("cluster") = wrap(cluster), 
                      Named("centroids") = wrap(centroids));
}

Analicemos un poco el programa:

  • La línea #include <RcppArmadillo.h> le indica al compilador de C++ que importe la cabecera que nos permitirá usar todas las herramientas para la interfaz con R.
  • using namespaces Rcpp nos permite utilizar todas las funciones de Rcpp sin tener que colocar Rcpp::NombreFuncion cada vez que la llamamos. También puede ser útil colocar using namespaces arma para ahorrarse los arma::.
  • // [[Rcpp::depends(RcppArmadillo)]] y // [[Rcpp::export]] le dicen a Rcpp que tiene que ocupar Armadillo y que debe exportar la función a R.

Luego viene la función en si. Tanto Rcpp como Armadillo poseen clases que permiten manejar vectores, matrices, listas, etc. En este caso, es importante notar que creamos una función que recibe y entrega objetos de Rcpp. La función recibe una matriz numérica (NumericMatrix, objeto de Rcpp), un entero y un double y entrega una lista de R. Los objetos de Rcpp son, a grandes rasgos, matrices (Matrix), vectores (Vector), listas (List) y Dataframes (Dataframe); mientras que los dos primeros pueden ser numéricos (Numeric), enteros (Integer), literales (String) o booleanos (Logical). Así, una matriz numérica corresponde a un objeto NumericMatrix, mientras que un vector booleano será un LogicalVector. Cada objeto tiene diversos métodos que se le pueden aplicar, pero en general las operaciones matemáticas funcionan correctamente. Para una descripción detallada, este tutorial es un buen punto de partida.

Dentro de la función, las primeras líneas inicializan los objetos que usaremos en la iteración. Por ejemplo, arma::mat x(xa.begin(), xa.nrow(), xa.ncol(), false); indica que crearemos una matriz x, objeto de Armadillo (los objetos más usados son vec, para vector numérico y mat para matrices numéricas) reutilizando la memoria que apunta a la matriz numérica de Rcpp xa, nuestro input. Los demás objetos se inicializan a partir de los tipos que usaremos.

El siguiente paso es el loop principal, que sigue la misma estructura que el loop de R, pero ahora utilizando las funciones propias de Armadillo o Rcpp, dependiendo de cuáles sean los objetos sobre los que se las estamos aplicando. Por ejemplo, el operador % corresponde a la multiplicación término a término entre matrices de Armadillo. Finalmente, devolvemos una lista de R con el cluster (ahora como una matriz con índices) y los centroides. La función wrap() permite transformar objetos de Armadillo en objetos de Rcpp fácilmente.

Para poder utilizar la función en R en nuestra sesión corremos el comando Rcpp::sourceCpp("kmeansCpp.cpp"). Con ello la función kmeansCpp() está disponible para su uso.

Comparando resultados

Comparemos los códigos, para ver cuánto ganamos por traspasar la función a C++. Utilizaré el paquete rbenchmark que permite rápidamente comparar distintas funciones en cuanto a su tiempo de ejecución. El código de comparación completo queda entonces:

library(RcppArmadillo)
library(rbenchmark)
Rcpp::sourceCpp("kmeansCpp.cpp")
x <- rbind(matrix(rnorm(1000 * 2, mean = 0), nrow = 1000, ncol = 2), 
           matrix(rnorm(1000 * 2, mean = 5), nrow = 1000, ncol = 2))
benchmark(kmeans(x,2), kmeansCpp(x,2), kmeans.R(x,2), order = "elapsed",
          columns = c("test", "replications", "elapsed", "relative", 
                      "user.self"))
# Resultados:
#              test replications elapsed relative user.self
# 2 kmeansCpp(x, 2)          100   0.040    1.000     0.036
# 1    kmeans(x, 2)          100   0.059    1.475     0.056
# 3  kmeans.R(x, 2)          100   7.775  194.375     7.752

Nuestra implementación en C++ corre en 0.036 segundos, 194 veces más rápido (!!!) que la implementación en R puro y corre además casi un 50% más rápido que la implementación en C que trae incorporado R. El resultado habla por si solo, vale la pena traspasar códigos con muchas iteraciones a C++, pues es más eficiente tanto en uso de memoria como en velocidad. La diferencia con kmeans de R es engañosa, pues no realizamos ningún control de errores ni depuraciones a las entradas que probablemente sí realiza esta función, por lo que las ganancias se deben a eso probablemente.

Así que ya saben, si quieren que su código vuele, bajen un poco de nivel y utilicen código en C++. R ofrece todas las herramientas para lograrlo de forma sencilla.

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.