[Rcpp-devel] Speed of RCppEigen Cholesky decomposition on sparse matrix (Serguei Sokol)
Serguei Sokol
serguei.sokol at gmail.com
Wed Nov 28 11:05:10 CET 2018
Le 27/11/2018 à 17:57, Dmitriy Selivanov a écrit :
> But adding 0 to a sparse matrix is expensive operation. It doesn't look
> fair to include it to benchmark.
It is true that adding 0 comes at some cost. But qualifying it expensive
or not is kind of subjective opinion. Let see how much does it cost:
system.time(replicate(10, C+0.))
#utilisateur système écoulé
# 0.017 0.030 0.047
So, in my opinion 0.017 s is negligible time laps compared to 5 or even
3 s needed for 10 matrix decompositions. And someone else could say
"it's important".
Anyhow, the goal was to show that Matrix::chol does only 1 decomposition
not 10) while Eigen did all 10 decompositions, and not to thoroughly
compare performances of two methods. The bias detected in the original
test was much higher than the bias induced by adding 0.
Moreover, if adding 0 would be really a problem, one could easily
exclude it from time accounting:
sum(sapply(1:10, function(i) {C=C+0.; system.time(chol(C))[1]}))
#[1] 4.488
sum(sapply(1:10, function(i) {C=C+0.; system.time(CholSparse(C))[1]}))
#[1] 3.341
and once more (to have an idea about time variability)
sum(sapply(1:10, function(i) {C=C+0.; system.time(CholSparse(C))[1]}))
#[1] 3.481
The conclusion remains the same. We see also that time variability is
almost 10 fold higher than time needed to add 0. So I conclude that
adding 0 did not perturbed so much the fairness of the test.
Best,
Serguei.
>
> вт, 27 нояб. 2018 г., 20:07
> rcpp-devel-request at lists.r-forge.r-project.org
> <mailto: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
> <mailto: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
> <mailto: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
> <mailto: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
> <mailto:serguei.sokol at gmail.com>>
> To: "Hoffman, Gabriel" <gabriel.hoffman at mssm.edu
> <mailto:gabriel.hoffman at mssm.edu>>,
> "rcpp-devel at lists.r-forge.r-project.org
> <mailto:rcpp-devel at lists.r-forge.r-project.org>"
> <rcpp-devel at lists.r-forge.r-project.org
> <mailto: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
> <mailto: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
> <mailto: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
> <mailto:barthriley at comcast.net>>
> To: " rcpp-devel at lists.r-forge.r-project.org
> <mailto:rcpp-devel at lists.r-forge.r-project.org>"
> <rcpp-devel at lists.r-forge.r-project.org
> <mailto: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 <mailto: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
> <mailto:iucar at fedoraproject.org>>
> To: barthriley at comcast.net <mailto:barthriley at comcast.net>
> Cc: "<rcpp-devel at lists.r-forge.r-project.org
> <mailto:rcpp-devel at lists.r-forge.r-project.org>>"
> <rcpp-devel at lists.r-forge.r-project.org
> <mailto: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
> <mailto:CALEXWq3tTXhGPidXiswCLfvhB4-f2RDAmn12q%2BCGeMYaNs8VLg 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
> <mailto: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
> <mailto:serguei.sokol at gmail.com>>
> To: rcpp-devel at lists.r-forge.r-project.org
> <mailto: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
> <mailto: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
> <mailto: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
> <mailto:barthriley at comcast.net>>
> To: I?aki Ucar <iucar at fedoraproject.org
> <mailto:iucar at fedoraproject.org>>, "
> rcpp-devel at lists.r-forge.r-project.org
> <mailto:rcpp-devel at lists.r-forge.r-project.org>"
> <rcpp-devel at lists.r-forge.r-project.org
> <mailto: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 <mailto: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
> <mailto: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
> *******************************************
>
>
>
> _______________________________________________
> 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
>
More information about the Rcpp-devel
mailing list