<div dir="ltr">Hi Hao,<div><br></div><div> from what I understand Armadillo doesn't have a forward or back-solver.</div><div><br></div><div>In any case it's pretty easy to create one, for example you can have a look here:</div>
<div><br></div><div><a href="https://github.com/mfasiolo/mvnfast/blob/master/src/mahaCpp.cpp">https://github.com/mfasiolo/mvnfast/blob/master/src/mahaCpp.cpp</a><br></div><div><br></div><div>Have a nice weekend,</div><div>
<br></div><div>Matteo</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Sat, Jun 28, 2014 at 12:25 AM, Hao Ye <span dir="ltr"><<a href="mailto:hye@ucsd.edu" target="_blank">hye@ucsd.edu</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Peter,<br>
<br>
I think Soren had a similar issue in February (about Armadillo not having an optimized solver for triangular matrices):<br>
<a href="http://lists.r-forge.r-project.org/pipermail/rcpp-devel/2014-February/007147.html" target="_blank">http://lists.r-forge.r-project.org/pipermail/rcpp-devel/2014-February/007147.html</a><br>
<br>
I don't think there was ever a reply though...<br>
<br>
In contrast, the notes for backsolve() show:<br>
> This is a wrapper for the level-3 BLAS routine dtrsm.<br>
<br>
which *is* specific for triangular matrices. So I'm not surprised that backsolve() is faster than RcppArmadillo's solve() in your situation.<br>
<br>
Best,<br>
--<br>
Hao Ye<br>
<a href="mailto:hye@ucsd.edu">hye@ucsd.edu</a><br>
<div><div class="h5"><br>
On Jun 27, 2014, at 4:14 PM, Peter Rossi <<a href="mailto:perossichi@gmail.com">perossichi@gmail.com</a>> wrote:<br>
<br>
> Folks-<br>
><br>
> In my package, bayesm, I use backsolve() to invert upper-triangular<br>
> arrays. I am in the process of converting my package to rccp-arma.<br>
><br>
> I thought I would test the analogous operation in arma by declaring<br>
> the matrix as upper-triangular and inverting using solve().<br>
><br>
> To my surprise, pure R code was considerably faster than<br>
> rcpp-armadillo for 50 x 50 and larger matrices.<br>
><br>
> I attach an rmd and cpp files necessary to run this benchmark.<br>
><br>
> I assume I'm doing something terribly wrong or naive in my use of<br>
> rcpp-armadillo and would appreciate any thoughts on what causes these<br>
> unexpected results.<br>
><br>
> p<br>
><br>
> Here are the results for those who do not want to compile the rmd file:<br>
><br>
> 10 by 10 uppertriangular matrix;<br>
><br>
> # 10 by 10<br>
> set.seed(777)<br>
> k <- 10<br>
> m <- k + 10<br>
> A <- matrix(rnorm(k*m), m, k)<br>
> R <- chol(t(A)%*%A)<br>
><br>
> rep <- 500<br>
> microbenchmark(backsolve_R(rep, R),backsolve_rcpp(rep, R),times=10)<br>
><br>
> ## Unit: microseconds<br>
> ## expr min lq median uq max neval<br>
> ## backsolve_R(rep, R) 3781.3 4002.4 4131.0 4958.8 6146.3 10<br>
> ## backsolve_rcpp(rep, R) 809.1 812.5 825.1 851.4 889.7 10<br>
><br>
> 50 by 50 uppertriangular matrix<br>
><br>
> # 50 by 50<br>
> set.seed(777)<br>
> k <- 50<br>
> m <- k + 10<br>
> A <- matrix(rnorm(k*m), m, k)<br>
> R <- chol(t(A)%*%A)<br>
><br>
> rep <- 500<br>
> microbenchmark(backsolve_R(rep, R),backsolve_rcpp(rep, R),times=10)<br>
><br>
> ## Unit: milliseconds<br>
> ## expr min lq median uq max neval<br>
> ## backsolve_R(rep, R) 18.92 19.40 19.97 20.94 44.22 10<br>
> ## backsolve_rcpp(rep, R) 36.67 36.88 37.13 37.28 37.41 10<br>
><br>
> 100 by 100 uppertriangular matrix<br>
><br>
> # 100 by 100<br>
> set.seed(777)<br>
> k <- 100<br>
> m <- k+10<br>
> A <- matrix(rnorm(k*m), m, k)<br>
> R <- chol(t(A)%*%A)<br>
><br>
> rep <- 500<br>
> microbenchmark(backsolve_R(rep, R),backsolve_rcpp(rep,R),times=10)<br>
><br>
> ## Unit: milliseconds<br>
> ## expr min lq median uq max neval<br>
> ## backsolve_R(rep, R) 101.5 102.4 103.2 104.4 126.8 10<br>
> ## backsolve_rcpp(rep, R) 251.6 251.8 251.8 252.0 254.0 10<br>
><br>
><br>
><br>
> --<br>
> Peter E. Rossi<br>
</div></div>> <backsolve_test.cpp><backsolve_test.Rmd>_______________________________________________<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>
_______________________________________________<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>
</blockquote></div><br></div>