4

Minimizing Tile Re-ordering Problem:

Suppose I had the following symmetric 9x9 matrix, N^2 interactions between N particles:

(1,2) (2,9) (4,5) (4,6) (5,8) (7,8),

These are symmetric interactions, so it implicitly implies that there exists:

(2,1) (9,2) (5,4) (6,4) (8,5) (8,7),

In my problem, suppose they are arranged in matrix form, where only the upper triangle is shown:

t       0     1     2     (tiles)
  #   1 2 3 4 5 6 7 8 9   
  1 [ 0 1 0 0 0 0 0 0 0 ]
0 2 [ x 0 0 0 0 0 0 0 1 ]
  3 [ x x 0 0 0 0 0 0 0 ]
  4 [ x x x 0 1 1 0 0 0 ] 
1 5 [ x x x x 0 0 0 1 0 ]
  6 [ x x x x x 0 0 0 0 ]
  7 [ x x x x x x 0 1 0 ]
2 8 [ x x x x x x x 0 0 ]
  9 [ x x x x x x x x 0 ] (x's denote symmetric pair)

I have some operation that's computed in 3x3 tiles, and any 3x3 that contains at least a single 1 must be computed entirely. The above example requires at least 5 tiles: (0,0), (0,2), (1,1), (1,2), (2,2)

However, if I swap the 3rd and 9th columns (and along with the rows since its a symmetric matrix) by permutating my input:

t       0     1     2
  #   1 2 9 4 5 6 7 8 3 
  1 [ 0 1 0 0 0 0 0 0 0 ]
0 2 [ x 0 1 0 0 0 0 0 0 ]
  9 [ x x 0 0 0 0 0 0 0 ]
  4 [ x x x 0 1 1 0 0 0 ] 
1 5 [ x x x x 0 0 0 1 0 ]
  6 [ x x x x x 0 0 0 0 ]
  7 [ x x x x x x 0 1 0 ]
2 8 [ x x x x x x x 0 0 ]
  3 [ x x x x x x x x 0 ] (x's denote symmetric pair)

Now I only need to compute 4 tiles: (0,0), (1,1), (1,2), (2,2).

The General Problem:

Given an NxN sparse matrix, finding an re-ordering to minimize the number of TxT tiles that must be computed. Suppose that N is a multiple of T. An optimal, but unfeasible, solution can be found by trying out the N! permutations of the input ordering.

For heuristics, I've tried bandwidth minimization routines (such as Reverse CutHill McKee), Tim Davis' AMD routines, so far to no avail. I don't think diagonalization is the right approach here.

Here's a sample starting matrix:

http://proteneer.com/misc/out2.dat

Hilbert Curve: Hilbert Curve

RCM: RCM

Morton Curve: Morton Curve

proteneer
  • 572
  • 5
  • 17
  • Just a first thought (haven't had my first coffee yet, so don't expect much) - if a greedy reordering algorithm isn't optimal, simulated annealing often gets good results. That is, for each possible "move" (of a column/row), how much improvement do you get? – Rafe Sep 28 '12 at 00:07
  • In general, it is very difficult (and nearly impossible) for a single column/row swap to improve the score, where the score is the # of TxT tiles to be computed. In realistic situations, most tiles already have at least two 1's. A single swap likely will never be able to remove a tile. A better scoring function would certainly help. – proteneer Sep 28 '12 at 18:02
  • When you say "However, if I swap the 3rd and 9th columns (and along with the rows since its a symmetric matrix) by permutating my input." Do you mean that any given row can be swapped with any other row? Does it matter if a column in either row to be swapped has already been swapped? Also, In your example when column 3 was swapped with column 9 - the x values do not appear to have been swapped, and zeros were then specified for the x values that were in column 3 and now are in column 9. So, if column 8 were to be swapped with column 1, would that be allowed? – Bob Bryan Sep 28 '12 at 20:43
  • @BobBryan A better way to illustrate the above is to consider a neighbour list of tuples: (1,2) (2,9) (4,5) (4,6) (5,8) (7,8), these are symmetric, so it implicitly implies that there exists a (2,1) (9,2) (5,4), .. etc. Both matrices represent this neighourlist. – proteneer Sep 28 '12 at 21:48
  • @proteneer: Can you be more specific when you say morton curve? Is this the z order morton curve or just the 4 copies of the hilbert curve stichting together to a morton curve name after the american mathemacien? – Micromega Sep 28 '12 at 22:12
  • @Chiyou the Z-Order curve generated by interleaving 3 sets of bits – proteneer Sep 28 '12 at 22:15
  • @proteneer: Thanks. Do you think the real morton curve can help you? There is some difference to the original hilbert curve. The morton curve makes a loop (when you reorder the points)? – Micromega Sep 28 '12 at 22:18
  • @proteneer: Yes, the moore curve. 3d is difficult but when you know how it isn't. It's about the grey code. – Micromega Sep 28 '12 at 22:29
  • Did you arrive at an acceptable solution? I would suggest trying a octtree approach for the fun of it. I do not have a 3D implementation, but you can try if it works in 2D using [a matlab mex function](http://sourceforge.net/projects/milamin/files/). You need to download the mutils-0.2 package. Let me know if you need any help with it. – angainor Oct 08 '12 at 11:04
  • Unfortunately no, so far nothing still beats the Hilbert curve - I'm going to turn to graph theory for some additional insights. The problem is actually reducible to a variant of a k-way partitioning problem. – proteneer Oct 14 '12 at 06:40

2 Answers2

3

There are several well-known options you can try (some of them you have, but still):

  • (Reverse) Cuthill-McKee reduced the matrix bandwidth, keeping the entries close to the diagonal.
  • Approximage Minimum Degree - a light-weight fill-reducing reordering.
  • fill-reducing reordering for sparse LU/LL' decomposition (METIS, SCOTCH) - quite computationally heavy.
  • space filling curve reordering (something in these lines)
  • quad-trees for 2D or oct-trees for 3D problems - you assign the particles to quads/octants and later number them according to the quad/octant id, similar to space filling curves in a sense.
  • Self Avoiding Walk is used on structured grids to traverse the grid points in such order that all points are only visited once
  • a lot of research in blocking of the sparse matrix entries has been done in the context of Sparse Matrix-Vector multiplication. Many of the researchers have tried to find good reordering for that purpose (I do not have the perfect overview on that subject, but have a look at e.g. this paper)

All of those tend to find structure in your matrix and in some sense group the non-zero entries. Since you say you deal with particles, it means that your connectivity graph is in some sense 'local' because of spatial locality of the particle interactions. In this case these methods should be of good use.

Of course, they do not provide the exact solution to the problem :) But they are commonly used in exactly such cases because they yield very good reorderings in practice. I wonder what do you mean by saying the methods you tried failed? Do you expect to find the optimum solution? Surely, they improve the situation compared to a random matrix ordering.

Edit Let me briefly go through a few pictures. I have created a 3D structured cartesian mesh composed of 20-node brick elements. I matched the size of the mesh so that it is similar to yours (~1000 nodes). Also, number of non-zero entries per row are not too far off (51-81 in my case, 59-81 in your case, both however have very different distributions) The pictures below show RCM and METIS reorderings for non-periodic mesh (left), and for mesh with complete x-y-z periodicity (right):RCM reordering

Next picture shows the same matrix reordered using METIS and fill-reducing reordering

metis reordering

The difference is striking - bad impact of periodicity is clear. Now your matrix reordered with RCM and METIS

enter image description here

WOW. You have a problem :) First of all, I think there is something wrong with your rcm, because mine looks different ;) Also, I am certain that you can not conclude anything general and meaningful about any reordering based on this particular matrix. This is because your system size is very small (less than roughly 10x10x10 points), and you seem to have relatively long-range interactions between your particles. Hence, introducing periodicity into such small system has a much stronger bad effect on reordering than is seen in my structured case.

I would start the search for a good reordering by turning off periodicity. Once you have a reordering that satisfies you, introduce periodic interactions. In the system you showed there is almost nothing but periodicity: because it is very smal and because your interactions are fairly long-range, at least compared to my mesh. In much larger systems periodicity will have a smaller effect on the center of the model.

Smaller, but still negative. Maybe you could change your approach to periodicity? Instead of including periodic connectivities explicitly in the matrix, construct and reorder a matrix without those and introduce explicit equations binding the periodic particles together, e.g.:

V_particle1 = V_particle100

or in other words

V_particle1 - V_particle100 = 0

and add those equations at the end of your matrix. This method is called the Lagrange multipliers. Here is how it looks for my system

enter image description here

You keep the reordering of the non-periodic system and the periodic connectivities are localized in a block at the end of the matrix. Of course, you can use it for any other reorderings.

The next idea is you start with a reordered non-periodic system and explicitly eliminate matrix rows for the periodic nodes by adding them into the rows they are mapped onto. You should of course also eliminate the columns.

Whether you can use these depends on what you do with your matrix. Lagrange multiplier for example introduce 0 on the diagonal - not all solvers like that..

Anyway, this is very interesting research. I think that because of the specifics of your problem (as I understand it - irregularly placed particles in 3D, with fairly long-range interactions) make it very difficult to group the matrix entries. But I am very curious what you end up doing. Please let me know!

angainor
  • 11,760
  • 2
  • 36
  • 56
  • I originally intended to make this as black-box as possible, without the introduction of a coordinate system. The main reason is that the space we're dealing with can be periodic, (3-dimensional flat toruses, topologically speaking.) It is AFAIK, unknown if there exists a space-filling curve that can preserve locality in periodic spaces. The hilbert curve (by simply assuming that the space is non-periodic) is what we currently use, but we end up with boundary problems. CHMK, AMD, SAW, Morton, all perform significantly worse than the Hilbert Curve. – proteneer Sep 28 '12 at 18:06
  • hm I'm not too familiar with how Lagrange multipliers work. I'd have to read up on it. The latter option seems a bit difficult, when you say "number the rest", what exactly are the rest? Do you consider all 26 periodic images of the center one? – proteneer Sep 28 '12 at 22:06
  • @proteneer I see now that my comment regarding the numbering is confusing - I removed it. It was late Friday night yesterday ;) I'll try to write something up later tonight. The general point is that reordering a matrix that explicitly includes connectivities between periodic nodes yields much worse reorderings than for the same system, just not periodic. Hence it might be beneficial for you to introduce periodicity in another way (like Lagrange multipliers, or row-column manipulation on already reordered matrix). But that also depends on what you later do with the matrix... – angainor Sep 29 '12 at 09:20
  • phenomenal reply. Yes the example posted is a bit of a silly example. It is in-fact a simulation done with ~900 atoms, in a periodic box of 3 angstroms with a cutoff radius of 0.8. If I understand the Lagrange multiplier method, it seems to be basically making explicit copies of all atoms whose periodic copy has an interaction with some particular atom, and decoupling it (and in addition, introducing new indices for these replacement copies). Now, there are technical limitations for increasing the size of the matrix. In addition, I'm curious to see the tile count. – proteneer Sep 30 '12 at 01:03
  • I will also tryout METIS later on tonight to see how it performs on a much more realistic system (24k atoms, 6.223 angstroms, 0.8 angstroms). I can send you some more info (that's not directly relevant to this topic), if you shoot me an e-mail: proteneer at gmail dot com – proteneer Sep 30 '12 at 01:05
  • angainor - I tried METIS on a much larger system, and unfortunately it still does not beat the Hilbert curve. Do you know if METIS is sensitive on input? ie. if I give it a Hilbert Curve to start with, versus I give it a Morton curve to start with. – proteneer Oct 01 '12 at 19:51
  • @proteneer It is not sensitive to input - it essentially does a nested dissection recursively dividing the graph into partitions. I really think your best hope is to treat the periodicity differently.. But that is a very interesting research question you have, not really a technical problem to solve :) – angainor Oct 01 '12 at 20:05
  • my METIS doesn't quite look like yours, (using a different cutoff) http://i.imgur.com/iNMqV.png - I used a smaller cutoff. any idea why it seems to have artifacts in the sparse regions? What settings are you using? are you calling from MATLAB? – proteneer Oct 01 '12 at 23:45
  • @proteneer what do you mean by cutoff? What is white and black in your figure? Usually white is zeros, non-white is non-zeros on this type of plots. I use MATLAB metis mex file compiled from SuiteSparse. – angainor Oct 02 '12 at 08:01
  • two points must be closer than _cutoff_ to have an interaction (or to have a 1). I'm using inverted colors - black - 0, white - 1. I ended up also using metismex and got much better results. Strange... hmm. – proteneer Oct 02 '12 at 21:24
  • @proteneer Did you use matlab this time? I told you I think your RCM version of the original problem looks strange to me. You have non-zeros outside of the main 'band' - that should not be so. Might it be you have some problem with permutation? Anyway, by cutoff you effectively change the number of non-zeros in your model. I would definitely play with that if you do have any freedom in this respect. – angainor Oct 03 '12 at 10:29
  • @angainor: proteneer seems happy with your answer but what do you think of using a peano curve for a 3x3 or 9x9 matrix? A z order or hilbert curve is good for power of 2 grids – Micromega Oct 07 '12 at 09:54
  • @Chiyou The problem is that proteneer does not have a regular grid. He has an irregularly placed cloud of points.. This curve will not work in that case. You could first create a quad-tree structure, which can then be numbered in different ways, I guess. – angainor Oct 08 '12 at 10:59
  • @angainor: My idea is that peano would fit it better. It usese trees with 9-leafs. – Micromega Oct 08 '12 at 11:08
  • @Chiyou I am not saying it won't. But the problem is not as simple. In his problem the 'points' have long-range interactions, so not only the immediate neighbors are in contact. The density of those interactions also differ through throughout the space. I do not know what is the best approach. One has to try them all I guess. – angainor Oct 08 '12 at 11:12
0

You can look for a data structure like kd-tree, R-tree, quadtree or a space filling curve. Especially a space filling curve can help because it reduce the dimension and also reorder the tiles and thus can add some new information to the grid. With a 9x9 grid it's probably good to look into peano curves. The z order morton curve is better for power of 2 grids.

Micromega
  • 12,486
  • 7
  • 35
  • 72