<div dir="ltr"><div dir="ltr"><br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, 18 Oct 2024 at 18:39, Robin Liu <<a href="mailto:robin28liu@gmail.com">robin28liu@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr">Dear all,<div><br></div><div>Computing the Cholesky factor of a matrix is on the order of n^3 while inverting a triangular system is n^2. However, it appears that these two operations take roughly the same time using armadillo.</div><div><br></div><div>Unit: milliseconds <br> ticker mean sd min max neval<br> chol 80.81 3.001 75.84 88.11 20<br> inv-chol 84.89 4.449 80.78 97.30 20<br> inv-direct 274.39 10.787 260.85 303.20 20<br></div><div><br></div><div>Inverting the triangular Cholesky factor is fast relative to direct inversion, but is of the same order as computing the factor itself. Is this to be expected?</div></div></div></blockquote><div><br></div><div>Maybe inverting a triangular matrix is not as optimized as the Cholesky decomposition in OpenBLAS? We can only guess here. This question should be addressed upstream to the OpenBLAS maintainers. </div><div><br></div><div>Best,</div><div>Iñaki</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr"><div>I am running with <br>BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 <br>LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/<a href="http://libopenblasp-r0.3.20.so" target="_blank">libopenblasp-r0.3.20.so</a>; LAPACK version 3.10.0 <br></div><div>along with OPENBLAS_NUM_THREADS=1.<br></div><br><div>I admit my question is not specific to Rcpp or RcppArmadillo, but I thought someone on this list could advise. </div><div><br></div><div>My code follows:</div><div>```</div><div>library(Rcpp)<br>library(RcppArmadillo)<br><br>src <-<br> r"(#include <RcppArmadillo.h><br>#include "/usr/local/lib/R/site-library/RcppClock/include/RcppClock.h"<br>using namespace arma;<br><br>// [[Rcpp:depends(RcppClock)]]<br>// [[Rcpp::depends(RcppArmadillo)]]<br>// [[Rcpp::export]]<br>void chol_and_inv(const arma::mat &V) {<br> Rcpp::Clock clock;<br> mat R, Rinv, Vinv;<br> for (int i = 0; i < 20; i++) {<br> clock.tick("chol");<br> R = chol(V);<br> clock.tock("chol");<br> clock.tick("inv-chol");<br> Rinv = inv(trimatu(R));<br> clock.tock("inv-chol");<br> clock.tick("inv-direct");<br> Vinv = inv(V);<br> clock.tock("inv-direct");<br> }<br> clock.stop("timer");<br> return;<br>})"<br><br># Create large PD matrix<br>size <- 2000<br>set.seed(123)<br>m <- matrix(rnorm(size^2), size, size)<br>m <- tcrossprod(m)<br>diag(m) <- rowSums(abs(m)) + 2<br><br>sourceCpp(code = src)<br>result <- chol_and_inv(m)<br>library(RcppClock)<br>timer</div><div>```</div><div>Best wishes,</div></div>
</div>
_______________________________________________<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" rel="noreferrer" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
</blockquote></div><br clear="all"><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><b>Iñaki Úcar</b><br>Assistant Professor of Statistics<div>Director of the <a href="https://www.uc3m.es/master/computational-social-science" target="_blank">Master in Computational Social Science</a></div><div><br></div><div>Department of Statistics | Big Data Institute<br><div>Universidad Carlos III de Madrid</div><div>Av. de la Universidad 30, 28911 Leganés, Spain<br>Office: 7.3.J25, Tel: +34 916248804</div></div></div></div></div>