# 10.7 SUDOKU¶

SUDOKU is a famous simple yet mind-blowing game. The objective is to fill a 9×9 grid with digits so that each column, each row, and each of the nine 3×3 sub-grids that compose the grid (also called boxes, blocks, regions, or sub-squares) contains all of the digits from 1 to 9. For more information see http://en.wikipedia.org/wiki/Sudoku. Here is a simple example:

In a more general setting we are given a grid of dimension $$n\times n$$, with $$n=m^2, m\in \natnumb$$. Each cell $$(i,j)$$ must be filled with an integer $$y_{ij} \in [1,n]$$. Along each row and each column there must be no repetitions. No repetitions are allowed also in each sub-grid with corners $$\{(mt,ml),\ (m(t+1)-1,m(l+1)-1)\}$$, for $$t,l=0,\ldots,m-1$$ (we index cells from $$(0,0)$$).

In general, each SUDOKU instance comes with a set $$F$$ of predetermined values which:

• reduce the complexity of the game by removing symmetries and guiding the initial moves of the player;
• ensure that there will be a unique solution.

We represent the set $$F$$ as list of triplets $$(i,j,v)$$, meaning that the cell $$(i,j)$$ contains the value $$v$$.

Note that SUDOKU is a feasibility problem. A typical Integer Programming formulation is straightforward: let $$x_{ijk}$$ be a binary variable that takes value $$1$$ if $$k$$ is written in cell $$(i,j)$$. Then we look for a feasible solution of a system of constraints given below.

SUDOKU is a typical assignment problem. Its constraints are commonly found in optimization problems concerning scheduling or resource allocation. SUDOKU has also been a nice problem to fiddle with for many researchers in the optimization community. Indeed, its simple structure and the easy way in which the results can be tested make it a perfect test problem.

We will approach SUDOKU as a standard integer linear program, and we will show how easily and elegantly it can be implemented in Fusion.

Mathematical Formulation

In this section we formulate SUDOKU as a mixed-integer linear optimization problem. Let’s introduce a binary variable $$x_{ijk}$$ that takes value $$1$$ if $$k$$ is written in the cell $$(i,j)$$, or $$0$$ otherwise. We first ask that for each cell exactly one digit is selected:

(1)$\sum_{k=0}^{n-1} x_{ijk} = 1, \qquad i,j=0,\ldots,n-1.$

Similar constraints can be used to force each digit to appear only once in each row or column:

(2)$\begin{split}\begin{array}{lr} \sum_{i=0}^{n-1} x_{ijk} = 1, & j,k=0,\ldots,n-1,\\ \sum_{j=0}^{n-1} x_{ijk} = 1, & i,k=0,\ldots,n-1. \end{array}\end{split}$

To force a digit to appear only once in each sub-grid we can use the following

(3)$\sum_{i=0}^{m-1} \sum_{j=0}^{m-1} x_{(i+tm)(j+tl)k} = 1 \qquad k=0,\ldots,n-1 \text{ and } t,l=0,\ldots,m-1$

If a cell $$(i,j)$$ has a predetermined value, i.e. $$(i,j,k)\in F$$ then we set

$x_{ijk} = 1.$

Summarizing, and considering that there is no objective function to minimize, the optimization model for the SUDOKU problem takes the form

(4)$\begin{split}\begin{array}{lr} \min 0 &\\ \mbox{s.t.} &\\ \sum_{i=0}^{n-1} x_{ijk} = 1, & j,k=0,\ldots,n-1,\\ \sum_{j=0}^{n-1} x_{ijk} = 1, & i,k=0,\ldots,n-1,\\ \sum_{k=0}^{n-1} x_{ijk} = 1, & i,j=0,\ldots,n-1,\\ \sum_{i=0}^{m-1} \sum_{j=0}^{m-1} x_{(i+tm)(j+tl)k} = 1, & k=0,\ldots,n-1 \text{ and } \\ & t,l=0,\ldots,m-1,\\ &\\ x_{ijk} = 1, & \forall (i,j,k)\in F. \end{array}\end{split}$

Implementation with Fusion

The implementation in Fusion is straightforward. First, we represent the variable $$x$$ using a three dimensional Fusion variable:

  Variable::t X = M->variable("X", new_array_ptr<int, 1>({n, n, n}), Domain::binary());


Then we can define constraints (1) and (2) simply using the Expr.sum operator, that allows to sum the elements of an expression (in this case of the variable itself) along arbitrary dimensions. The code reads:

  //each value only once per dimension
for (int d = 0; d < m; d++)
M->constraint( Expr::sum(X, d), Domain::equalsTo(1.) );


The last set of constraints (3) , i.e. the sum over block, needs a little more effort: we must loop over all blocks and select the proper slice:

  //each number must appear only once in a block
for (int k = 0; k < n ; k++)
for (int i = 0; i < m ; i++)
for (int j = 0; j < m ; j++)
M->constraint( Expr::sum( X->slice( new_array_ptr<int, 1>({i * m, j * m, k}),
new_array_ptr<int, 1>({(i + 1)*m, (j + 1)*m, k + 1}) ) ),
Domain::equalsTo(1.) );


To set the triplets given in the set $$F$$ we can use the Variable.pick method that returns a one dimensional view of an arbitrary set of elements of the variable.

  auto fixed = std::shared_ptr< ndarray<int, 2> >( new ndarray<int, 2>( shape(nfixed, 3) ) );

for (int i = 0; i < nfixed; i++)
for (int d = 0; d < m; d++)
(*fixed)(i, d) =  (*hr_fixed)(i, d) - 1;

M->constraint( X->pick( fixed ) , Domain::equalsTo(1.0) ) ;


SUDOKU: the complete example code.

The complete code for the SUDOKU problem is shown in Listing 32.

Listing 32 Fusion implementation to solve SUDOKU. Click here to download.
#include <iostream>
#include <sstream>
#include <cmath>

#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

void print_solution(int n, Variable::t X)
{
using namespace std;

cout << "\n";
int m( std::sqrt(n) );
for (int i = 0; i < n; i++)
{
stringstream ss;

for (int j = 0; j < n; j++)
{
if (j % m == 0) ss << " |";

for (int k = 0; k < n; k++)
{
auto x = X->index(new_array_ptr<int, 1>({i, j, k}))->level();
if ( (*x)[0] > 0.5 )
{
ss << " " << (k + 1);
break;
}
}
}
cout << ss.str() << " |";

cout << "\n";
if ((i + 1) % m == 0)
cout << "\n";
}
}

int main(int argc, char ** argv)
{

int m = 3;
int n = m * m;

//fixed cells in human readable (i.e. 1-based) format
auto hr_fixed = new_array_ptr<int, 2>(
{ {1, 5, 4},
{2, 2, 5}, {2, 3, 8}, {2, 6, 3},
{3, 2, 1}, {3, 4, 2}, {3, 5, 8}, {3, 7, 9},
{4, 2, 7}, {4, 3, 3}, {4, 4, 1}, {4, 7, 8}, {4, 8, 4},
{6, 2, 4}, {6, 3, 1}, {6, 6, 9}, {6, 7, 2}, {6, 8, 7},
{7, 3, 4}, {7, 5, 6}, {7, 6, 5}, {7, 8, 8},
{8, 4, 4}, {8, 7, 1}, {8, 8, 6},
{9, 5, 9}
}
);

int nfixed = hr_fixed->size() / m;

Model::t M = new Model("SUDOKU"); auto _M = finally([&]() { M->dispose(); });

M->setLogHandler([ = ](const std::string & msg) { std::cout << msg << std::flush; });

Variable::t X = M->variable("X", new_array_ptr<int, 1>({n, n, n}), Domain::binary());

//each value only once per dimension
for (int d = 0; d < m; d++)
M->constraint( Expr::sum(X, d), Domain::equalsTo(1.) );

//each number must appear only once in a block
for (int k = 0; k < n ; k++)
for (int i = 0; i < m ; i++)
for (int j = 0; j < m ; j++)
M->constraint( Expr::sum( X->slice( new_array_ptr<int, 1>({i * m, j * m, k}),
new_array_ptr<int, 1>({(i + 1)*m, (j + 1)*m, k + 1}) ) ),
Domain::equalsTo(1.) );

auto fixed = std::shared_ptr< ndarray<int, 2> >( new ndarray<int, 2>( shape(nfixed, 3) ) );

for (int i = 0; i < nfixed; i++)
for (int d = 0; d < m; d++)
(*fixed)(i, d) =  (*hr_fixed)(i, d) - 1;

M->constraint( X->pick( fixed ) , Domain::equalsTo(1.0) ) ;

M->solve();

//print the solution, if any...
if ( M->getPrimalSolutionStatus() == SolutionStatus::Optimal ||
M->getPrimalSolutionStatus() == SolutionStatus::NearOptimal)
print_solution(n, X);
else
std::cout << "No solution found!\n";

return 0;
}


The problem instance corresponding to Fig. 8 is hard-coded for the sake of simplicity. It will produce the following output

Problem
Name                   : SUDOKU
Objective sense        : min
Type                   : LO (linear optimization problem)
Constraints            : 350
Cones                  : 0
Scalar variables       : 1000
Matrix variables       : 0
Integer variables      : 729

Optimizer started.
Mixed integer optimizer started.
Presolve started.
Presolve terminated. Time = 0.00
Presolved problem: 0 variables, 0 constraints, 0 non-zeros
Presolved problem: 0 general integer, 0 binary, 0 continuous
Clique table size: 0
BRANCHES RELAXS   ACT_NDS  DEPTH    BEST_INT_OBJ         BEST_RELAX_OBJ       REL_GAP(%)  TIME
0        1        0        0        0.0000000000e+00     0.0000000000e+00     0.00e+00    0.0
An optimal solution satisfying the relative gap tolerance of 1.00e-02(%) has been located.
The relative gap is 0.00e+00(%).
An optimal solution satisfying the absolute gap tolerance of 0.00e+00 has been located.
The absolute gap is 0.00e+00.

Objective of best integer solution : 0.000000000000e+00
Best objective bound               : -0.000000000000e+00
Construct solution objective       : Not employed
Construct solution # roundings     : 0
User objective cut value           : 0
Number of cuts generated           : 0
Number of branches                 : 0
Number of relaxations solved       : 1
Number of interior point iterations: 0
Number of simplex iterations       : 0
Time spend presolving the root     : 0.00
Time spend in the heuristic        : 0.00
Time spend in the sub optimizers   : 0.00
Time spend optimizing the root   : 0.00
Mixed integer optimizer terminated. Time: 0.02

Optimizer terminated. Time: 0.02

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

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

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