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

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.

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.

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:

- Each row contains the digits from $1$ to $9$ without overlap.
- Each column contains the digits $1$ to $9$ without overlap.
- 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.

$$ \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{subgrid}}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. $$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.

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$.

In [ ]:

```
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)
```

In [ ]:

```
from amplify import BinaryPoly, gen_symbols
q = gen_symbols(BinaryPoly, 9, 9, 9)
```

`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:

In [ ]:

```
q[0][0]
```

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.

In [ ]:

```
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
```

In [ ]:

```
q[0][0]
```

In [ ]:

```
q[0][1]
```

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.

In [ ]:

```
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)
]
```

In [ ]:

```
# (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)
]
```

In [ ]:

```
# (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)
]
```

In [ ]:

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

In [ ]:

```
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)
```

`constraints`

as a logical model class by giving it to
`BinaryQuadraticModel`

, and we can get the result by giving it to the
`Solve.solver`

method of the previously set `solver`

.

In [ ]:

```
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
```

`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.

In [ ]:

```
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.

In [ ]:

```
print_sudoku(answer)
```

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$ mass Sudoku as an example.

In [ ]:

```
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)
```

In [ ]:

```
# 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)])
```

In [ ]:

```
print_sudoku(initial)
```

In [ ]:

```
print_sudoku(answer)
```