[Rcpp-devel] Speed of RCppEigen Cholesky decomposition on sparse matrix

Le 26/11/2018 ? 18:23, Hoffman, Gabriel a ?crit :

I am developing a statistical model and I have a prototype working in R
code.  I make extensive use of sparse matrices, so the R code is pretty
fast, but hoped that using RCppEigen to evaluate the log-likelihood
function could avoid a lot of memory copying and be substantially
faster.  However, in a simple  example I am seeing that RCppEigen is
3-5x slower than standard R code for cholesky decomposition of a sparse
matrix.  This is the case on R 3.5.1 using RcppEigen_0.3.3.4.0 on both
OS X and CentOS 6.9.

Since this simple operation is so much slower it doesn't seem like
using RCppEigen is worth it in this case.  Is this an issue with BLAS,
some libraries or compiler options, or is R code really the fastest
option?

After few checks, it seems to be a test issue. Matrix package stores the
decomposition somewhere in attributes of the submitted matrix. So the
the repetitive calls requiring chol() decomposition are not really doing
the job. Instead, previously stored result is reused.

You can easily convince yourself by "modifying" the matrix C (and thus
invalidating previous decomposition) by doing something like "C + 0." :

system.time(replicate(10, chol( C )))
#utilisateur     syst?me      ?coul?
#      0.459       0.011       0.471
system.time(replicate(10, chol( C+0. )))
#utilisateur     syst?me      ?coul?
#      5.365       0.060       5.425
system.time(replicate(10, CholSparse( C+0. )))
#utilisateur     syst?me      ?coul?
#      3.627       0.035       3.665

On my machine, I have almost identical timing for CholSparse() with or
without C modification:

system.time(replicate(10, CholSparse( C )))
#utilisateur     syst?me      ?coul?
#      3.283       0.004       3.289
which proves that Eigen doesn't store the decomposition for future reuse
and redo the decomposition at each call on the same matrix.

Best,
Serguei.