This is really a brilliant observation!
Best, Yixuan On Tue, Nov 27, 2018 at 9:34 AM Serguei Sokol <serguei.so...@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. > _______________________________________________ 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
