Trying to solve Sudoku with cvxpy

Question

I'm trying to solve a Sudoku with cvxpy optimization package. I'm really new to both optimization and cvxpy.

The constraints are:

1. all the values are between 1 to 9
2. sum of all rows = 45
3. sum of all columns = 45
4. sum of all squares = 45
5. the given numbers (I'm trying to solve a 17-clues Sudoku).

So, here is my code:

from cvxpy import *
import numpy as np

x = Variable((9, 9), integer=True)
obj = Minimize(sum(x)) #whatever, if the constrains are fulfilled it will be fine 
const = [x >= 1, #all values should be >= 1
         x <= 9, #all values should be <= 9
         sum(x, axis=0) == 45,  # sum of all rows should be 45
         sum(x, axis=1) == 45,  # sum of all cols should be 45
         sum(x[0:3, 0:3]) == 45, sum(x[0:3, 3:6]) == 45, #sum of all squares should be 45
         sum(x[0:3, 6:9]) == 45,
         sum(x[3:6, 0:3]) == 45, sum(x[3:6, 3:6]) == 45,
         sum(x[3:6, 6:9]) == 45,
         sum(x[6:9, 0:3]) == 45, sum(x[6:9, 3:6]) == 45,
         sum(x[6:9, 6:9]) == 45,
         x[0, 7] == 7, #the values themselves
         x[0, 8] == 1,
         x[1, 1] == 6,
         x[1, 4] == 3,
         x[2, 4] == 2,
         x[3, 0] == 7,
         x[3, 4] == 6,
         x[3, 6] == 3,
         x[4, 0] == 4,
         x[4, 6] == 2,
         x[5, 0] == 1,
         x[5, 3] == 4,
         x[6, 3] == 7,
         x[6, 5] == 5,
         x[6, 7] == 8,
         x[7, 1] == 2,
         x[8, 3] == 1]

prob = Problem(objective=obj, constraints=const)
prob.solve()

What I get from cvxpy is

prob.status
Out[2]: 'infeasible_inaccurate'

This is for sure a valid Sudoku, and I checked the numbers million times.

Why am I not getting the answer?

Any help would be appreciated!

Answer 1

This is an ECOS_BB problem which you are using by default. It is not a reliable integer programming solver and I suggest not to use it.

Other recommendation: do not use import *. It is much better to use import cvxpy as cp to avoid confusion with other functions with the same name. Also, numpy is not needed here by the way.

The following script returns a feasible solution with GUROBI (you can also use GLPK if you do not have a GUROBI license):

import cvxpy as cp

x = cp.Variable((9, 9), integer=True)

# whatever, if the constrains are fulfilled it will be fine
objective = cp.Minimize(cp.sum(x))
constraints = [x >= 1,  # all values should be >= 1
               x <= 9,  # all values should be <= 9
               cp.sum(x, axis=0) == 45,  # sum of all rows should be 45
               cp.sum(x, axis=1) == 45,  # sum of all cols should be 45
               # sum of all squares should be 45
               cp.sum(x[0:3, 0:3]) == 45, cp.sum(x[0:3, 3:6]) == 45,
               cp.sum(x[0:3, 6:9]) == 45,
               cp.sum(x[3:6, 0:3]) == 45, cp.sum(x[3:6, 3:6]) == 45,
               cp.sum(x[3:6, 6:9]) == 45,
               cp.sum(x[6:9, 0:3]) == 45, cp.sum(x[6:9, 3:6]) == 45,
               cp.sum(x[6:9, 6:9]) == 45,
               x[0, 7] == 7,  # the values themselves
               x[0, 8] == 1,
               x[1, 1] == 6,
               x[1, 4] == 3,
               x[2, 4] == 2,
               x[3, 0] == 7,
               x[3, 4] == 6,
               x[3, 6] == 3,
               x[4, 0] == 4,
               x[4, 6] == 2,
               x[5, 0] == 1,
               x[5, 3] == 4,
               x[6, 3] == 7,
               x[6, 5] == 5,
               x[6, 7] == 8,
               x[7, 1] == 2,
               x[8, 3] == 1]

prob = cp.Problem(objective, constraints)
prob.solve(solver=cp.GUROBI)

print(x.value)

That's the output

In [2]: run sudoku.py
[[1. 6. 1. 4. 7. 9. 9. 7. 1.]
 [6. 6. 1. 1. 3. 9. 9. 9. 1.]
 [8. 7. 9. 1. 2. 9. 1. 7. 1.]
 [7. 7. 1. 9. 6. 1. 3. 2. 9.]
 [4. 9. 5. 9. 5. 1. 2. 1. 9.]
 [1. 2. 9. 4. 9. 1. 9. 1. 9.]
 [8. 1. 1. 7. 8. 5. 2. 8. 5.]
 [9. 2. 9. 9. 4. 1. 1. 1. 9.]
 [1. 5. 9. 1. 1. 9. 9. 9. 1.]]

Answer 2

The documentation (http://cvxr.com/cvx/doc/solver.html) suggests that the root cause may be a close to feasible solution that looks infeasible due to floating point inaccuracy above tolerance, which squares with all the equality constraints.

Stepping back, though, that set of constraints isn't enough to specify a valid Sudoku solution. A better formulation would be to have 0-1 variables indexed by (row, column, large square, digit) and constraints that, for each row/column/square, row/digit, column/digit, square/digit combination, the sum of matching variables is equal to 1.