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-requ...@lists.r-forge.r-project.org
<mailto:rcpp-devel-requ...@lists.r-forge.r-project.org>:
Send Rcpp-devel mailing list submissions to
rcpp-devel@lists.r-forge.r-project.org
<mailto: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
<mailto: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
<mailto: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
<mailto:serguei.so...@gmail.com>>
To: "Hoffman, Gabriel" <gabriel.hoff...@mssm.edu
<mailto:gabriel.hoff...@mssm.edu>>,
"rcpp-devel@lists.r-forge.r-project.org
<mailto:rcpp-devel@lists.r-forge.r-project.org>"
<rcpp-devel@lists.r-forge.r-project.org
<mailto: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
<mailto: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
<mailto: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
<mailto:barthri...@comcast.net>>
To: " rcpp-devel@lists.r-forge.r-project.org
<mailto:rcpp-devel@lists.r-forge.r-project.org>"
<rcpp-devel@lists.r-forge.r-project.org
<mailto: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 <mailto: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
<mailto:iu...@fedoraproject.org>>
To: barthri...@comcast.net <mailto:barthri...@comcast.net>
Cc: "<rcpp-devel@lists.r-forge.r-project.org
<mailto:rcpp-devel@lists.r-forge.r-project.org>>"
<rcpp-devel@lists.r-forge.r-project.org
<mailto:rcpp-devel@lists.r-forge.r-project.org>>
Subject: Re: [Rcpp-devel] Problems with Rcout
Message-ID:
<calexwq3ttxhgpidxiswclfvhb4-f2rdamn12q+cgemyans8...@mail.gmail.com
<mailto:calexwq3ttxhgpidxiswclfvhb4-f2rdamn12q%2bcgemyans8...@mail.gmail.com>>
Content-Type: text/plain; charset="UTF-8"
On Tue, 27 Nov 2018 at 16:35, Barth Riley <barthri...@comcast.net
<mailto: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
<mailto:serguei.so...@gmail.com>>
To: rcpp-devel@lists.r-forge.r-project.org
<mailto:rcpp-devel@lists.r-forge.r-project.org>
Subject: Re: [Rcpp-devel] Problems with Rcout
Message-ID: <ccc55a85-79cb-c3c6-14e8-222223abe...@gmail.com
<mailto: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
<mailto: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
<mailto:barthri...@comcast.net>>
To: I?aki Ucar <iu...@fedoraproject.org
<mailto:iu...@fedoraproject.org>>, "
rcpp-devel@lists.r-forge.r-project.org
<mailto:rcpp-devel@lists.r-forge.r-project.org>"
<rcpp-devel@lists.r-forge.r-project.org
<mailto: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 <mailto: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
<mailto: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
_______________________________________________
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
