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!

Salomone

Suggestions: (1) add the verbose flag and inspect the solver log (2) choose a different solver (cvxpy can talk to different solvers) (3) Sudoku puzzles are usually modeled in a different fashion (the all-different constraints are difficult with your setup) – Sphagnum

I had never seen an integer model. I had seen binary models, e.g., Solve Sudoku puzzles via Integer Programming (MATLAB) and A Sudoku problem formulated as an LP (Python). – Ultraconservative

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

Odalisque

Thanks! Just to be accurate the solver is GLPK_MI (not GLPK). Happily accepted! – Salomone

This isn't a valid sudoku. The original post's constraints were incorrect. Each value can only appear once per row, column and submatrix. So the sum() ==45 is incorrect – Blamed

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.

Definitely

I don't care if at this stage I won't get a valid Sudoku, I will add constrains if necessary. All I'm after right now is getting a feasible solution the fulfills my current constraints. – Salomone

@Binyamin then dump the variable settings and see what's going on. – Definitely

it returns a matrix of only 5s, which is valid. If I add just one variable (for example x[0,0]==7), it's not feasible anymore... – Salomone

@Binyamin does it throw an exception if you dump the values anyway? If so, try increasing the tolerance. – Definitely

It doesn't throw an exception. – Salomone

Answer 3

Sum == 45 is not a good constraint since it does not prevent duplicate values. You don't need either numpy or cvxpy to solve this as a simple backtrack algorithm can solve it.

count=0

def findnext():
    global x
    for row in range(0,9):
        for col in range(0, 9):
            if x[row][col] == 0:
                return (row, col)
    return (-1,-1)

def colok(row, col):
    global x
    for c in range(0,9):
        if c != col and x[row][c] == x[row][col]:
            # print "x[%d][%d] == x[%d][%d]" % (row,c,row,col)
            return False
    return True

def rowok(row, col):
    global x
    for r in range(0,9):
        if r != row and x[r][col] == x[row][col]:
            return False
    return True

def sqok(row, col):
    global x
    sqy = int(row/3)*3
    sqx = int(col/3)*3
    for r in range(0,3):
        for c in range(0, 3):
            if (sqx+c != col or sqy+r != row) and x[sqy+r][sqx+c] == x[row][col]:
                print "x[%d][%d] == x[%d][%d] = %d" % (row,c,row,col,x[row][col])
                return False

    return True

def backtrack():
    global x
    global count
    (row, col) = findnext()
    if row < 0:
        return True
    for v in range(1,10):
        x[row][col] = v
        if rowok(row, col) and colok(row,col) and sqok(row,col):
            count += 1
            if backtrack():
                return True
    x[row][col] = 0
    return False

Results in:

ending after calling backtrack 98248 times
2 8 3 6 5 4 9 7 1  
5 6 1 9 3 7 4 2 8  
9 7 4 8 2 1 5 6 3  
7 5 8 2 6 9 3 1 4  
4 3 6 5 1 8 2 9 7  
1 9 2 4 7 3 8 5 6  
3 1 9 7 4 5 6 8 2  
8 2 7 3 9 6 1 4 5  
6 4 5 1 8 2 7 3 9  

Knocker

Yes but I want to use cvx to solve it. I don't care about solving Sudoku, I care about learning optimization and Sudoku is just a use-case. – Salomone

@BinyaminEven how can you learn by not getting the right solution? – Agna