[Rcpp-devel] Speed of RCppEigen Cholesky decomposition on sparse matrix (Serguei Sokol)
Dmitriy Selivanov
selivanov.dmitriy at gmail.com
Tue Nov 27 17:57:07 CET 2018
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-request at lists.r-forge.r-project.org:
> Send Rcpp-devel mailing list submissions to
> rcpp-devel at 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-request at lists.r-forge.r-project.org
>
> You can reach the person managing the list at
> rcpp-devel-owner at 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.sokol at gmail.com>
> To: "Hoffman, Gabriel" <gabriel.hoffman at mssm.edu>,
> "rcpp-devel at lists.r-forge.r-project.org"
> <rcpp-devel at lists.r-forge.r-project.org>
> Subject: Re: [Rcpp-devel] Speed of RCppEigen Cholesky decomposition on
> sparse matrix
> Message-ID: <91e163f1-c169-217f-739b-2a379c63c35c at 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 at 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 <barthriley at comcast.net>
> To: " rcpp-devel at lists.r-forge.r-project.org"
> <rcpp-devel at lists.r-forge.r-project.org>
> Subject: [Rcpp-devel] Problems with Rcout
> Message-ID:
> <
> mailman.18867.1543334818.2604.rcpp-devel at 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 <iucar at fedoraproject.org>
> To: barthriley at comcast.net
> Cc: "<rcpp-devel at lists.r-forge.r-project.org>"
> <rcpp-devel at lists.r-forge.r-project.org>
> Subject: Re: [Rcpp-devel] Problems with Rcout
> Message-ID:
> <
> CALEXWq3tTXhGPidXiswCLfvhB4-f2RDAmn12q+CGeMYaNs8VLg at mail.gmail.com>
> Content-Type: text/plain; charset="UTF-8"
>
> On Tue, 27 Nov 2018 at 16:35, Barth Riley <barthriley at 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.sokol at gmail.com>
> To: rcpp-devel at lists.r-forge.r-project.org
> Subject: Re: [Rcpp-devel] Problems with Rcout
> Message-ID: <ccc55a85-79cb-c3c6-14e8-222223abe2c5 at 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 at 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 <barthriley at comcast.net>
> To: I?aki Ucar <iucar at fedoraproject.org>, "
> rcpp-devel at lists.r-forge.r-project.org"
> <rcpp-devel at lists.r-forge.r-project.org>
> Subject: Re: [Rcpp-devel] Problems with Rcout
> Message-ID:
> <
> mailman.18868.1543334818.2604.rcpp-devel at 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 at 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
> *******************************************
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20181127/20f79db6/attachment-0001.html>
More information about the Rcpp-devel
mailing list