This post is the third in a series of posts about my computational approach to the Zarankiewicz problem.

A quick summary of where we left off:

• We're looking for the maximum possible number of $1$s in a grid with $m$ rows and $n$ columns, with the constraint that no four $1$s form a rectangle. This is called the Zarankiewicz problem.
• Rather than try every possible grid, we're going to build up grids row-by-row, and reject partial grids which definitely can't be completed (because they already contain a rectangle).

A grid can only contain $0$s and $1$s, which computers are good at manipulating, so we'll represent each row as a binary number; a grid will be an array of rows.

In the previous post we discussed how partial grids form a tree structure, with completed solutions at the ends of each branch. Let's elaborate on that a bit more. Each node in the tree represents either a partial or complete grid. There are a very large number of possible grids, so it's hard to think about the whole tree at once. It's also unnecessary; the power of algorithms is that they solve problems by accumulating the effects of small, manageable steps, and at each step we only need to consider a small part of the tree.

What does part of the tree look like?

 $\vdots$ $\bbox[border:1px solid black]{\begin{matrix} 011 \\ \bbox[2pt,#8080FF]{101} \\ * \\ * \end{matrix}}$ $\swarrow$ $\swarrow$ $\swarrow$ $\swarrow$ $\searrow$ $\searrow$ $\searrow$ $\searrow$ $\bbox[border:1px solid black]{\begin{matrix} 011 \\ 101 \\ \bbox[2pt,#80FF80]{000} \\ * \end{matrix}}$ $\bbox[border:1px solid black]{\begin{matrix} 011 \\ 101 \\ \bbox[2pt,#80FF80]{001} \\ * \end{matrix}}$ $\bbox[border:1px solid black]{\begin{matrix} 011 \\ 101 \\ \bbox[2pt,#80FF80]{010} \\ * \end{matrix}}$ $\bbox[border:1px solid black]{\begin{matrix} 011 \\ 101 \\ \bbox[2pt,#FF8080]{011} \\ * \end{matrix}}$ $\bbox[border:1px solid black]{\begin{matrix} 011 \\ 101 \\ \bbox[2pt,#80FF80]{100} \\ * \end{matrix}}$ $\bbox[border:1px solid black]{\begin{matrix} 011 \\ 101 \\ \bbox[2pt,#FF8080]{101} \\ * \end{matrix}}$ $\bbox[border:1px solid black]{\begin{matrix} 011 \\ 101 \\ \bbox[2pt,#80FF80]{110} \\ * \end{matrix}}$ $\bbox[border:1px solid black]{\begin{matrix} 011 \\ 101 \\ \bbox[2pt,#FF8080]{111} \\ * \end{matrix}}$ $\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$

This picture contains enough information to think about one step in the algorithm. This is what happens when we "visit" the current node, which has just been reached by adding the blue row, $\bbox[2pt,#8080FF]{101}$:

1. Since the partial grid at this node doesn't contain any rectangle so far, we can continue.
2. Descend to the next layer of the tree, by inserting $\bbox[2pt,#80FF80]{000}$.
3. "Visit" the node with the partial grid we just created.
4. Once we're done with that, advance to the next sibling in the tree, and do step 3 again.
5. Once we're done with $\bbox[2pt,#FF8080]{111}$, we can't advance any further, so we've finished visiting the blue node. Delete the child, and go back up to the previous layer in the tree.

This is a recursive algorithm: "visiting" a node requires "visiting" its children.

Note that three of the new partial solutions, in red, have rectangles; these won't need to be explored any further. There are some more special cases to handle: what happens when we reach the bottom of the tree (i.e. there are no more rows to add), and what happens when we "go back up one" from the top of the tree (i.e. there are no rows to delete). These occur when we find a solution, and when the algorithm terminates.

The fact that we can sometimes reject nodes without having to visit all of their children is what distinguishes this algorithm from brute force. However, it's not quite the same as the algorithm for solving (weak) combination locks: in that problem, once a digit worked, we knew it was correct, so we could simply continue solving by moving onto the next digit. That's a greedy algorithm. For this problem, there are two subtle but important differences:

1. Completing a grid doesn't mean we've finished, because other grids might turn out to be better solutions, with more $1$s.
2. We don't immediately know whether a partial grid should be rejected or not: a grid might be rejected immediately (e.g. because it already contains a rectangle), but it also might be rejected later due to all of its children being rejected. This is like a dead end in a maze: a corridor is a dead end if it's blocked right in front of you, but it could also be a "dead end" in the sense that each way forward eventually leads to a dead end itself.

Both are addressed by the instruction to "go back up" in step 5; like solving a maze, once we we've finished exploring everywhere we can get to from this node, we must backtrack and explore somewhere else. We've designed a backtracking algorithm.

We'll write in Java. I've put the full code in a git repository, but let's look at the important part:

// solve using a backtracking algorithm
public int solve() {
grid = new Grid(width, height);
// count how many nodes we visit
nodesVisited = 0;
// this method does all of the work.
visitCurrentNode();
// after visiting the root node, we're done!
return bestWeight;
}

// visit one node in the search tree
private void visitCurrentNode() {
// we are visiting another node, so increment this counter
nodesVisited++;

if(rejectNode()) { return; }

if(grid.isComplete()) {
// we reached the bottom of the search tree, so the grid is complete.
int newWeight = grid.weight();
// check if this solution is better than what we had before
if(newWeight > bestWeight) {
bestWeight = newWeight;
// uncomment this line to see the solutions found.
//System.out.println(grid);
}
return;
}

// move down the tree to visit this node's children
descend();

// visit every child node
while(true) {
// visit this child node
visitCurrentNode();
// stop if there are no more siblings to visit
if(!canAdvance()) { break; }
// otherwise move on to the next sibling
}

// after visiting all of this node's children, backtrack
ascend();
}


Notice that this is written fairly abstractly: this makes it easier to understand the code, and also to improve the algorithm later.*

To fill in some of the missing details, we're representing a row as a Row object, which represents the current state of the row as a binary number; for now, descend() adds a new Row like $00000$ (represented by the number 0) to gridadvance() increments this number, and ascend() deletes it once we've tried all of its possibilities.

The only remaining non-trivial detail is the code to check whether the grid has a rectangle. To do this, we take advantage of the binary number representation in the Row class:

// determines whether any rectangle occurs across these two rows
boolean hasRectangleWith(Row other) {
int t = this.asNumber & other.asNumber;
return (t&(t-1)) != 0;
}


A rectangle occurs when two rows have two or more $1$ bits in common in their binary numbers. The bitwise operation & finds the common $1$ bits. The next line determines whether the result t has two or more $1$ bits; it's non-obvious why this works, but rather than explain it here, I think it's a good puzzle to try yourself.**

At last, we have an algorithm! It's certainly not efficient, but it beats brute force:***

By one measure, we only beat brute force by $1$: our backtracking algorithm can solve a $7{\times}7$ grid before the search tree explodes beyond feasibility, whereas brute force can get up to $6{\times}6$. By another measure, we beat brute force by a factor of about $400$: for a $7{\times}7$ grid, our algorithm visits about $1/400$th the number of nodes that brute force would. This will be a recurring theme as we continue: we can improve the algorithm by orders of magnitude, again and again, but it only gets us a little further each time.

In general, there's not going to be much point drawing graphs like this; every algorithm is going to cling to the $n$ axis until it explodes. If we squash the vertical axis logarithmically, we can see more detail:

This is comparing by the number of nodes visited. The efficiency of the algorithm is mostly determined by how the number of nodes visited grows with the size of the grid, but of course this doesn't correspond exactly with running time. It might still be worth visiting more nodes if we can visit each node faster. That said, the time taken to visit each node is not going to make a big difference in what inputs are feasible or infeasible; what matters is how long we can stave off the explosion. We're not trying to squeeze every millisecond out of the algorithm.

Where to go from here? There are many avenues for improvement, and we've identified a few already.

Since we want to find good solutions as soon as possible, perhaps we should start with a row like $11111$, and make advance() decrement this number. We're trying to maximise the number of $1$s, so perhaps we should try bigger numbers first, as they tend to have more $1$s in binary.

We are rejecting partial grids which have rectangles; however, we could also reject partial grids if we can know they won't lead to an optimal solution. To do this, we need an upper bound on how many $1$s the completion of a partial grid can have, and if this is too low, the node can be rejected by the rejectNode() method.

How about symmetry? We're free to permute the rows and columns in a solution, while preserving the rectangle-free property and the number of $1$s in the grid. A large proportion of the nodes are symmetrical to other nodes which we'll already have visited; if we can identify which nodes, we don't need to visit them.

We'll explore these options in the next post in this series.

*Of course, since I'm writing these posts long after writing the algorithms, I can say with the benefit of hindsight that we will indeed be able to improve substantially on the algorithms described here.

**Even if you won't find any use for them, bit manipulation hacks like this are at least good logic puzzles. Of course, because this is Java we could have written Integer.bitCount(t) > 1; but that's no fun.

***We can implement the brute force algorithm by making rejectNode() always return false for partial grids; this forces it to visit the whole search tree. That said, brute force is far too slow and I really don't want to wait around for it; instead, we'll use the exact formula$$\sum_{i=0}^{m} 2^{in} = \frac{2^{(m+1)n} - 1}{2^{n} - 1}$$to calculate exactly how many nodes the brute force algorithm would visit; $2^{in}$ is the number of nodes in layer $i$ of the tree.