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

[Yixuan Qiu]

This is really a brilliant observation!


Best,
Yixuan


On Tue, Nov 27, 2018 at 9:34 AM Serguei Sokol <serguei.sokol at gmail.com> wrote:
>
> 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.
>