<div><div dir="auto">If I recall correctly, Eigen does not use the BLAS, but has reprogrammed the operations in C++. If you want to take advantage of a fast system BLAS, try Armadillo, or in this case RcppArmadillo. </div></div><div dir="auto"><br></div><div dir="auto">Thanks,</div><div dir="auto"><br></div><div dir="auto">Avi</div><div><br><div class="gmail_quote"><div dir="ltr">On Mon, Nov 26, 2018 at 12:23 PM Hoffman, Gabriel <<a href="mailto:gabriel.hoffman@mssm.edu">gabriel.hoffman@mssm.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">



<div style="word-wrap:break-word;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="m_-7268223502459530928MAC_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>
</div>

_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" rel="noreferrer" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a></blockquote></div></div>-- <br><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature">Sent from Gmail Mobile</div>