Databases
StackOverflow2013
Postsdbo
POdmvnorm MVN density - RcppArmadillo implementation slower than R package including a bit of Fortran
 

Note that there are some explanatory texts on larger screens.

plurals
  1. POdmvnorm MVN density - RcppArmadillo implementation slower than R package including a bit of Fortran
    primarykey
    data
    text
    <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>
    singulars
    1. PO
      singulars
      1. PTAnswer
    1. This table or related slice is empty.
    1. PTQuestion
    1. USInferrator
    1. USInferrator
    plurals
    1. PL
      singulars
      1. LTLinked
    1. PL
      singulars
      1. LTLinked
    1. This table or related slice is empty.
    1. PO
      singulars
      1. PTAnswer
    1. VO
      singulars
      1. POdmvnorm MVN density - RcppArmadillo implementation slower than R package including a bit of Fortran
      1. This table or related slice is empty.
      1. VTUpMod
    2. VO
      singulars
      1. POdmvnorm MVN density - RcppArmadillo implementation slower than R package including a bit of Fortran
      1. This table or related slice is empty.
      1. VTUpMod
    1. COI don't understand what your question is
      singulars
      1. POdmvnorm MVN density - RcppArmadillo implementation slower than R package including a bit of Fortran
      1. USDavid Marx
    2. COlook at the code of `mvtnorm::dmvnorm`. It uses `eigen`, which is calling LAPACK code for symmetric eigenvalues ... it may be hard to get faster for this very standard piece of computational linear algebra ...
      singulars
      1. POdmvnorm MVN density - RcppArmadillo implementation slower than R package including a bit of Fortran
      1. USBen Bolker
    3. COwell, but Armadillo is supposed to be fast with linear algebra as well ... I think it uses Lapack. Is there so much overhead when calling an Rcpp function? That would surprise me. Besides, the mahalanobis distance computation should be faster in Rcpp.
      singulars
      1. POdmvnorm MVN density - RcppArmadillo implementation slower than R package including a bit of Fortran
      1. USInferrator
 

Querying!

 
Guidance
A row detail

Detail views are divided into sections. All the information in the data section comes from columns in the selected row. The other sections display data from other, related rows.

Related data can be related in a to-one or a to-many fashion. Captions of data related in a to-many fashion link to a list view showing a filtered view of the table.

Try moving around until you find a non-empty to-many entry and click on the label to get to one. You can move back to the root by clicking on the database name in the header.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload