Take advantage of the sparsity of b in AX=b

Question

There is a lot of info about how to use the sparsity pattern of A in order to solve $Ax = b$. However I can't find much about using the sparsity pattern of b. Let me take a concrete example:

Let us assume A is a large $N \times N$ matrix ($N > 10^6$) with a known sparsity pattern. Now let us assume I want to solve repeatedly linear systems of the form $Ax = b$ in which $b$ is a column vector with only one non-zero entry. (this seems to be a similar problem as computing some given columns of the matrix inverse $A^{-1}$).

I know I can either through direct or iterative solvers compute the solution of the linear system and at the same time take advantage of the known sparsity structure of A (which can be precomputed if needed), however what about $b$ ? Doing this feels like I am going to spend as much computing time as solving $Ax=c$ where c is a full vector, but here only one element of b is non-zero, is there a way to take advantage of that specifically ?

For example, I can think of a way if using a direct LU decomposition: if $b = e_N$ , the unit column vector with the last element equal to 1, for example (and we can always apply a permutation on A to make it that way I believe), then the forward substitution $L(Ux) = e_N$ consists indeed of only one operation. I assume that there might be other "tricks" like this out there that I am not aware of, hopefully included in some packages (I personnally like using petsc), so any link towards a reference on that topic would be greatly appreciated.

Answer 1

First let's rule out iterative methods like Krylov subspace methods. If you form products $Ax$, $A^2x$, etc. and they have a predictable sparsity pattern, then your matrix $A$ has some reducible property to it (it has some block that is decoupled from the rest, or something like that). Generally, these products will not maintain sparsity as the iteration progresses, and computing residuals like $b-Ax$ will destroy the sparsity of intermediates. No luck here.

For direct methods involving a matrix factorization, you are presumably already taking advantage of sparsity in the $O(N^3)$ factorization phase. As you mentioned, the sparsity pattern of $b$ can lead to an obvious reduction in the cost of the back-substitution (it is $O(kN)$ where $k$ is the number of nonzeros in $b$. This is already optimal, so it doesn't look like you can do better here either.

Now, as you pointed out, you are essentially just interested in a few selected columns of $A^{-1}$. The only algorithms I know of for doing this, like SelInv, require a matrix factorization, so it's no better than the direct method.

Answer 2

It is a fundamental intuition that the inverse of any "interesting" matrix is dense. This means that even if $b$ is sparse, $x = A^{-1} b$ will be dense. Thus in the best case (with a direct solver), you could get a speedup of 2x by using sparsity in the solve, but factorization is no faster.

With an iterative solver, the preconditioner either does its job to spread information globally as fast as possible (destroying sparsity immediately) or the iteration will converge very slowly with the vector still becoming dense before convergence. Even without considering implementation performance, this would provide very little benefit.

If you only need a few entries of the inverse, or perhaps the diagonal, algorithms like SelInv intelligently carry around the factors from a multifrontal solve instead of computing an explicit dense representation of the inverse. (This idea is old, but it has only recently been restated in common language.) There are also probing methods to estimate properties like the trace of the inverse.

Answer 3

Decomposing $A$ with $N>10^6$ can be even a problem for sparse direct solvers like UMFPACK. If you want to take an advantage out of the sparsity of $b$ you have to compute $LU=A$ and then compute the reachability set of the one non zero entry of $b$ over the graph defined by $L$. In this way you are able to solve $Ly=b$ and then compute the reachability set of the non zero elements of $y$ over the graph defined by $U$ and solve $Ux=y$ in this way. A easily understandable way computing these sets is described in the book "Direct Solution of Sparse Linear Systems" by Tim Davis.