[Rcpp-devel] RcppEigen: Avoiding accessing elements with coeff(i, j) in sparse matrix

Dirk Eddelbuettel edd at debian.org
Tue Mar 25 13:16:25 CET 2014


On 25 March 2014 at 11:42, Søren Højsgaard wrote:
| Dear all, 
| I have a large sparse adjacency matrix X for an undirected graph. For a subset 'idx' of the vertices I want to find out if the subgraph defined by this subset is complete (i.e. has edges between all variables). So in R, one would check if X[idx,idx] has only ones outside the diagonal, but this is slow and I want to do it with Rcpp. Here is one attempt where I simply add up the elements above (or below) the diagonal, and to access the elements of X I use coeff which is relatively slow (because X is a sparse matrix).
| 
| #include <RcppEigen.h>
| typedef Eigen::MappedSparseMatrix<double> MSpMat;
| typedef Eigen::SparseVector<double> SpVec;
| typedef SpVec::InnerIterator InIter;
| 
| //[[Rcpp::export]]
| double foo (MSpMat X, Eigen::VectorXi idx){
|   SpVec sidx = idx.sparseView();
|   double out=0;
|   int i2, j2;
|   for (InIter i_(sidx); i_; ++i_){
|     i2 = i_.index(); 
|     for (InIter j_(sidx); j_; ++j_){
|       j2 = j_.index(); 
|       if (i2>j2)
| 	out += X.coeff( i2, j2);
|     }
|   }
|   return out;
| }
| 
| /*** R
| library(Matrix)
| M1 <- new("dgCMatrix"
|     , i = c(1L, 2L, 3L, 0L, 2L, 3L, 0L, 1L, 3L, 0L, 1L, 2L, 4L, 5L, 3L, 5L, 3L, 4L)
|     , p = c(0L, 3L, 6L, 9L, 14L, 16L, 18L)
|     , Dim = c(6L, 6L)
|     , Dimnames = list(c("a", "b", "c", "d", "e", "f"), c("a", "b", "c", "d", "e", "f"))
|     , x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
|     , factors = list()
| )
| M1 
| 
| foo(M1, c(2,3,4))
| 
| */
| Can anyone see a better way of doing this? 

Not really from the top of my head. You do need to access all those elements,
and the traversal in C++ ought to / will be faster than doing it in R.
 
| I was thinking about whether the sparse matrix iterators could be a solution, but I can't a hold on it. I was also thinking about using sparse matrices in RcppArmadillo, but I guess the problem (speed) is the same (haven't tried!).

Sparse matrices are complicated. I heard from Conrad and Ryan (who, over at
mlpack.org is driving some sparse matrix adoption for Armadillo) that they
put in a bunch of work to get a sparse solver going (via an external library)
only to find out that for common problem sizes it was slower than the dense
solver :-/

So while Rcpp may give you some tools, you may still have to do some
algorithmic work.

Dirk
 
| Any advice will be greatly appreciated.
| Cheers
| Søren
| 
| 
| 
| _______________________________________________
| 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

-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com


More information about the Rcpp-devel mailing list