Sudoku

This tutorial explains how to solve a Sudoku puzzle using an Ising machine with Amplify.

About Sudoku

Sudoku is a puzzle that fills a \(9\times9\) grid with the numbers from \(1\) to \(9\) accordingly to the following rules:

  • Place one number from \(1\) to \(9\) to each cell.

  • No identical numbers may be placed in each row and column, and in each of the nine \(3\times3\) subgrids.

First, we are given initial clues, where 17 or more numbers are placed (it has been proven that no solution can be found with 16 or less clues). The game is played by determining the numbers that go into the unfilled cells accordingly to the rules above. When the game’s difficulty level is low, it is relatively easy to determine which numbers to place. As the difficulty of the game increases, and it requires some experience to solve the puzzles.

https://upload.wikimedia.org/wikipedia/commons/thumb/e/e0/Sudoku_Puzzle_by_L2G-20050714_standardized_layout.svg/361px-Sudoku_Puzzle_by_L2G-20050714_standardized_layout.svg.png

Reference: Wikipedia Sudoku

Various algorithms have been devised to solve Sudoku puzzles using computers, such as depth-first search, probabilistic methods, constraint satisfaction problems, exact cover problems, etc. By following these algorithms, Sudoku puzzles can be solved mechanically.

In this tutorial, we will show how to solve Sudoku puzzles using an Ising machine. Formulating the problem from the above Sudoku rules, we define the corresponding function, where finding its minimum value is equivalent to solving the Sudoku puzzle. In the case of Sudoku, only constraint conditions are needed, and no cost function appears.

Now let us take a look at how the code for solving a Sudoku puzzle is implemented using Amplify.

Formulation of Constraints

How to Express Sudoku Constraints on Annealing Machines

First, we consider how to express the constraint conditions that correspond to the Sudoku rules, using binary variables whose values are \(\{0, 1\}\).

The rules to be expressed are the followings:

  1. Each row contains the digits from \(1\) to \(9\) without overlap.

  2. Each column contains the digits \(1\) to \(9\) without overlap.

  3. Each \(3\times3\) subgrid contains the digits \(1\) to \(9\) without overlap.

Without imposing the above rules 1-3, \(9\) possible numbers can be placed in each of \(9\times9=81\) cells. We thus consider assigning \(9\) binary variables to each of the cells, so that \(9\times9\times9=729\) variables are prepared in total. One way to describe this situation is that we stack \(9\) layers of \(9\times9=81\) grids, so that each layer corresponds to the numbers from \(1\) to \(9\), respectively. For the \(9\times9=81\) grid in each layer, we let \(i\) and \(j\) be the indices expressing the rows and columns, respectively, with \(i,j=0,\cdots,8\), corresponding to \(1, \cdots, 9\) th row and column. We also assign the layer index \(k=0,\dots,8\), so that the layer specified by the index corresponds to numbers \(1,\cdots,9\), respectively.

From this observation, we can express the binary variables by \(q_{i,j,k}\), where \(i,j,k=0,\cdots,8\) correspond to the indices for rows, columns, and layers, respectively. For example, \(q_{2,4,6}=1\) means that number \(7\) is placed in the cell at the 3rd row and 5th column, and \(q_{2,4,6}=0\) means that \(7\) is not to be placed here.

Using these variables, the constraints 1-3 can be written as the following one-hot constraints, respectively.

\[\begin{split}\left\{ \begin{align} &\begin{split} (a) \quad &\sum_{j=0}^8 q_{i,j,k}=1 \end{split}\\ &\begin{split} (b) \quad &\sum_{i=0}^8 q_{i,j,k}=1 \end{split}\\ &\begin{split} (c) \quad &\sum_{i,j\in 3\times3\,\text{subgrids}}q_{i,j,k}=1 \end{split}\\ &\begin{split} (d) \quad &\sum_{k=0}^8 q_{i,j,k}=1 \end{split} \end{align} \right.\end{split}\]

The constraints \((a)\), \((b)\) and \((c)\) correspond to rules 1, 2 and 3, respectively. The \((d)\) constraint corresponds to the basic condition that each cell can only contain one number.

Initialization

In Sudoku, clues are given as an initial placement with 17 or more numbers already placed. Here we use the following initial placement which is considered a difficult problem. In the following notation, the cells with no filled numbers are marked with \(0\).

import numpy as np

# Indicate initial placement as a list
# Reference: http://www.sudokugame.org/archive/printable.php?nd=4&y=2018&m=03&d=1
initial = np.array(
    [
        [2, 0, 5, 1, 3, 0, 0, 0, 4],
        [0, 0, 0, 0, 4, 8, 0, 0, 0],
        [0, 0, 0, 0, 0, 7, 0, 2, 0],
        [0, 3, 8, 5, 0, 0, 0, 9, 2],
        [0, 0, 0, 0, 9, 0, 7, 0, 0],
        [0, 0, 0, 0, 0, 0, 4, 5, 0],
        [8, 6, 0, 9, 7, 0, 0, 0, 0],
        [9, 5, 0, 0, 0, 0, 0, 3, 1],
        [0, 0, 4, 0, 0, 0, 0, 0, 0],
    ]
)

# Display Sudoku grid with numbers
def print_sudoku(sudoku):
    print("\n\n")
    for i in range(len(sudoku)):
        line = ""
        if i == 3 or i == 6:
            print("---------------------")
        for j in range(len(sudoku[i])):
            if j == 3 or j == 6:
                line += "| "
            line += str(sudoku[i][j]) + " "
        print(line)
>>> print_sudoku(initial)
2 0 5 | 1 3 0 | 0 0 4
0 0 0 | 0 4 8 | 0 0 0
0 0 0 | 0 0 7 | 0 2 0
---------------------
0 3 8 | 5 0 0 | 0 9 2
0 0 0 | 0 9 0 | 7 0 0
0 0 0 | 0 0 0 | 4 5 0
---------------------
8 6 0 | 9 7 0 | 0 0 0
9 5 0 | 0 0 0 | 0 3 1
0 0 4 | 0 0 0 | 0 0 0

Creating Constraints

Define Variables and Reflect Initial Placement

First, we prepare variables using gen_symbols() and BinaryPoly() provided by Amplify.

from amplify import BinaryPoly, gen_symbols

q = gen_symbols(BinaryPoly, 9, 9, 9)

This gives a three-dimensional array with \(9^3=729\) variables. Each of 9, 9, 9 represents the number of elements in the rows, columns, and number layers, and their indices are i, j, and k, respectively. Each element can be accessed with q[i][j][k]. For example, the 9 variables for row number \(0\) and column number \(0\) can be displayed as follows:

>>> q[0][0]
[q_0, q_1, q_2, q_3, q_4, q_5, q_6, q_7, q_8]

From the initial placement stored in initial, we can narrow down the unknown variables that can be determined accordingly to the rules of Sudoku. For example, the cell indicated by i=1 and j=5 already has a number \(8\) (initial[1][5]=8), so the variable in the layer k=7 corresponding to the number \(8\) can be specified by q[1][5][7]=1. In this specific case, we can place some numbers accordingly to the rules in the following way.

Rules 1 and 2 allow us to further determine the values of the variables q[i][5][7]=0 (\(i\neq1\)), q[1][j][7]=0 (\(j\neq5\)), and q[i][5][7]=0 (\(j\neq5\)), since the row and column containing the cell indicated by i=1 and j=5 in the layer k=7 cannot contain the same number as \(8\). This corresponds to imposing the constraints \((a)\) and \((b)\).

Also, by rule 3, the \(3\times3\) subgrid, which this number \(8\) above belongs to, cannot contain \(8\) any more, so in \((i,j)\in\{(0,3),(0,4),(0,5),(1,3),(1,4),(2,3),(2,4),(2,5)\}\), we can let q[i][j][7]=0. We can imposed the constraint \((c)\) in this manner.

Furthermore, we impose the constraint \((d)\) to place only one of the numbers in a cell where the number is fixed. In the above example, we have q[1][5][k]=0 ( \(k\neq7\) ). By performing similar operations for all the cells given from the initial layout, we can narrow down the required variables so that we can perform the calculations with fewer number of variables.

for i, j in zip(*np.where(initial != 0)):  # Obtain the indices of non-zero rows, columns, and layers
    k = initial[i, j] - 1  # Note that -1 is used to convert from value to index

    for m in range(9):
        q[i][m][k] = BinaryPoly(0)  # Constraint(a)
        q[m][j][k] = BinaryPoly(0)  # Constraint(b)
        q[i][j][m] = BinaryPoly(0)  # Constraint(d)
        q[3 * (i // 3) + m // 3][3 * (j // 3) + m % 3][k] = BinaryPoly(0)  # Constraint(c)

    q[i][j][k] = BinaryPoly(1)  # Set the value of the variable to 1

We have now prepared the initial settings. As an example, if the 9 variables with the row number \(0\) and the column number \(0\) are displayed, the second element is fixed as \(1\), that is, the number \(2\) is assigned to this cell.

>>> q[0][0]
[0, 1, 0, 0, 0, 0, 0, 0, 0]

Similarly, if the 9 variables with the row number \(0\) and column number \(1\) are displayed, it is seen that the numbers \(1,2,3,4,5,6\) are not candidates.

>>> q[0][1]
[0, 0, 0, 0, 0, 0, q_15, q_16, q_17]

Setting Constraints

Next, we show how the constraint parts are implemented. The one-hot constraint in \((a)\)-\((d)\) can be expressed using Amplify’s equal_to() function.

We start by defining the constraint condition corresponding to \((a)\) that each row of cannot contain the same number. Since the sum of the variables for all columns specified by row i and layer k is given by \(\sum_{j=0}^{8}q_{i,j,k}\), the constraint can be obtained by setting this sum to \(1\) as follows.

from amplify.constraint import equal_to

# (a): Constraint that each row cannot contain the same number
row_constraints = [
    equal_to(sum([q[i][j][k] for j in range(9)]), 1) for i in range(9) for k in range(9)
]

Similarly, the constraint \((b)\) that each column has cannot contain the same number and \((d)\) that only one number can be placed in a given cell can be expressed as follows:

# (b): Constraint that each column cannot contain the same number
col_constraints = [
    equal_to(sum([q[i][j][k] for i in range(9)]), 1) for j in range(9) for k in range(9)
]

# (d): Constraint that only one number can be placed in a cell
num_constraints = [
    equal_to(sum([q[i][j][k] for k in range(9)]), 1) for i in range(9) for j in range(9)
]

Finally, we express the constraint \((c)\) that each \(3\times3\) subgrid cannot contain the same number. We take the sum of the variables in each \(3\times3\) subgrid for all layers and impose the one-hot constraint as follows:

# (c): Constraint that 3x3 subgrids cannot contain the same number
block_constraints = [
    equal_to(sum([q[i + m // 3][j + m % 3][k] for m in range(9)]), 1)
    for i in range(0, 9, 3)
    for j in range(0, 9, 3)
    for k in range(9)
]

We have now obtained all the constraints, and we can add all these constraints together to obtain the total constraint object.

constraints = (
    sum(row_constraints)
    + sum(col_constraints)
    + sum(num_constraints)
    + sum(block_constraints)
)

This completes the formulation of the Sudoku problem. If the Ising machine can find the combination of variables that satisfies all the constraints, this allows us to determine the numbers to be placed for all the cells, giving the solution for Sudoku.

Running the Ising Machine

We will now prepare to run the Ising machine using the constraints we created in the previous section. First, create a client for the Ising machine and set parameters such as timeout. We then create the solver from the configured client.

from amplify import Solver
from amplify.client import FixstarsClient

client = FixstarsClient()
client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
client.parameters.timeout = 1000  # Timeout is 1 second

solver = Solver(client)

We formulate the constraint constraints as a logical model class by giving it to BinaryQuadraticModel(), and we can get the result by giving it to the solver() method of the previously set solver.

from amplify import BinaryQuadraticModel

model = BinaryQuadraticModel(constraints)
result = solver.solve(model)
if len(result) == 0:
    raise RuntimeError("Some of the constraints are not satisfied.")

values = result[0].values

The result of the execution is stored in values. By giving the variable array q and values to the decode_solution() function, the variables are evaluated with the result values. Then, by searching for k such that q[i][j][k] = 1 for all i and j, we can get the number k + 1 as the solution for the indicated cell.

from amplify import decode_solution

q_values = decode_solution(q, values)
answer = np.array([np.where(np.array(q_values[i]) != 0)[1] + 1 for i in range(9)])

Finally, the print_sudoku function is used to output the solution.

>>> print_sudoku(answer)
2 8 5 | 1 3 9 | 6 7 4
6 7 3 | 2 4 8 | 5 1 9
4 1 9 | 6 5 7 | 3 2 8
---------------------
7 3 8 | 5 6 4 | 1 9 2
5 4 2 | 3 9 1 | 7 8 6
1 9 6 | 7 8 2 | 4 5 3
---------------------
8 6 1 | 9 7 3 | 2 4 5
9 5 7 | 4 2 6 | 8 3 1
3 2 4 | 8 1 5 | 9 6 7

Generalization of Sudoku

So far we have dealt with a \(9 \times 9\) grid with \(3 \times 3\) subgrids, but the Ising machine can easily handle Sudoku puzzles with extended problem sizes, such as \(16 \times 16\) and \(25 \times 25\). We let the number of cells in a Sudoku puzzle be \(N\times N,\,(N\in\mathbb{Z})\) and the corresponding subgrids be in the form of \(n\times n\) (where \(N=n^2,\,(n\in\mathbb{Z})\)). For example, a basic \(9\times9\) Sudoku would have \(N=9\) and \(n=3\).

We generalize the previous code using \(N\) and \(n\) and solve the \(16\times16\) grid Sudoku as an example.

from amplify import (
    BinaryPoly,
    BinaryQuadraticModel,
    gen_symbols,
    Solver,
    decode_solution,
)
from amplify.constraint import equal_to
import numpy as np
import math

# Display a sudoku grid with numbers
def print_sudoku(sudoku):
    print("\n\n")
    N = len(sudoku)
    n = int(math.sqrt(N))
    width = len(str(N))
    for i in range(len(sudoku)):
        line = ""
        if i % n == 0 and i != 0:
            print("-" * ((width + 1) * n * n + 2 * (n - 1)))
        for j in range(len(sudoku[i])):
            if j % n == 0 and j != 0:
                line += "| "
            line += str(sudoku[i][j]).rjust(width) + " "
        print(line)

# Construct variables and constraints from initial configuration and return them
def get_model(initial):
    N = len(initial)
    n = int(math.sqrt(N))
    q = gen_symbols(BinaryPoly, N, N, N)

    for i, j in zip(*np.where(initial != 0)):
        k = initial[i, j] - 1

        for m in range(N):
            q[i][m][k] = BinaryPoly(0)  # constraint (a)
            q[m][j][k] = BinaryPoly(0)  # constraint (b)
            q[i][j][m] = BinaryPoly(0)  # constraint (d)
            q[n * (i // n) + m // n][n * (j // n) + m % n][k] = BinaryPoly(0)  # constraint (c)

        q[i][j][k] = BinaryPoly(1)  # Substitute 1

    # (a): Constraints that each row cannot contain the same number
    row_constraints = [
        equal_to(sum(q[i][j][k] for j in range(N)), 1) for i in range(N) for k in range(N)
    ]

    # (b): Constraints that each column cannot contain the same number
    col_constraints = [
        equal_to(sum(q[i][j][k] for i in range(N)), 1) for j in range(N) for k in range(N)
    ]

    # (d): Constraints that only one number can be in a cell
    num_constraints = [
        equal_to(sum(q[i][j][k] for k in range(N)), 1) for i in range(N) for j in range(N)
    ]

    # (c): Constraints that a nxn subgrid cannot contain the same number
    block_constraints = [
        equal_to(sum([q[i + m // n][j + m % n][k] for m in range(N)]), 1)
        for i in range(0, N, n)
        for j in range(0, N, n)
        for k in range(N)
    ]

    constraints = (
        sum(row_constraints)
        + sum(col_constraints)
        + sum(num_constraints)
        + sum(block_constraints)
    )

    return q, BinaryQuadraticModel(constraints)
# Initial configuration for n = 4 (N = 16)
# Reference: https://www.free-sudoku-puzzle.com/puzzle_fours/solve/3/238
initial = np.array(
    [
        [0, 0, 0, 7, 0, 0, 0, 1, 5, 0, 3, 16, 4, 0, 15, 0],
        [0, 11, 0, 0, 0, 0, 5, 0, 0, 2, 12, 6, 0, 0, 7, 14],
        [4, 0, 0, 0, 7, 8, 9, 0, 11, 0, 1, 15, 0, 0, 10, 0],
        [10, 0, 0, 0, 0, 0, 0, 15, 13, 0, 9, 7, 8, 0, 0, 1],
        [13, 0, 0, 16, 15, 0, 4, 9, 0, 0, 14, 0, 11, 0, 1, 0],
        [8, 0, 5, 0, 0, 0, 10, 0, 0, 0, 0, 0, 15, 0, 0, 0],
        [0, 0, 0, 11, 0, 0, 0, 8, 16, 7, 0, 9, 0, 0, 0, 0],
        [0, 0, 0, 14, 0, 0, 3, 0, 4, 0, 0, 5, 13, 0, 0, 0],
        [0, 0, 0, 0, 3, 0, 0, 14, 0, 0, 4, 0, 9, 12, 8, 15],
        [0, 0, 0, 0, 0, 1, 7, 10, 0, 15, 8, 11, 0, 0, 0, 0],
        [0, 0, 0, 0, 11, 12, 0, 0, 0, 0, 16, 0, 3, 5, 0, 0],
        [0, 0, 0, 0, 0, 16, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 6, 16, 0, 15, 1, 5, 0, 14, 2, 0, 0],
        [0, 0, 0, 3, 0, 0, 0, 0, 9, 0, 0, 14, 0, 1, 0, 4],
        [2, 0, 12, 0, 0, 0, 0, 0, 0, 16, 13, 0, 6, 0, 3, 5],
        [1, 0, 0, 0, 0, 15, 0, 0, 2, 11, 6, 12, 7, 9, 0, 10],
    ]
)

# Obtain variables and model
q, model = get_model(initial)
from amplify import Solver, decode_solution
from amplify.client import FixstarsClient

client = FixstarsClient()
client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
client.parameters.timeout = 10000  # Timeout is 10 seconds

solver = Solver(client)
result = solver.solve(model)
if len(result) == 0:
    raise RuntimeError("Some of the constraints are not satisfied.")

values = result[0].values

q_values = decode_solution(q, values)
N = len(q)
answer = np.array([np.where(np.array(q_values[i]) != 0)[1] + 1 for i in range(N)])
>>> print_sudoku(initial)
 0  0  0  7 |  0  0  0  1 |  5  0  3 16 |  4  0 15  0
 0 11  0  0 |  0  0  5  0 |  0  2 12  6 |  0  0  7 14
 4  0  0  0 |  7  8  9  0 | 11  0  1 15 |  0  0 10  0
10  0  0  0 |  0  0  0 15 | 13  0  9  7 |  8  0  0  1
------------------------------------------------------
13  0  0 16 | 15  0  4  9 |  0  0 14  0 | 11  0  1  0
 8  0  5  0 |  0  0 10  0 |  0  0  0  0 | 15  0  0  0
 0  0  0 11 |  0  0  0  8 | 16  7  0  9 |  0  0  0  0
 0  0  0 14 |  0  0  3  0 |  4  0  0  5 | 13  0  0  0
------------------------------------------------------
 0  0  0  0 |  3  0  0 14 |  0  0  4  0 |  9 12  8 15
 0  0  0  0 |  0  1  7 10 |  0 15  8 11 |  0  0  0  0
 0  0  0  0 | 11 12  0  0 |  0  0 16  0 |  3  5  0  0
 0  0  0  0 |  0 16  0  0 |  0  5  0  0 |  0  0  0  0
------------------------------------------------------
 0  0  0  0 |  0  6 16  0 | 15  1  5  0 | 14  2  0  0
 0  0  0  3 |  0  0  0  0 |  9  0  0 14 |  0  1  0  4
 2  0 12  0 |  0  0  0  0 |  0 16 13  0 |  6  0  3  5
 1  0  0  0 |  0 15  0  0 |  2 11  6 12 |  7  9  0 10
>>> print_sudoku(answer)
12  2  8  7 | 14 11  6  1 |  5 10  3 16 |  4 13 15  9
 9 11  1 15 | 10  4  5 13 |  8  2 12  6 | 16  3  7 14
 4  3 16 13 |  7  8  9  2 | 11 14  1 15 |  5  6 10 12
10  5 14  6 | 16  3 12 15 | 13  4  9  7 |  8 11  2  1
------------------------------------------------------
13 12 10 16 | 15  5  4  9 |  6  3 14  8 | 11  7  1  2
 8  7  5  9 |  6 14 10 16 |  1 13 11  2 | 15  4 12  3
 3  4  2 11 | 12 13  1  8 | 16  7 15  9 | 10 14  5  6
 6  1 15 14 |  2  7  3 11 |  4 12 10  5 | 13  8  9 16
------------------------------------------------------
 5 16 11  1 |  3  2 13 14 |  7  6  4 10 |  9 12  8 15
14  6  9 12 |  5  1  7 10 |  3 15  8 11 |  2 16  4 13
15 10 13  2 | 11 12  8  4 | 14  9 16  1 |  3  5  6  7
 7  8  3  4 |  9 16 15  6 | 12  5  2 13 |  1 10 14 11
------------------------------------------------------
11  9  7 10 |  4  6 16 12 | 15  1  5  3 | 14  2 13  8
16 15  6  3 | 13 10  2  5 |  9  8  7 14 | 12  1 11  4
 2 14 12  8 |  1  9 11  7 | 10 16 13  4 |  6 15  3  5
 1 13  4  5 |  8 15 14  3 |  2 11  6 12 |  7  9 16 10