Solving Sudoku with Integer Linear Programming

Introduction

If you've ever solved a Sudoku by hand you know the drill — scan rows, columns, boxes, pencil in candidates, eliminate, repeat. It gets the job done, but on hard puzzles it's painfully slow and you inevitably mess something up.

I wanted to see if there was a more… mechanical way to do it. Turns out, you can describe every single rule of Sudoku as a math equation, toss the whole thing at an optimization solver, and get the answer back in milliseconds. No guessing, no erasing.

The technique is called Integer Linear Programming (ILP). It's the same math behind airline scheduling, warehouse logistics, and chip design — except here we're using it to crack a newspaper puzzle. Feels a bit like bringing a cannon to a knife fight, but it works beautifully.

1. Sudoku Rules — Quick Refresher

A standard 9×9 Sudoku grid has three rules:

1. Row: each row contains digits 1–9 exactly once
2. Column: each column contains digits 1–9 exactly once
3. Box: each 3×3 sub-grid contains digits 1–9 exactly once

Some cells come pre-filled (the "givens"). You fill the rest.

What's nice is that these three rules translate almost word-for-word into linear constraints.

2. The ILP Formulation

Decision Variables

We introduce binary variables:

where:

• $r$ = row (0–8)
• $c$ = column (0–8)
• $v$ = value (1–9)

$x_{r,c,v} = 1$ means "cell $(r, c)$ holds digit $v$."

So that's $9 \times 9 \times 9 = 729$ binary variables total. Sounds like a lot, but it's tiny for a solver.

Constraints

1. One value per cell:

In plain English — for every cell, exactly one of the nine digit-slots is turned on.

2. Each digit once per row:

3. Each digit once per column:

4. Each digit once per 3×3 box:

5. Pre-filled cells:

If the puzzle tells us cell $(r, c)$ is $v^*$, we pin it:

Objective Function

There isn't one — Sudoku is a feasibility problem, not an optimization problem. We just need any assignment that passes all the constraints. In practice solvers want something, so we set a dummy:

Summary

For a modern solver, this is trivially small.

3. Python Implementation

Here's a complete solver using Google's OR-Tools (free, open-source). I've been using it for various projects and it handles this kind of thing effortlessly:

from ortools.sat.python import cp_model

def solve_sudoku(grid):
    """
    Solve a 9x9 Sudoku puzzle using Integer Linear Programming.
    grid: 9x9 list of ints (0 = empty cell)
    Returns the solved 9x9 grid, or None if unsolvable.
    """
    model = cp_model.CpModel()

    # Decision variables: x[r][c][v] = 1 if cell (r,c) holds digit v+1
    x = [[[model.NewBoolVar(f'x_{r}_{c}_{v}')
            for v in range(9)]
           for c in range(9)]
          for r in range(9)]

    # Constraint 1: each cell has exactly one digit
    for r in range(9):
        for c in range(9):
            model.Add(sum(x[r][c][v] for v in range(9)) == 1)

    # Constraint 2: each digit appears once per row
    for r in range(9):
        for v in range(9):
            model.Add(sum(x[r][c][v] for c in range(9)) == 1)

    # Constraint 3: each digit appears once per column
    for c in range(9):
        for v in range(9):
            model.Add(sum(x[r][c][v] for r in range(9)) == 1)

    # Constraint 4: each digit appears once per 3x3 box
    for box_r in range(3):
        for box_c in range(3):
            for v in range(9):
                model.Add(sum(
                    x[r][c][v]
                    for r in range(box_r * 3, box_r * 3 + 3)
                    for c in range(box_c * 3, box_c * 3 + 3)
                ) == 1)

    # Constraint 5: fix the givens
    for r in range(9):
        for c in range(9):
            if grid[r][c] != 0:
                model.Add(x[r][c][grid[r][c] - 1] == 1)

    # Solve
    solver = cp_model.CpSolver()
    status = solver.Solve(model)

    if status in (cp_model.OPTIMAL, cp_model.FEASIBLE):
        solution = []
        for r in range(9):
            row = []
            for c in range(9):
                for v in range(9):
                    if solver.Value(x[r][c][v]) == 1:
                        row.append(v + 1)
                        break
            solution.append(row)
        return solution

    return None  # no solution exists

Install OR-Tools:

pip install ortools

4. Walking Through an Example

Let's throw a real puzzle at it. Here's one you might recognise from textbooks (0 means empty):

 5 3 .  | . 7 .  | . . .
 6 . .  | 1 9 5  | . . .
 . 9 8  | . . .  | . 6 .
---------+---------+--------
 8 . .  | . 6 .  | . . 3
 4 . .  | 8 . 3  | . . 1
 7 . .  | . 2 .  | . . 6
---------+---------+--------
 . 6 .  | . . .  | 2 8 .
 . . .  | 4 1 9  | . . 5
 . . .  | . 8 .  | . 7 9

Feed it in:

puzzle = [
    [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],
]

solution = solve_sudoku(puzzle)

for row in solution:
    print(" ".join(str(d) for d in row))

Output:

5 3 4 6 7 8 9 1 2
6 7 2 1 9 5 3 4 8
1 9 8 3 4 2 5 6 7
8 5 9 7 6 1 4 2 3
4 2 6 8 5 3 7 9 1
7 1 3 9 2 4 8 5 6
9 6 1 5 3 7 2 8 4
2 8 7 4 1 9 6 3 5
3 4 5 2 8 6 1 7 9

Done in under 10 ms. Every row, column, and 3×3 box checks out.

What's happening inside the solver?

Roughly:

1. It takes in the 729 variables and 324+ constraints
2. Constraint propagation narrows things down — similar to what you'd do mentally, eliminating impossibilites
3. When propagation stalls, the solver branches: picks an undecided variable, tries both values
4. If it hits a contradiction, it backtracks
5. Repeat until it finds a valid assignment (or proves there's none)

Modern solvers also throw in clause learning and cutting planes to speed things up further. For a problem this size it barely breaks a sweat.

5. How Does ILP Compare?

There are plenty of ways to solve Sudoku. Here's a rough comparison:

Where ILP really shines is extensibility. Want to bolt on extra rules? Just add more constraints. No new algorithm needed:

• Diagonal Sudoku — digits 1–9 on both main diagonals:

# main diagonal
for v in range(9):
    model.Add(sum(x[i][i][v] for i in range(9)) == 1)

# anti-diagonal
for v in range(9):
    model.Add(sum(x[i][8 - i][v] for i in range(9)) == 1)

• Anti-King Sudoku — identical digits can't touch diagonally
• Killer Sudoku — cages with sum constraints
• Thermo Sudoku — digits increase along thermometer paths

Each variant is just a few extra lines. The solver doesn't care.

6. Beyond Sudoku

The modeling pattern we used here isn't specific to puzzles. The exact same approach works on real industrial problems:

The recipe is always:

1. Define binary decision variables
2. Express your rules as linear equalities / inequalities
3. Optionally add an objective to minimize or maximize
4. Let the solver handle the rest

Try It Yourself

Theory only gets you so far. Use the interactive solver below to enter any Sudoku and watch it get cracked in your browser — no server, no install. Pick one of the example puzzles or type your own.

Wrapping Up

Sudoku turns out to be a great on-ramp to Integer Linear Programming. The rules are intuitive, the formulation is clean (729 variables, 324 constraints), and the solver nails it in milliseconds. And the best part — the exact same framework scales to problems that matter in industry.

Next time you're staring at a hard grid, you can either keep erasing… or just model it and let the solver do the work.

Full code above is ready to copy-paste. Grab ortools and try it on your own puzzles.