[Rcpp-devel] RcppEigen: Avoiding accessing elements with coeff(i, j) in sparse matrix
Douglas Bates
bates at stat.wisc.edu
Tue Mar 25 15:48:54 CET 2014
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/7885576f/attachment.html>
More information about the Rcpp-devel
mailing list