But adding 0 to a sparse matrix is expensive operation. It doesn't look fair to include it to benchmark.
вт, 27 нояб. 2018 г., 20:07 rcpp-devel-requ...@lists.r-forge.r-project.org: > Send Rcpp-devel mailing list submissions to > rcpp-devel@lists.r-forge.r-project.org > > To subscribe or unsubscribe via the World Wide Web, visit > > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel > > or, via email, send a message with subject or body 'help' to > rcpp-devel-requ...@lists.r-forge.r-project.org > > You can reach the person managing the list at > rcpp-devel-ow...@lists.r-forge.r-project.org > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of Rcpp-devel digest..." > > > Today's Topics: > > 1. Re: Speed of RCppEigen Cholesky decomposition on sparse > matrix (Serguei Sokol) > 2. Problems with Rcout (Barth Riley) > 3. Re: Problems with Rcout (I?aki Ucar) > 4. Re: Problems with Rcout (Serguei Sokol) > 5. Re: Problems with Rcout (Barth Riley) > > > ---------------------------------------------------------------------- > > Message: 1 > Date: Tue, 27 Nov 2018 15:33:55 +0100 > From: Serguei Sokol <serguei.so...@gmail.com> > To: "Hoffman, Gabriel" <gabriel.hoff...@mssm.edu>, > "rcpp-devel@lists.r-forge.r-project.org" > <rcpp-devel@lists.r-forge.r-project.org> > Subject: Re: [Rcpp-devel] Speed of RCppEigen Cholesky decomposition on > sparse matrix > Message-ID: <91e163f1-c169-217f-739b-2a379c63c...@gmail.com> > Content-Type: text/plain; charset=utf-8; format=flowed > > 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. > > > > > > > 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 > > > > > > ------------------------------ > > Message: 2 > Date: Tue, 27 Nov 2018 09:35:10 -0600 > From: Barth Riley <barthri...@comcast.net> > To: " rcpp-devel@lists.r-forge.r-project.org" > <rcpp-devel@lists.r-forge.r-project.org> > Subject: [Rcpp-devel] Problems with Rcout > Message-ID: > < > mailman.18867.1543334818.2604.rcpp-de...@lists.r-forge.r-project.org> > Content-Type: text/plain; charset="utf-8" > > Dear Rcppers > > I am encountering a problem with Rcout. Even with basic string output, > when I run the function containing the call to Rcout, no output is > generated. Here is a simple example of what I?m trying to do: > > // [[Rcpp::export]] > void testFunc() { > Rcpp::Rcout << "testFunc begins" << std:endl; > . . . > } > > My code is part of a package I?m developing, using R 3.51 and Rcpp > 0.12.19. The Rcpp code compiles without a problem. > > Thanks > > Barth > > > -------------- next part -------------- > An HTML attachment was scrubbed... > URL: < > http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20181127/6adbedfe/attachment.html > > > > ------------------------------ > > Message: 3 > Date: Tue, 27 Nov 2018 16:50:13 +0100 > From: I?aki Ucar <iu...@fedoraproject.org> > To: barthri...@comcast.net > Cc: "<rcpp-devel@lists.r-forge.r-project.org>" > <rcpp-devel@lists.r-forge.r-project.org> > Subject: Re: [Rcpp-devel] Problems with Rcout > Message-ID: > < > calexwq3ttxhgpidxiswclfvhb4-f2rdamn12q+cgemyans8...@mail.gmail.com> > Content-Type: text/plain; charset="UTF-8" > > On Tue, 27 Nov 2018 at 16:35, Barth Riley <barthri...@comcast.net> wrote: > > > > Dear Rcppers > > > > I am encountering a problem with Rcout. Even with basic string output, > when I run the function containing the call to Rcout, no output is > generated. Here is a simple example of what I?m trying to do: > > > > // [[Rcpp::export]] > > void testFunc() { > > Rcpp::Rcout << "testFunc begins" << std:endl; > > . . . > > } > > Note that it should be "std::endl", with double colon. Assuming this > is a transcription error, you'll have to provide more context (and, > ideally, some kind of reproducible example), because this works just > fine. > > I?aki > > > ------------------------------ > > Message: 4 > Date: Tue, 27 Nov 2018 16:51:28 +0100 > From: Serguei Sokol <serguei.so...@gmail.com> > To: rcpp-devel@lists.r-forge.r-project.org > Subject: Re: [Rcpp-devel] Problems with Rcout > Message-ID: <ccc55a85-79cb-c3c6-14e8-222223abe...@gmail.com> > Content-Type: text/plain; charset=utf-8; format=flowed > > Le 27/11/2018 ? 16:35, Barth Riley a ?crit?: > > Dear Rcppers > > > > I am encountering a problem with Rcout. Even with basic string output, > > when I run the function containing the call to Rcout, no output is > > generated. Here is a simple example of what I?m trying to do: > > > > // [[Rcpp::export]] > > > > void testFunc() { > > > > ?????????? Rcpp::Rcout << "testFunc begins" << std:endl; > > > > ?????????? . . . > > > > } > > > This example works for me: > > library(Rcpp) > sourceCpp(code="#include <Rcpp.h>\n// [[Rcpp::export]]\nvoid testFunc() > {\nRcpp::Rcout << \"testFunc begins\" << std::endl;\n}") > testFunc() > #testFunc begins > > May be in your session you have redirected stdout elsewhere? > > Serguei. > > > My code is part of a package I?m developing, using R 3.51 and Rcpp > > 0.12.19. The Rcpp code compiles without a problem. > > > > Thanks > > > > Barth > > > > > > > > _______________________________________________ > > 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 > > > > > > ------------------------------ > > Message: 5 > Date: Tue, 27 Nov 2018 10:06:47 -0600 > From: Barth Riley <barthri...@comcast.net> > To: I?aki Ucar <iu...@fedoraproject.org>, " > rcpp-devel@lists.r-forge.r-project.org" > <rcpp-devel@lists.r-forge.r-project.org> > Subject: Re: [Rcpp-devel] Problems with Rcout > Message-ID: > < > mailman.18868.1543334818.2604.rcpp-de...@lists.r-forge.r-project.org> > Content-Type: text/plain; charset="utf-8" > > Here is a more complete example. Note that I want to output text strings > for debugging purposes as the code for treatAsVector = true is never > executed. > > Barth > > NumericVector getValidCount(Rcpp::NumericMatrix m, > bool treatAsVector) { > > Rcpp::Rcout << "getValidCount BEGINS" << std::endl; > > int N = m.cols(); > NumericVector u, vec; > NumericVector count (N); > > if(!treatAsVector) { > Rcpp::Rcout << "Treating as matrix" << std::endl; > for(int i = 0; i < N; i++) { > vec = m(_,i); > vec = vec[!Rcpp::is_na(vec)]; > u = Rcpp::unique(vec); > count[i] = u.length(); > } > } else { > Rcpp::Rcout << "treating as vector" << std::endl; > vec = as<NumericVector>(m); > vec = vec[!Rcpp::is_na(vec)]; > u = Rcpp::unique(vec); > > count.fill(u.length()); > } > > return count; > } > -------------- next part -------------- > An HTML attachment was scrubbed... > URL: < > http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20181127/607ae61f/attachment.html > > > > ------------------------------ > > Subject: Digest Footer > > _______________________________________________ > 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 > > ------------------------------ > > End of Rcpp-devel Digest, Vol 109, Issue 11 > ******************************************* >
_______________________________________________ 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