That's useful to know about what svd_econ() exactly does. I will give that a shot and report back.<br><br>Has
anyone looked at integrating something like SLEPc, Anasazi(via
Trilinos) or ARPACK++ into rcpp? These would be some really cool tools
to have available. <br>
<br><a href="http://www.grycap.upv.es/slepc/description/summary.htm" target="_blank">http://www.grycap.upv.es/slepc/description/summary.htm</a><br><br><a href="http://trilinos.sandia.gov/packages/anasazi/" target="_blank">http://trilinos.sandia.gov/packages/anasazi/</a><br>
<br><a href="http://www.ime.unicamp.br/%7Echico/arpack++/" target="_blank">http://www.ime.unicamp.br/~chico/arpack++/</a><br><br>Someone actually wrote a wrapper for ARPACK++ for Eigen<br><br><a href="https://github.com/beam2d/arpaca/blob/master/README.md" target="_blank">https://github.com/beam2d/arpaca/blob/master/README.md</a><div class="yj6qo ajU">
<div id=":ll" class="ajR" tabindex="0"><img class="ajT" src="images/cleardot.gif"></div></div><br><br><div class="gmail_quote">On Thu, Jun 14, 2012 at 9:04 AM, Conrad Sand <span dir="ltr"><<a href="mailto:conradsand@gmail.com" target="_blank">conradsand@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="im"><p>On Jun 15, 2012 12:11 AM, "Dirk Eddelbuettel" <<a href="mailto:edd@debian.org" target="_blank">edd@debian.org</a>> wrote:<br>
> Thanks for that earlier hint re 'thin' and 'full' SVDs.</p>
</div><p>armadillo has the standard svd() and the thin version too: svd_econ().</p><div class="im">
<p>> Conrad, any interest in switching to dgesdd?</p>
</div><p>yes, but probably not simply changing svd() directly to dgedd(). lapack gave the function a different name for a reason: the results may be slightly different. if armadillo starts giving different results all of a sudden, there would be a lot of displeased people. a big no no, given that armadillo is used for critical stuff.</p>
<p>I'll probably add an option to svd() to optionally use dgedd().</p>
<p>I've done something very similar for eig_sym() in armadillo 3.2, where an alternative faster algorithm for eigen decomposition can be optionally used.<br><br></p><div class="HOEnZb"><div class="h5">
<p>> Dirk, at useR and across the room from Doug<br>
><br>
> | > On Wed, Jun 13, 2012 at 4:07 PM, Douglas Bates <<a href="mailto:bates@stat.wisc.edu" target="_blank">bates@stat.wisc.edu</a>> wrote:<br>
> | >><br>
> | >> On Wed, Jun 13, 2012 at 5:16 PM, Dirk Eddelbuettel <<a href="mailto:edd@debian.org" target="_blank">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<br>
> | >> > 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<br>
> | >> > session<br>
> | >> > | | uses. So MKL, OpenBlas, ... are all options. Eigen actually has its<br>
> | >> > 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?<br>
> | >> > 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<br>
> | >> > 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>
> | >> 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>
> | >><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" target="_blank">edd@debian.org</a>><br>
> | >> > 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<br>
> | >> > few days<br>
> | >> > | and<br>
> | >> > | | don't know which library to continue using. RcppEigen seems<br>
> | >> > really slick,<br>
> | >> > | but<br>
> | >> > | | appears to be lacking some of the decompositions I want and<br>
> | >> > isn't nearly<br>
> | >> > | as<br>
> | >> > | | fast to code. RcppArmadillo seems about as fast, easier to code<br>
> | >> > 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<br>
> | >> > find<br>
> | >> > | easier to get my head around. Doug, however, moved from<br>
> | >> > [Rcpp]Armadillo<br>
> | >> > | to<br>
> | >> > | [Rcpp]Eigen as it has some things he needs. Eigen should have a<br>
> | >> > "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<br>
> | >> > went as far<br>
> | >> > | as building 'robustHD' using RcppArmadillo with a drop-in for<br>
> | >> > RcppEigen (in<br>
> | >> > | package 'sparseLTSEigen'; both package names from memmory and I<br>
> | >> > may have<br>
> | >> > | mistyped). He reported a performance gain of around 25% for his<br>
> | >> > 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<br>
> | >> > RcppEigen<br>
> | >> > | to<br>
> | >> > | | call LAPACK functions? Want some other decomposition methods,<br>
> | >> > dont like<br>
> | >> > | the<br>
> | >> > | | JacobiSVD method in Eigen.<br>
> | >> > |<br>
> | >> > | You need to differentiate between the Eigen and Armadillo docs<br>
> | >> > _for their<br>
> | >> > | libraries_ and what happens when you access the Rcpp* variant from<br>
> | >> > R.<br>
> | >> > |<br>
> | >> > | RcppArmadillo will use the very same LAPACK and BLAS libs your R<br>
> | >> > session<br>
> | >> > | uses. So MKL, OpenBlas, ... are all options. Eigen actually has<br>
> | >> > 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>
> | >> > | | _______________________________________________<br>
> | >> > | | Rcpp-devel mailing list<br>
> | >> > | | <a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank">Rcpp-devel@lists.r-forge.r-project.org</a><br>
> | >> > | |<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" target="_blank">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" target="_blank">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" target="_blank">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>
> | ><br>
><br>
> --<br>
> Dirk Eddelbuettel | <a href="mailto:edd@debian.org" target="_blank">edd@debian.org</a> | <a href="http://dirk.eddelbuettel.com" target="_blank">http://dirk.eddelbuettel.com</a><br>
</p>
</div></div></blockquote></div><br>