Doesn't svd in R by default compute D, U and V?<br><br><a href="http://stat.ethz.ch/R-manual/R-patched/library/base/html/svd.html">http://stat.ethz.ch/R-manual/R-patched/library/base/html/svd.html</a><br><br><div class="gmail_quote">
On Wed, Jun 13, 2012 at 4:07 PM, Douglas Bates <span dir="ltr"><<a href="mailto:bates@stat.wisc.edu" target="_blank">bates@stat.wisc.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="HOEnZb"><div class="h5">On Wed, Jun 13, 2012 at 5:16 PM, Dirk Eddelbuettel <<a href="mailto:edd@debian.org">edd@debian.org</a>> wrote:<br>
><br>
> On 13 June 2012 at 15:05, Julian Smith wrote:<br>
> | I agree that RcppEigen is a little bit faster, but ease of use is important to<br>
> | me, so I feel like RcppArmadillo might win out in my application.<br>
><br>
> Yup, that my personal view too.<br>
><br>
> | | RcppArmadillo will use the very same LAPACK and BLAS libs your R session<br>
> | | uses. So MKL, OpenBlas, ... are all options.  Eigen actually has its own<br>
> | code<br>
> | | outperforming LAPACK, so it doesn't  as much there.<br>
> |<br>
> | Why do you think R outperforms RcppArmadillo in this example below? Anyway to<br>
> | speed this up?<br>
><br>
> That is odd. "I guess it shouldn't." I shall take another look -- as I<br>
> understand it both should go to the same underlying Lapack routine.  I may<br>
> have to consult with Conrad on this.<br>
><br>
> Thanks for posting a full and reproducible example!<br>
><br>
> Dirk<br>
><br>
> | require(RcppArmadillo)<br>
> | require(inline)<br>
> |<br>
> | arma.code <- '<br>
> |   using namespace arma;<br>
> |   NumericMatrix Xr(Xs);<br>
> |   int n = Xr.nrow(), k = Xr.ncol();<br>
> |   mat X(Xr.begin(), n, k, false);<br>
> |   mat U;<br>
> |   vec s;<br>
> |   mat V;<br>
> |   svd(U, s, V, X);<br>
> |   return wrap(s);<br>
> | '<br>
<br>
</div></div>Because the arma code is evaluating the singular vectors (U and V) as<br>
well as the singular values (S) whereas the R code is only evaluating<br>
the singular values.  There is considerably more effort required to<br>
evaluate the singular vectors in addition to the singular values.<br>
<div class="HOEnZb"><div class="h5"><br>
> | rcppsvd <- cxxfunction(signature(Xs="numeric"),<br>
> |                         arma.code,<br>
> |                         plugin="RcppArmadillo")<br>
> |<br>
> | A<-matrix(rnorm(5000^2), 5000)<br>
> |<br>
> | > system.time(rcppsvd(A))<br>
> |     user   system  elapsed<br>
> | 1992.406    4.862 1988.737<br>
> |<br>
> | > system.time(svd(A))<br>
> |    user  system elapsed<br>
> | 652.496   2.641 652.614<br>
> |<br>
> | On Wed, Jun 13, 2012 at 11:43 AM, Dirk Eddelbuettel <<a href="mailto:edd@debian.org">edd@debian.org</a>> wrote:<br>
> |<br>
> |<br>
> |     On 13 June 2012 at 10:57, Julian Smith wrote:<br>
> |     | I've been toying with both RcppArmadillo and RcppEigen the past few days<br>
> |     and<br>
> |     | don't know which library to continue using. RcppEigen seems really slick,<br>
> |     but<br>
> |     | appears to be lacking some of the decompositions I want and isn't nearly<br>
> |     as<br>
> |     | fast to code. RcppArmadillo seems about as fast, easier to code up etc.<br>
> |     What<br>
> |     | are some of the advantages/disadvantages of both?<br>
> |<br>
> |     That's pretty close.  I have been a fan of [Rcpp]Armadillo which I find<br>
> |     easier to get my head around.  Doug, however, moved from [Rcpp]Armadillo<br>
> |     to<br>
> |     [Rcpp]Eigen as it has some things he needs.  Eigen should have a "larger"<br>
> |     API<br>
> |     than Armadillo, but I find the code and docs harder to navigate.<br>
> |<br>
> |     And you should find Eigen to be a little faster. Andreas Alfons went as far<br>
> |     as building 'robustHD' using RcppArmadillo with a drop-in for RcppEigen (in<br>
> |     package 'sparseLTSEigen'; both package names from memmory and I may have<br>
> |     mistyped).  He reported a performance gain of around 25% for his problem<br>
> |     sets.  On the 'fastLm' benchmark, we find the fast Eigen-based<br>
> |     decompositions<br>
> |     to be much faster than Armadillo.<br>
> |<br>
> |     | Can you call LAPACK or BLAS from either? Is there a wrapper in RcppEigen<br>
> |     to<br>
> |     | call LAPACK functions? Want some other decomposition methods, dont like<br>
> |     the<br>
> |     | JacobiSVD method in Eigen.<br>
> |<br>
> |     You need to differentiate between the Eigen and Armadillo docs _for their<br>
> |     libraries_ and what happens when you access the Rcpp* variant from R.<br>
> |<br>
> |     RcppArmadillo will use the very same LAPACK and BLAS libs your R session<br>
> |     uses. So MKL, OpenBlas, ... are all options.  Eigen actually has its own<br>
> |     code<br>
> |     outperforming LAPACK, so it doesn't  as much there.<br>
> |<br>
> |     Hope this helps,   Dirk (at useR!)<br>
> |<br>
> |     |<br>
> |     | ----------------------------------------------------------------------<br>
> |     | _______________________________________________<br>
> |     | Rcpp-devel mailing list<br>
> |     | <a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
> |     | <a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
> |     --<br>
> |     Dirk Eddelbuettel | <a href="mailto:edd@debian.org">edd@debian.org</a> | <a href="http://dirk.eddelbuettel.com" target="_blank">http://dirk.eddelbuettel.com</a><br>
> |<br>
> |<br>
><br>
> --<br>
> Dirk Eddelbuettel | <a href="mailto:edd@debian.org">edd@debian.org</a> | <a href="http://dirk.eddelbuettel.com" target="_blank">http://dirk.eddelbuettel.com</a><br>
> _______________________________________________<br>
> Rcpp-devel mailing list<br>
> <a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
> <a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
</div></div></blockquote></div><br>