5
$\begingroup$

Consider a square $N \times N$ block-matrix $\mathbf{A}$, where each $n \times n$ block $\mathbf{A}_{ii}$ is either a dense block or a zero-block. So, $N$ denotes the number of blocks, $n$ denotes the size of each square block.

Random example for illustration: $$ \mathbf{A}=\left[ \begin{array}{ccccc} \mathbf{A}_{11} & \mathbf{0} & \mathbf{A}_{13} & \mathbf{A}_{14} & \mathbf{0} \\ \mathbf{A}_{21} & \mathbf{A}_{22} & \mathbf{0} & \mathbf{0} & \mathbf{A}_{25} \\ \mathbf{0} & \mathbf{0} & \mathbf{A}_{33} & \mathbf{A}_{34} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{A}_{43} & \mathbf{A}_{44} & \mathbf{0} \\ \mathbf{A}_{51} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{A}_{55} \\ \end{array} \right] $$

Now, one wants to compute an LU decomposition of this block matrix with minimal fill-in.

I am generally aware of the algorithms for reducing the fill-in for general sparse matrices. Say, Matlab reference gives an overview of several methods and their effect on Cholesky factorization. However, these algorithms certainly do not give an optimal reordering to give an optimal structure for the factorization in preserving zeros. This is done for a very good reason: for a large sparse matrix, finding an optimal reordering is very expensive, unneeded, and is unlikely to give substantial benefits compared to the common relatively fast techniques.

However, I wonder if that can be slightly different for the factorization of a block-matrix. In this case, $N$ is an order of magnitudes smaller, and one might be willing to use very expensive [graph?] algorithms to find this optimal fill-in for a given block pattern.

I wonder:

  • what family of algorithms should I look at?
  • is there a canonical reference on this aspect of reducing fill-in? (my search is constantly swamped by results on reducing fill-in of simple sparse matrices not the block-wise case with its own peculiarities – as much more expensive algorithm can become feasible)?
  • is it even a tangible problem for relatively small $N$?

Notes:

  • In this example, all blocks have the same $n \times n$ size for simplicity, and I am interested in it. I might be interested in the example where each block is still square, but the block-sizes differ.
  • For my applications, the number of blocks $N$ can range from 1 (trivial) to, say, 100, and $n$ is on the order of thousands (not that $n$ matters at all).
  • I am willing to spend significant precomputation time (say 10$\times t_\text{LU}$) to analyze the pattern and come up with the best reordering
$\endgroup$
  • $\begingroup$ You say: "reducing fill-in of simple sparse matrices not the block-wise case with its own peculiarities". It seems to me that as far as determining a fill-reducing ordering, it is irrelevant whether the block size is 100x100 or 1x1 (e.g. scalar). Why do you see the problem differently? $\endgroup$ – Bill Greene Dec 11 '19 at 0:37
  • $\begingroup$ @BillGreene you are totally right. $n$, block-size is totally irrelevant (especially when all blocks have the same size $n \times n$) $\endgroup$ – Anton Menshov Dec 11 '19 at 0:38
  • $\begingroup$ I wanted to highlight that $n$ is usually in the order of tens of thousands when the reordering algorithms for reducing fill in are applied (and $N=1$). I am interested in the case when $N$ (num blocks, the dimension that matters) is relatively small, but larger than 1 -> block matrix. $\endgroup$ – Anton Menshov Dec 11 '19 at 0:40
5
$\begingroup$

Unfortunately, optimal reordering is NP-complete, so all the algorithms in this space are kinda heuristic. At first glance my instinct would be to try some variant of minimum degree, mainly because the bisection-based algorithms (nested-dissection / recursive-spectral-bisection) would seem not really applicable.. their heuristics are kinda rooted in the geometry/sparse connectivity of a 2D/3D mesh.

The particular example is interesting in that the pattern is not symmetric, and most of the reordering codes in wide use are meant for symmetric matrices or nearly-symmetric perturbations. I think KLU is a noteworthy exception, a sparse solver that is not so tailored for symmetry/FEM (it's for circuits). Perhaps take a look at how it reorders matrices?

For extremely small numbers of blocks (less than 7 or so?) there's always the possibility of trying every permutation, performing symbolic factorization on the pattern, then just picking the best? I don't think I would recommend this since it won't scale, but if you have known upper bounds on N, it could be an easy thing to try.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ that's a very interesting algorithm I was not aware of, and this documentation is a very interesting read on every level! $\endgroup$ – Anton Menshov Dec 15 '19 at 22:00
1
$\begingroup$

I don't know if this will be helpful to you (since it is a late answer and you already accepted the other answer) but there is a decent literature review of bandwidth minimizing algorithms in https://hal.archives-ouvertes.fr/hal-01166658/document .

I would think the matrix as a graph and collapse each block into a single node, i.e. if two nodes of the directed graph -corresponding to the non-symmetric sparse matrix- are connected via a bidirected edge then collapse them. That would end up reducing the graph size significantly enough that most bandwidth minimizing algorithms would be feasible. I probably would use Reverse Cuthill-McKee ( MATLAB implementation https://www.mathworks.com/help/matlab/ref/symrcm.html ) but that is mostly due to its availability.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for the overview, I haven't seen this one. However, these types of algorithms do not provide optimal reordering. Which is certainly a good decision for a large sparse matrix, but maybe not for a block-matrix question, where more expensive and less scalable algorithms can be used. $\endgroup$ – Anton Menshov Mar 11 at 14:54
  • $\begingroup$ They cite Caprara A. and J. J. Salazar, 2005 as an optimal labeling algorithm for medium size graphs. I don't know if there is a public implementation of their 'branch and bound' method. $\endgroup$ – Abdullah Ali Sivas Mar 12 at 15:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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