Solve banded linear system with large bandwidth but sparse interior band structure

linear algebramatrix decompositionpartial differential equationssparse matrices

Assume the linear system $Ax = b$, where $A$ is a $N \times N$ banded matrix with lower and upper bandwidth $l$, and $N >> l >> 1$. $A$ has the following structure:

All entries of $A$ are zero, except for the following: the diagonal of $A$, the first super-diagnonal band and the first sub-diagonal band are non-zero; the $(l-2)$nd, $(l-1)$th and the $l$-th super-diagonal and sub-diagonal bands are also nonzero. Therefore, $A$ contains 9 non-zero bands.

The system can be solved using a standard band-matrix $LU$ decomposition, the complexity of the solve then scales as $N\times l^2$ for $l << N$. However, this method does not exploit the fact that all bands between the first off-diagonals and the $(l-2)$nd off-diagonals are identically zero.

Such a system arises in the discretization of second-order, two-dimensional differential equations with mixed derivatives, using a three-point compact stencil at interior nodes and a two-point stencil at boundary nodes.

Besides bandwidth-reducing algorithms, does a faster method exist to solve such a system, that makes use of the sparse band structure described here?

Best Answer

Note that standard banded $LU$ factorization for a matrix of bandwidth $\ell$ has time complexity $O(N \ell^2)$.${}^*$ In your question, you specifically address discretized differential equations.

There has been considerable work on iterative methods for such discretizations. It is a truly vast field, but very popular is a Krylov subspace method preconditioned by, say, algebraic multigrid. For many PDEs, these can provably produce an approximate solution to a specified tolerance in $O(N)$ time.

Let me now address the question you asked. Given an arbitrary banded matrix with bandwidth $\ell$ with only $9$ nonzero diagonals as specified in your answer, does there exist a method faster than banded $LU$ factorization? If the answer to this question is no, then it would be very hard to show this, as you must demonstrate that every method which solves the problem must have a certain time complexity. There is some very sophisticated and interesting work in theoretical computer science which is addressing this problem, which I must admit I don't fully understand which might be able to give very precise "lower bound" statements on how fast problems of this type might be solved.

Here is a specific result which you may or may not be aware of which might temper your expectations about what time complexity you can hope to achieve. Informally, the result might be stated as follows:

For any matrix $A$ fitting your description with $\ell \approx \sqrt{N}$, there is no reordering of the rows and columns of the matrix $A$ such that an $LU$ factorization of the reordered $A$ can be computed in time faster than $O(N^{3/2})$. Moreover, this time complexity is obtained by reordering $A$ according to the nested dissection ordering.

This result is more or less a consequence of the results of this paper. The case $\ell \approx \sqrt{N}$ is important because this corresponds to the discretization of a 2D PDE on a square mesh. The result does not rule out some other algorithm solving the problem in faster than $O(N^{3/2})$ time, but it does rule out a class of algorithms which are "$LU$ factorization-like". My hunch is that $N^{3/2}$ is as good as you can get for this problem without some additional assumption (like $A$ being diagonally dominant, which is often the case for discretized PDEs).

I would look into numerical methods for linear systems arising specifically from discretized PDEs, for which there is additional structure (like perhaps diagonal dominance) which can open up certain algorithmic approaches that won't work on a fully general matrix $A$ with the sparsity pattern you describe.


${}^*$ To see why this must be the case, note that if $A$ is dense and thus $\ell = N-1$, then a complexity of $O(N\ell)$ would mean we could solve every system of linear equations in $O(N^2)$ operations, which has not been shown to be possible. $O(N\ell^2)$ recovers the standard $O(N^3)$ for Gaussian elimination.

Related Question