Library that performs sparse matrix-vector and matrix-transpose-vector multiplication

Question

My main interest is sparse matrix-vector and matrix-transpose-vector multiplication of the form y=y+AA'x. Is there any library that performs y=y+AA'x efficiently? I have investigated SPARSKIT, MKL, ARPACK,... but I was not able to find the related routine. Isn't the y=y+AA'x operation widely used?

Answer 1

The matrix product $B = AA^T$ is generally faster to apply as $A (A^T x)$ even when the product matrix $B$ is already available. Only for peculiar graphs with extremely low expansion factor would the product matrix be faster.

PDE graphs in two or more dimensions have higher expansion factors. For example, if $A$ is a 9-point differencing operation on a 2D rectangular grid, then $A^2$ (or $A A^T$) is a 25-point operation; applying the 9-point operation twice is faster. In 3D, the analog is a 27-point $A$ and 125-point $A^2$.

Note the worst case: if $A$ contains a single dense column, then $A A^T$ is dense.

All sparse linear algebra libraries have support for applying $A$ and $A^T$.

Answer 2

Usually you would do one of two things:

1. Compute $z = Ax$, then $y = y+Az$;
2. If you anticipate needing to do this operation many times, you might compute the matrix $B = AA'$ and then evaluate $y = y+Bx$, but it will be much less sparse than the original matrix $A$.

I don't know of any libraries that will compute $AA'x$ directly, and if there are any that do so, they probably use (1) behind the scenes.

SPARSKIT has routines for both $A\cdot x$ and $A'\cdot x$, they're called amux and atmux respectively; these are for matrices in the compressed sparse row format, which is described in the SPARSKIT documentation.

Do you have a language preference? If you're just experimenting and don't need your program to be super fast, you could try the Scipy sparse matrix modules in python, which would be less of a hassle than using Fortran. These will also have the functionality you're looking for.

EDIT: Since you want fast and you mentioned GPU computing, you might want to consider the CUSP library. The ellpack matrix format is particularly amenable to usage on GPU machines, although there are exceptional cases.

Answer 3

If a C++ library is acceptable, I suggest trying Eigen.

http://eigen.tuxfamily.org/index.php?title=Main_Page

The performance is good and it is easy to use. In particular, your expression above can be written like this:

Eigen::SparseMatrix<double> a(n,n);
Eigen::VectorXd y(n), x(n);
y += a*a.transpose()*x;

Answer 4

If you want good parallel performance in shared memory, you can try out the Compressed Sparse Blocks (CSB) code: http://gauss.cs.ucsb.edu/~aydin/csb/csb.html

It is specifically designed to give equally good parallel performance for SpMV and SpMV_T (the transposed matrix case) without explicitly transposing the matrix. It also has provable parallelism guarantees regardless of the nonzero distribution (it can be as skewed as it gets)