Sudoku Solver

Prepared by: William Bartel (wbartel@nd.edu, 2024) and Andrew Brooks (abrooks7@nd.edu, 2024). Adapted from material by Prof. Alexander Dowling and Elvis Eugene at the University of Notre Dame

Learning Objectives¶

Formulating a sudoku puzzle as a mixed-integer feasibility problem
Reading a committed dataset of randomly generated puzzles
Visualizing a puzzle’s starting position alongside its solution
Comparing two MILP solvers, and seeing how the model size scales

Imports¶

# This code cell installs packages on Colab

import sys

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper

    helper.easy_install()
else:
    sys.path.insert(0, "../")
    import helper

helper.set_plotting_style()

# Import Libraries
import math
import time

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pyomo.environ as pyo
from matplotlib.patches import Rectangle

# This notebook draws no random numbers: the puzzles are read from committed
# CSV files. Generating those files IS random, and that code is seeded --
# see scripts/generate_sudoku_puzzles.py in the course repository.

Defining Sudoku Puzzles¶

Sudoku puzzles are one of the most widely known logical puzzles. The rules of a sudoku puzzle are simple. For an $n \times n$ puzzle, each row, each column, and each subsquare must have exactly one of each number 1 through $n$ . One can look at the numbers given in the puzzle and use logic to conclude which squares must contain which numbers. An example puzzle is shown below.

A 9 by 9 sudoku puzzle on the left, with about half of its cells blank, and the same puzzle completed on the right, where the filled-in digits are printed in red and the digits that were given are printed in black.

Formulating an Optimization Problem for $n\times n$ Puzzle¶

Sets:

Rows $R=\{1,2, ...,n\}$
Columns $C=\{1,2, ...,n\}$
Values $V=\{1,2, ...,n\}$
Subgrids $K=\{1,2, ...,n\}$

Parameters:

$s_{v,i,j} =$ is value $v$ present in box located at $(i,j)$ from unsolved puzzle for $i\in R, j\in C, v \in V$ . (Binary)
$SG_{i} =$ coordinate pairs for the boxes in each subgrid for $i\in K$

Variables:

$y_{v,i,j} =$ is value $v$ present in box located at $(i,j)$ when creating solved puzzle for $i\in R, j\in C, v \in V$ . (Binary)

Objective:

$\min 0$

Constraints:

\sum_{j\in C}y_{v,i,j}=1, \forall i\in R, v\in V

(1)

\sum_{i\in R}y_{v,i,j}=1, \forall j\in C, v\in V

(2)

\sum_{(i,j)\in SG[\ l \ ]}y_{v,i,j}=1, \forall l\in K, v\in V

(3)

\sum_{v\in V}y_{v,i,j}=1, \forall i\in R, j\in C

(4)

y_{v,i,j} * s_{v,i,j} = s_{v,i,j}, \forall i\in R, j\in C, v\in V

(5)

Any feasible completion is acceptable, so the objective is the constant 0. A puzzle can have no completion, one completion or several; a zero objective does not establish uniqueness. For standard square subgrids, use $n=a^2$ with integer $a$ , and index each $a\times a$ subgrid by its $(r,c)$ pairs.

The constraints are described as follows:

There can only be one instance of each possible value in each row.
There can only be one instance of each possible value in each column.
There can only be one instance of each possible value in each subgrid.
There can only be one value in each square.
Each value given by the unsolved puzzle must be the same in the solved puzzle.

Rows, columns and a subgrid in the nine by nine model.

The indices held fixed and summed in three exactly-one families.

Explore the LP relaxation¶

For the Lecture 3 study activity, replace $y\in\{0,1\}$ with $0\le y\le1$ . Before solving, predict whether the returned point will be integral or fractional. Then solve and check every binary component against an integrality tolerance and every row against a feasibility tolerance.

Reference: BGW (1997), Appendix A, §A.3.2, pp. 758, 760–761, for integer relaxation and bound quality. For the Sudoku formulation, Bynum et al., Pyomo, 3rd ed., §5.4.2, pp. 76–79.

$A one-hot cell and a uniform fractional cell each sum to one.$

# Define class 'Puzzle' to store an unsolved or solved sudoku board
# Stores size of puzzle, grid values, and subgrid coordinates


class Puzzle:
    """One sudoku board: its size, its cells, and its subgrid geometry.

    Arguments:
        board: an n-by-n nested list (or array) of ints, 1..n for a filled
            cell and 0 for a blank one
    """

    def __init__(self, board):

        # n = number of rows = number of columns
        n = len(board[0])
        # side length of one subgrid
        sub = int(math.sqrt(n))

        self.size = n
        self.grid = [list(int(v) for v in row) for row in board]

        # s[(v,r,c)] = 1 if the cell in position (r,c) has value v
        # s[(v,r,c)] = 0 otherwise
        # Note: rows and columns begin at 1 and end at n
        self.params = {
            (v, r, c): int(self.grid[r - 1][c - 1] == v)
            for v in range(1, n + 1)
            for r in range(1, n + 1)
            for c in range(1, n + 1)
        }

        # subgrids[i] = collection of ordered pairs (r,c) such that (r,c) is in
        # the i-th subgrid. Subgrids are sqrt(n) by sqrt(n) and are indexed top
        # to bottom, then left to right.
        self.subgrids = {
            i: [
                (r, c)
                for r in range(1 + sub * ((i - 1) % sub), 1 + sub * ((i - 1) % sub + 1))
                for c in range(
                    1 + sub * ((i - 1) // sub),
                    1 + sub * ((i - 1) // sub + 1),
                )
            ]
            for i in range(1, n + 1)
        }

The Sudoku Datasets¶

The puzzles used below were generated ahead of time and are committed to the course repository as notebooks/data/sudoku_puzzles_4.csv, ..._9.csv, ..._16.csv, ..._25.csv and ..._36.csv. Each file holds 100 puzzles of one size, one puzzle per row, flattened row-major, with 0 for a blank cell.

The generator itself is not in this notebook. It lives in scripts/generate_sudoku_puzzles.py, it is seeded, and re-running it reproduces the committed CSV files byte for byte:

python scripts/generate_sudoku_puzzles.py            # rewrite the CSV files
python scripts/generate_sudoku_puzzles.py --check    # verify, write nothing

How it works, in one paragraph. It starts from a stored solved board for each of $n=4,9,16,25,36$ and scrambles it with moves that preserve the sudoku property: interchanging two rows that pass through the same band of subgrids, interchanging two whole bands of rows, and the two column analogues. Any composition of those maps a solved board to a different solved board. It then blanks out a size-dependent number of cells to make the puzzle.

For instance, swap_cols_within_band interchanges two columns that pass through the same subgrids.

A 9 by 9 solved sudoku with two columns inside the middle band of subgrids highlighted, and the same board on the right after those two columns have been interchanged; the result is still a valid solved sudoku.

Reading the Puzzles¶

load_puzzles reads one CSV file and returns every puzzle in it as an array of shape (number of puzzles, n, n). On Google Colab the relative path below is rewritten to a raw.githubusercontent.com URL when the notebook is published, so the data downloads automatically and nothing else has to change.

def load_puzzles(size):
    """Read every puzzle of one size from the committed dataset.

    Arguments:
        size: puzzle side length n, one of 4, 9, 16, 25, 36

    Returns:
        an integer array of shape (number of puzzles, size, size); 0 is a blank
    """
    table = pd.read_csv(f"https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/data/sudoku_puzzles_{size}.csv", header=None)
    return table.to_numpy(dtype=int).reshape(-1, size, size)


# Read the first 9 by 9 puzzle in the dataset
puzzles_9 = load_puzzles(9)
print(f"{len(puzzles_9)} puzzles of size 9 x 9")
print(puzzles_9[0])

100 puzzles of size 9 x 9
[[3 2 0 0 8 1 0 0 7]
 [1 0 0 7 2 6 3 0 0]
 [6 8 0 3 0 5 9 0 0]
 [2 6 8 4 0 3 1 7 9]
 [4 0 0 0 0 0 0 3 0]
 [5 0 3 6 1 9 0 0 4]
 [8 0 5 2 0 0 7 9 1]
 [9 4 2 0 0 7 8 5 6]
 [7 1 6 0 0 8 0 0 3]]

Sudoku Optimization Solver¶

Below is the optimization formulation to solve a sudoku puzzle. The formulation was described earlier. Following the Pyomo style guide, the model is assembled by a function that takes data and returns a fresh ConcreteModel, and the solve is a separate step.

def build_sudoku_model(puzzle):
    """Build the MILP feasibility model for one sudoku puzzle.

    Arguments:
        puzzle: a Puzzle holding the starting position

    Returns:
        a Pyomo ConcreteModel
    """
    m = pyo.ConcreteModel()

    n = puzzle.size

    # Index sets. All four run 1..n; they are declared separately because they
    # mean different things, and m.pprint() then says which is which.
    m.VALUES = pyo.RangeSet(1, n)
    m.ROWS = pyo.RangeSet(1, n)
    m.COLS = pyo.RangeSet(1, n)
    m.SUBGRIDS = pyo.RangeSet(1, n)

    # Is value v printed in cell (r,c) of the starting puzzle? s_{v,i,j} [-]
    m.given = pyo.Param(
        m.VALUES, m.ROWS, m.COLS, initialize=puzzle.params, domain=pyo.Binary
    )

    # Does cell (r,c) hold value v in the completed puzzle? y_{v,i,j} [-]
    m.y = pyo.Var(m.VALUES, m.ROWS, m.COLS, domain=pyo.Binary, initialize=0)

    # The DECLARATION ORDER of the five constraints below is deliberate: it
    # fixes the row order of the LP file handed to the solver, and CBC's
    # branch-and-bound is sensitive to it. Reordering these (subgrid before
    # cell) leaves the model mathematically identical but made the hardest
    # 36 x 36 puzzle in the dataset take over ten minutes instead of about
    # one. Measured 2026-08-25; the LP file is byte-identical in this order.

    @m.Constraint(m.VALUES, m.ROWS)
    def one_per_row(b, v, r):
        """Value v appears exactly once in row r."""
        return sum(b.y[v, r, c] for c in b.COLS) == 1

    @m.Constraint(m.VALUES, m.COLS)
    def one_per_col(b, v, c):
        """Value v appears exactly once in column c."""
        return sum(b.y[v, r, c] for r in b.ROWS) == 1

    @m.Constraint(m.ROWS, m.COLS)
    def one_value_per_cell(b, r, c):
        """Cell (r,c) holds exactly one value."""
        return sum(b.y[v, r, c] for v in b.VALUES) == 1

    @m.Constraint(m.VALUES, m.ROWS, m.COLS)
    def match_given(b, v, r, c):
        """A value printed in the starting puzzle carries into the solution."""
        return b.y[v, r, c] * b.given[v, r, c] == b.given[v, r, c]

    @m.Constraint(m.VALUES, m.SUBGRIDS)
    def one_per_subgrid(b, v, k):
        """Value v appears exactly once in subgrid k."""
        return sum(b.y[v, r, c] for (r, c) in puzzle.subgrids[k]) == 1

    # A sudoku is a feasibility problem: any completion is as good as any other
    m.obj = pyo.Objective(expr=0)

    return m


def solve_sudoku(puzzle, solver="appsi_highs", tee=False):
    """Build and solve the sudoku MILP.

    Arguments:
        puzzle: a Puzzle holding the starting position
        solver: a SolverFactory name; this notebook uses "appsi_highs" and "cbc"
        tee: show the solver log

    Returns:
        the solved Pyomo ConcreteModel
    """
    m = build_sudoku_model(puzzle)

    opt = pyo.SolverFactory(solver)
    results = opt.solve(m, tee=tee)

    # Never touch pyo.value() before checking the termination condition
    assert pyo.check_optimal_termination(results), (
        f"Solve failed: status={results.solver.status}, "
        f"termination={results.solver.termination_condition}"
    )

    return m


def extract_solution(model, puzzle):
    """Pull the completed grid out of a solved model.

    Separating extraction from plotting is the house convention: the solve is
    expensive, the picture is not, so the numbers come out once and the figure
    can be redrawn as often as you like.

    Arguments:
        model: a solved model from solve_sudoku
        puzzle: the Puzzle it was built from

    Returns:
        an n-by-n nested list of ints, the completed grid
    """
    n = puzzle.size
    grid = [[0] * n for _ in range(n)]
    for v, r, c in model.y:
        # HiGHS can return 0.9999... or -0.0 for a binary; compare with a
        # tolerance, never with == 1
        if pyo.value(model.y[v, r, c]) >= 0.5:
            grid[r - 1][c - 1] = v
    return grid

Visualizing a Puzzle and Its Solution¶

plot_sudoku draws one board. Given a solution as well, it draws the completed board with the starting pieces distinguished from the cells the solver filled in --- which is the interesting part of the picture, because it shows how much of the puzzle the solver actually had to do.

The distinction is carried three ways, not one: the given cells sit on a shaded background, their digits are set in bold, and they are black where the solved digits are vermillion. That redundancy is the course figure rule --- a distinction may be carried by color, but never by color alone, because students print these handouts on black-and-white printers. See figures/README.md.

# Okabe-Ito vermillion, the same color the house style cycles second.
# Color is the third cue here, not the only one: see the docstring.
SOLVED_COLOR = "#D55E00"
GIVEN_SHADE = "0.87"


def plot_sudoku(puzzle, solution=None, ax=None, title=None):
    """Draw a sudoku board, distinguishing the givens from the solved cells.

    Arguments:
        puzzle: a Puzzle holding the starting position
        solution: a completed grid from extract_solution, or None to draw only
            the starting position
        ax: an Axes to draw on, or None to create a new figure
        title: axes title

    Returns:
        the Axes drawn on

    Givens and solved cells are told apart three ways so the figure survives
    grayscale printing: the givens are shaded, bold and black; the solved cells
    are unshaded, regular weight and vermillion.
    """
    n = puzzle.size
    sub = int(math.sqrt(n))

    if ax is None:
        side = 1.0 + 0.32 * n
        _, ax = plt.subplots(figsize=(side, side))

    # Cell (r,c), 1-indexed from the top left, occupies the unit square with
    # lower-left corner (c-1, n-r) so that row 1 is drawn at the top.
    for r in range(1, n + 1):
        for c in range(1, n + 1):
            given = puzzle.grid[r - 1][c - 1]
            if given:
                ax.add_patch(
                    Rectangle(
                        (c - 1, n - r),
                        1,
                        1,
                        facecolor=GIVEN_SHADE,
                        edgecolor="none",
                        zorder=0,
                    )
                )

    fontsize = max(4.5, 95.0 / n)
    for r in range(1, n + 1):
        for c in range(1, n + 1):
            given = puzzle.grid[r - 1][c - 1]
            value = given if given else (solution[r - 1][c - 1] if solution else 0)
            if not value:
                continue
            ax.text(
                c - 0.5,
                n - r + 0.5,
                str(value),
                ha="center",
                va="center",
                fontsize=fontsize,
                fontweight="bold" if given else "normal",
                color="black" if given else SOLVED_COLOR,
                zorder=2,
            )

    # Grid lines: heavy on the subgrid boundaries, light inside them
    for i in range(n + 1):
        lw = 2.0 if i % sub == 0 else 0.6
        ax.axhline(i, color="black", linewidth=lw, zorder=1)
        ax.axvline(i, color="black", linewidth=lw, zorder=1)

    ax.set_xlim(0, n)
    ax.set_ylim(0, n)
    ax.set_aspect("equal")
    ax.axis("off")
    if title:
        ax.set_title(title, fontsize=12)
    return ax

A Worked Example¶

The 9 by 9 board below is the one drawn at the top of this notebook. We keep it as a named constant rather than reading it from the dataset precisely because it is a worked example: it has to stay the same board as the figure. Every other puzzle in this notebook comes from the CSV files.

# The 9 by 9 puzzle drawn in the figure at the top of this notebook.
# A worked example, deliberately hardcoded so it cannot drift from that figure.
EXAMPLE_PUZZLE_9 = [
    [5, 3, 0, 0, 7, 0, 0, 0, 0],
    [6, 0, 0, 1, 9, 5, 0, 0, 0],
    [0, 9, 8, 0, 0, 0, 0, 6, 0],
    [8, 0, 0, 0, 6, 0, 0, 0, 3],
    [4, 0, 0, 8, 0, 3, 0, 0, 1],
    [7, 0, 0, 0, 2, 0, 0, 0, 6],
    [0, 6, 0, 0, 0, 0, 2, 8, 0],
    [0, 0, 0, 4, 1, 9, 0, 0, 5],
    [0, 0, 0, 0, 8, 0, 0, 7, 9],
]

example = Puzzle(EXAMPLE_PUZZLE_9)
model = solve_sudoku(example)
solution = extract_solution(model, example)

fig, axes = plt.subplots(1, 2, figsize=(9, 4.6))
plot_sudoku(example, ax=axes[0], title="Starting position")
plot_sudoku(example, solution=solution, ax=axes[1], title="Solved")
fig.suptitle("9 x 9 sudoku: shaded bold digits are given, plain digits are solved")
plt.show()

# The same two pictures for a 4 by 4 puzzle read from the dataset

small = Puzzle(load_puzzles(4)[0])
small_solution = extract_solution(solve_sudoku(small), small)

fig, axes = plt.subplots(1, 2, figsize=(7, 3.8))
plot_sudoku(small, ax=axes[0], title="Starting position")
plot_sudoku(small, solution=small_solution, ax=axes[1], title="Solved")
plt.show()

Count the formulation actually declared¶

For $n=a^2$ , the four families contain $n^2$ rows each: cell $(r,c)$ , value-row $(v,r)$ , value-column $(v,c)$ , and value-subgrid $(v,k)$ . The model declares $n^3$ binaries in $[0,1]$ and $4n^2$ exactly-one equalities. These rows are not all independent: summing cell rows in one row of the grid gives the same total as summing its value-row rows. Do not call $n^3-4n^2$ a continuous degree-of-freedom count for this integer feasibility problem.

Given value $v$ in cell $(r,c)$ can enter as a fixing $y_{v,r,c}=1$ , as in Bynum et al. (2021), p. 79. The other values in that cell then become zero through nonnegativity and the cell equation. With $G$ givens, the fixed-variable implementation retains $4n^2$ declared constraint rows and fixes $G$ variables.

The canonical builder in this notebook intentionally retains its existing different implementation: a binary parameter $s_{v,r,c}$ and $n^3$ rows $s_{v,r,c}y_{v,r,c}=s_{v,r,c}$ . Exactly $G$ rows impose a given; the remaining $n^3-G$ are $0=0$ . These equations are linear because $s$ is data, not a variable. For $n=9$ , the code therefore declares 729 binaries and 1,053 equalities, including the 324 exactly-one rows. For the worked board $G=30$ , so 699 matching rows are trivial. Presolve may remove those rows and other redundancies; declared counts and solver-reported counts need not coincide.

Both representations have the same feasible $y$ values; they are not the same declarations. The diagnostic below builds the fixed-variable alternative by deleting match_given and applying .fix(1) at the givens, keeping all other rows unchanged.

Subgrid labels are conventions. The handout and Pyomo book number left-to-right, then down; Puzzle.subgrids in this notebook numbers top-to-bottom, then right. For $n=9$ , the handout’s top-middle $S_2$ is puzzle.subgrids[4]; the notebook’s block 2 is middle-left. The collection of nine blocks and its feasible set are identical. The established code order is retained for reproducibility of the benchmark. Each constraint still sums over the actual $(r,c)$ pairs in its block.

Equivalent fixing and data-matching representations for a given Sudoku value.

In the handout, $K=\{1,\ldots,9\}$ indexes blocks, $S_k\subset R\times C$ is the coordinate-pair set for block $k$ , and $\mathcal S=\{S_k:k\in K\}$ is the collection of those sets. The quantifiers for a block rule are $v\in V$ and $k\in K$ ; one does not index $S_k$ by an element of the collection that is itself a set. This matches the notebook’s block-index/coordinate-map construction despite the different label order.

# Compare equivalent ways to enter the same 30 givens.
matching = build_sudoku_model(example)
fixed_givens = matching.clone()
fixed_givens.del_component(fixed_givens.match_given)
for (v, r, c), given in example.params.items():
    if given:
        fixed_givens.y[v, r, c].fix(1)
for label, candidate in [("matching rows", matching), ("fixed variables", fixed_givens)]:
    result = pyo.SolverFactory("appsi_highs").solve(candidate)
    assert pyo.check_optimal_termination(result)
    rows = list(candidate.component_data_objects(pyo.Constraint))
    variables = list(candidate.component_data_objects(pyo.Var))
    print(f"{label}: variables={len(variables)}, rows={len(rows)}, "
          f"fixed={sum(v.fixed for v in variables)}, objective={pyo.value(candidate.obj):.0f}")
    assert extract_solution(candidate, example) == solution

# A solver return on one board is evidence about that return, not all boards.
for label, board in [("30-given board", example), ("empty 9x9", Puzzle([[0] * 9 for _ in range(9)]))]:
    relaxed = build_sudoku_model(board)
    pyo.TransformationFactory("core.relax_integer_vars").apply_to(relaxed)
    result = pyo.SolverFactory("appsi_highs").solve(relaxed)
    assert pyo.check_optimal_termination(result)
    fractions = [abs(pyo.value(v) - round(pyo.value(v))) for v in relaxed.y.values()]
    violation = max(abs(pyo.value(c.body) - pyo.value(c.lower))
                    for c in relaxed.component_data_objects(pyo.Constraint))
    print(f"LP {label}: fractional={sum(d > 1e-7 for d in fractions)}, "
          f"max distance={max(fractions):.6f}, max row residual={violation:.2e}, "
          f"objective={pyo.value(relaxed.obj):.0f}")

# A directly checked fractional optimum of the empty-board relaxation.
uniform = build_sudoku_model(Puzzle([[0] * 9 for _ in range(9)]))
pyo.TransformationFactory("core.relax_integer_vars").apply_to(uniform)
for var in uniform.y.values():
    var.set_value(1 / 9)
residual = max(abs(pyo.value(c.body) - pyo.value(c.lower))
               for c in uniform.component_data_objects(pyo.Constraint))
assert residual < 1e-12
print(f"Uniform empty board: each variable=1/9, max row residual={residual:.2e}, objective=0")

matching rows: variables=729, rows=1053, fixed=0, objective=0

fixed variables: variables=729, rows=324, fixed=30, objective=0

LP 30-given board: fractional=0, max distance=0.000000, max row residual=0.00e+00, objective=0

LP empty 9x9: fractional=238, max distance=0.498629, max row residual=5.77e-15, objective=0
Uniform empty board: each variable=1/9, max row residual=0.00e+00, objective=0

Computing Average Sudoku Solve Times¶

For a given positive integer $a\geq2$ , let $n=a^2$ . We are curious how the complexity of solving an $n\times n$ sudoku scales. It is fair to hypothesize that the complexity scales at a superlinear rate, given that the number of cells in an $n\times n$ sudoku is $n^2=a^4$ , and the number of binary variables --- one per (value, row, column) triple --- is $n^3$ .

Our goal is to test two solvers (HiGHS and CBC) to determine which solver, if any, is fastest at solving sudokus. The following program reads the first 10 puzzles of each size from the committed CSV files. The number read can be modified (iterations below), but it is set to 10 by default for a quick computation time.

We test the average solve time for an $n\times n$ sudoku for $a=2,3,4,5,6$ , and we record the number of binary variables in each model at the same time so the two can be plotted together.

# Computes the average time to build and solve a sudoku of each size

# Set iterations to a number from 1 to the number of puzzles in the CSV file
# (100 by default). Determines how many puzzles go into each average.
# A full run of this cell takes roughly 3-4 minutes, almost all of it CBC at
# n = 36. Drop 36 from `sizes` for a much faster run.
iterations = 10

sizes = [4, 9, 16, 25, 36]

# DELIBERATE EXCEPTION to the course default. Everywhere else in this repo we
# call HiGHS ("appsi_highs") for LPs and MILPs -- see org/pyomo-style-guide.md,
# section 7a. This cell keeps "cbc" because BENCHMARKING THE TWO IS THE POINT of
# this section; the conclusion below is drawn from the comparison. Do not
# "standardise" this list.
solvers = ["appsi_highs", "cbc"]

# Read every dataset once, rather than once per solve
datasets = {size: load_puzzles(size) for size in sizes}

# Number of binary variables in the model, counted from Pyomo rather than
# assumed. The prose below claims this is n**3; the assertion checks it.
binary_counts = []
for size in sizes:
    m = build_sudoku_model(Puzzle(datasets[size][0]))
    n_binary = sum(1 for v in m.component_data_objects(pyo.Var) if v.is_binary())
    assert n_binary == size**3, f"{size}: {n_binary} binaries, expected {size**3}"
    binary_counts.append(n_binary)

print("size  binary variables")
for size, n_binary in zip(sizes, binary_counts):
    print(f"{size:4d}  {n_binary:16,d}")
print()

# Average wall-clock time to build and solve, per solver and size
average_times = {}

for solver in solvers:
    averages = []
    print(f"For the solver {solver}:")
    for size in sizes:
        times = []
        for i in range(iterations):
            puzzle = Puzzle(datasets[size][i])

            start_time = time.time()
            solve_sudoku(puzzle, solver=solver)
            times.append(time.time() - start_time)

        averages.append(sum(times) / len(times))
        print(f"  Average time for {size} by {size} sudoku: {averages[-1]:.5f} s")
    average_times[solver] = averages
    print()

size  binary variables
   4                64
   9               729
  16             4,096
  25            15,625
  36            46,656

For the solver appsi_highs:

  Average time for 4 by 4 sudoku: 0.00529 s

  Average time for 9 by 9 sudoku: 0.02751 s

  Average time for 16 by 16 sudoku: 0.14794 s

  Average time for 25 by 25 sudoku: 0.52916 s

  Average time for 36 by 36 sudoku: 1.51825 s

For the solver cbc:

  Average time for 4 by 4 sudoku: 0.07930 s

  Average time for 9 by 9 sudoku: 0.09627 s

  Average time for 16 by 16 sudoku: 0.20839 s

  Average time for 25 by 25 sudoku: 0.53489 s

  Average time for 36 by 36 sudoku: 17.03840 s

Solve Time and Model Size¶

The figure below puts both quantities on one pair of logarithmic axes: average solve time on the left, number of binary variables on the right. On log-log axes a power law $y \propto n^p$ is a straight line of slope $p$ , so the binary count --- exactly $n^3$ , asserted against Pyomo in the cell above --- must plot as a perfectly straight line of slope 3, and the two time curves can be read against it.

# Least-squares slopes on log-log axes: y ~ n**p plots as a line of slope p
log_sizes = np.log(sizes)
slopes = {s: np.polyfit(log_sizes, np.log(average_times[s]), 1)[0] for s in solvers}
slope_binary = np.polyfit(log_sizes, np.log(binary_counts), 1)[0]

print(f"log-log slope, binary variables : {slope_binary:.4f}  (n**3 predicts 3)")
for s in solvers:
    print(f"log-log slope, {s:12s}: {slopes[s]:.2f}")

fig, ax_time = plt.subplots(figsize=(7.2, 5))

# Two solvers on the left axis. Every series carries a marker AND a linestyle,
# so the figure is still readable printed in black and white.
ax_time.loglog(
    sizes,
    average_times["appsi_highs"],
    color="#0072B2",
    marker="o",
    linestyle="-",
    label="HiGHS, solve time",
)
ax_time.loglog(
    sizes,
    average_times["cbc"],
    color="#D55E00",
    marker="s",
    linestyle="--",
    label="CBC, solve time",
)
ax_time.set_xlabel("Puzzle side length $n$")
ax_time.set_ylabel("Average solve time [s]")

# Model size on the right axis, also logarithmic
ax_binary = ax_time.twinx()
ax_binary.loglog(
    sizes,
    binary_counts,
    color="black",
    marker="^",
    linestyle=":",
    label="binary variables ($n^3$)",
)
ax_binary.set_ylabel("Binary variables [-]")

ax_time.set_xticks(sizes)
ax_time.set_xticklabels([str(s) for s in sizes])
ax_time.minorticks_off()

handles = ax_time.get_lines() + ax_binary.get_lines()
ax_time.legend(handles, [h.get_label() for h in handles], loc="upper left", fontsize=10)
plt.show()

log-log slope, binary variables : 3.0000  (n**3 predicts 3)
log-log slope, appsi_highs : 2.61
log-log slope, cbc         : 2.09

The binary count is exactly $n^3$ ; on log-log axes its slope is 3. The stored run returned HiGHS times of approximately 0.0053, 0.0275, 0.1479, 0.5292 and 1.5183 seconds, versus CBC times 0.0793, 0.0963, 0.2084, 0.5349 and 17.0384 seconds for $n=4,9,16,25,36$ . The empirical log-log slopes in that run were 2.61 and 2.09. They are not asymptotic complexity bounds, and neither establishes growth faster than the cubic variable count.

The interfaces differ: appsi_highs builds the model in memory through highspy, while cbc launches an executable and exchanges files. This can contribute to small-instance overhead. These end-to-end timings do not isolate interface cost, model construction, presolve or search, so they do not establish which component caused the observed gap.

The published benchmark was measured on August 25, 2026, on an Apple silicon macOS laptop, Python 3.13, Pyomo 6.10, HiGHS through highspy, and CBC 2.10 from idaes get-extensions. The full cell took about 3 minutes 27 seconds there. Results depend on machine, load, solver version, settings and the chosen puzzles. The dataset and row ordering are held fixed. One of the first four $36\times36$ CBC instances took 62.4 seconds, versus 1.2, 0.8 and 1.0 seconds for the next three; a few hard instances can dominate an average. Re-run and report your own conditions before comparing performance.

Conclusion¶

The notebook supplies a complete model, a reproducible puzzle dataset, solved-board pictures and a measured comparison of two solver interfaces. For side length $n=a^2$ , the declared binary count is $n^3$ . The finite timing sweep describes the tested instances and implementations; it does not prove an asymptotic running-time law.

The seeded generator preserves the Sudoku property of a completed board before blanking cells. It does not test whether the remaining givens have a unique completion. An integer no-good cut can exclude one completed assignment, after which a second solve checks whether another completion exists; Bynum et al., Pyomo, 3rd ed., pp. 78–81 develops that enumeration. Solver tolerances and termination must be checked before declaring uniqueness. More demanding generators, larger boards or other solvers are possible extensions; they require new measurements rather than extrapolation of the current timings.