- こちらのマハラノビス距離の繰り返し計算でスピードアップの様子を見てみよう
- cxxfunction()関数を使っている例をsourceCpp()用に書き換えてみる
- データ型指定はRcpp namespace形式らしい
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>
using namespace Rcpp ;
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);
}
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