[Rcpp-devel] Sparse matrix operations

[Douglas Bates]

On Fri, Jan 31, 2014 at 1:14 PM, French, Joshua
<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
> http://math.ucdenver.edu/~jfrench/
> Ph:  303-315-1709  Fax:  303-315-1701
>
> _______________________________________________
> 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/20140131/a410704d/attachment.html>