On Fri, 18 Oct 2024 at 18:39, Robin Liu <robin28...@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@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
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
