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

[Iñaki Ucar]

On Fri, 18 Oct 2024 at 18:39, Robin Liu <robin28liu at gmail.com> wrote:

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

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.

Best,
Iñaki


> 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,
> _______________________________________________
> Rcpp-devel mailing list
> Rcpp-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>


-- 
*Iñaki Úcar*
Assistant Professor of Statistics
Director of the Master in Computational Social Science
<https://www.uc3m.es/master/computational-social-science>

Department of Statistics | Big Data Institute
Universidad Carlos III de Madrid
Av. de la Universidad 30, 28911 Leganés, Spain
Office: 7.3.J25, Tel: +34 916248804
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20241020/8eb9f444/attachment.htm>