繰り返し距離:RcppArmadilloを使ってみる

  • こちらのマハラノビス距離の繰り返し計算でスピードアップの様子を見てみよう
  • cxxfunction()関数を使っている例をsourceCpp()用に書き換えてみる
    • データ型指定はRcpp namespace形式らしい
## Using RcppArmadillo
code <- 'arma::mat x = Rcpp::as<arma::mat>(X);
arma::mat mu = Rcpp::as<arma::mat>(Mu);
arma::mat sigma = Rcpp::as<arma::mat>(Sigma);
int n = x.n_rows;
arma::vec md(n);
for (int i=0; i<n; i++){
arma::mat x_i = x.row(i) - mu;
arma::mat Y = arma::solve( sigma, arma::trans(x_i) );
md(i) = arma::as_scalar(x_i * Y);
}
return wrap(md);'
mahalanobis_RcppArma <- cxxfunction( signature(X="numeric",
Mu="numeric", Sigma="numeric"), code, plugin="RcppArmadillo")
#include <RcppArmadilloExtensions/sample.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp ;

// [[Rcpp::export]]

NumericVector mahalanobis_RcppArma(NumericMatrix X,NumericMatrix Mu, NumericMatrix Sigma){
	arma::mat x = Rcpp::as<arma::mat>(X);
	arma::mat mu = Rcpp::as<arma::mat>(Mu);
	arma::mat sigma = Rcpp::as<arma::mat>(Sigma);
	int n = x.n_rows;
	arma::vec md(n);
	for (int i=0; i<n; i++){
		arma::mat x_i = x.row(i) - mu;
		arma::mat Y = arma::solve( sigma, arma::trans(x_i) );
		md(i) = arma::as_scalar(x_i * Y);
	}
	return wrap(md);
}
## Using R
mahalanobis_R <- function(X, Mu, Sigma){
md <- rep(NA, nrow(X) )
for(i in 1:nrow(X) )
md[i] <- mahalanobis(X[i,], center=Mu, cov=Sigma)
return(md)
}
X <- matrix(runif(10000),100,100)
Mu <- matrix(runif(100),1,100)
Sigma <- matrix(runif(10000),100,100)

benchmark(
mahalanobis_R(X, c(Mu), Sigma),
mahalanobis_RcppArma(X, Mu, Sigma),
columns=c("test", "replications","elapsed", "relative"),
replications=100,
order="relative")
                                test replications elapsed relative
2 mahalanobis_RcppArma(X, Mu, Sigma)          100   10.11    1.000
1     mahalanobis_R(X, c(Mu), Sigma)          100   49.02    4.849