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

[Douglas Bates]

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> wrote:

> On Tue, Mar 25, 2014 at 6:42 AM, Søren Højsgaard <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
>> 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/fa900e16/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: foo.cpp
Type: text/x-c++src
Size: 1586 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140325/fa900e16/attachment-0001.cpp>