how to sort a matrix into a block matrix (zero blocks), preserving naming order on ties using row/col Sums?

Question

I have an adjacency matrix:

myV = c(0,0,0,0,0,0,0,0,0,0,
        0,0,0,0,0,0,0,0,0,0,
        0,0,0,0,0,0,0,0,0,0,
        0,1,1,0,0,0,0,0,0,0,
        1,1,0,0,0,0,0,0,0,0,
        0,0,0,0,0,0,0,0,0,0,
        1,1,0,0,1,1,0,0,0,0,
        1,0,0,1,0,1,0,0,0,0,
        0,0,1,0,0,0,0,0,0,0,
        0,0,0,0,0,1,0,1,0,0);

myA = matrix(myV, nrow=10, byrow=TRUE);
  rownames(myA) = colnames(myA) = paste0("P.",1:10);
myA;

That I want to sort jointly by rows and columns.

I can sort independently by rows, then columns ... or columns, then rows.

# row sort
rs = rowSums(myA);
rs.s = sort(rs);
rs.idx = as.numeric( gsub("P.","",names(rs.s)) );

myA.s = myA[rs.idx,];
myA.s;

Outputs:

     P.1 P.2 P.3 P.4 P.5 P.6 P.7 P.8 P.9 P.10
P.1    0   0   0   0   0   0   0   0   0    0
P.2    0   0   0   0   0   0   0   0   0    0
P.3    0   0   0   0   0   0   0   0   0    0
P.6    0   0   0   0   0   0   0   0   0    0
P.9    0   0   1   0   0   0   0   0   0    0
P.4    0   1   1   0   0   0   0   0   0    0
P.5    1   1   0   0   0   0   0   0   0    0
P.10   0   0   0   0   0   1   0   1   0    0
P.8    1   0   0   1   0   1   0   0   0    0
P.7    1   1   0   0   1   1   0   0   0    0

And

cs = colSums(myA.s);
cs.s = sort(cs);
cs.idx = as.numeric( gsub("P.","",names(cs.s)) );

myA.ss = myA.s[,cs.idx];
myA.ss;

Outputs:

     P.7 P.9 P.10 P.4 P.5 P.8 P.3 P.1 P.2 P.6
P.1    0   0    0   0   0   0   0   0   0   0
P.2    0   0    0   0   0   0   0   0   0   0
P.3    0   0    0   0   0   0   0   0   0   0
P.6    0   0    0   0   0   0   0   0   0   0
P.9    0   0    0   0   0   0   1   0   0   0
P.4    0   0    0   0   0   0   1   0   1   0
P.5    0   0    0   0   0   0   0   1   1   0
P.10   0   0    0   0   0   1   0   0   0   1
P.8    0   0    0   1   0   0   0   1   0   1
P.7    0   0    0   0   1   0   0   1   1   1

As an adjacency matrix, that is not exactly what I want ... The rows and cols don't match... e.g., P.1 as a row and P.1 as a column needs to be in the same location. Adjacency matrix. Symmetry.

So how do I sort simultaneously (jointly), I would like to be able to block it into partions of zero sub-matrices.

In this example, I am looking for an answer like 1,2,3,6 (rowSum=0, do I care about ties if I need symmetry?) ... 4,5,8 (colSum != 0 && rowSum !=0, natural order for symmetry?) ... 7,9,10 (colSum=0, natural order again, or use rowSum to sort on ties).

The goal is to place the most dense "block" in the bottom-left hand corner of the matrix ... And the current correct answer 1,2,3,6 ... 4,5,8 ... 7,9,10 may not accomplish that. Or maybe it does given the symmetry constraint.

Obviously I am looking for a function or code for a generalized solution and offer this example to demonstrate my need. I would prefer a solution using base R that accomplishes this block sorting so the "top" and the "right" has sub-grids of zero matrices, and the adjacency symmetry is maintained. A general c-based solution I could port, if given the correct pseudo-code logic.

How to accomplish this?

Something like: myA[order,order]; ... how to build order?

Answer 1

TL DR

You can solve the problem by implementing a scoring method, combined with a hill climbing algorithm.

Summary

A scoring method takes an adjacency matrix as input, and returns a scalar value that indicates the quality of the matrix. The scoring method can be positive (higher scores are better), or negative (lower scores are better). I'll show a negative scoring method, based on a penalty that must be paid to place a 1 at a position in the matrix.

The hill climbing algorithm is based on simple trial and error. Given a matrix and its score, we attempt to improve the matrix by swapping rows. If a swap results in a better score, keep it. Otherwise, try a different swap.

When rows are swapped, the corresponding columns need to be swapped as well to maintain the adjacency matrix symmetry. So there's only one degree of freedom, which is the row order. With 10 rows, as in the example, there are only 10! = 3628800 possible row orders, so a brute force algorithm can try them all, and find the one that gives the lowest score.

Results

The adjacency matrix in the question scores 204. A brute force search for the best row order takes 2.3 seconds to find the matrix shown below, with a score of 143.

0 0 0 0   0 0 0   0 0 0 
0 0 0 0   0 0 0   0 0 0 
0 0 0 0   0 0 0   0 0 0 
0 0 0 0   0 0 0   0 0 0 

0 1 0 1   0 0 0   0 0 0 
1 1 0 0   0 0 0   0 0 0 
1 0 1 0   1 0 0   0 0 0 

0 0 0 1   0 0 0   0 0 0 
0 0 1 0   0 0 1   0 0 0 
1 1 1 0   0 1 0   0 0 0 

Hill climbing finds a very similar matrix (with score 143) in 5 milliseconds.

0 0 0 0   0 0 0   0 0 0 
0 0 0 0   0 0 0   0 0 0 
0 0 0 0   0 0 0   0 0 0 
0 0 0 0   0 0 0   0 0 0 

1 0 0 1   0 0 0   0 0 0 
1 1 0 0   0 0 0   0 0 0 
0 1 1 0   1 0 0   0 0 0 

0 0 0 1   0 0 0   0 0 0 
0 0 1 0   0 0 1   0 0 0 
1 1 1 0   0 1 0   0 0 0 

Computing the score

The penalty matrix I used looks like this:

penalty[10][10] =
{
 45  46  48  51  55  60  66  73  81  90 
 36  37  39  42  46  51  57  64  72  81 
 28  29  31  34  38  43  49  56  64  73 
 21  22  24  27  31  36  42  49  57  66 
 15  16  18  21  25  30  36  43  51  60 
 10  11  13  16  20  25  31  38  46  55 
  6   7   9  12  16  21  27  34  42  51 
  3   4   6   9  13  18  24  31  39  48 
  1   2   4   7  11  16  22  29  37  46 
  0   1   3   6  10  15  21  28  36  45 
}

Interpretation: the penalty for putting a 1 in the lower left of the adjacency matrix is 0. But putting a 1 at the upper right has a penalty of 90. In general, the penalty for a 1 at adjacency[row][col] is penalty[row][col]. The total score is the sum of the penalties. Lower scores are better.

To generate a penalty matrix of size NxN, start with the bottom row, and the first column. The values in the bottom row are penalty[N-1][i] = sum(0..i). The first column has the same values from bottom to top. All the other values are just the sum of the row and column penalties, i.e. penalty[row][col] = penalty[N-1][col] + penalty[row][0]. The largest number in the penalty matrix is N*(N-1), e.g. 10*9=90 for the 10x10 matrix.

Climbing the hill

The score can be improved by rearranging rows in the adjacency matrix. We can model this with a row order array. Initially the array elements are sequential starting at zero, i.e. order[10] = {0,1,2, ... ,9}.

The simplest method to improve the row order is a "bubble sort"¹. This means repeatedly scanning the order array, swapping adjacent elements if the swap would improve the score. The bubble sort is finished when there are no more beneficial swaps between adjacent elements.

The bubble sort method is good enough to solve the example in the question. However, it's possible that the score can still be improved even though swapping adjacent rows doesn't help. After the bubble sort stops, we need to look at non-adjacent swaps. There are Nc2 = N*(N-1) / 2 pairs that could be swapped. Each of those should be checked. If an improvement is found, then it's back to the bubble sort.

The next level is to take all Nc3 triplets, and try rotating left and right.

Worst case, finding the optimal solution requires trying all N! possible row orders. But an implementation that

• bubble sorts
• swaps non-adjacent pairs, and
• rotates triplets

should get pretty close to the optimum in a reasonable amount of time.

Pseudo-code for the hill climbing algorithm

swapped = true
while (swapped)

   // bubble sort
   while (swapped)
      swapped = false
      for each pair of adjacent elements
         if swapping the elements improves the score
            swap the elements
            swapped = true

   // swapping non-adjacent pairs
   for all possible pairs
      if swapping the pair improves the score
         swap the pair
         swapped = true
         break

   // rotating triplets
   if not swapped
      for all possible triplets
         if rotating the triplet (left or right) improves the score
            rotate the triplet
            swapped = true
            break

---

^{1) A real bubble sort needs total order to work. We don't have that with this problem. The code to sort the order array looks a lot like a bubble sort, but it's not guaranteed to find the final answer.}