<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
</head>
<body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space; color: rgb(0, 0, 0); font-size: 14px; font-family: Calibri, sans-serif;">
<div>
<div style="font-family: -webkit-standard;">I am developing a statistical model and I have a prototype working in R</div>
<div style="font-family: -webkit-standard;">code. I make extensive use of sparse matrices, so the R code is pretty</div>
<div style="font-family: -webkit-standard;">fast, but hoped that using RCppEigen to evaluate the log-likelihood</div>
<div style="font-family: -webkit-standard;">function could avoid a lot of memory copying and be substantially</div>
<div style="font-family: -webkit-standard;">faster. However, in a simple example I am seeing that RCppEigen is</div>
<div style="font-family: -webkit-standard;">3-5x slower than standard R code for cholesky decomposition of a sparse</div>
<div style="font-family: -webkit-standard;">matrix. This is the case on R 3.5.1 using RcppEigen_0.3.3.4.0 on both</div>
<div style="font-family: -webkit-standard;">OS X and CentOS 6.9.</div>
<div style="font-family: -webkit-standard;"><br>
</div>
<div style="font-family: -webkit-standard;">Since this simple operation is so much slower it doesn't seem like</div>
<div style="font-family: -webkit-standard;">using RCppEigen is worth it in this case. Is this an issue with BLAS,</div>
<div style="font-family: -webkit-standard;">some libraries or compiler options, or is R code really the fastest</div>
<div style="font-family: -webkit-standard;">option?</div>
</div>
<div style="font-family: -webkit-standard;"><br>
</div>
<div style="font-family: -webkit-standard;">
<div><br>
</div>
<div>library(Matrix)</div>
<div>library(inline)</div>
<div><br>
</div>
<div># construct sparse matrix</div>
<div>#########################</div>
<div><br>
</div>
<div># construct a matrix C that is N x N with S total entries</div>
<div># set C = crossprod(X)</div>
<div>N = 100000</div>
<div>S = 1000000</div>
<div>i = sample(1:1000, S, replace=TRUE)</div>
<div>j = sample(1:1000, S, replace=TRUE)</div>
<div>values = runif(S, 0, .3)</div>
<div>X = sparseMatrix(i=i, j=j, x = values, symmetric=FALSE )</div>
<div><br>
</div>
<div>C = as(crossprod(X), "dgCMatrix")</div>
<div><br>
</div>
<div># check sparsity fraction</div>
<div>S / N^2</div>
<div><br>
</div>
<div># define RCppEigen code</div>
<div>CholeskyCppSparse<-' </div>
<div>using Rcpp::as; </div>
<div>using Eigen::Map;</div>
<div>using Eigen::SparseMatrix;</div>
<div>using Eigen::MappedSparseMatrix;</div>
<div>using Eigen::SimplicialLLT;</div>
<div><br>
</div>
<div>// get data into RcppEigen</div>
<div>const MappedSparseMatrix<double> Sigma(as<MappedSparseMatrix<double> >(Sigma_in));</div>
<div><br>
</div>
<div>// compute Cholesky</div>
<div>typedef SimplicialLLT<SparseMatrix<double> > SpChol;</div>
<div>const SpChol Ch(Sigma);</div>
<div>'</div>
<div><br>
</div>
<div>CholSparse <- cxxfunction(signature(Sigma_in = "dgCMatrix"), CholeskyCppSparse, plugin = "RcppEigen")</div>
<div><br>
</div>
<div>
<div># compare times</div>
<div>system.time(replicate(10, chol( C )))</div>
<div># output:</div>
<div># user system elapsed</div>
<div># 0.341 0.014 0.355</div>
<div><br>
</div>
<div>system.time(replicate(10, CholSparse( C )))</div>
<div># output:</div>
<div># user system elapsed</div>
<div># 1.639 0.046 1.687</div>
<div><br>
</div>
<blockquote id="MAC_OUTLOOK_ATTRIBUTION_BLOCKQUOTE" style="padding: 0px 0px 0px 5px; margin: 0px 0px 0px 5px;">
sessionInfo()</blockquote>
<div>R version 3.5.1 (2018-07-02)</div>
<div>Platform: x86_64-apple-darwin15.6.0 (64-bit)</div>
<div>Running under: macOS 10.14</div>
<div><br>
</div>
<div>Matrix products: default</div>
<div>BLAS:</div>
<div>/Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib</div>
<div>LAPACK:</div>
<div>/Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib</div>
<div><br>
</div>
<div>locale:</div>
<div>[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8</div>
<div><br>
</div>
<div>attached base packages:</div>
<div>[1] stats graphics grDevices datasets utils methods base</div>
<div><br>
</div>
<div>other attached packages:</div>
<div>[1] inline_0.3.15 Matrix_1.2-15</div>
<div><br>
</div>
<div>loaded via a namespace (and not attached):</div>
<div>[1] compiler_3.5.1 RcppEigen_0.3.3.4.0 Rcpp_1.0.0</div>
<div>[4] grid_3.5.1 lattice_0.20-38</div>
<div><br>
</div>
<div>Changing the size of the matrix and the number of entries does not</div>
<div>change the relative times much</div>
</div>
<div><br>
</div>
<div>Thanks, </div>
<div>- Gabriel</div>
</div>
</body>
</html>