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