I guess BLAS plays a minor role here. The operation here is a sparse Cholesky decomposition, and if I'm correct the Matrix package in R uses SuiteSparse/CHOLMOD to do the factorization, while Eigen uses its own implementation. The underlying numerical algorithms may be different.
Best, Yixuan On Mon, Nov 26, 2018 at 12:23 PM Hoffman, Gabriel <gabriel.hoff...@mssm.edu> wrote: > 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? > > > library(Matrix) > library(inline) > > # construct sparse matrix > ######################### > > # construct a matrix C that is N x N with S total entries > # set C = crossprod(X) > N = 100000 > S = 1000000 > i = sample(1:1000, S, replace=TRUE) > j = sample(1:1000, S, replace=TRUE) > values = runif(S, 0, .3) > X = sparseMatrix(i=i, j=j, x = values, symmetric=FALSE ) > > C = as(crossprod(X), "dgCMatrix") > > # check sparsity fraction > S / N^2 > > # define RCppEigen code > CholeskyCppSparse<-' > using Rcpp::as; > using Eigen::Map; > using Eigen::SparseMatrix; > using Eigen::MappedSparseMatrix; > using Eigen::SimplicialLLT; > > // get data into RcppEigen > const MappedSparseMatrix<double> Sigma(as<MappedSparseMatrix<double> > >(Sigma_in)); > > // compute Cholesky > typedef SimplicialLLT<SparseMatrix<double> > SpChol; > const SpChol Ch(Sigma); > ' > > CholSparse <- cxxfunction(signature(Sigma_in = "dgCMatrix"), > CholeskyCppSparse, plugin = "RcppEigen") > > # compare times > system.time(replicate(10, chol( C ))) > # output: > # user system elapsed > # 0.341 0.014 0.355 > > system.time(replicate(10, CholSparse( C ))) > # output: > # user system elapsed > # 1.639 0.046 1.687 > > sessionInfo() > > R version 3.5.1 (2018-07-02) > Platform: x86_64-apple-darwin15.6.0 (64-bit) > Running under: macOS 10.14 > > Matrix products: default > BLAS: > /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib > LAPACK: > /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib > > locale: > [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 > > attached base packages: > [1] stats graphics grDevices datasets utils methods base > > other attached packages: > [1] inline_0.3.15 Matrix_1.2-15 > > loaded via a namespace (and not attached): > [1] compiler_3.5.1 RcppEigen_0.3.3.4.0 Rcpp_1.0.0 > [4] grid_3.5.1 lattice_0.20-38 > > Changing the size of the matrix and the number of entries does not > change the relative times much > > Thanks, > - Gabriel > _______________________________________________ > 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
_______________________________________________ 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
