On 18 October 2024 at 09:39, Robin Liu 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? 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.
The point of OpenBLAS is to be multithreaded. Why do you constrain it to one
thread?
| I admit my question is not specific to Rcpp or RcppArmadillo, but I thought
| someone on this list could advise.
What is the reference point? "Slow" is relative to what here? Did you try to
vary the matrix dimension too and see if the change in time is aligned with
your expectations?
In any event, Rcpp "is a messenger here" as is RcppArmadillo as is Armadillo
as the actual operations are plain BLAS operations and this list may no be
the most focussed forum to discuss the topic.
Dirk
| 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
--
dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org
_______________________________________________
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
