<div dir="auto"><div>But adding 0 to a sparse matrix is expensive operation. It doesn't look fair to include it to benchmark.<br><br><div class="gmail_quote"><div dir="ltr">вт, 27 нояб. 2018 г., 20:07 <a href="mailto:rcpp-devel-request@lists.r-forge.r-project.org">rcpp-devel-request@lists.r-forge.r-project.org</a>:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Send Rcpp-devel mailing list submissions to<br>
<a href="mailto:rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">rcpp-devel@lists.r-forge.r-project.org</a><br>
<br>
To subscribe or unsubscribe via the World Wide Web, visit<br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" rel="noreferrer noreferrer" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
<br>
or, via email, send a message with subject or body 'help' to<br>
<a href="mailto:rcpp-devel-request@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">rcpp-devel-request@lists.r-forge.r-project.org</a><br>
<br>
You can reach the person managing the list at<br>
<a href="mailto:rcpp-devel-owner@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">rcpp-devel-owner@lists.r-forge.r-project.org</a><br>
<br>
When replying, please edit your Subject line so it is more specific<br>
than "Re: Contents of Rcpp-devel digest..."<br>
<br>
<br>
Today's Topics:<br>
<br>
1. Re: Speed of RCppEigen Cholesky decomposition on sparse<br>
matrix (Serguei Sokol)<br>
2. Problems with Rcout (Barth Riley)<br>
3. Re: Problems with Rcout (I?aki Ucar)<br>
4. Re: Problems with Rcout (Serguei Sokol)<br>
5. Re: Problems with Rcout (Barth Riley)<br>
<br>
<br>
----------------------------------------------------------------------<br>
<br>
Message: 1<br>
Date: Tue, 27 Nov 2018 15:33:55 +0100<br>
From: Serguei Sokol <<a href="mailto:serguei.sokol@gmail.com" target="_blank" rel="noreferrer">serguei.sokol@gmail.com</a>><br>
To: "Hoffman, Gabriel" <<a href="mailto:gabriel.hoffman@mssm.edu" target="_blank" rel="noreferrer">gabriel.hoffman@mssm.edu</a>>,<br>
"<a href="mailto:rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">rcpp-devel@lists.r-forge.r-project.org</a>"<br>
<<a href="mailto:rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">rcpp-devel@lists.r-forge.r-project.org</a>><br>
Subject: Re: [Rcpp-devel] Speed of RCppEigen Cholesky decomposition on<br>
sparse matrix<br>
Message-ID: <<a href="mailto:91e163f1-c169-217f-739b-2a379c63c35c@gmail.com" target="_blank" rel="noreferrer">91e163f1-c169-217f-739b-2a379c63c35c@gmail.com</a>><br>
Content-Type: text/plain; charset=utf-8; format=flowed<br>
<br>
Le 26/11/2018 ? 18:23, Hoffman, Gabriel a ?crit?:<br>
> I am developing a statistical model and I have a prototype working in R<br>
> code.??I make extensive use of sparse matrices, so the R code is pretty<br>
> fast, but hoped that using RCppEigen to evaluate the log-likelihood<br>
> function could avoid a lot of memory copying and be substantially<br>
> faster.??However, in a simple??example I am seeing that RCppEigen is<br>
> 3-5x slower than standard R code for cholesky decomposition of a sparse<br>
> matrix.??This is the case on R 3.5.1 using RcppEigen_0.3.3.4.0 on both<br>
> OS X and CentOS 6.9.<br>
> <br>
> Since this simple operation is so much slower it doesn't seem like<br>
> using RCppEigen is worth it in this case.??Is this an issue with BLAS,<br>
> some libraries or compiler options, or is R code really the fastest<br>
> option?<br>
After few checks, it seems to be a test issue. Matrix package stores the <br>
decomposition somewhere in attributes of the submitted matrix. So the <br>
the repetitive calls requiring chol() decomposition are not really doing <br>
the job. Instead, previously stored result is reused.<br>
<br>
You can easily convince yourself by "modifying" the matrix C (and thus <br>
invalidating previous decomposition) by doing something like "C + 0." :<br>
<br>
system.time(replicate(10, chol( C )))<br>
#utilisateur syst?me ?coul?<br>
# 0.459 0.011 0.471<br>
system.time(replicate(10, chol( C+0. )))<br>
#utilisateur syst?me ?coul?<br>
# 5.365 0.060 5.425<br>
system.time(replicate(10, CholSparse( C+0. )))<br>
#utilisateur syst?me ?coul?<br>
# 3.627 0.035 3.665<br>
<br>
On my machine, I have almost identical timing for CholSparse() with or <br>
without C modification:<br>
<br>
system.time(replicate(10, CholSparse( C )))<br>
#utilisateur syst?me ?coul?<br>
# 3.283 0.004 3.289<br>
which proves that Eigen doesn't store the decomposition for future reuse <br>
and redo the decomposition at each call on the same matrix.<br>
<br>
Best,<br>
Serguei.<br>
<br>
> <br>
> <br>
> library(Matrix)<br>
> library(inline)<br>
> <br>
> # construct sparse matrix<br>
> #########################<br>
> <br>
> # construct a matrix C that is N x N with S total entries<br>
> # set C = crossprod(X)<br>
> N = 100000<br>
> S = 1000000<br>
> i = sample(1:1000, S, replace=TRUE)<br>
> j = sample(1:1000, S, replace=TRUE)<br>
> values = runif(S, 0, .3)<br>
> X = sparseMatrix(i=i, j=j, x = values, symmetric=FALSE )<br>
> <br>
> C = as(crossprod(X), "dgCMatrix")<br>
> <br>
> # check sparsity fraction<br>
> S / N^2<br>
> <br>
> # define RCppEigen code<br>
> CholeskyCppSparse<-'<br>
> using Rcpp::as;<br>
> using Eigen::Map;<br>
> using Eigen::SparseMatrix;<br>
> using Eigen::MappedSparseMatrix;<br>
> using Eigen::SimplicialLLT;<br>
> <br>
> // get data into RcppEigen<br>
> const MappedSparseMatrix<double> Sigma(as<MappedSparseMatrix<double> <br>
> >(Sigma_in));<br>
> <br>
> // compute Cholesky<br>
> typedef SimplicialLLT<SparseMatrix<double> > SpChol;<br>
> const SpChol Ch(Sigma);<br>
> '<br>
> <br>
> CholSparse <- cxxfunction(signature(Sigma_in = "dgCMatrix"), <br>
> CholeskyCppSparse, plugin = "RcppEigen")<br>
> <br>
> # compare times<br>
> system.time(replicate(10, chol( C )))<br>
> # output:<br>
> #?? user??system elapsed<br>
> #??0.341?? 0.014?? 0.355<br>
> <br>
> system.time(replicate(10, CholSparse( C )))<br>
> # output:<br>
> #?? user??system elapsed<br>
> # 1.639?? 0.046?? 1.687<br>
> <br>
> sessionInfo()<br>
> <br>
> R version 3.5.1 (2018-07-02)<br>
> Platform: x86_64-apple-darwin15.6.0 (64-bit)<br>
> Running under: macOS??10.14<br>
> <br>
> Matrix products: default<br>
> BLAS:<br>
> /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib<br>
> LAPACK:<br>
> /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib<br>
> <br>
> locale:<br>
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8<br>
> <br>
> attached base packages:<br>
> [1] stats???? graphics??grDevices datasets??utils???? methods?? base<br>
> <br>
> other attached packages:<br>
> [1] inline_0.3.15 Matrix_1.2-15<br>
> <br>
> loaded via a namespace (and not attached):<br>
> [1] compiler_3.5.1??????RcppEigen_0.3.3.4.0 Rcpp_1.0.0<br>
> [4] grid_3.5.1??????????lattice_0.20-38<br>
> <br>
> Changing the size of the matrix and the number of entries does not<br>
> change the relative times much<br>
> <br>
> Thanks,<br>
> - Gabriel<br>
> <br>
> <br>
> _______________________________________________<br>
> Rcpp-devel mailing list<br>
> <a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">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 noreferrer" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
> <br>
<br>
<br>
<br>
------------------------------<br>
<br>
Message: 2<br>
Date: Tue, 27 Nov 2018 09:35:10 -0600<br>
From: Barth Riley <<a href="mailto:barthriley@comcast.net" target="_blank" rel="noreferrer">barthriley@comcast.net</a>><br>
To: " <a href="mailto:rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">rcpp-devel@lists.r-forge.r-project.org</a>"<br>
<<a href="mailto:rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">rcpp-devel@lists.r-forge.r-project.org</a>><br>
Subject: [Rcpp-devel] Problems with Rcout<br>
Message-ID:<br>
<<a href="mailto:mailman.18867.1543334818.2604.rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">mailman.18867.1543334818.2604.rcpp-devel@lists.r-forge.r-project.org</a>><br>
Content-Type: text/plain; charset="utf-8"<br>
<br>
Dear Rcppers<br>
<br>
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:<br>
<br>
// [[Rcpp::export]]<br>
void testFunc() {<br>
Rcpp::Rcout << "testFunc begins" << std:endl;<br>
. . . <br>
}<br>
<br>
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.<br>
<br>
Thanks <br>
<br>
Barth<br>
<br>
<br>
-------------- next part --------------<br>
An HTML attachment was scrubbed...<br>
URL: <<a href="http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20181127/6adbedfe/attachment.html" rel="noreferrer noreferrer" target="_blank">http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20181127/6adbedfe/attachment.html</a>><br>
<br>
------------------------------<br>
<br>
Message: 3<br>
Date: Tue, 27 Nov 2018 16:50:13 +0100<br>
From: I?aki Ucar <<a href="mailto:iucar@fedoraproject.org" target="_blank" rel="noreferrer">iucar@fedoraproject.org</a>><br>
To: <a href="mailto:barthriley@comcast.net" target="_blank" rel="noreferrer">barthriley@comcast.net</a><br>
Cc: "<<a href="mailto:rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">rcpp-devel@lists.r-forge.r-project.org</a>>"<br>
<<a href="mailto:rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">rcpp-devel@lists.r-forge.r-project.org</a>><br>
Subject: Re: [Rcpp-devel] Problems with Rcout<br>
Message-ID:<br>
<<a href="mailto:CALEXWq3tTXhGPidXiswCLfvhB4-f2RDAmn12q%2BCGeMYaNs8VLg@mail.gmail.com" target="_blank" rel="noreferrer">CALEXWq3tTXhGPidXiswCLfvhB4-f2RDAmn12q+CGeMYaNs8VLg@mail.gmail.com</a>><br>
Content-Type: text/plain; charset="UTF-8"<br>
<br>
On Tue, 27 Nov 2018 at 16:35, Barth Riley <<a href="mailto:barthriley@comcast.net" target="_blank" rel="noreferrer">barthriley@comcast.net</a>> wrote:<br>
><br>
> Dear Rcppers<br>
><br>
> 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:<br>
><br>
> // [[Rcpp::export]]<br>
> void testFunc() {<br>
> Rcpp::Rcout << "testFunc begins" << std:endl;<br>
> . . .<br>
> }<br>
<br>
Note that it should be "std::endl", with double colon. Assuming this<br>
is a transcription error, you'll have to provide more context (and,<br>
ideally, some kind of reproducible example), because this works just<br>
fine.<br>
<br>
I?aki<br>
<br>
<br>
------------------------------<br>
<br>
Message: 4<br>
Date: Tue, 27 Nov 2018 16:51:28 +0100<br>
From: Serguei Sokol <<a href="mailto:serguei.sokol@gmail.com" target="_blank" rel="noreferrer">serguei.sokol@gmail.com</a>><br>
To: <a href="mailto:rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">rcpp-devel@lists.r-forge.r-project.org</a><br>
Subject: Re: [Rcpp-devel] Problems with Rcout<br>
Message-ID: <<a href="mailto:ccc55a85-79cb-c3c6-14e8-222223abe2c5@gmail.com" target="_blank" rel="noreferrer">ccc55a85-79cb-c3c6-14e8-222223abe2c5@gmail.com</a>><br>
Content-Type: text/plain; charset=utf-8; format=flowed<br>
<br>
Le 27/11/2018 ? 16:35, Barth Riley a ?crit?:<br>
> Dear Rcppers<br>
> <br>
> I am encountering a problem with Rcout. Even with basic string output, <br>
> when I run the function containing the call to Rcout, no output is <br>
> generated. Here is a simple example of what I?m trying to do:<br>
> <br>
> // [[Rcpp::export]]<br>
> <br>
> void testFunc() {<br>
> <br>
> ?????????? Rcpp::Rcout << "testFunc begins" << std:endl;<br>
> <br>
> ?????????? . . .<br>
> <br>
> }<br>
> <br>
This example works for me:<br>
<br>
library(Rcpp)<br>
sourceCpp(code="#include <Rcpp.h>\n// [[Rcpp::export]]\nvoid testFunc() <br>
{\nRcpp::Rcout << \"testFunc begins\" << std::endl;\n}")<br>
testFunc()<br>
#testFunc begins<br>
<br>
May be in your session you have redirected stdout elsewhere?<br>
<br>
Serguei.<br>
<br>
> My code is part of a package I?m developing, using R 3.51 and Rcpp <br>
> 0.12.19. The Rcpp code compiles without a problem.<br>
> <br>
> Thanks<br>
> <br>
> Barth<br>
> <br>
> <br>
> <br>
> _______________________________________________<br>
> Rcpp-devel mailing list<br>
> <a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">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 noreferrer" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
> <br>
<br>
<br>
<br>
------------------------------<br>
<br>
Message: 5<br>
Date: Tue, 27 Nov 2018 10:06:47 -0600<br>
From: Barth Riley <<a href="mailto:barthriley@comcast.net" target="_blank" rel="noreferrer">barthriley@comcast.net</a>><br>
To: I?aki Ucar <<a href="mailto:iucar@fedoraproject.org" target="_blank" rel="noreferrer">iucar@fedoraproject.org</a>>, "<br>
<a href="mailto:rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">rcpp-devel@lists.r-forge.r-project.org</a>"<br>
<<a href="mailto:rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">rcpp-devel@lists.r-forge.r-project.org</a>><br>
Subject: Re: [Rcpp-devel] Problems with Rcout<br>
Message-ID:<br>
<<a href="mailto:mailman.18868.1543334818.2604.rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">mailman.18868.1543334818.2604.rcpp-devel@lists.r-forge.r-project.org</a>><br>
Content-Type: text/plain; charset="utf-8"<br>
<br>
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.<br>
<br>
Barth<br>
<br>
NumericVector getValidCount(Rcpp::NumericMatrix m, <br>
bool treatAsVector) {<br>
<br>
Rcpp::Rcout << "getValidCount BEGINS" << std::endl;<br>
<br>
int N = m.cols();<br>
NumericVector u, vec;<br>
NumericVector count (N); <br>
<br>
if(!treatAsVector) {<br>
Rcpp::Rcout << "Treating as matrix" << std::endl;<br>
for(int i = 0; i < N; i++) {<br>
vec = m(_,i);<br>
vec = vec[!Rcpp::is_na(vec)];<br>
u = Rcpp::unique(vec);<br>
count[i] = u.length();<br>
} <br>
} else {<br>
Rcpp::Rcout << "treating as vector" << std::endl;<br>
vec = as<NumericVector>(m);<br>
vec = vec[!Rcpp::is_na(vec)];<br>
u = Rcpp::unique(vec);<br>
<br>
count.fill(u.length());<br>
}<br>
<br>
return count;<br>
}<br>
-------------- next part --------------<br>
An HTML attachment was scrubbed...<br>
URL: <<a href="http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20181127/607ae61f/attachment.html" rel="noreferrer noreferrer" target="_blank">http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20181127/607ae61f/attachment.html</a>><br>
<br>
------------------------------<br>
<br>
Subject: Digest Footer<br>
<br>
_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank" rel="noreferrer">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 noreferrer" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
<br>
------------------------------<br>
<br>
End of Rcpp-devel Digest, Vol 109, Issue 11<br>
*******************************************<br>
</blockquote></div></div></div>