LeetCode Sudoku Solver with Scipy Integer Programming

Published on 2023-04-02

There’s a catch though: It’s not available yet.

Preface

I’ve never been an avid user of LeetCode until lately. People occasionally ask me about some of the interesting questions from there, to which I reply with some thoughts, without taking a deep look into what else the website had to offer. Things changed recently, since I’ve started to apply for internship this summer, and the hiring process may involve some pre-interview coding tests.

As I was powering through the list of problems, the 37. Sudoku Solver caught my eyes. While it was labeled as a hard problem, all it takes to solve the problem is a simple backtracking search over the possible combinations.

Interestingly, under the Stochastic search / optimization methods heading of the Wikipedia entry for Sudoku solving algorithms, it’s stated that the problem of solving Sudoku can be reduced to an integer linear programming problem. The Wikipedia entry lead to me implementing it with scipy.optimize.milp. Sadly though, my code ended up not running on LeetCode due to their version of Scipy being slightly out of date¹, but they’re probably gonna update it in the future.

Encoding Sudoku as Integer Programming Problem

Bartlett et al. [1] provides an integer programming formulation for the Sudoku problem. The key step is to map the puzzle on a $9 \times 9$ matrix to a puzzle on a $9 \times 9 \times 9$ matrix called $x$, indexed by $(v, r, c)$ tuples of value, row, and column. The element at $x_{vrc}$ is $1$ iff the number at $(r, c)$ is $v + 1$ on the original puzzle².

This representation for the puzzle is used for other combinatorial optimization methods, such as using SAT solvers, as well. See Sudoku as a SAT Problem [2] for a SAT formulation that encodes the puzzle in the same way.

The integer programming formulation is for the encoding is shown below. Note that the objective function is zero, as the goal is to find any solution satisfying the constraints. The constraints 2-4 are symmetric.

1. The bounds — Every element is either $0$ or $1$.

   $x_{vrc} \in \{ 0, 1 \},\; \forall v, r, c$

2. Every position of the original puzzle is filled.

   $\sum_{v=0}^{8} x_{vrc} = 1,\; \forall c, r \in [0, 8]$

3. Every column has exactly one $v + 1$.

   $\sum_{r=0}^8 x_{vrc} = 1,\; \forall v, c \in [0, 8]$

4. Every row has exactly one $v + 1$.

   $\sum_{c=0}^8 x_{vrc} = 1,\; \forall v, r \in [0, 8]$

5. Every 3-by-3 block has exactly one $v + 1$.

   $\sum_{r=3p}^{3p+2} \sum_{c=3q}^{3q+2} x_{vrc},\; \forall p, q \in [0, 2]$

6. The preset values on the puzzle.

   $x_{vrc} = 1$, for each $v + 1$ preset at position $(r, c)$ on the puzzle.

The formulation maps pretty much directly to the API provided by scipy.optimize.milp. The bounds can be specified using Bounds(lb=np.zeros(729), ub=np.ones(729)). Each constraint can be specified using LinearConstraint(A, lb=1, ub=1), where A is an array of length $729$ where the elements at index $(v,r,c)$ is $1$ iff $x_{vrc}$ appears on the LHS of the constraint.

It’s Not Fast

It’s not fast. Well, actually it’s fast, but only for the hardest puzzles.

I’ve made a quick benchmark comparing the total execution time of my MILP solution against some of the most voted solutions on LeetCode.

• Python Easiest Recursive Solution is the top voted Python solution. It’s implemented using simple backtracking search, with no additional bookkeeping other than the board itself.
• 🐍 97% faster || Clean & Concise || Well-Explained 📌📌 ³ is the second voted Python solution. As with the first solution, it’s also implemented using backtracking search, but it makes use of additional sets to keep track of the constraints.

The set of LeetCode’s official test cases contains only 6 instances of the puzzle, and they are subjectively easy as I was able to hand solve an instance in less than twenty minutes. Therefore, another set of harder instances was brought into place as a fun experiment as well.

The execution times for the LeetCode test cases are shown in this figure. While the Scipy MILP solution runs slightly faster (25%) than the simple DFS solution, the more optimized version of DFS crashes the MILP solution by a 8x speedup.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import ticker

exec_secs = np.load('./results/leetcode_exec_secs.npy').tolist()

plt.style.use('bmh')

fig, ax = plt.subplots()

labels = ['ILP', 'DFS', 'Faster DFS']
ax.set_xticks([1, 2, 3])
ax.set_xticklabels(labels)
ax.set_axis_on()
ax.set_ylabel('Execution Time (secs)')

ax.grid(visible=True, which='major')
ax.violinplot(exec_secs, showmeans=True)

plt.savefig("./leetcode_bench.svg")

The execution times for the harder cases are shown in this figure. I’ve plotted them in logarithmic scale, because otherwise the outlier values (up to 7 minutes!) make the rest of the data points degenerate into tiny horizontal bars at the bottom.

It’s interesting to see that for these hard puzzles, the MILP solution runs notably faster than the DFS solutions. The MILP solution also has less dispersed execution times compared to the others.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import ticker

exec_secs = np.load('./results/hard_exec_secs.npy')
mean_exec_secs = exec_secs.mean(axis=1)
log_exec_secs = np.log10(exec_secs).tolist()

plt.style.use('bmh')

fig, ax = plt.subplots()

labels = ['ILP', 'DFS', 'Faster DFS']
ax.set_xticks([1, 2, 3])
ax.set_xticklabels(labels)
ax.set_axis_on()
ax.set_ylabel('Execution Time (secs, log10 scale)')
ax.yaxis.set_major_formatter(ticker.StrMethodFormatter('$\mathdefault{{10^{{{x:.0f}}}}}$'))

ax.grid(visible=True, which='major')
r = ax.violinplot(log_exec_secs, showmeans=True)

# Fix mean bars for log10 scale
paths = r['cmeans'].get_paths()
for i in range(3):
  vertices = paths[i].vertices
  vertices[:, 1] = np.log10(mean_exec_secs[i])
  paths[i].vertices = vertices

plt.savefig("./hard_bench.svg")

Code Listing

from typing import List
from itertools import product
from scipy.optimize import milp, LinearConstraint, Bounds
import numpy as np

class IlpSolution():
    def solveSudoku(self, board: List[List[str]]) -> None:
        constraints = create_constraints(board)
        milp_output = milp(
            np.zeros(729),
            bounds=Bounds(lb=np.zeros(729), ub=np.ones(729)),
            constraints=constraints,
            # 1 means that the element must be an integer
            integrality=np.full(729, 1),
        )
        return set_board_to_solution(board, milp_output.x)


def set_board_to_solution(board: List[List[str]], solution: np.ndarray):
    for r, c in product(range(9), repeat=2):
        for v in range(9):
            if solution[index(v, r, c)] == 1.0:
                board[r][c] = str(v + 1)
                break


def create_constraints(board: List[List[str]]) -> List[LinearConstraint]:
    rule_constraints = create_rule_constraints()
    instance_constraints = create_instance_constraints(board)
    return rule_constraints + instance_constraints


def create_instance_constraints(board: List[List[str]]) -> List[LinearConstraint]:
    constraints = []
    for r, c in product(range(9), repeat=2):
        ch = board[r][c]
        if ch == '.':
            continue
        v = int(ch) - 1
        A = np.zeros(729)
        A[index(v, r, c)] = 1
        constraints.append(LinearConstraint(A, lb=1, ub=1))

    return constraints


def create_rule_constraints() -> List[LinearConstraint]:
    constraints = []
    for r, c in product(range(9), repeat=2):
        A = np.zeros(729)
        for v in range(9):
            A[index(v, r, c)] = 1
        constraints.append(LinearConstraint(A, lb=1, ub=1))

    for v, r in product(range(9), repeat=2):
        A = np.zeros(729)
        for c in range(9):
            A[index(v, r, c)] = 1
        constraints.append(LinearConstraint(A, lb=1, ub=1))

    for v, c in product(range(9), repeat=2):
        A = np.zeros(729)
        for r in range(9):
            A[index(v, r, c)] = 1
        constraints.append(LinearConstraint(A, lb=1, ub=1))

    for x, y, v in product([0, 3, 6], [0, 3, 6], range(9)):
        A = np.zeros(729)
        for r, c in product(range(x, x + 3), range(y, y + 3)):
            A[index(v, r, c)] = 1
        constraints.append(LinearConstraint(A, lb=1, ub=1))

    return constraints


def index(v: int, r: int, c: int) -> int:
    return v * 81 + r * 9 + c