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?
Here is my example:
library(Matrix)
library(inline)
# construct sparse matrix
#########################
# construct a matrix C that is N x X with S total entries
N = 10000
S = 1000000
i = sample(1:1000, S, replace=TRUE)
j = sample(1:1000, S, replace=TRUE)
idx = i >= j
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
Thanks,
- Gabriel
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
