En este post mostrar茅 algunos trucos que he ido aprendiendo al utilizar Random Forests (RF), una de mis t茅cnicas de clasificaci贸n favoritas, para mejorar el poder de clasificaci贸n que poseen. Este m茅todo, desarrollado por Leo Breiman en el a帽o 2001 tiene sus or铆genes en esta publicaci贸n en Machine Learning y se ha convertido, especialmente en los 煤ltimos cinco o seis a帽os, en una de las preferidas por cient铆ficos y practitioners. La t茅cnica es tremendamente poderosa, es muy f谩cil de usar (en comparaci贸n a las redes neuronales por ejemplo), y entrega bastante informaci贸n sobre los resultados al usuario.

La l贸gica de los RF es generar muchos 谩rboles de decisi贸n sobre los datos de entrada, utilizando una cantidad reducida de las variables que se encuentran disponibles. La t茅cnica explora as铆 diferentes subsegmentos del espacio de entrada, y aprende patrones muy complejos aprovechando la diversidad de los clasificadores que entrena. Una prueba de su 茅xito se ha visto en las competencias de Kaggle: es una de las dos t茅cnicas m谩s importantes que identifican los administradores del sitio.

En este ejemplo entrenar茅 un RF para un dataset del repositorio de la UCI, y luego mostrar茅 algunos trucos: c贸mo acelerar el entrenamiento utilizando programaci贸n paralela para aprovechar mejor los recursos disponibles y c贸mo graficar una curva ROC utilizando

ggplot2

. En un post futuro mostrar茅 c贸mo balancear las muestras para solucionar los problemas de desbalance y finalmente c贸mo graficar la importancia de las variables y analizar estos resultados.

Entrenamiento Secuencial

Primero importemos los datos. Para este ejemplo usar茅 la base de datos de Bank Marketing, de Moro et. al. Esta base de datos posee 41.188 instancias, 20 variables predictoras, y tiene como objetivo predecir si el cliente toma o no un dep贸sito a plazo. Usar茅 R para descargar directamente el archivo, descomprimirlo, y cargar los datos. Al final agrego una variable adicional para los test (muestra holdout).

# Descargar el archivo. Se debe usar "method=curl" si la p谩gina usa https.
url <- "https://archive.ics.uci.edu/ml/machine-learning-databases/00222/bank-additional.zip";
download.file(url,"bank-additional.zip", mode="wb",method="curl");
# Descomprimir el archivo y cargar la tabla.
unzip("bank-additional.zip");
bank.data <- read.table("bank-additional-full.csv", sep=";", header = TRUE);
bank.data$if.test <- rbinom(nrow(bank.data), 1, 0.7) # Agrega variable para test

Con esto el archivo est谩 cargado en la variable

bank.data

y ya se puede utilizar. Entrenemos un primer Random Forest y veamos sus resultados. Para entrenar los modelos est谩 el paquete

randomForest

de Breiman mismo, y adem谩s utilizar茅 los paquetes

verification

,

rbenchmark

,

e1071

y

caret

para analizar los resultados. Carguemos los paquetes o instal茅moslos si es necesario.

# library(..., logical.return = TRUE) devuelve falso si el paquete no existe.
if(!library(randomForest, logical.return = TRUE)) {
  install.packages('randomForest');
}
if(!library(verification, logical.return = TRUE)) {
  install.packages('verification');
}
if(!library(caret, logical.return = TRUE)) {
  install.packages('caret');
}
if(!library(e1071, logical.return = TRUE)) {
  install.packages('e1071');
  library(e1071)
}
if(!library(rbenchmark, logical.return = TRUE)) {
  install.packages('rbenchmark');
  library(rbenchmark)
}

Ahora podemos correr un primer Random Forest.

# Crea f贸rmula con todas las variables. Con los paste evito escribir de m谩s.
vars.bank <- colnames(bank.data)
formula.rf <- as.formula(paste0("y ~", paste(vars.bank[1:20], collapse = "+")))

# Funci贸n para facilitar futuras ejecuciones.
rfNaive <- function(x){randomForest(formula.rf, data = x, 
                                     subset = (if.test == 0), ntree = 500,
                                     importance = TRUE)
}
# Entrena sin ninguna precauci贸n. Corre usando 650MB en RAM y en
# 36 segundos en un procesador core i7-3770k.
rf.naive <- rfNaive(bank.data)

# Calcula AUC y Crea matriz  de confusi贸n con corte en 0.5.
prob.rfNaive <- predict(rf.naive, newdata = bank.data[bank.data$if.test == 1 , ],
                        type = "prob")[, 2]
roc.area(as.numeric(bank.data$y[bank.data$if.test == 1] == "yes"), prob.rfNaive)
confusionMatrix(round(prob.rfNaive),  as.numeric(bank.data$y[bank.data$if.test == 1] == "yes"))
# AUC = 0.9426409
# Matriz de confusi贸n:
#               Reference
# Prediction     0     1
#         0   24711  1546
#         1   939  1703

Entrenando en Paralelo

Los Random Forest son altamente paralelizables pues cada 谩rbol puede entrenarse por si solo, y al final del entrenamiento basta con juntar cada uno de los 谩rboles que se entrenaron. El paquete

randomForest

toma ventaja de eso y entrega directamente una manera de combinar resultados cuando se entrenan bosques en paralelo.

Para paralelizar en R es necesario registrar un backend que es el que se encarga de comunicarse con los procesadores, y adem谩s usar un paquete que permita paralelizar el c贸digo. Hay distintos paquetes con distintos niveles de profundidad en cu谩nto a control (y cuanta dificultad) entregan al usuario. Lo que me ha resultado m谩s f谩cil es:

  • Usar como backend ya sea el paquete doMC (Linux) o el paquete doParallel(Windows). Lamentablemente no existe (que yo sepa) un paquete que funcione en ambos sistemas.
  • Usar el paquete foreach para realizar la paralelizaci贸n.

Las opciones de estos paquetes son amplias, y escapan el alcance de este post, as铆 que me enfocar茅 en lo necesario para correr randomForests. Para todos los detalles de paralelizar R se puede comenzar por la vi帽eta de foreach.

驴Qu茅 es paralelizar? La gran mayor铆a de los computadores modernos traen m谩s de un procesador (n煤cleo) incorporado para correr procesos. La ventaja de tener m谩s de un procesador radica en que el usuario usa m煤ltiples programas y servicios en su equipo en cualquier momento, y el tener m谩s de un n煤cleo permite poder realizar varias actividades de forma simult谩nea. Nosotros podemos tomar ventaja de esto para generar c贸digos m谩s eficientes, que utilicen todo el poder computacional que tenemos. Para paralelizar un Random Forest la estrategia a seguir ser谩 simplemente decirle al programa que entrene una cierta cantidad de 谩rboles en cada n煤cleo, y despu茅s combine los resultados en un 煤nico bosque aleatorio.

La computaci贸n paralela tiene una propiedad que es importante tener presente: dependiendo de cu谩n eficiente es la implementaci贸n del c贸digo, puede multiplicar la cantidad de memoria utilizada por la cantidad de procesadores que se utilizan. As铆, si una aplicaci贸n utiliza 500MB en RAM en un n煤cleo, perfectamente puede utilizar 2GB si se utilizan cuatro, por ejemplo. Por lo mismo, mucho cuidado al correr este ejemplo, ajusten los valores de los par谩metros acorde a la memoria RAM que tengan disponibles.

Para paralelizar nuestro Random Forest, lo primero es determinar cu谩ntos n煤cleos poseemos y cu谩ntos 谩rboles queremos entrenar. En este ejemplo estamos entrenando 500 谩rboles, y en mi PC tengo ocho n煤cleos y suficiente RAM, as铆 que los usar茅 todos. El n煤mero de procesadores se puede obtener googleando el tipo de procesador que tienen, o en Windows con el backend parallel pueden correr

detectCores()

en R, o en Linux el backend doMC utiliza la mitad de los n煤cleos disponibles de forma predeterminada.

Con eso listo ahora podemos paralelizar. Le diremos a R que entrene sub-bosques, y despu茅s los combine. El c贸digo es el siguiente:

cores <- 4 # Depende de cada computador.
ntrees <- 500 # Decisi贸n del usuario.
library(doMC)
registerDoMC(cores) # Registra el backend. Solo en Linux.
# En Windows:
# library(doParallel)
# registerDoParallel(cores = cores)

rfParallel <- function(x){
  foreach(ntree.iter = rep(ceiling(ntrees / cores), cores), 
          .combine = combine, .packages = "randomForest") %dopar% {
            randomForest(formula.rf, data = x, ntree = ntree.iter,
                         do.trace = FALSE, subset = (if.test == 0),
                         keep.forest = TRUE, importance = TRUE)
          }
}

rf.parallel <- rfParallel(bank.data)

# Calcular medidas efectividad
prob.rfParallel <- predict(rf.parallel, newdata = bank.data[bank.data$if.test == 1 , ],
                        type = "prob")[, 2]
roc.area(as.numeric(bank.data$y[bank.data$if.test == 1] == "yes"), prob.rfParallel)
confusionMatrix(round(prob.rfParallel),  
                as.numeric(bank.data$y[bank.data$if.test == 1] == "yes"))
# AUC = 0.9439
# Matriz de confusi贸n:
#               Reference
# Prediction     0     1
#       0      24665  1584
#       1        905  1673
# Accuracy: 91.37%
# Balanced Accuracy: 73.91%

Como podemos ver, los resultados son equivalentes entre los dos modelos. Existe una peque帽a ca铆da en la capacidad discriminante, pero esto es s贸lo por razones aleatorias.

Los par谩metros de la funci贸n

foreach

son los importantes:

  •  ntree.iter = rep(ceiling(ntrees / cores), cores)

    : Indica que cada iteraci贸n de foreach entrene un total de 500 divido por la cantidad de n煤cleos que disponen.

  • .combine = combine

    : Indica a

    foreach

    que utilice la funci贸n 芦combine禄 del paquete

    randomForest

    para combinar los resultados. Esta funci贸n permite unir bosques aleatorios en bosques m谩s grandes.

  • .packages='randomForest'

    : Indica a

    foreach

    que, en cada n煤cleo donde correr谩 la funci贸n, cargue el paquete

    randomForest

    .

  • %dopar%

    : Indica a

    foreach

    que corra en paralelo utilizando el backend que se encuentre registrado (

    doMC

    en este caso).

Ahora comparemos el tiempo de corrida. Repetir茅 tres veces el entrenamiento utilizando el excelente paquete

rBenchmark

.

benchmark(rfNaive(bank.data), rfParallel(bank.data), 
          order = "elapsed", replications = 3,
          columns = c("test", "replications", "elapsed", "relative",
                      "user.self"))

#         test           replications elapsed relative user.self
#2 rfParallel(bank.data)       3       77.760   1.00     2.204
#1    rfNaive(bank.data)       3      148.557   1.91   138.380

Este ejemplo lo corr铆 en mi laptop con dos n煤cleos reales y cuatro virtuales. La programaci贸n en paralelo corre casi dos veces m谩s r谩pido que el c贸digo secuencial. Esto es en general lo que ocurre: existe una ganacia que es casi lineal con el n煤mero de n煤cleos f铆sicos que se disponen.

 

Consejo 2: Graficar las curvas ROC con ggplot2

Para finalizar este post mostrar茅 algo simple para quienes reportan los resultados de modelos. Las curvas ROC de R en general se ven bastante mal, incluso en el mismo paquete

verification

indican que existen mejores alternativas para graficar.

El paquete m谩s potente para gr谩ficos en R corresponde a ggplot2, pero utilizarlos es bastante complejo. El paquete utiliza el concepto de est茅ticas, o variables que se asignan a distintas partes del gr谩fico. Toma un poco de tiempo acostumbrarse a esto, pero una vez que se logra los resultados son muy buenos.

Para graficar una curva ROC, es necesario primero que nada de disponer un data.frame tal que muestre los puntos de la curva ROC y a qui茅n se refiere. Luego creamos el gr谩fico, asignamos la variable que identifica al modelo tanto a la est茅tica de grupos como a la de colores, y luego le damos los valores correspondientes. El c贸digo es el siguiente:

# Crea estructura con curva ROC
roc.rfParallel <- roc.plot(as.numeric(bank.data$y[bank.data$if.test == 1] == "yes"), 
                           prob.rfParallel)
roc.rfParallel <- as.data.frame(roc.rfParallel$plot.data)
roc.rfNaive <- roc.plot(as.numeric(bank.data$y[bank.data$if.test == 1] == "yes"), 
                           prob.rfNaive)
roc.rfNaive <- as.data.frame(roc.rfNaive$plot.data)

temp <- data.frame(fpr = roc.rfParallel$V3, var = 'En Paralelo', value = roc.rfParallel$V2)
temp <- rbind(temp, data.frame(fpr = roc.rfNaive$V3, var = 'Secuencial', value = roc.rfNaive$V2))

# Leyenda del gr谩fico
legend.roc <- rep("",2)
legend.roc[1] <- paste0("En Paralelo: AUC = ", round(with(bank.data, roc.area(as.numeric(y[if.test == 1] == "yes"),
                                                                              prob.rfParallel))$A,3))
legend.roc[2] <- paste0("Secuencial: AUC = ", round(with(bank.data, roc.area(as.numeric(y[if.test == 1] == "yes"),
                                                                                 prob.rfNaive))$A,3))

breaks <- c("En Paralelo", "Secuencial")

# Gr谩fico
library(ggplot2)
library(grid)

p <- (ggplot(temp, aes(fpr, value, group = var, color = var)) 
      + geom_line(size = 0.5, aes(linetype = var))+ theme_bw()
      + scale_linetype_manual(name = 'leyenda', values=c("dotted", "solid"), 
                              labels = legend.roc, breaks = breaks)
      + scale_color_discrete(name = 'leyenda', breaks = breaks, labels = legend.roc)
      + scale_x_continuous("False Positive Rate (1-Specificity)")
      + scale_y_continuous("True Positive Rate (Sensitivity)")
      + coord_fixed()
      + theme(legend.justification=c(1,0), legend.position=c(1,0),
              legend.title = element_blank(), 
              legend.text = element_text(size = 14),
              legend.key.width = unit(3, "cm"),
              legend.key = element_blank()
      )
)
plot(p)

El resultado es:

 

Curva ROC Random Forests

Curva ROC Random Forests

Como se observa, ambos modelos son id茅nticos. Esto puede cambiar en una comparaci贸n entre distintos tipos de modelos.

Este post mostr贸 tanto como entrenar los bosques aleatorios en paralelo y como generar curvas ROC con

ggplot2

. Cualquier comentario, por twitter o directamente ac谩 en el blog. 隆Gracias por leer!

Actualizaci贸n: Una versi贸n anterior de este post ten铆a algunos paquetes faltantes en el c贸digo. 隆Gracias Pablo!

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.