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

[Søren Højsgaard]

Doug: Thanks a lot; I see the logic now! For my specific application the second argument is a sparse vector, but that modification should be straight forward with more iterators (and they are definitely sorted).

To the list: Does it have enough general interest to justify writing it up for the gallery?

Cheers
Søren





From: dmbates at gmail.com [mailto:dmbates at gmail.com] On Behalf Of Douglas Bates
Sent: 25. marts 2014 16:51
To: Søren Højsgaard
Cc: rcpp-devel at lists.r-forge.r-project.org
Subject: Re: [Rcpp-devel] RcppEigen: Avoiding accessing elements with coeff(i, j) in sparse matrix

Two other small points.  I am assuming that the set of indices in the second argument is sorted.  Also, it would be best to declare the function as

bool foo(const MSpMat X, const Eigen::Map<Eigen::VectorXi> idx) as in the enclosed

On Tue, Mar 25, 2014 at 10:45 AM, Douglas Bates <bates at stat.wisc.edu<mailto:bates at stat.wisc.edu>> wrote:
The enclosed works on this example.  The logic is to check each column in the index set for all the elements of the index set, except the one on the diagonal, being in the set of row indices. I have added some diagnostic output to demonstrate the flow.



On Tue, Mar 25, 2014 at 9:48 AM, Douglas Bates <bates at stat.wisc.edu<mailto:bates at stat.wisc.edu>> wrote:
On Tue, Mar 25, 2014 at 6:42 AM, Søren Højsgaard <sorenh at math.aau.dk<mailto:sorenh at math.aau.dk>> 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?

I was thinking about whether the sparse matrix iterators could be a solution, but I can't a hold on it.

Sparse matrix inner and outer iterators are definitely the way to go.  I'll write a version using them.

I was also thinking about using sparse matrices in RcppArmadillo, but I guess the problem (speed) is the same (haven't tried!).

Any advice will be greatly appreciated.
Cheers
Søren



_______________________________________________
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


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140325/ca0a0bc0/attachment.html>