linear programming – Book Proofs

Tiling squares

This week’s Fiddler is about tiling a square with smaller squares.

Suppose you have infinitely many 3-by-3 cm tiles and infinitely many 5-by-5 cm tiles. You want to use some of these tiles to precisely cover a square whose side length is a whole number of centimeters. Tiles may not overlap, and they must completely cover the larger square, without jutting beyond its borders. What is the smallest side length this larger square can have, such that it can be precisely covered using at least one 3-by-3 tile and at least one 5-by-5 tile?

Extra credit:
This time, you have an infinite supply of square tiles for each odd whole number side length (as measured in centimeters) greater than 1 cm. In other words, you have infinitely many 3-by-3 cm tiles, infinitely many 5-by-5 cm tiles, infinitely many 7-by-7 cm tiles, and so on. You want to use one or more of these tiles to precisely cover a square whose side length is $N$ cm, where $N$ is an integer. Once again, tiles may not overlap, and they must completely cover the larger square without jutting beyond its borders. What is the largest integer N for which this task is not possible?

My solution:
[Show Solution]

This problem is about tiling squares with smaller squares, which is a famous problem known as “squaring the square”. There is a nice wikipedia article about it, and also a website I found with a comprehensive database of different tilings. There are certain properties of interest when it comes to squaring the square:

simple tiling: it contains no sub-tilings
perfect tiling: all sub-squares are of different size
symmetric: there are certain rotational or reflection symmetries present in the final shape.

However, none of this is going to help us, since we don’t care about these properties for the purpose of this problem!

In our case, we can can repeat tiles (not perfect), we can have sub-tilings (not simple) and we do not enforce any special symmetry. Moreover, we are restricted in the kinds of tiles we can use. To tackle this problem, I formulated the problem as a mixed integer program using a column enumeration approach. This is the same approach I used to solve a similar past Riddler problem from 2017, so if you’re interested in how I did the modeling, be sure to read that previous post!

First problem: minimum size

I re-used the code from the 2017 post linked above and made two changes:

restricted the available tiles to 3×3 and 5×5 only.
added a constraint that requires using at least one of each type of tile.

I looped through possible total sizes starting from 4×4 and stopped once the integer program found a tiling. The smallest size is $N=18$ and one possible tiling is shown below.

The code made short work of this problem, testing all $N\leq 18$ in a couple seconds.

Second problem: maximum size

The second problem is a bit trickier, since we are asked for the largest square that cannot be tiled. So we’ll need to do a bit of deductive work before we turn to computation. If there exists a tiling of an $N\times N$ square using only odd-sized square tiles, I will say that $N$ is “tileable”. Now some observations:

Any odd $N$ is tileable (use a single large tile of size $N\times N$).
If $M$ is tileable, then $N = k M$ is tileable for any $k$. We can simply start with the tiling for the $M\times M$ square, and then repeat it in a $k\times k$ pattern to obtain the $N\times N$ tiling.
By items 1 and 2, if $N$ is not tileable, it must be a power of $2$.
By item 2, if $N = 2^k$ is tileable, so is any larger power of $2$.

In other words, it suffices to find the smallest tileable power of $2$, and all larger $N$ will also be tileable. Using the code again, this time allowing for tiles of all odd sizes, I was able to verify that $2$, $4$, $8$, and $16$ are not tileable, but $32$ is tileable!!! Here is one possible tiling, which uses tiles of size $3$, $5$, $7$, and $11$.

Therefore, we conclude that $N=16$ is the largest integer for which an $N\times N$ square cannot be tiled using smaller squares with odd sidelengths greater than $1$.

Computational note: I used the Gurobi solver to solve this problem, and it took about 16 seconds on my laptop. I also tried Mosek, but it was slower, coming in at about 40 seconds.

A contradiction

How is it possible that the smallest tileable $N$ is $18$ but the largest non-tileable $N$ is 16?

These are actually different problems! For the first problem, we were only allowed to use tiles of size 3×3 and 5×5. In this case, 18 is the smallest tileable number, but not all subsequent numbers are tileable! It turns out 19, 22, 23, 26, 28, 29 (and more…) are not tileable. In the second problem, we had more tiles at our disposal, so naturally there must be fewer non-tileable $N$’s in this case.

Ostomachion coloring

The following problems appeared in The Riddler. And it features an interesting combination of an ancient game with the four-color theorem.

The famous four-color theorem states, essentially, that you can color in the regions of any map using at most four colors in such a way that no neighboring regions share a color. A computer-based proof of the theorem was offered in 1976.

Some 2,200 years earlier, the legendary Greek mathematician Archimedes described something called an Ostomachion (shown below). It’s a group of pieces, similar to tangrams, that divides a 12-by-12 square into 14 regions. The object is to rearrange the pieces into interesting shapes, such as a Tyrannosaurus rex. It’s often called the oldest known mathematical puzzle.

Your challenge today: Color in the regions of the Ostomachion square with four colors such that each color shades an equal area. (That is, each color needs to shade 36 square units.) The coloring must also satisfy the constraint that no adjacent regions are the same color.

Extra credit: How many solutions to this challenge are there?

Here are the details of how I solved the problem:
[Show Solution]

I’d like to start by pointing out that Hector Pefo also posted a nice solution to this problem. There were also diagrams of solutions posted on Twitter, such as those from @DimEarly and from @chattinthehatt. The solution I’ll describe here is similar to Hector’s, but I’ll explain how to model it as an integer linear program and solve it that way. Of course both solutions produce the same answers (because they’re both correct!) but it’s nice to see different approaches. So without further ado…

Integer linear program

Let’s suppose there are $N$ shapes and we’d like to use $C$ colors. Define:
\[
x_{ij} = \begin{cases}1&\text{if shape }i\text{ is colored using color }j\\0&\text{otherwise}\end{cases}
\]The first constraint is that each shape must be assigned to exactly one color. We can express this algebraically as follows:
\[
\sum_{j=1}^C x_{ij} = 1\qquad\text{for }i=1,\dots,N
\]Now suppose that shape $j$ has area $a_j$. The constraint that each color has an equal share of the area means that the sum of the areas of the shapes of each color must be $A_\text{tot}/C$, where $A_\text{tot}$ is the total area of the Ostomachion. Algebraically, we can write this as:
\[
\sum_{i=1}^N a_i x_{ij} = A_\text{tot}/C\qquad\text{for }j=1,\dots,C
\]Finally, we must deal with the edge constraints. We say that $(p,q)$ is an “edge” if shape $p$ shares an edge with shape $q$. Shapes that share an edge must be different colors. We can encode this algebraically as:
\[
x_{pj} + x_{qj} \le 1\qquad\text{for all edges }(p,q)\text{ and for }j=1,\dots,C
\]And that’s it! We described the feasible set of solutions as the set of matrices $x \in \{0,1\}^{N\times C}$ that satisfy the three sets of linear constraints above. This sort of optimization model is a modified version of a Knapsack problem. More generally, it’s an Integer Linear Program (ILP). Such problems are known to be NP-hard, which roughly means that in the worst-case, there is no efficient way to solve them short of trying all possible colorings. However, modern ILP solvers are quite sophisticated and can usually solve these types of problems very efficiently.

Finding all solutions

When using an ILP model, most commercial solvers are designed to find a single solution. For this problem, we are interested in finding all solutions. One way to get around the limitation is to use Barvinok’s algorithm. One can also use variants of branch-and-cut. Yet another way, which is not as systematic but much easier to explain and implement, is to use a randomized approach.

Since all our variables are binary (0 or 1), the convex hull of the feasible set cannot contain any interior solution points. All solutions will be on the boundary. This means that any solution can be found by specifying the right objective function for the optimization. As stated, our optimization is just a feasibility problem and there is no objective specified. For this problem, we start by letting $r_{ij}$ be random numbers (we can draw them i.i.d. from a standard normal distribution, for example). Then, we simply add the objective function:
\[
\text{Minimize}\qquad \sum_{i=1}^N\sum_{j=1}^C r_{ij} x_{ij}
\]and solve the problem. If we do this repeatedly with different realizations of the $\{r_{ij}\}$, we will eventually enumerate all possible solutions of the optimization problem.

To jump straight to the solution (and pretty pictures!) see below.
[Show Solution]

Squaring the square

This Riddler puzzle is about tiling a square using smaller squares.

You are handed a piece of paper containing the 13-by-13 square shown below, and you must divide it into some smaller square pieces. If you are only allowed to cut along the lines, what is the smallest number of squares you can divide this larger square into? (You could, for example, divide it into one 12-by-12 square and 25 one-by-one squares for a total of 26 squares, but you can do much better.)

Here is how I solved the problem:
[Show Solution]

This is a challenging puzzle and there is unfortunately no clever trick one can use to solve it quickly. If we consider the general problem of tiling an $N\times N$ grid using smaller squares, we can make several observations:

The smallest answer we could ever hope for is $4$. If $N=2k$ is even, we can achieve the minimum by splitting the square into four $k\times k$ squares. It’s clear that this solution cannot be improved.
If $N = 3k$ is a multiple of $3$, we can split it into six pieces: one $2k\times 2k$ square and five $k\times k$ squares. I believe this also cannot be improved.
When $N$ has many divisors, it becomes easier to tile it. The most challenging cases are therefore the cases where $N$ is prime. So a larger grid doesn’t necessarily mean a more difficult problem! We’ll see examples of this later.

I’ll solve this problem by converting it into an integer linear program (specifically, a variant of a knapsack problem), which can be solved much more efficiently than by brute-force checking the astronomical number of possible tilings.

I’ll illustrate the approach for the $5\times 5$ case. We’ll represent each possible square as a $5\times 5$ matrix of $0$’s and $1$’s. For example, we can use $25$ different $1\times 1$ tiles:
\[
\scriptsize \begin{bmatrix}{\color{red}{\color{red}1}}&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix},
\begin{bmatrix}0&{\color{red}1}&0&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix},
\begin{bmatrix}0&0&{\color{red}1}&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix},\,\dots\,,
\begin{bmatrix}0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&{\color{red}1}\end{bmatrix}
\]We can also use $16$ different $2\times 2$ tiles:
\[
\scriptsize \begin{bmatrix}{\color{red}1}&{\color{red}1}&0&0&0\\{\color{red}1}&{\color{red}1}&0&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix},
\begin{bmatrix}0&{\color{red}1}&{\color{red}1}&0&0\\0&{\color{red}1}&{\color{red}1}&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix},
\begin{bmatrix}0&0&{\color{red}1}&{\color{red}1}&0\\0&0&{\color{red}1}&{\color{red}1}&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix},\,\dots\,,
\begin{bmatrix}0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&{\color{red}1}&{\color{red}1}\\0&0&0&{\color{red}1}&{\color{red}1}\end{bmatrix}
\]Continuing in this fashion, we have $9$ different $3\times 3$ tiles and $4$ different $4\times 4$ tiles. So in total, we have $25+16+9+4=54$ different tiles. The name of the game is to pick the smallest possible subset of these matrices that sums to the matrix of all $1$’s. We can represent this compactly by vectorizing each of the $54$ matrices above into $25\times 1$ columns, and then stacking the columns side by side to form a matrix $A \in \{0,1\}^{25\times 54}$. Our goal is then to solve the problem:
\begin{align}
\underset{x \in \{0,1\}^{54}}{\text{minimize}}\qquad& \sum_{i=1}^{54} x_i \\
\text{subject to:}\qquad & Ax = 1
\end{align}Here, each $x_i$ is a binary variable that selects whether we’ll be using the $i^\text{th}$ column of $A$ or not. It turns out the optimal solution for this $5\times 5$ case uses $8$ squares; or the sum:
\begin{multline}
\scriptsize \begin{bmatrix}0&0&0&0&0\\0&0&0&0&0\\{\color{red}1}&{\color{red}1}&{\color{red}1}&0&0\\{\color{red}1}&{\color{red}1}&{\color{red}1}&0&0\\{\color{red}1}&{\color{red}1}&{\color{red}1}&0&0\end{bmatrix}
+\begin{bmatrix}{\color{red}1}&{\color{red}1}&0&0&0\\{\color{red}1}&{\color{red}1}&0&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix}
+\begin{bmatrix}0&0&0&{\color{red}1}&{\color{red}1}\\0&0&0&{\color{red}1}&{\color{red}1}\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix}
+\begin{bmatrix}0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&{\color{red}1}&{\color{red}1}\\0&0&0&{\color{red}1}&{\color{red}1}\end{bmatrix}\\
\scriptsize+\begin{bmatrix}0&0&{\color{red}1}&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix}
+\begin{bmatrix}0&0&0&0&0\\0&0&{\color{red}1}&0&0\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix}
+\begin{bmatrix}0&0&0&0&0\\0&0&0&0&0\\0&0&0&{\color{red}1}&0\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix}
+\begin{bmatrix}0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&{\color{red}1}\\0&0&0&0&0\\0&0&0&0&0\end{bmatrix}
\end{multline}This general strategy of enumerating choices and then selecting among them doesn’t scale all that well because as $N$ gets large $A$ grows pretty quickly. In fact, $A$ will be $N^2 \times \tfrac{1}{6}(N-1)(2N^2+5N+6)$. However, in these cases we can use a heuristic such as column generation to quickly obtain approximate answers that can be iteratively and efficiently refined.

I solved the above integer linear program using Julia with JuMP and the Mosek solver. Read ahead for the results!

And here is the tl;dr, just the solutions!
[Show Solution]

Timing a stoplight just right

This Riddler is about how to perfectly time a stoplight, something we’ve all had to deal with!

You are driving your car on a perfectly flat, straight road. You are the only one on the road and you can see anything ahead of you perfectly. At time t=0, you are at Point A, cruising along at a speed of 100 kilometers per hour, which is the speed limit for the whole road. You want to reach Point C, exactly 4 kilometers ahead, in the shortest time possible. But, at Point B, 2 kilometers ahead of you, there is a traffic light.

At time t=0, the light is green, but you don’t know how long it has been green. You do know that at the beginning of each second, there is a 1 percent chance that the light will turn yellow. Once it turns yellow, it remains yellow for 5 seconds and then turns red for 20 seconds. Your car can accelerate or decelerate at a maximum rate of 2 meters per second-squared. You must always drive at or below the speed limit. You can pass through the intersection when the traffic light is yellow, but not when it is red.

What is the best strategy to reach your destination as soon as possible?

Here is my solution:
[Show Solution]

Let’s give names to the variables so that the unit conversions don’t become cumbersome. Define:
\begin{align}
x_\text{light} &= \text{position of the traffic light (2 km)} \\
x_\text{end} &= \text{position of the finish line (4 km)} \\
v_\text{max} &= \text{speed limit (100 km/h)} \\
a_\text{max} &= \text{maximum acceleration (2 m/sec}^2\text{)} \\
\tau_y &= \text{time that a yellow light lasts (5 sec)} \\
\tau_r &= \text{time that a red light lasts (20 sec)}
\end{align}

It is implicitly assumed in the problem that we can never run a red light. So even though we don’t know when the light will turn yellow, we can’t take any chances; we must make sure that even the unlikeliest most ill-timed yellow light can never cause us to run a red.

There are many different approaches to solving this problem but since I’m currently teaching a class on optimization modeling, I figured I would go that route. The idea behind a modeling approach is to imagine that the position, velocity, and acceleration of the car are decision variables that must satisfy several constraints.

Basic variables and constraints

Let’s discretize the time interval and suppose we make decisions at every $\Delta t$ seconds. So if $\Delta t = 1$, we make decisions every second. Then, define the following decision variables:
\begin{align}
x_t & \ge 0 && \text{position at time }t\\
0 \le v_t &\le v_\text{max} && \text{speed at time }t\\
-a_\text{max} \le a_t &\le a_\text{max} && \text{acceleration at time }t
\end{align}where the time indices range over $t=0,1,\dots,T$ where $T$ is sufficiently large. Of course, we cannot choose positions and velocities arbitrarily; they must be consistent. By approximating the acceleration as constant over each $\Delta t$ step, we may write that:
\begin{align}
v_{t+1} &= v_t + a_t \Delta t && \text{(change in speed with constant accel.)}\\
x_{t+1} &= x_t + v_t\Delta t + \tfrac{1}{2}a_t \Delta t^2 && \text{(change in pos. with constant accel.)}
\end{align}We also have initial conditions for position and speed:
\[
x_0 = 0 \qquad\text{and}\qquad v_0 = v_\text{max}
\]Note that we constrained the speed to be nonnegative. It turns out this doesn’t matter and there is nothing to be gained by allowing backward movement. We get the same answer either way!

Contingency planning

Since we can never risk running a red light, we must ensure that at every instant, one of two things must happen. Either:

We can reach $x_\text{light}$ within $\tau_y$ sec. (we run the yellow) or:
We never pass $x_\text{light}$ within $\tau_y+\tau_r$ sec. (we can avoid running a red).

These constraints can be encoded algebraically. For each $t$, We’ll define a new set of positions, velocities, and accelerations that last for $\tau_y+\tau_r$ seconds. These are contingency variables that allow us to ensure nothing can go wrong if the light turns yellow at time $t$. So we define the variables:
\begin{align}
\hat x_{t,k} & \ge 0 && \text{(contingency position)}\\
0 \le \hat v_{t,k} &\le v_\text{max} && \text{(contingency velocity)}\\
-a_\text{max} \le \hat a_{t,k} &\le a_\text{max} && \text{(contingency acceleration)}
\end{align}These new contingency variables are defined over $0\le t \le T$ and $0 \le k \le \tau_y+\tau_r$. Just like the ordinary position and velocity, they must satisfy contingency constraints:
\begin{align}
\hat v_{t,k+1} &= \hat v_{t,k} + \hat a_{t,k} \Delta t && \text{(change in speed with const. accel.)}\\
\hat x_{t,k+1} &= \hat x_{t,k} + \hat v_{t,k} \Delta t + \tfrac{1}{2}\hat a_{t,k} \Delta t^2 && \text{(change in pos. with const. accel.)}
\end{align}We also have initial conditions for position and speed:
\[
\hat x_{t,0} = x_t \qquad\text{and}\qquad \hat v_{t,0} = v_t
\]Now that we have the contingency trajectories set up, we must encode our binary constraint. Either we can run the yellow, or we are safe from running the red. Let’s define binary variables:
\begin{align}
z^y_t \in \{0,1\} & \qquad \text{(can we run the yellow if light turns yellow at }t\text{?)}\\
z^r_t \in \{0,1\} & \qquad \text{(are we safe from running the red if light turns yellow at }t\text{?)}
\end{align}The constraints that we can run the yellow and not run a red are:
\begin{align}
\hat x_{t,\tau_y} &\ge x_\text{light} z^y_t &&\quad\text{for all }t\text{, and}\\
\hat x_{t,\tau_y+\tau_r} &\le x_\text{light} + (1-z^r_t)x_\text{end} &&\quad\text{for all }t\text{, respectively.}
\end{align}together with the logic constraint $z^y_t + z^r_t \ge 1$ for all $t$, which ensures that we can run the yellow OR we are safe from running the red. And that’s it!

Objective function

The work above describes the feasible set of position, velocity, and acceleration profiles that guarantee we never break the law. Among these profiles, we should choose the one that gets us to the finish line as fast as possible on average. I struggled to find a simple objective function that captures this notion precisely. As a compromise, I simply minimized the duration of the trip assuming the light never turns red. This should be a decent approximation to the true solution because yellow lights are relatively rare. Most of the time when the light turns yellow, it doesn’t actually affect us; it’s only a concern if the light turns yellow when we are close to it, which won’t happen all that often.

To characterize the finish time, I defined the binary variables:
\[
z^e_t \in \{0,1\} \qquad \text{(are we finished at time }t\text{?)}
\]together with the constraint:
\[
x_t \ge x_\text{end} z^e_t\qquad\text{for all }t.
\]Now, our binary variable will be $1$ whenever we have crossed the finish line. So by maximizing $\sum_{t=0}^T z^e_t$, we get to the end as soon as possible.

Solution

From a modeling standpoint, all our constraints are linear, which is nice. We do have some binary variables, which makes this a mixed-integer linear program (MILP). Nevertheless, the problem is relatively easy to solve using a standard solver. For this problem, I coded the model in Julia together with the JuMP package and the Gurobi solver. Some details:

Solving with $\Delta t = 1$. After presolve, the problem had: 3492 constraints, 6216 continuous variables, and 108 binary variables. Solving took about 0.3 seconds on a laptop.
Solving with $\Delta t = 0.5$. After presolve, the problem had: 15499 constraints, 29602 continuous variables, and 216 binary variables. Solving took about 1.3 seconds on a laptop.
Solving with $\Delta t = 0.25$. After presolve, the problem had: 63925 constraints, 125250 continuous variables, and 436 binary variables. Solving took about 12 seconds on a laptop.

No matter the time discretization, all solutions looked the same. Here is a plot of the solution for $\Delta t = 0.25$, zoomed in near the point where the car approaches the traffic light.

The plot looks constant outside this range. So, the optimal thing to do is maintain a speed of 100 km/h for 65 seconds, then decelerate as hard as possible for 2.5 seconds, then accelerate as hard as possible for 2.5 seconds, and then continue on at 100 km/h. The reason for the dip in velocity is that we would otherwise be unable to stop in time if the light were to turn yellow, so we must be risk-averse. After 67.5 seconds, we can safely accelerate, knowing that if the light turns yellow we will make it within the next 5 seconds.

If you’re interested in seeing my code in full detail, you can view/download my notebook here.

Baby poker

Another game theory problem from the Riddler. This game is a simplified version of poker, but captures some interesting behaviors!

Baby poker is played by two players, each holding a single die in a cup. The game starts with each player anteing \$1. Then both shake their die, roll it, and look at their own die only. Player A can then either “call,” in which case both dice are shown and the player with the higher number wins the \$2 on the table, or Player A can “raise,” betting one more dollar. If A raises, then B has the option to either “call” by matching A’s second dollar, after which the higher number wins the \$4 on the table, or B can “fold,” in which case A wins but B is out only his original \$1. No other plays are made, and if the dice match, a called pot is split equally.

What is the optimal strategy for each player? Under those strategies, how much is a game of baby poker worth to Player A? In other words, how much should A pay B beforehand to make it a fair game?

If you’re interested in the derivation (and maybe learning about some game theory along the way), you can read my full solution here:
[Show Solution]

Our first task will be to figure out the payout of the game as a function of the players’ strategies. For each possible dice roll, the Player A must decide whether to call ($0$) or raise ($1$). Similarly, Player B must decide whether to call ($0$) or fold ($1$). An example strategy for Player A is:
\[
p = \begin{bmatrix} 0 & 0 & 0 & 0 & 1 & 1 \end{bmatrix}
\]This corresponds to Player A calling if they roll $\{1,2,3,4\}$ and raising if they roll $\{5,6\}.$ This is an example of a pure strategy because it completely determines what the player does as a function of the information available. We can encode each pure strategy as a six-digit binary number. There are $2^6 = 64$ possible pure strategies for each player and we can number them from $0$ to $63$ based on the decimal representation. For example, strategy $44$ corresponds to $101100$ in binary, i.e. strategy $p = \begin{bmatrix}1 & 0 & 1 & 1 & 0 & 0\end{bmatrix}$. We can also consider mixed strategies, where actions are probabilistic. An example of a mixed strategy for Player A is:
\[
p = \begin{bmatrix} 0.3 & 0 & 1 & 1 & 0 & 0 \end{bmatrix}
\]This strategy is identical to the previous example, except if Player A rolls a $1$, they will raise $30$% of the time and call $70$% of the time. Mixed strategies allows players to play randomly if they so desire! Pure strategies are special cases of mixed strategies in which the probabilities are either $0$ or $1$ (deterministic).

Let’s call $0 \le p_i \le 1$ the strategy of Player A if they roll $i \in \{1,\dots,6\}.$ Likewise, let’s call $0 \le q_j \le 1$ the strategy of Player B if they roll $j \in \{1,\dots,6\}.$ Define the $6\times 6$ matrix $E$ as:
\[
E_{ij} = \begin{cases}
1 &\text{if }i > j \\
0 &\text{if }i = j \\
-1 &\text{if }i < j \end{cases} \]This matrix tells us how much Player A wins if they call with an $i$ and the opponent rolls $j.$ Player A raises with probability $p_i$ and calls with probability $1-p_i.$ Similarly, Player B folds with probability $q_j$ and calls with probability $1-q_j.$ Since each pair $(i,j)$ is equally likely, the expected winnings for Player A are: \[ W = \frac{1}{36} \sum_{i=1}^6 \sum_{j=1}^6 \bigl( (1-p_i)E_{ij} + p_i\left( 2(1-q_j) E_{ij} + q_j \right) \bigr) \]The last term looks the way it does because if Player B folds, Player A wins $+1$ with certainty. If Player B calls, then we use the $E_{ij}$ matrix to determine the payoff, but with a factor of 2 because the stakes were doubled.

Pure strategy equilibria

Define $P_{ab}$ to be the winnings of Player A from the $W$ formula above when players adopt pure strategies $a,b \in \{0,1,\dots,63\}$ respectively. This yields the $64 \times 64$ matrix $P$ corresponding to the winnings for all possible combinations of pure strategies. Here is what $P$ looks like:

Player A is trying to maximize $P_{ab}$ while Player B is trying to minimize it (it’s a zero-sum game, so minimizing the winnings of one player is equivalent to maximizing the losses of the other player). Put another way, Player A selects a row of the matrix $P$ (corresponding to a choice among 64 pure strategies), while Player B selects a column (again, corresponding to a choice among pure strategies). The number in the chosen cell determines how much Player A will win. A Nash Equilibrium is a cell $P_{ab}$ that satisfies
\[
P_{\bar a b} \le P_{ab} \le P_{a \bar b}\qquad \text{for all }\bar a,\bar b \in \{1,\dots,64\}
\]In other words, $P_{ab}$ is the largest element in its column (so Player A has no incentive to pick a different row) and it’s also the smallest element in its row (so Player B has no incentive to pick a different column).

A natural question is to ask is: is there an optimal pair of pure strategies? This amounts to searching the above $P$ matrix for an element that is simultaneously the largest in its column and smallest in its row. A quick search reveals that no such element exists. In other words, there are no pure equilibria for this game!

Mixed strategy equilibria

Every mixed strategy can be expressed as a convex combination of pure strategies. For example:
\[
\begin{bmatrix}1\\\tfrac{1}{3}\\0\\0\\\tfrac{3}{4}\\1\end{bmatrix} =
\frac{1}{3}\begin{bmatrix}1\\1\\0\\0\\1\\1\end{bmatrix} +
\frac{5}{12}\begin{bmatrix}1\\0\\0\\0\\1\\1\end{bmatrix} +
\frac{1}{4}\begin{bmatrix}1\\0\\0\\0\\0\\1\end{bmatrix}
\]This follows from the fact that each mixed strategies lies in the convex hull of all pure strategies. Also, since the winnings $W$ are bilinear in $p$ and $q$, the winnings for a convex combination of strategies is the convex combination of the corresponding winnings.

Now let $u \in [0,1]^{64}$ and $v \in [0,1]^{64}$ be the coefficients of the convex combinations of the pure strategies for Players 1 and 2 respectively. The winnings of Player A can be written compactly as $u^\textsf{T} P v$. Player A would like to select $u$ in order to maximize this quantity, while Player B would like to select $v$ in order to minimize it. All this is subject to the constraint that $0 \le u_i \le 1$ for all $i$ and $u_1 + \dots + u_{64} = 1$, and similarly for $v$. Note that in the special case where $u$ and $v$ contain a single “1” and the rest is “0”, we recover the case of seeking pure strategies.

Nash famously proved that every game with finitely many players and actions has a mixed strategy equilibrium. So all that’s left to do is find it! The set of Nash equilibria can be expressed as the optimal set of a linear program. In particular, the optimal strategies $u,v$ are solutions to the dual linear programs:
\[
\begin{aligned}
\underset{u,\mu}{\max}\quad& \mu \\
\text{s.t.}\quad& u\ge 0,\quad \mathbf{1}^\textsf{T} u = 1 \\
& P^\textsf{T} u \ge \mu \mathbf{1}
\end{aligned}
\qquad
\begin{aligned}
= \qquad
\underset{v,t}{\min}\quad& t \\
\text{s.t.}\quad& v\ge 0,\quad \mathbf{1}^\textsf{T} v = 1 \\
& P v \le t \mathbf{1}
\end{aligned}
\]The first program (Player A) has a unique solution, which is the mixture:
\[
p_\text{opt} =
\frac{1}{3}\begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 1 \\ 1 \end{bmatrix} +
\frac{2}{3}\begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \\ 1 \\ 1 \end{bmatrix} =
\begin{bmatrix} \tfrac{2}{3} \\ 0 \\ 0 \\ 0 \\ 1 \\ 1 \end{bmatrix}
\]The second program (Player B) the solution is not unique, and is rather complicated. As with any linear program, the solution set is a polytope. In this case, the polytope has 15 vertices and lives in a 7-dimensional subspace of $\mathbb{R}^{64}$ (neat, huh?). Projected back down to $\mathbb{R}^6$ where $q$ lives, every optimal solution occupies a 3-dimensional subspace. Specifically,
\[
q_\text{opt} \sim
\begin{bmatrix} 1 \\ x \\ y \\ z \\ 0 \\ 0 \end{bmatrix}
\]where $0 \le x,y \le 1$ and $0 \le z \le \tfrac{2}{3}$ and $x+y+z =\tfrac{4}{3}$.

Each optimal solution has the same optimal value (i.e. the optimal value of $\mu$ or $t$ in the linear programs for $u$ and $v$). The optimal value is $\tfrac{5}{54} \approx 0.0926$. So on average, Player A comes out ahead to the tune of about 9.26 cents when both players are playing optimally.

This alternate solution was proposed by a commenter named Chris. Same answer, but a simpler argument!
[Show Solution]

The solution proceeds much like the previous solution, but this time, we rewrite the cost as a quadratic form in the vectors $p$ and $q$:
\begin{align}
W &= \frac{1}{36} \sum_{i=1}^6 \sum_{j=1}^6 \bigl( (1-p_i)E_{ij} + p_i\left( 2(1-q_j) E_{ij} + q_j \right) \bigr) \\
&= \frac{1}{36} \left( p^\textsf{T}(E \mathbf{1}) + p^\textsf{T}(\textbf{1}\textbf{1}^\textsf{T}-2E)q \right) \\
&= p^\textsf{T} b + p^\textsf{T} A q
\end{align}where $b = \tfrac{1}{36}E\mathbf{1}$, and $A = \tfrac{1}{36}\left(\textbf{1}\textbf{1}^\textsf{T}-2E\right)$, and $\textbf{1}$ is the vector of all ones. We can now directly write linear programs for $p$ and $q$. The first player is trying to maximize $W$ under the assumption that the second player is trying to minimize $W$. In other words, we should solve:
\begin{align}
&\phantom{=}\underset{\textbf{0} \le p \le \textbf{1}}{\max}\quad\underset{\textbf{0} \le q \le \textbf{1}}{\min}\quad p^\textsf{T} b + p^\textsf{T} A q\\
&= \underset{\textbf{0} \le p \le \textbf{1}}{\max} \left( p^\textsf{T} b + \sum_{j=1}^6 \underset{0\le q_j \le 1}{\min} q_j (A^\textsf{T}p)_j \right) \\
&= \underset{\textbf{0} \le p \le \textbf{1}}{\max} \left( p^\textsf{T} b + \sum_{j=1}^6 \min\left( 0, (A^\textsf{T}p)_j \right) \right) \\
&= \left\{\begin{aligned}
\underset{p,t}{\text{maximize}}&\quad p^\textsf{T} b + t^\textsf{T} \textbf{1} \\
\text{subject to:}&\quad \textbf{0} \le p \le \textbf{1} \\
& t \le \textbf{0},\quad t \le A^\textsf{T} p
\end{aligned}\right.
\end{align}Likewise, the second player is trying to minimize $W$ under the assumption that the first player is trying to maximize $W$. We can proceed in a similar manner to above and we obtain:
\begin{align}
&\phantom{=}\underset{\textbf{0} \le q \le \textbf{1}}{\min}\quad\underset{\textbf{0} \le p \le \textbf{1}}{\max}\quad p^\textsf{T} b + p^\textsf{T} A q\\
&= \underset{\textbf{0} \le q \le \textbf{1}}{\min} \left( \sum_{i=1}^6 \max_{0\le p_i \le 1} p_i(Aq+b)_i \right) \\
&= \underset{\textbf{0} \le q \le \textbf{1}}{\min} \left( \sum_{i=1}^6 \max\left( 0, (Aq+b)_i \right) \right) \\
&= \left\{\begin{aligned}
\underset{q,\mu}{\text{minimize}}&\quad \mu^\textsf{T} \textbf{1} \\
\text{subject to:}&\quad \textbf{0} \le q \le \textbf{1} \\
& \mu \ge \textbf{0},\quad \mu \ge A q + b
\end{aligned}\right.
\end{align}Solving these two linear programs yields the same optimal $p$ and $q$ vectors as in the previous approach, but this time our linear programs have far fewer variables! The objective values of both LPs match as well and are both equal to $\tfrac{5}{54} \approx 0.926$, as before. In fact, we can check algebraically that the two LPs are dual to one another. Since both are clearly feasible (easy to check that the feasible set is compact and nonempty), strong duality holds and the two LPs must have the same optimal value. This is in essence Nash’s celebrated result; that the minimax and maximin formulations are equivalent!

If you’d just like to know the answer along with a brief explanation, here is the tl;dr version:
[Show Solution]

Rig the election with math!

This Riddler problem is about gerrymandering. How to redraw borders to sway a vote one way or another.

Below is a rough approximation of Colorado’s voter preferences, based on county-level results from the 2012 presidential election, in a 14-by-10 grid. Colorado has seven districts, so each would have 20 voters in this model. In each district, the party with the most votes will win. The districts must be non-overlapping and contiguous (that is, each square in a district must share an edge with at least one other square in the district). What is the most districts that the Red Party could win? What about the Blue Party? (Assume ties within a district are considered wins for the party of your choice.)

Two boards are provided, a test 5×5 grid and the larger 14×10 grid:

Here is a first solution, using some simple logic and intuition:
[Show Solution]

And here is a computational approach, using integer programming:
[Show Solution]

One way to solve this problem numerically is to formulate it as an integer programming problem. The benefit is that it will always find an optimal solution eventually, but the downside is that it may take awhile. My code solved the 5-by-5 case in a matter of seconds, but the 14-by-10 case computed for over 12 hours with no solution, so I gave up. Of course, it may be possible to solve the larger case using this method, but it would require either a more efficient model or a faster computer!

The idea is to think of the grid cells as nodes of a directed graph. We place an edge between two nodes (one in each direction) if the cells are adjacent. Number the cells $1,2,\dots,C$ and number the edges $1,2,\dots,E$. Define the adjacency matrix $A \in \{0,1\}^{C\times C}$, the incidence matrix $B\in \{0,1\}^{C\times E}$, and the vector $v\in\{0,1\}^C$ as follows:
\begin{align}
A_{ij} &= \begin{cases}
1&\text{if there is an edge from cell }i\text{ to cell }j\\
0&\text{otherwise}
\end{cases}\\[2mm]
B_{ij} &= \begin{cases}
1&\text{if edge }j\text{ starts at cell }i\\
-1&\text{if edge }j\text{ ends at cell }i\\
0&\text{otherwise}
\end{cases}\\[2mm]
v_i &= \begin{cases}
1&\text{if cell }i\text{ is blue}\\
0&\text{if cell }i\text{ is red}
\end{cases}
\end{align}Suppose our districts are numbered $1,2,\dots,D$. The first set of variables in our model are $x\in\{0,1\}^{C\times D}$ and $w\in\{0,1\}^D$ and they represent:
\begin{align}
x_{ij} &= \begin{cases}
1&\text{if cell }i\text{ belongs to district }j\\
0&\text{otherwise}
\end{cases}\\[2mm]
w_j &= \begin{cases}
1&\text{if district }j\text{ is a winner}\\
0&\text{otherwise}
\end{cases}
\end{align}We then have several intuitive constraints which we can encode algebraically as linear constraint for our model:

Each cell belongs to exactly one district:
\[\sum_{j=1}^D x_{ij} = 1\quad\text{ for }i=1,\dots,C\]
Each district contains exactly $S$ cells:
\[\sum_{i=1}^C x_{ij} = S\quad\text{ for }j=1,\dots,D\]
A district wins if and only if there are at least $W$ votes:
\[W w_j \le \sum_{i=1}^C x_{ij} v_i \le (S+1)w_j + (W-1)\quad\text{ for }j=1,\dots,D\]

Unfortunately, these constraints alone don’t force the districts to be contiguous. For this, we must make use of the graph structure. The idea is to set it up as a flow problem. Imagine two new nodes, a universal source and a universal sink. There are edges from the source to every cell, and edges from every cell to the sink. Then imagine a flow that propagates along the edges. The source provides $S$ units of flow (at most one unit of flow per source edge), the sink must accept $S$ units of flow (and it must all occur on a single sink edge), and each cell in the graph conserves flow. These combined constraints ensure that selected cells are path-connected in the graph, which means they are contiguous in the grid. The variables we’ll use are $f \in \mathbb{Z}^{E\times D}$ and $f^\text{out} \in \{0,1\}^{C\times D}$, defined as:
\begin{align}
f_{ij} &= \begin{cases}
1&\text{if edge }i\text{ is used in the flow for district }j\\
0&\text{otherwise}
\end{cases}\\[2mm]
f^\text{out}_{ij} &= \begin{cases}
1&\text{if cell }i\text{ contains the flow to the sink for district }j\\
0&\text{otherwise}
\end{cases}
\end{align}We then have the constraints:

Each edge has a flow capacity of $S$:
\[
0 \le f_{kj} \le S\quad\text{ for }k=1,\dots,E\text{ and }j=1,\dots,D
\]
At most one sink edge is used to return the flow:
\[
\sum_{i=1}^C f^\text{out}_{ij} = 1\quad\text{ for }j=1,\dots,D
\]
Flow is conserved at each cell (including flow involving source and sink)
\[
\sum_{k=1}^E B_{ik}f_{kj} + S f_{ij} = x_{ij}\quad\text{ for }i=1,\dots,C\text{ and }j=1,\dots,D
\]
If an edge is used, the associated cells must be selected
\[
x_{ij} \ge \alpha \sum_{k=1}^E |B_{ik}| f_{kj}\quad\text{ for }i=1,\dots,C\text{ and }j=1,\dots,D
\]

In the final constraint, $\alpha$ is a constant that must be chosen sufficiently small to be less than 1 for any amount of flow. Since there are at most 6 active edge (4 internal, 1 source, 1 sink), it suffices to choose $\alpha < \frac{1}{5S+1}$.

Finally, our objective function is to minimize (or maximize) the total number of winning districts. In other words, $w_1+\dots+w_D$. That’s it!

I used IJulia, JuMP, and Gurobi to solve the integer programs, and on my laptop it took on the order of seconds for the 5×5 case. Unfortunately, the 14×10 case ran for over 12 hours without terminating. So it seems my approach might not be viable for large versions of this problem. The solution my code found for the 5×5 was:
If you’re interested in seeing my code in full detail, you can view/download the notebook here.

Non-intersecting chessboard paths

This Riddler classic puzzle is about finding non-intersecting paths on a chessboard:

First, how long is the longest path a knight can travel on a standard 8-by-8 board without letting the path intersect itself?

Second, there are unorthodox chess pieces that don’t exist in the standard game, which are known as fairy chess pieces. What are the longest nonintersecting paths that can be taken by the camel (which moves like a knight, except 3 squares by 1 square), the zebra (3 by 2), and the giraffe (4 by 1)?

This is a very challenging problem, and there doesn’t appear to be any way to solve it via some clever observation or simplification. Of course, we can try to come up with ever longer tours by hand, but we’ll never know for sure that we have found the longest one.

Much like the recent Pokemon Go problem, we must resort to computational means to obtain a solution. In this case, however, the problem is “small enough” that we can find exact solutions!

Here are some optimal tours:
[Show Solution]

If you’re interested in the details of how I found my solutions, read on:
[Show Solution]

Even though we’re resorting to a numerical solution, we must still be careful about how we proceed. There is an astronomical number of possible paths on a chessboard, so simply enumerating and checking each path is out of the question! Instead, we’ll think about the problem as a visiting nodes on a graph. Here, each node is a cell and edges connect pairs of cells between which a knight can move.

Let’s define the following matrices:

Adjacency matrix: Suppose there are $N$ total cells and we enumerate them $1,2,\dots,N$. $A \in \{0,1\}^{N\times N}$ encodes which cells are reachable from which other cells. $A_{ij} = 1$ if you can move from cell $i$ to cell $j$.
Incidence matrix: Suppose there are $M$ total edges (nonzero entries in $A$) and we enumerate them $1,2,\dots,M$. $B \in \{0,1\}^{N\times M}$ encodes which cells are connected by each edge. $B_{ij} = 1$ if edge $j$ involves cell $i$.
Conflict matrix: The matrix $C\in\{0,1\}^{M\times M}$ indicates which edges cross. $C_{ij} = 1$ if edges $i$ and $j$, excluding the case where the edges share one cell in common.

It’s a bit tedious to figure out what $A$, $B$, and $C$ are, but it’s relatively straightforward to do in code by looping through all nodes or edges as necessary. We now define the following variables for our optimization model:
\begin{align*}
x_1,\dots,x_M &: \text{ variables indicating which edges are used} \\
z_1,\dots,z_N &: \text{ variables indicating which cells are used} \\
w_1,\dots,w_N &: \text{ variables indicating endpoint cells}
\end{align*}We then have three main constraints:
\begin{align*}
\sum_{i=1}^N w_i &\le 2 & & \text{(we have at most two endpoints)}\\
B x &= 2z-w & & \text{(two edges per cell, one per endpoint)}\\
C x &\le H(1-x) & & \text{(no self-intersecting paths)}
\end{align*}
In the last constraint, $H$ is a constant that upper-bounds the number of paths that can intersect any given path. Generally, we can set $H$ to be equal to the maximum row sum of $C$. In the case of a knight tour, $H = 9$. This is an example of the M-trick in integer programming.

Of course, our goal here is to maximize $\sum_{i=1}^M x_i$, the total length of the path. As formulated, this is an integer program that is a version of the Traveling Salesman Problem (TSP) with additional constraints. Unfortunately, the formulation above doesn’t always give a valid solution, since it doesn’t explicitly forbid disjoint paths, e.g. paths with separate closed loops. These are known as subtours.

To get rid of subtours, I used a standard heuristic:

Solve the problem as-is
See if the solution has any subtours. If not, we’re done!
If $\{x_{k_1},\dots,x_{k_\ell}\}$ is a subtour, add the constraint: $\sum_{i=1}^\ell x_{k_i} \le \ell-1$ to explicitly forbid this subtour from occurring in the solution.
Re-solve the problem and return to step 2.

This heuristic always eventually yields an optimal solution, but it may take awhile. For this problem, I only had to re-solve about 10 times before finding a subtour-free (and hence optimal) path.

I used IJulia, JuMP, and Gurobi to solve the integer programs, and on my laptop it took about 30 seconds for 6×6, 1 minute for 7×7 and 5-10 minutes for 8×8. These times only account for a single solve. Subtour eliminations required up to 10 solves per case to find an optimal path.

If you’re interested in seeing my code in full detail, you can view/download my notebook here.

I should also point out that this is a well-studied problem. Some references include: the encyclopedia of integer sequences, non-crossing knight tours, and non-intersecting paths.

Pokémon Go Efficiency

This Riddler puzzle is about a topic near and dear to many hearts: Pokémon!

Your neighborhood park is full of Pokéstops — places where you can restock on Pokéballs to, yes, catch more Pokémon! You are at one of them right now and want to visit them all. The Pokéstops are located at points whose (x, y) coordinates are integers on a fixed coordinate system in the park.

For any given pair of Pokéstops in your park, it is possible to walk from one to the other along a path that always goes from one Pokéstop to another Pokéstop adjacent to it. (Two Pokéstops are considered adjacent if they are at points that are exactly 1 unit apart. For example, Pokéstops at (3, 4) and (4, 4) would be considered adjacent.)

You’re an ambitious and efficient Pokémon trainer, who is also a bit of a homebody: You wish to visit each Pokéstop and return to where you started, while traveling the shortest possible total distance. In this open park, it is possible to walk in a straight line from any point to any other point — you’re not confined to the coordinate system’s grid. It turns out that this is a really hard problem, so you seek an approximate solution.

If there are N Pokéstops in total, find the upper and lower bounds on the total length of the optimal walk. (Your objective is to find bounds whose ratio is as close to 1 as possible.)

Advanced extra credit: For solvers who prefer a numerical question with this theme, suppose that the Pokéstops are located at every point with coordinates (x, y), where x and y are relatively prime positive integers less than or equal to 1,000. Find upper and lower bounds for the length of the optimal walk, again seeking bounds whose ratio is as close to 1 as possible.

The problem of visiting a set of locations while minimizing total distance traveled is known as a Traveling Salesman Problem (TSP), and it is indeed a famous and notoriously difficult problem in computer science. That being said, bounding the solution to a particular TSP instance can be easy if we take advantage of its structure.

Here is my solution to the first part:
[Show Solution]

To clarify, we know there are $N$ Pokéstops but not how they are arranged. The question asks us to find the smallest and largest possible value for the minimum tour that visits each Pokéstop and returns to the starting point.

Upper bound

The worst-case arrangement of Pokéstops is to have them lie in a line, which leads to a path length of $2(N-1)$. Proving that this is indeed the worst possible arrangements of stops is a bit more involved, so read on if you’re interested.

Let $T$ be a minimum spanning tree for the set of Pokéstops. Any tree with $N$ nodes has $N-1$ edges, so $T$ has $N-1$ edges. It also follows from the adjacency property of Pokéstops that each edge of $T$ must have length $1$. One way to visit all the nodes is by using a depth-first tree traversal:

Imagine all edges of $T$ are initially colored black.
Pick a node at random to be the current node.
Find a black edge emanating from the current node. Paint that edge blue and follow it to the next node, which becomes the new current node. If there are no more black edges, follow a blue edge instead and paint it red.
Repeat until all edges are painted red.

This process clearly visits all the nodes and visits each edge exactly twice, for a total tour length of $2(N-1)$. This construction works for any arrangement of Pokéstops. Of course this is merely an upper bound; the optimal path may be shorter. By arranging the $N$ stops in a line, we attain our upper bound of $2(N-1)$, so it must be a tight bound.

Lower bound

Since the Pokéstops lie at integer coordinates, the shortest possible distance from one to the next is $1$. So a path visiting $N$ Pokéstops and then returning to the start has length at least $N$. Does there exist an arrangement of Pokéstops that achieves this minimum? If $N$ is even, the answer is yes; arrange the stops in a $2\times N$ grid:
If the number of stops is odd, however, it’s impossible to arrange them in a way that creates a tour of length $N$. To see why, imagine that the stops are on a checkerboard. Each minimum-length move takes you from black square to a white square or vice versa. Therefore, we can’t return to where we started after an odd number of moves. The second-shortest move has length $\sqrt{2}$, so the best we can do is:
Therefore, a lower bound for the length of the optimal walk is given by:

\[
\text{Lower Bound} = \begin{cases}
0 & \text{if }N = 1 \\
N & \text{if }N\text{ is even} \\
N-1+\sqrt{2} & \text{if }N\text{ is odd and }N \ge 3
\end{cases}
\]

Putting it all together

So in the end, the smallest possible ratio of bounds assuming there are at least two Pokéstops is given by:

$\displaystyle
\text{Bound Ratio} = \begin{cases}
2-\tfrac{2}{N} & \text{if }N\text{ is even} \\
2-\tfrac{2\sqrt{2}}{N-1+\sqrt{2}} & \text{if }N\text{ is odd}
\end{cases}
$

The smallest bound ratio that works for all $N$ is $2$, so this is the ratio we should use if we don’t know how many Pokéstops there are.

Here is my solution to the second part:
[Show Solution]

For the second part of the Riddle, we are given a very specific arrangement of Pokéstops: the set of points $1 \le i,j \le M$ such that $\mathrm{gcd}(i,j) = 1$ (greatest common divisor is $1$). The problem asks for $M=1000$, but we’ll keep $M$ as a parameter for now.

Upper bounds

Any tour that visits all the Pokéstop and returns to the starting point gives us an upper bound. With this many nodes, it’s impossible to check every tour, so we’ll compare different heuristics for choosing admissible paths and see how they compare.

Trivial upper bound: An obvious way to hit all the stops is to simply visit every point with integer coordinates in the $M\times M$ square. We know how to do this optimally from the first part of the problem; the optimal tour has length $M^2$ if $M$ is even and $M^2+\sqrt{2}-1$ if $M$ is odd. See the illustration below for example constructions:
So the trivial bound for the case of interest is 1,000,000. We can definitely do better!

Sierpiński bound: One way to find a reasonable tour is to use a space-filling curve, which maps every point in the square to a point on a curve. This gives us a way to order the stops based on where they lie on the curve. We then visit the stops in that order. This heuristic was pioneered by Prof. John Bartholdi III and others, and does very well in many cases.
The Sierpiński curve is often used here, and it’s what I used to find an upper bound. Here are some examples of solutions for different $M$ values.
These paths have lengths 313.9 and 1828.8, respectively. Using this method, our upper bound improves to 815,234 for the $M=1000$ case, and this was computed in 0.35 seconds.

Columnwise bound: The Sierpiński approach takes advantage of the pseudouniform distribution of points, but it doesn’t leverage the fact that points often appear in lines. For example, the rows and columns corresponding to prime coordinates are always full. Another heuristic, suggested by Diarmuid Early in a comment, is to proceed in a lawnmower-like pattern, but two columns at a time. So as we move up or down a pair of columns, we pick the next nearest point from the pair of available points at each step. Here are examples for the same grid sizes.
These paths are slightly shorter than the Sierpiński ones, with lengths of 308.3 and 1721.8, respectively. Using this method, our upper bound improves to 735,394 for the $M=1000$ case, and was computed in 2 seconds.

More on upper bounds: Other notable heuristics for upper-bounding TSP problems include the nearest neighbor algorithm and the Christofides algorithm. These don’t work well for our instance because they require computing all pairwise distances, a prohibitively expensive task when we have roughly 600,000 Pokéstops. The Sierpiński and columnwise algorithms have complexity $\mathcal{O}(n)$, where $n$ is the number of stops, whereas computing pairwise distances is $\mathcal{O}(n^2)$. This is a significant difference because $n$ is so large — if an $\mathcal{O}(n)$ algorithm takes 1 second then $\mathcal{O}(n^2)$ will take roughly 7 days!

Lower bounds

Trivial lower bound: One possibility is to count how many total stops there are, and assume we can find a path consisting entirely of adjacent stops. This isn’t possible of course, but it yields a lower bound. It’s relatively easy to write computer code to count this quantity, or we can use the fact that the density of stops in the $M \times M$ square approaches $6/\pi^2 \approx 0.6079$ as $M\to\infty$ (source). For the case $M=1000$, the lower bound turns out to be 608,383. Lots of room for improvement!

Averaging bound: Another observation, again courtesy of Diarmuid Early, is that each Pokéstop has two edges emanating from it. So if we assume each stop is connected to its two nearest stops, we will get another lower bound on the shortest possible path. Computing this quantity for the case $M=1000$, we obtain a much-improved 643,811.

Held-Karp relaxation: This final bounding method is significantly more advanced. Finding a set of edges (i.e. assigning each edge a label 0 or 1) that yields a valid tour requires two things:

each node has exactly two incident edges labeled 1
there are no disjoint cycles (it’s a connected path)

What makes this problem so hard is the second constraint. There are exponentially many cycles, so it’s not feasible to rule each one out. If we relax the problem by removing the disjoint cycles constraint, and allow the label on each edge to be any number between 0 and 1, the problem turns into a linear program (LP), which can be solved very efficiently. Successive refinements can then be made to improve the solution by ruling out particular cycles. This is known as the Held-Karp LP relaxation. Unfortunately, this still requires computing all pairwise distances, which we already decided was too costly (see above). The LP also has $\mathcal{O}(n^2)$ variables, since we need a variable for each edge. So I made two simplifying assumptions:

Each node may only be connected to one of its $k$ nearest neighbors. This is a relaxation to the GTSP, (traveling salesmsan problem on a graph). This means the number of variables, constraints, and distances we need to compute are now $\mathcal{O}(kn)$.
The final path will be symmetric ($x=y$). This allows us to use half as many variables and constraints and greatly speeds up computation.

To find an appropriate $k$, I solved the relaxed LP with $k$ values of 4, 8, 12, 16, 20. The solutions were 695251, 694956, 694954, 694952, 694952. So it seems there is no use in adding any more edges. For those that are interested, I implemented the solution in Julia using JuMP with Gurobi. On my laptop, the LPs took 30–90 seconds to solve.

Putting everything together

Here is a plot comparing all the bounds above.
So ultimately, our best bound for the case $M=1000$ is

\[
694952 \le L^\star \le 735394
\]

for an approximation ratio of about $1.06$.

Tag: linear programming