<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? 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>