[Rcpp-devel] Speed of RCppEigen Cholesky decomposition on sparse matrix (Serguei Sokol)

Dmitriy Selivanov

2018-11-27 16:57:07 UTC

But adding 0 to a sparse matrix is expensive operation. It doesn't look
fair to include it to benchmark.

Send Rcpp-devel mailing list submissions to
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
You can reach the person managing the list at
When replying, please edit your Subject line so it is more specific
than "Re: Contents of Rcpp-devel digest..."
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
Subject: Re: [Rcpp-devel] Speed of RCppEigen Cholesky decomposition on
sparse matrix
Content-Type: text/plain; charset=utf-8; format=flowed

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
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
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.

(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 )))
#?? user??system elapsed
#??0.341?? 0.014?? 0.355
system.time(replicate(10, CholSparse( C )))
#?? 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

/Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
/Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
[1] stats???? graphics??grDevices datasets??utils???? methods?? base
[1] inline_0.3.15 Matrix_1.2-15
[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
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel

------------------------------
Message: 2
Date: Tue, 27 Nov 2018 09:35:10 -0600
Subject: [Rcpp-devel] Problems with Rcout
<
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
// [[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
Subject: Re: [Rcpp-devel] Problems with Rcout
<
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

// [[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
Subject: Re: [Rcpp-devel] Problems with Rcout
Content-Type: text/plain; charset=utf-8; format=flowed

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
// [[Rcpp::export]]
void testFunc() {
?????????? Rcpp::Rcout << "testFunc begins" << std:endl;
?????????? . . .
}

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
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel

------------------------------
Message: 5
Date: Tue, 27 Nov 2018 10:06:47 -0600
Subject: Re: [Rcpp-devel] Problems with Rcout
<
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
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
------------------------------
End of Rcpp-devel Digest, Vol 109, Issue 11
*******************************************

Serguei Sokol

2018-11-28 10:05:10 UTC

Permalink

Post by Dmitriy Selivanov
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.

Post by Dmitriy Selivanov
вт, 27 нояб. 2018 г., 20:07
Send Rcpp-devel mailing list submissions to
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
You can reach the person managing the list at
When replying, please edit your Subject line so it is more specific
than "Re: Contents of Rcpp-devel digest..."
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
Subject: Re: [Rcpp-devel] Speed of RCppEigen Cholesky decomposition on
sparse matrix
Content-Type: text/plain; charset=utf-8; format=flowed

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?

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 )))
#?? user??system elapsed
#??0.341?? 0.014?? 0.355
system.time(replicate(10, CholSparse( C )))
#?? 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

/Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
/Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
------------------------------
Message: 2
Date: Tue, 27 Nov 2018 09:35:10 -0600
Subject: [Rcpp-devel] Problems with Rcout
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
// [[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...
<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
Subject: Re: [Rcpp-devel] Problems with Rcout
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

// [[Rcpp::export]]
void testFunc() {
Rcpp::Rcout << "testFunc begins" << std:endl;
. . .
}

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
// [[Rcpp::export]]
void testFunc() {
?????????? Rcpp::Rcout << "testFunc begins" << std:endl;
?????????? . . .
}

https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
------------------------------
Message: 5
Date: Tue, 27 Nov 2018 10:06:47 -0600
Subject: Re: [Rcpp-devel] Problems with Rcout
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...
<http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20181127/607ae61f/attachment.html>
------------------------------
Subject: Digest Footer
_______________________________________________
Rcpp-devel mailing list
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
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel