[Rcpp-devel] Inverting a triangular matrix appears slow using Armadillo and openblas

Robin Liu robin28liu at gmail.com
Fri Oct 18 18:39:03 CEST 2024


Dear all,

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.

Unit: milliseconds
     ticker   mean     sd    min    max neval
       chol  80.81  3.001  75.84  88.11    20
   inv-chol  84.89  4.449  80.78  97.30    20
 inv-direct 274.39 10.787 260.85 303.20    20

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
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;
 LAPACK version 3.10.0
along with OPENBLAS_NUM_THREADS=1.

I admit my question is not specific to Rcpp or RcppArmadillo, but I thought
someone on this list could advise.

My code follows:
```
library(Rcpp)
library(RcppArmadillo)

src <-
  r"(#include <RcppArmadillo.h>
#include "/usr/local/lib/R/site-library/RcppClock/include/RcppClock.h"
using namespace arma;

// [[Rcpp:depends(RcppClock)]]
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
void chol_and_inv(const arma::mat &V) {
  Rcpp::Clock clock;
  mat R, Rinv, Vinv;
  for (int i = 0; i < 20; i++) {
    clock.tick("chol");
    R = chol(V);
    clock.tock("chol");
    clock.tick("inv-chol");
    Rinv = inv(trimatu(R));
    clock.tock("inv-chol");
    clock.tick("inv-direct");
    Vinv = inv(V);
    clock.tock("inv-direct");
  }
  clock.stop("timer");
  return;
})"

# Create large PD matrix
size <- 2000
set.seed(123)
m <- matrix(rnorm(size^2), size, size)
m <- tcrossprod(m)
diag(m) <- rowSums(abs(m)) + 2

sourceCpp(code = src)
result <- chol_and_inv(m)
library(RcppClock)
timer
```
Best wishes,
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20241018/dd1fb7f0/attachment.htm>


More information about the Rcpp-devel mailing list