18
$\begingroup$

I'm trying to solve a 2D Poisson equation by finite differences. In the process, I obtain a sparse matrix with only $5$ variables in each equation. For example, if the variables were $U$, then the discretization would yield:

$$U_{i-1,j} + U_{i+1,j} -4U_{i,j} + U_{i,j-1} + U_{i,j+1} = f_{i,j}$$

I know that I can solve this system by an iterative method, but the thought occurred to me that if I ordered the variables appropriately, I might be able to obtain a banded matrix which could be solved by a direct method (i.e., Gaussian elimination w/o pivoting). Is this possible? Are there any strategies for doing this for other, perhaps less structured sparse systems?

Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ Something like Cuthill-McKee, then? $\endgroup$
    – J. M.
    Commented Jan 18, 2012 at 17:46
  • $\begingroup$ Interesting... i've never heard of the Cuthill-McKee algorithm before! :) $\endgroup$
    – Paul
    Commented Jan 18, 2012 at 18:18
  • 1
    $\begingroup$ There's also a Reverse Cuthill-McKee as well. $\endgroup$
    – Geoff Oxberry
    Commented Jan 19, 2012 at 4:38
  • 1
    $\begingroup$ I hope it is clear from the answers, but you do not want to use a banded solver for this problem, nor choose an ordering that minimizes bandwidth. Perhaps the question or the chosen answer can be edited to make this clear, otherwise I fear that this myth will be perpetuated. I gave a visual comparison and compared fill in scicomp.stackexchange.com/a/880/119. $\endgroup$
    – Jed Brown
    Commented Jan 19, 2012 at 13:40
  • $\begingroup$ @JedBrown: Actually, I'm not quite working with a poisson problem, per se... My problem has a similar structure to the poisson problem... The indicies of the variables (i's and j's) are exactly the same, and the matrix is diagonally dominant with the off-diagonal entries (within the same row) add to exactly the sum of the diagonal entry. $\endgroup$
    – Paul
    Commented Jan 19, 2012 at 14:09

6 Answers 6

Reset to default
13
$\begingroup$

This is a well-studied problem in the field of sparse-direct solvers. I highly recommend reading Joseph Liu's overview of the multifrontal method in order to get a better idea of how reorderings and supernodes effect fill-in and solution time.

Nested dissection is an extremely common way to generate the reordering, and essentially consists of recursive graph partitioning. MeTiS is the de facto standard for graph partitioning, and you can read about some of the ideas behind it here. Another commonly used package is SCOTCH, and Chaco is also important, as its authors introduced multi-level graph partitioning, which is also the fundamental idea behind MeTiS.

George and Liu showed in their classic book that 2d sparse-direct solutions only require $O(n^{3/2})$ work and $O(n \log n)$ memory, while 3d sparse-direct requires $O(n^2)$ work and $O(n^{4/3})$ memory.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Do you have a citation for the George and Liu reference? $\endgroup$
    – Paul
    Commented Jan 18, 2012 at 18:49
  • $\begingroup$ Added; I was about to get out of the car when I first submitted it. I know that there exists a freely available version of the book online somewhere (Jed knows where it is), but I could not find it. $\endgroup$
    – Jack Poulson
    Commented Jan 18, 2012 at 19:04
  • $\begingroup$ I updated the link to point to the PDF of the book instead of the book review. $\endgroup$
    – Jed Brown
    Commented Jan 18, 2012 at 20:56
  • $\begingroup$ @JedBrown That was a great reference! Thanks so much! :) $\endgroup$
    – Paul
    Commented Jan 19, 2012 at 5:34
  • 1
    $\begingroup$ @Alexander Everyone attributes the 3D bound to George and Liu, though I don't know if they explicitly point it out in the book. It is obvious from the theory however. The minimal vertex separator for a $n = m\times m\times m$ grid is $n^{2/3} = m\times m$. The dense matrix associated with that supernode has $(n^{2/3})^2 = n^{4/3}$ entries and requires $(n^{2/3})^3 = n^2$ operations to factor. The logarithmic term in the 2D case is more subtle and is treated in Chapter 8 on Nested Dissection, which achieves the lower bound. $\endgroup$
    – Jed Brown
    Commented Aug 27, 2012 at 13:59
6
$\begingroup$

Cuthill-McKee is the de facto standard for what you want to do. If you wanted to play with this method, there's an easy-to-use implementation of the algorithm (and its reverse) in the Boost Graph Library (BGL), and the documentation contains examples how to use it.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ actually reverse Cuhill-McKee; it usually gives less fill. But a nested dissection ordering is far superior to a low bandwidth ordering. $\endgroup$
    – Arnold Neumaier
    Commented Mar 29, 2012 at 19:46
4
$\begingroup$

Speaking of multifrontal methods, Tim Davis, who works on multifrontal methods for LU factorization (UMFPACK) has a number of routines that will reorder matrices to minimize fill-in. You can find them as here as part of SuiteSparse. SuiteSparse uses MeTiS.

One other thing to note: In some problems, you can be clever about ordering variables so that you get banded, or close to banded, patterns, which can save you the trouble (and the CPU time) of calling these algorithms. However, this clever reordering requires insight on your part and is nowhere near as general as the graph-theory-based reordering algorithms people have mentioned in their answers here.

Improve this answer
$\endgroup$
1
  • $\begingroup$ You're welcome, Paul. If you like it, vote it up. $\endgroup$
    – Geoff Oxberry
    Commented Jan 19, 2012 at 5:31
3
$\begingroup$

There's an algorithm called ADI (Alternating Direction Implicit) in applied math circles and Split-operator in physics circles that does basically what you describe. It's an iterative method, and it follows this basic procedure:

  1. For every value of $y$ , relax in the $x$-direction. This matrix should be tridiagonal, so it can be solved directly in relatively little time.

  2. For every value of $x$ , relax in the $y$-direction. Again, this should be pretty quick.

  3. Repeat 1 and 2 until the error is as small as you want it to be.

I don't know the formal complexity of this algorithm, but I've found it to converge in fewer iterations than things like Jacobi and Gauss-Seidel every time I've used it.

Improve this answer
$\endgroup$
4
  • 2
    $\begingroup$ If you decide to go the operator splitting route, something you'll want to be careful about is that operator splitting is known to result in errors in steady state solutions in some cases. (One of my labmates has developed a way to overcome this difficulty, but I don't believe he's published it yet.) Also, operator splitting is known to cause numerical errors. There are well-established ways to estimate these errors a posteriori; Don Estep has done excellent work in that area. $\endgroup$
    – Geoff Oxberry
    Commented Jan 19, 2012 at 4:37
  • $\begingroup$ @GeoffOxberry It sounds like you are referring to a different splitting. You can use ADI in a fully implicit scheme which has no splitting error because it actually solves the system. There are also IMEX methods that rigorously control splitting errors. $\endgroup$
    – Jed Brown
    Commented Jan 19, 2012 at 5:14
  • $\begingroup$ @JedBrown: I was talking about Godunov and Strang splitting, which can yield similar tridiagonal matrices if you split up the $x$ and $y$ terms. My bad. (I keep telling people I want to learn PDE methods in a postdoc...) $\endgroup$
    – Geoff Oxberry
    Commented Jan 19, 2012 at 5:24
  • $\begingroup$ I've never heard of Godunov and Strang splitting. I tend to split my operator with the Baker-Campbell-Hausdorf formula. Is that the same thing? $\endgroup$
    – Dan
    Commented Jan 19, 2012 at 6:20
2
$\begingroup$

There are many matrix reordering algorithms, and they seek a permutation array that maps from the original coordinates to the reordered coordinates.

Behrisch et al. (2016) gave a nice review of different reordering algorithms, with an analysis of their usage, complexity, objective and applications. The objective is not always to minimize matrix bandwidth, although it lists many that do, including the Reverse Cuthill McKee algorithm. I recommend that to minimize bandwidth, and have used the Python implementation in Scipy.

Improve this answer
$\endgroup$
1
$\begingroup$

While all the present answers are valid solutions to the practical problem, technically the answer to the question in your title (how to reorder variable to minimize bandwidth) is "it's an NP-complete problem". Article: https://link.springer.com/article/10.1007/BF02280884 .

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.