The Cleaning Robot Puzzle: A Lower Bound from a Linear Program

Part 3’s fast route cleans a furnished 40×40 room in 1,327 commands, and all Part 3 could prove is that no route needs fewer than 1,037. Between the two lie 290 commands nobody could account for.

This is Part 8 of the cleaning robot puzzle, and the first of the parts that take the robot’s optimisation problems to the tools operations research really uses. Part 3 asked for the fewest commands that clean the whole floor, built a fast solver, an exact one that gives up on big rooms, and a lower bound from three simple facts. It ended on an open question: “On big furnished rooms the gap to the lower bound stays wide. The bound is weak there, and nobody knows how far from the best the fast routes really are.”

This part attacks that gap from below. It writes the question as a linear program, whose minimum is a lower bound no route can beat, and on that 40×40 room it proves 1,122: 85 commands higher than Part 3 could, which rules out 29% of the gap. On smaller rooms it often proves that the fast route is the best there is. Part 9 will come from above and close what is left. The new algorithms:

• modelling: a route written as a handful of numbers per field and per pair of fields, with rules they must obey;
• the simplex method, from scratch, in exact fractions, with Bland’s rule so it can never go round in circles;
• duality: prices on the rules that prove no answer is cheaper, and that turn out to rediscover Part 3’s checkerboard argument on their own;
• cutting planes: a program with more rules than there are atoms in the universe, solved by adding a rule only when a minimum cut, found with Edmonds and Karp’s algorithm, shows the current answer breaks it;
• HiGHS, the solver behind SciPy’s linprog, doing the same loop on the big rooms.

The simplex method and the cutting-plane loop here use only Python’s standard library, like everything in Parts 1 to 7. The big rooms need a real solver, and that part of the code uses HiGHS through its own Python package, highspy. Every file behind all the parts is listed in the code’s README.

The question again

The rules are Part 1’s: the robot starts on its field, moves one field at a time, must visit every floor field, and may finish anywhere. Part 3 asked for the fewest commands that do it, and gave three facts that no route can beat:

• One new field per command. With N floor fields, at least N − 1 commands.
• Checkerboard colours. Every command changes colour, so with S fields of the start’s colour and O of the other, at least 2S − 2 commands and at least 2O − 1.
• Dead ends. Every dead end but the last must be stepped out of again.

The fast solver’s routes are upper bounds: they exist, so the best is no worse. The three facts are lower bounds: the best is no better. When the two meet, the fast route is proven best. On Part 3’s 40×40 rooms with furniture they didn’t meet, and the exact search, which tries every (field, cleaned set) state, runs out of memory long before a room that size.

A route as numbers

A route is a walk. Close it with an imaginary free step from wherever it ends back to the start, and it becomes a closed walk, which has two properties any schoolchild who has tried to draw a figure without lifting the pen knows:

• every field is entered as often as it is left, so it touches an even number of steps;
• it is connected: every field is reachable from the start along its steps.

The converse holds too: a connected set of steps in which every field touches an even number of them can always be walked in one go. Euler stated it in his paper on the bridges of Königsberg, presented in 1735, without a proof; the first proof is Carl Hierholzer’s, written up from memory by colleagues and published in 1873, two years after his death (Hierholzer, 1873). So instead of searching for a route, count how many times it walks each pair of neighbouring fields, and ask for counts that are even at every field and connected. For each pair of neighbours $e$, let $x_e$ be the number of times the route walks between them, and for each field $v$, let $z_v$ be 1 if the route ends on $v$ and 0 if not. The free closing step joins the end to the start, and it costs nothing.

Here $\delta(S)$ is the set of neighbouring pairs with one field inside $S$ and one outside, and $x(\delta(S))$ the number of steps the route takes across that boundary. The start needs a parity rule of its own, since the closing step always touches it, and it is in the code. The connected rule says that the route crosses into every set of fields that doesn’t hold the start and comes out again, counting the closing step if the route ends inside.

Why at most 2 per pair? If a route walked between two fields three times, drop two of those walks: every field still touches an even number of steps, and nothing gets disconnected, because one walk still joins them. The same argument, about roads a salesman may drive more than once, is where the graphical travelling salesman problem starts (Cornuéjols, Fonlupt and Naddef, 1985). Every solution of this program is a route, and every shortest route is a solution, so its minimum is exactly Part 3’s answer.

Part 3’s facts are already inside

Drop parity and integrality and one might expect a weak bound. But look at the connected rule for a set of just one field $v$: $x(\delta(v)) + z_v \ge 2$. Add it up over every field of one colour. Every pair of neighbours has exactly one field of each colour, so the left side adds up to $\sum_e x_e$ plus the $z$ of that colour, at most 1, and the right side to twice the number of fields of that colour. That is Part 3’s checkerboard fact. A dead end has only one neighbour, so its rule says that one pair is walked at least $2 - z_v$ times: Part 3’s dead-end fact. The linear program knows everything Part 3 knew, and the test checks it: on 300 small rooms, its bound was never below Part 3’s.

What it knows beyond that comes from the rules for bigger sets. There is one for every set of fields, $2^{N-1} - 1$ of them: for a room of 1,000 fields, a number with 301 digits. They can’t all be written down. The way round that is the cutting-plane loop below. First, how to solve a linear program at all.

The simplex method

A linear program asks for the best point of a polyhedron, and the best point, when there is one, is always at a corner. George Dantzig’s simplex method walks from corner to corner of the polyhedron, along its edges, each step to a corner at least as good, until no edge leads anywhere better. Then that corner is the best, because a linear objective that doesn’t improve along any edge from a corner of a convex polyhedron doesn’t improve anywhere.

In more than two dimensions a corner is a choice of which variables are allowed to be nonzero, the basis, and moving along an edge swaps one variable into the basis and one out: a pivot. simplex.py keeps the rules in a table, the tableau, and does the pivots with Python’s Fraction, so there is no rounding anywhere:

def run(cost, phase, allowed):
    """Minimise cost.x from the current basis; returns False if unbounded."""
    while True:
        # reduced costs r_j = c_j - c_B . column j
        r = [cost[j] - sum(cost[basis[i]] * T[i][j] for i in range(m)) for j in range(cols)]
        enter = next((j for j in range(cols) if allowed(j) and j not in basis and r[j] < 0), None)
        if enter is None:
            return True
        best, leave = None, None
        for i in range(m):
            if T[i][enter] > 0:
                ratio = T[i][cols] / T[i][enter]
                if best is None or ratio < best or (ratio == best and basis[i] < basis[leave]):
                    best, leave = ratio, i
        if leave is None:
            return False
        out = basis[leave]
        pivot(leave, enter)
        ...                 # with trace=True, record the pivot (for the figure below)

The reduced cost $r_j$ of a variable outside the basis is how much the objective changes per unit of it brought in. If none is negative, no edge leads downhill and the corner is optimal. Otherwise one enters, and the ratio test picks the row that stops it first: bring it in any further and some basic variable would go negative.

Two details keep it honest:

• Phase I. The loop needs a corner to start from. When the rules include “at least” or “equal to”, the origin is not one, so the method first adds an artificial variable to each such rule and minimises their sum. If that can’t reach zero there is no answer at all; otherwise it ends at a genuine corner, and phase II starts from there.
• Bland’s rule. At a degenerate corner, where several rules meet more than they need to, a pivot can change the basis without moving, and a careless choice of pivots can go round in a circle forever. Always taking the lowest-numbered variable that improves, and breaking ties in the ratio test by the lowest-numbered one leaving, never cycles (Bland, 1977). The robot’s programs are highly degenerate, so this matters here.

The test runs 600 random small programs through it and through HiGHS: 127 with a best answer, 319 with no answer at all and 154 whose answers get better without end. Both agree on every one, on the verdict and on the value.

Prices that prove it

When the simplex method stops, how would anyone else know it found the minimum? It can hand over a proof, and the proof is as interesting as the answer.

Give every rule a price, a number $y_S \ge 0$ for each “at least 2” rule and a number $\mu$ for “the $z$s add up to 1”. Multiply each rule by its price and add them all up. If the prices are chosen so that no variable’s total coefficient exceeds its cost, which is 1 for each $x_e$ and 0 for each $z_v$, then every answer satisfies

because the left side is at least the priced sum of the rules’ left sides, which is at least the priced sum of their right sides. Finding the best such prices is itself a linear program, the dual, and the duality theorem says its maximum equals the original’s minimum. The simplex method finds both at once: the prices are sitting in the final tableau, in the columns of the rules’ slack and artificial variables.

simplex.check_certificate checks a claimed answer and claimed prices in exact fractions, without trusting whatever produced them: the answer obeys every rule, the prices have the right signs, no variable is overpriced, and the two totals are equal. The test checks the certificate of every one of the 127 random programs with an answer, and of every room’s final program.

The prices mean something. Here is the proof the program found for Part 3’s corridor, whose best route has 7 commands:

#########
#.....*.#
#########

It prices the rule $x(\delta(v)) + z_v \ge 2$ at 1 for four fields, the first, third, fifth and seventh from the left, and “the $z$s add up to 1” at −1. That gives 2 × 4 − 1 = 7. It is Part 3’s checkerboard argument, rediscovered: those four fields share a colour, each must be walked into and out of, and only the last can be the end. On Part 3’s 5×5 example, with the start in the middle, the certificate for 24 prices the twelve fields of the other colour at 1, plus one big set, every field except the start, also at 1, and the end at −2: 2 × 13 − 2 = 24.

On bigger rooms it finds proofs Part 3 never had. The figure below ends with one: a room where Part 3 could prove only 33 and the program proves 40 with 24 priced sets, 16 of them larger than one field, nested around the rooms and corridors the route has to enter and leave.

Adding rows until nothing is broken

The program has one rule for every set of fields without the start. For the 29-field room below that is $2^{28} - 1$, 268 million rules; for a 40×40 room, a number with hundreds of digits. The cutting-plane method never writes them all:

1. Start with a few rules: one for every field on its own, and one for every 2 × 2 block of fields.
2. Solve the program with the rules so far.
3. Look for a rule the answer breaks. If there is none, the answer obeys all of them, and its value is the program’s minimum.
4. Otherwise add the broken rules and go back to 2.

Step 3 is where the work is: among hundreds of millions of rules, find one that is broken, or prove none is. Think of $x_e$ as a capacity on the pair of fields $e$, and of $z_v$ as a capacity on the closing step from $v$ to the start. A rule $x(\delta(S)) + z(S) \ge 2$ is broken exactly when the total capacity across the boundary of $S$ is under 2: when $S$ is cut off from the start by a cut of capacity under 2. So step 3 is: for every field $v$, find the smallest cut between the start and $v$, and if it is under 2, the side of it that holds $v$ is a broken rule.

The smallest cut is found as a maximum flow. Push as much as possible from the start to $v$, respecting the capacities; Ford and Fulkerson proved in 1956 that the most that gets through equals the capacity of the smallest cut, and the fields the flow can still reach at the end are its start side. maxflow.py pushes along shortest paths, which Edmonds and Karp showed takes a number of steps polynomial in the size of the graph whatever the capacities are:

def min_cut(self, s, t, enough=None, eps=0):
    """Returns (flow value, set of nodes on s's side of a minimum cut). With `enough`, stops as
    soon as the flow reaches it and returns (value, None): the cut is at least that big."""
    head, arcs = self.head, self.arcs
    residual = self.capacity[:]
    flow = 0
    while True:
        came_by = [-1] * self.n                 # the arc each node was reached by
        came_by[s] = -2
        queue = deque([s])
        while queue and came_by[t] == -1:
            u = queue.popleft()
            for a in arcs[u]:
                v = head[a]
                if came_by[v] == -1 and residual[a] > eps:
                    came_by[v] = a
                    queue.append(v)
        if came_by[t] == -1:
            return flow, {v for v in range(self.n) if came_by[v] != -1}
        path, v = [], t
        while v != s:
            a = came_by[v]
            path.append(a)
            v = head[a ^ 1]                     # arc a ^ 1 runs the other way
        push = min(residual[a] for a in path)
        for a in path:
            residual[a] -= push
            residual[a ^ 1] += push
        flow += push
        if enough is not None and flow >= enough - eps:
            return flow, None

Pushing flow along an arc makes room to push it back along the reverse arc, which is what lets a later path undo an earlier bad choice. The enough argument stops as soon as 2 units get through: then no rule is broken for this field, and the exact size of the cut doesn’t matter.

Running a flow to every field every round is the slow part, so route_lp.broken_cuts tries something cheaper first. The fields reached only by steps the answer walks at least ¼, ½, ¾ or all the way, starting from the start, form a piece; every other piece is a candidate set, and its rule is checked directly. Only when none of those is broken does it run the flows. The 2 × 2 blocks in the first round are there for the same reason: a little loop round four fields, cut off from everything, is the commonest way the program cheats.

In the first round the program’s cheapest answer costs 29, which is Part 3’s first fact, N − 1: it walks every field twice in little loops that never meet the start. Each round the minimum cut finds some of those loops, their rules join the program, and the answer has to connect them to the start some other way. After 19 rounds nothing is broken, the bound is 40, and the fast route, 40 commands, is proven the best there is.

HiGHS on the big rooms

The exact simplex method is fine for rooms of a few dozen fields: on the 300 small rooms of the test it took 618 ms per room on average. The 40×40 rooms have a thousand fields and thousands of rules, and there a tableau of fractions is hopeless. The same loop runs on HiGHS, an open-source solver from the University of Edinburgh, and the default behind SciPy’s linprog since version 1.9. It took 4 ms per small room, on the same rooms, with the same answers.

It is the same method with every engineering trick. It runs the dual simplex method, which walks corners of the dual program (Huangfu and Hall, 2018); that suits a cutting-plane loop perfectly, because adding a rule leaves the previous corner feasible for the dual, so each round restarts from the last round’s basis instead of from nothing. It keeps the basis as a sparse factorisation rather than a dense table. And it works in floating point, which is why route_lp.py treats a rule as broken only by more than 10⁻⁷, and why the test checks that the exact and floating-point loops agree on every small room. For the big rooms the loop keeps one HiGHS model for the whole run and adds each round’s rules to it:

def add(self, cuts):
    np = self.np
    starts, index, values = [], [], []
    for S in cuts:
        row = self.model.cut_row(S)
        starts.append(len(index))
        index += row.keys()
        values += row.values()
    self.h.addRows(len(cuts), np.full(len(cuts), 2.0), np.full(len(cuts), np.inf), len(index),
                   np.array(starts, dtype=np.int32), np.array(index, dtype=np.int32),
                   np.array(values, dtype=float))
    self.cuts += cuts

What the bound buys

One room of each kind, from the report:

And across random furnished rooms:

The program never does worse than Part 3’s three facts, and on furnished rooms it does much better. On the 300 small rooms, where Part 3’s exact search can still find the true best, it was the true best every time. What it can’t say on big rooms is which of two things the remaining gap is: a fast route that could be shorter, or a bound that could be higher. Parity is the missing rule, and the linear program’s answers use it freely: an answer can walk a pair of fields half a time, or finish a fraction of the way through a field. Part 9 puts parity and whole numbers back, which makes it an integer program, and solves it exactly with branch and bound: the gap closes from both sides at once.

Testing Part 8

test_route_lp.py produces every number above:

What Part 8 teaches

• Write the question down as numbers and rules. A route is a count per pair of fields that is even, connected and small; once it is written that way, fifty years of tools apply.
• Relaxing gives bounds. Drop the hard rules and the minimum can only go down, so whatever the relaxation proves, the real problem obeys.
• An answer can carry its own proof. The simplex method’s prices certify the minimum, can be checked by anyone in exact arithmetic, and on the corridor turn out to be Part 3’s checkerboard argument.
• You don’t need every rule, only a way to find a broken one. Hundreds of millions of rules, and a minimum cut finds the one that matters.
• Teach with the small solver, measure with the big one. The exact simplex method shows what happens; HiGHS, checked against it on every small room, makes the big rooms possible.