title: dmvnorm MVN density - RcppArmadillo implementation slower than R package including a bit of Fortran
body: <p>The <strong>solution</strong> is now online in the <strong><a href="http://gallery.rcpp.org/articles/dmvnorm_arma/" rel="nofollow">Rcpp Gallery</a></strong></p>
  
  <hr>
  
  <p>I re-implemented dmvnorm from the mvtnorm package in RcppArmadillo. I somehow like Armadillo, but I guess it would also work in plain Rcpp. The approach from dmvnorm is based on the mahalanobis distance, so I have a function for that and then the multivariate normal density function.</p>
  
  <p>Let me show you my code:</p>
  
  <pre><code>#include &lt;RcppArmadillo.h&gt;
  #include &lt;Rcpp.h&gt;
  
  // [[Rcpp::depends("RcppArmadillo")]]
  
  // [[Rcpp::export]]
  arma::vec mahalanobis_arma( arma::mat x ,  arma::mat mu, arma::mat sigma ){
  
    int n = x.n_rows;
    arma::vec md(n);
      for (int i=0; i&lt;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 md;
  
      }
  
  
  
  // [[Rcpp::export]]
  arma::vec dmvnorm ( arma::mat x,  arma::mat mean,  arma::mat sigma, bool log){ 
  
  arma::vec distval = mahalanobis_arma(x,  mean, sigma);
  
      double logdet = sum(arma::log(arma::eig_sym(sigma)));
      double log2pi = 1.8378770664093454835606594728112352797227949472755668;
      arma::vec logretval = -( (x.n_cols * log2pi + logdet + distval)/2  ) ;
  
         if(log){ 
           return(logretval);
  
         }else { 
         return(exp(logretval));
           }
  }
  </code></pre>
  
  <p>So, and not to my big disappointment:</p>
  
  <p>simulate some data</p>
  
  <pre><code>sigma &lt;- matrix(c(4,2,2,3), ncol=2)
  x &lt;- rmvnorm(n=5000000, mean=c(1,2), sigma=sigma, method="chol")
  </code></pre>
  
  <p>and benchmark </p>
  
  <pre><code>system.time(mvtnorm::dmvnorm(x,t(1:2),.2+diag(2),F))
     user  system elapsed 
     0.05    0.02    0.06 
  
  system.time(dmvnorm(x,t(1:2),.2+diag(2),F))
     user  system elapsed 
     0.12    0.02    0.14 
  </code></pre>
  
  <p>No!!!!!! :-(</p>
  
  <p>[EDIT]</p>
  
  <p>The <strong>questions</strong> are:
  1) Why is the RcppArmadillo implementation slower than a plain R implementation?
  2) How do I create an Rcpp/RcppArmadillo implementation that beats the R implementation?</p>
  
  <p>[EDIT 2]</p>
  
  <p>I put in the mahalanobis_arma into the mvtnorm::dmvnorm function and it also slows down.</p>