[Rcpp-devel] Sparse matrix operations

[French, Joshua]

Soren and Doug,

Thanks for the info about the sparse matrices.  I'll give it a go and seen what happens.  If it works out, then perhaps I'll be able to provide a nice RcppGallery example.

Joshua
- -
Joshua French, Ph.D.
Assistant Professor
Department of Mathematical and Statistical Sciences
University of Colorado Denver
Joshua.French at ucdenver.edu
http://math.ucdenver.edu/~jfrench/
Ph:  303-556-6253  Fax:  303-556-8550

From: Søren Højsgaard <sorenh at math.aau.dk<mailto:sorenh at math.aau.dk>>
Date: Friday, January 31, 2014 at 4:22 PM
To: Douglas Bates <bates at stat.wisc.edu<mailto:bates at stat.wisc.edu>>, Joshua French <joshua.french at ucdenver.edu<mailto:joshua.french at ucdenver.edu>>
Cc: "rcpp-devel at lists.r-forge.r-project.org<mailto:rcpp-devel at lists.r-forge.r-project.org>" <rcpp-devel at r-forge.wu-wien.ac.at<mailto:rcpp-devel at r-forge.wu-wien.ac.at>>
Subject: RE: [Rcpp-devel] Sparse matrix operations

I’ll just add that I have used the sparse matrices in RcppEigen for representing graphs as adjacency matrices, as well as for implementing graph algorithms (with some help from Doug Bates; thanks Doug!). At the time when I implemented that I don’t think sparse matrices were available in RcppArmadillo, so I can’t tell which is easier.

A nice feature of the sparse matrices in RcppEigen is that there are iterators for those classes such that one can iterate over a sparse matrix efficiently by only visiting the “non-zero elements”; this is used in the gRbase package.

I find the documentation of Armadillo somewhat easier to read than that of Eigen (perhaps less abstract), but it is well worth spending the time digging into the Eigen documentation.

All the best
Søren






From: rcpp-devel-bounces at lists.r-forge.r-project.org<mailto:rcpp-devel-bounces at lists.r-forge.r-project.org> [mailto:rcpp-devel-bounces at lists.r-forge.r-project.org] On Behalf Of Douglas Bates
Sent: 31. januar 2014 21:45
To: French, Joshua
Cc: rcpp-devel at lists.r-forge.r-project.org<mailto:rcpp-devel at lists.r-forge.r-project.org>
Subject: Re: [Rcpp-devel] Sparse matrix operations

On Fri, Jan 31, 2014 at 1:14 PM, French, Joshua <JOSHUA.FRENCH at ucdenver.edu<mailto:JOSHUA.FRENCH at ucdenver.edu>> wrote:
Hello everyone,

For those of you who have used the sparse matrix capabilities of Armadillo/RcppArmadillo, how seamless are the operations in moving between sparse and dense matrices?  I know there are some tricks for getting sparseMatrix objects into SpMat objects (Dirk has a nice Rcpp Gallery post on this), but is it just as easy when it comes to doing Linear Algebra? My google ninja skills didn't seem to produce any helpful information.

E.g., if A is a dense matrix and B is a sparse matrix in Armadillo, it is pretty much seamless to do calculations like A + B, A*B, solve(B, A), etc.?

I'm not that familiar with the sparse matrix capabilities of Armadillo but I do know those of Eigen fairly well.  Operations on sparse matrices require more care and thought than do those on dense matrices.  It is possible to use "one size fits all" algorithms but they can be slow or inaccurate compared to more specialized approaches to operations like solve(B,A).  For dense matrices you could use an LU factorization or a QR factorization on a symmetric B and never notice the difference relative to using a Cholesky factorization.  With sparse matrices you could notice the difference.

Operations like A + B would not be terribly meaningful because the result would be dense if A is dense so you may as well extract a dense version of B and use that.  A * B is definitely available in Eigen and I imagine also available with Armadillo.  A general solve(B, A) could be written using a sparse LU but that may not be an effective way of solving such a system.  Most sparse matrix operations have a symbolic phase and a numeric phase.  In the symbolic phase the problem is analyzed to determine the positions of the nonzeros in the result and possibly auxiliary information such as a fill-reducing permutation.  Often if you need to solve several systems of the same form but with different numerical values you can save and update the factorization rather than starting from scratch each time.

This is why I prefer Eigen which has classes representing factorizations in addition to the one-shot methods.

I'm trying to decide whether it would be worth it to convert my code from R using sparseMatrix objects (Matrix package) to a C++ equivalent using RcppArmadillo.  On the same topic, would it be easier/more difficult if I used RcppEigen instead?

Thanks,

Joshua
--
Joshua French, Ph.D.
Assistant Professor
Department of Mathematical and Statistical Sciences
University of Colorado Denver
Joshua.French at ucdenver.edu<mailto:Joshua.French at ucdenver.edu>
http://math.ucdenver.edu/~jfrench/
Ph:  303-315-1709<tel:303-315-1709>  Fax:  303-315-1701<tel:303-315-1701>

_______________________________________________
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/20140203/eb9cb942/attachment.html>