Achieving Quantum Advantage at Solving Super Sudokus with D-Wave
Of Sudokus, Humans, Classical and Quantum Computers
Of Sudokus, Humans, Classical and Quantum Computers
Sudokus are puzzles some people enjoy solving manually. It’s a 9x9 (or (3*3)*(3*3)) puzzle, where in each of the squares, you are supposed to fill out a number from 1 to 9. Each row and column as well as each of the marked 3x3 sub-squares can only contain each number from 1 to 9 once.
The most trained humans solve the 9x9 version in a few minutes, by trial and error, or in computer science terms, by a backtracking algorithm.
Some people prefer to let computers solve sudokus. However, computers find traditional (3*3)*(3*3) Sudokus laughably easy. With a classical Mixed Integer Programming approach these are typically solved in about 15 ms. To increase the challenge, we considered larger (n*n)*(n*n) for n >=4, ‘super’ or ‘meta’ sudokus. When n gets bigger, (n*n)*(n*n) gets bigger at a much faster rate. So at n around 6, 7, 8 classical computers start to struggle, meaning they need hours or days to solve such a sudoku.
This is where, potentially, quantum computers could perform better.
Researching that question is the topic of this article.
What is Quantum Computation?
Thanks to the quantum phenomena (1) superposition (unlike classical bits that are either 1 or 0 at any time, a quantum bit or qubit can have two states at a time (yes, that’s weird!)) and (2) entanglement (multiple qubits at a distance can have guaranteed the same or the opposite state at the same time (‘spooky action at a distance’ said Einstein, yes, weird again!)), quantum computers can try multiple potential solutions at a time and so can be faster at it than classical computers that have to go through each, be it educated guess sequentially.
So far for the quantum computation theory. But has it been demonstrated in practice yet? You can read here and here that the answer to the question of ‘quantum supremacy’ more modestly called ‘quantum advantage’ is still subject to controversy.
Why Sudokus?
To pitch classical versus quantum computers on a concrete problem, we took the sudoku puzzle because (1) it is the most widely known optimisation problem around and (2) it is easy to scale up to arbitrary large and so arbitrary difficult versions so that an advantage on either system, classical or quantum, if at all present, can be blown up to arbitrarily large proportions.
The relevance of solving sudokus in an automated or even optimally fast way, to the world at large, with all its serious problems, may look futile. However, sudokus are just one instance of a large class of problems called ‘optimisation problems’ (more specifically, Mixed Integer Programming or ‘MIP’ problems). In that class, more practically relevant problems reside, like supply chain optimisation problems, train schedule optimisation problems, nurse scheduling problems, automated energy market trading problems and more. Today, solving some of these problems can take hours or days on classical computers and so, prevent their application in real time, as well as consume large amounts of electricity. Some think that, soon, or already today, quantum computers could/can save time and energy in solving these problems.
So let’s get to our challenge now.
1. We solved Super Sudokus on a Classical Computer
As mentioned, with a classical MIP approach, solving times rise quickly with n getting larger. The model used and the calculation times we obtained, are reported in this article from March 2020, for the Sudoku sizes n in [1, .. 7]. We redid the experiments, now for n in [1,..8], with an updated version of the MIP solver Gurobi, v9.1.0. The newly experienced solver times are given in Table 1 below.
The top half of Table 1 shows Sudokus without any squares having preset values. So the solver has maximal freedom, and any accidental classical versus quantum bias due to ‘specifically hard/easy square value choosing’ can be avoided.
In the MIP model, the solver has to decide on the 0 or 1 value of exactly n⁶ bits and also must take care that 4 types of constraints are satisfied, each of them occurring n⁴ times. (See below for the model variables and constraints that will make this clear.) This is why the time in the solve_model column increases quickly with increasing n.
The lower half of Table 1 reports calculation times for sudokus of the same sizes, but with some squares having pre-filled values. Compared with no preset values above, this is sometimes harder (e.g. for n=6) but sometimes easier (e.g. for n=7) and sometimes extremely hard. Yes, after 80 hours there was no result yet while for the case with no pre-filled squares case for n=8 it really only took 6s, no reporting mistake there!
2. We Solved Super Sudokus on Quantum Computers
We took a Quantum Computing course from IBM and one from D-Wave and the D-Wave implementation suite of hardware (and software) very clearly looked most promising, the main factor being that D-Wave chips using the quantum annealing model, contain up to 5000 qubits while the IBM, Google and other quantum computers which use the quantum gate model are at around 50 to 60 qubits.
We then looked around if people had attacked the Meta-Sudoku problem before us, and a few had.
2.0 Prior Work vs Our Work Here
2.0.1 Online Code From D-Wave: Works up to Sudokus of size (3*3)*(3*3)
D-Wave itself has published code online at https://github.com/dwave-examples/sudoku. We tried this out and it correctly solves a (3*3)*(3*3) sudoku case that is given in the supplied problem.txt file. The code is written in a way that depends on the parameter n, so it can theoretically solve any (n*n)*(n*n) sudoku. However, when we tried it out for n>3, it quite consistently returned invalid solutions.
2.0.2 Article explaining Quantum Supremacy compared to Classical Simulated Annealing on Sudokus of size (2*2)*(2*2)
At https://towardsdatascience.com/solving-sudoku-with-a-quantum-power-up-5bb4f64f3944, Julian Hatski described how quantum annealing on a D-Wave pure quantum (so QPU, and non-hybrid) machine is faster by about a factor of 5 than a certain classical simulated annealing algorithm on a (2*2)*(2*2) sudoku. He experienced, as did we, that the standard (3*3)*(3*3) sudoku cannot be automatically embedded yet on a QPU.
2.0.3. Article showing an Experiment with a 4x4 Sudoku
This Medium article by Eyan Fosythe explains how a Sudoku with 16 squares (so a (2*2)*(2*2) one) can be solved with a quantum computer based on the gate model. The solution returned is not always valid and scaling it up to larger sudokus is currently impossible due to lack of more qubits in quantum computers based on that model. Eyan mentions that the classical complexity of a sudoku with m undecided squares is of order O(N^m) and a non-trivial theoretical result that quantum computing has been shown to be of complexity order O(log(n)), so demonstrating exponential speedup. Note that it would not depend on m at all! Academic publications here and here show the same practical results for boards of at most 16 squares.
2.0.4 This article: Scaling up and Using Hybrid Solver, Comparing to MIP Gurobi solving times
We intend to use D-Wave systems for larger Sudokus, up to (8*8)*(8*8) and to compare solving times to what we got on classical computers with Mixed Integer Programming (MIP) solvers. We use Gurobi, which in our experiments, was the fastest available solver for this kind of problems. There may be algorithms dedicated to solve Sudokus that get to results quicker than Gurobi with its algorithms dedicated to solve the general class of MIP problems, but they require a lot of development time. We think comparing a general Quantum solver (with QUBO models) to a general classical solver (with MIP models) is the most fair one from the following viewpoints.
(1) The development time of classical MIP models and quantum QUBO algorithms is similar (if you use the newest D-Wave ocean libraries with add_variable, add_linear_constraint).
(2) The financial cost of use is similar (classical 1 PC, 10 minutes of ‘cloud’ quantum time during the D-Wave course at 700 USD for the course).
(3) So if the above (1) and (2) are similar and the time to solution is better for quantum than classical, then there is a very practical and concrete quantum advantage.
2.1 Converting MIP constraints to Objective Function terms of a QUBO: General Principle
We have to cast the Sudoku MIP model into its Quadratic Unconstrained Binary Optimisation (QUBO) equivalent model form, which is the format needed to be able to program a d-wave quantum annealing machine.
2.1.1 Performing the Conversion Algebra
The QUBO format allows only an objective (function) and, unlike a MIP model, does not directly allow constraints. Constraints can be converted into terms of the objective function via the well known principle of Lagrange multipliers. So a say linear constraint f(b)=c, with b a vector of some of the binary variables in the MIP problem, is converted to a term 𝜆 (f(b)-c)². 𝜆 is a weight and can be chosen higher than the value of other terms in the objective, when the minimisation of (f(b)-c)² (so consequently, the near-satisfaction of f(b)=c) is more important than the minimisation of other objective function terms.
Working out the expression (f(b)-c)², when f(b) is linear in each b in b is easy. To save time, effort and avoid mistakes, one can use Wolfram Alpha for that. For example, the constraint b1+b2+b3=1 converts to (b1+b2+b3-1)² and can be expanded by WolframAlpha by this link. This returns the expanded form of the expression: b1² + b2² +b3² + 2 b1 b2 + 2 b1 b3 + 2 b2 b3 -2 b1 -2 b2 -2 b3+ 1.
Now, because any binary value b squared is equal to the binary variable b itself (since 1²=1 and 0²=0), we can simplify this down to (2 b1 b2 + 2 b1 b3 + 2 b2 b3) +(-b1-b2-b3) + 1, and one can see that this expression is composed of quadratic terms in the binary variable cross terms plus linear terms in the binary variables.
For generality, note that also constraints that are inequalities i.o. equalities, so of the form f(b) ≤ c can be converted to f(b) + c’ * b’ = c with b’ a newly introduced binary variable, and can then be converted to an objective function term (f(b)+ c’ * b’ -c)², which can be worked out similarly to an expression consisting of linear terms and cross product terms and then be converted to a QUBO problem with one new and fresh variable b’.
So any linear constraint in a MIP problem can be converted into linear and quadratic terms in the objective function of its QUBO form.
2.1.2 Constructing the QUBO matrix
The quantum computer needs a QUBO in the form of a matrix. (This is in fact similar to the MPS format from the early days of MIP solver APIs.) If q is the amount of qubits used in the QUBO problem, the QUBO matrix is of dimension q*q. Its diagonal contains the coefficients of the linear terms bi and the upper right triangle contains the coefficients of the terms bi bj, where i<j. (In fact one can also supply a symmetric matrix, where the i>j elements will be added to the i<j matrix elements before the D-Wave machine consumes it.)
Now, with the understanding of this QUBO matrix encoding, loosely formulated, the physical workings of a QUBO based quantum annealing computer are as follows. Nature in general, and so also the physics of the quantum system strives for minimal energy, which corresponds to minimising the objective function. A positive QUBO coefficient wants the binary (+bi) or binary product (+bi * bj) with this coefficient to be zero since the whole objective represents energy of the encoded system and will be minimised, while a negative coefficient will push the thing following (bi or bi * bj) to be 1. Coefficients that are larger than others will push harder for their terms than these other terms are pushed. As such (1) individual terms (of shape bi, so on the diagonal) can be pushed up or down, and cross terms (bi * bj, off-diagonal) (representing constraints) can pushed up or down, and the relevant importance can be ‘weighted’.
More formally, the following expression representing the objective function is minimised.
‘I’ is the set of qubit indices. So #(I) is the number of qubits in the problem, which is n⁶. The bi are the qubits, the value of which has to be determined. The ci and cij are the fixed coefficients. The QUBO matrix consists of the ci and cij constants, representing the encoded problem.
Remember that the MIP objective and each constraint in the MIP problem is transformed, as described above, in an expression that contributes to the cii and cij values. Per expression, weights can be chosen.
The question can arise whether, per constraint, correct weights can always be found, and it seems that it has been proven that this is not the case. However, in practice, it seems that for quite some problems, reaching satisfaction of all constraints can be attained by some experimenting with these weights. For the sudokus for n≤7 without pre-filled squares, we found out that this works out fine with all weights (Lagrange parameter values) of all constraints set equal. For the cases with pre-filled squares, we had to set the weights of the fixed squares satisfaction constraints to a higher number than the weight of all the other constraints for a n in [1,2] on a QPU and do the opposite for n in [3,..7] on a hybrid D-Wave system. Figuring this out takes experimentation time and this is something that is not needed with a classical MIP solver where constraints don’t need weights specified by the user.
2.2 Converting MIP constraints to Objective Function terms of a QUBO: Sudoku Application
Summarized, we get for the MIP model for a #(R) rows * #(C) columns-sized Sudoku to be filled out with the #(N) numbers from N=[1..#N], where #N=#R=#C and D=#(S) where D²=#N:
The top 3 constraints together ensure that each of the numbers in N only occurs once per row and column. The fourth constraint enforces uniqueness of each number in N occurs only once per subsquare. The last constraint imposes the preset values per square, if there are any.
For the sub-board constraint, in a standard sudoku D=3 and S=[0,1,2].
For the preset squares, f(r,c) is a function that maps some (row,column) indices to the value that square should contain. Not all rows and squares will have preset values, since that would be too easy to solve. dom(f(r,c)) is the domain of f, so the set of (r,c) couples for which we have preset values. For the quantum version, the model only contains binaries and the last constraint is omitted. The integer values (n+1) to fill in the squares are (trivially) computed post-optimisation in a classical way.
Note that the indices of b here are still relatively simple formulae thanks to the fact of having offset 0 for rows and columns.
For the top 4 constraints, we can see that we have n⁴ (with n=3 for a classical 9*9 Sudoku) of each type because #(R)=#(C)=#(N)=#(S)*#(S)=n and for example for the first type, we have a number of #(R)*#(C)=n²*n² of them. So in total we have 4 n⁴ constraints, excluding the square-preset constraints.
We can indeed see that all the above constraints are convertible to QUBO form with the general principles described in section 2.1.
After applying these rules in a pretty mechanical way, which is reflected by the following Python code, we arrived at our QUBO model. We believe it the process described above applied to a sudoku will be clear from the comments in the code.
2.3 Python code constructing the QUBO
The below code contains all functions to solve Sudokus both classically and quantumly.
See the fully executable (and executed) notebook at https://gist.github.com/PeterSels/0ef625b0340598fde4f2f2b6f3a8e5d6
All cases executed split into 4 categories: (without/with preset squares) * (with Gurobi/with D-Wave). To keep code size and execution times small and modular, we put each of these 4 cases in a separate notebook that includes the above notebook via the command %run QuantumSudokuLib.py.
(1) No Pre-filled Squares, Classical Gurobi:
For the full gist of the test runs, so including both code and test result output, see https://gist.github.com/PeterSels/139f8477133c4604dcac8399e231788e .
We could run all the n in [1,..8] experiments here.
(2) No Pre-filled Squares, Quantum D-Wave
See https://gist.github.com/PeterSels/8364118a8f3d1275a69eafe290ff1d96 for the notebook gist of the test runs.
For this case, we could do our experiments for n in [1..5], but not redo our n in [6, 7] experiments we ran earlier. We ran out of quantum time on our account.
But we rely on our output log of Python code (prior to plugging it in a notebook). We refer to these 3 logs here. The code in the notebook should regenerate similar results if sufficient quantum time is available.
Our logs for n=6, no preset squares, quantum
We published the whole log here on GitHub so that you can verify that it is correct. The summary is this:
(ocean) $ time python helper_functions.py sudokuDim6.json dwave
...
is_correct check done.
Sudoku correctly solved. :)
tests: ['solver'] done
real 22m4.531s
user 0m20.908s
sys 0m2.179s
(ocean) $So one can see from “real 22m4.531s”, that this took just over 22 minutes for n=6. The “Sudoku correctly solved. :)” indicates that the solution returned is checked (on a classical computer) and indeed valid.
We did the experiment again, again for n=6 and achieved a valid solution again in 22 minutes 50 seconds. The full log is present here on GitHub, with summary:
(ocean) $ time python helper_functions.py sudokuDim6.json dwave
...
is_correct check done.
Sudoku correctly solved. :)
tests: ['solver'] done
real 22m20.279s
user 0m23.123s
sys 0m2.560s
(ocean) $Our log for n=7, no preset squares, quantum
The full log is here on GitHub, with summary:
(ocean) $ time python helper_functions.py sudokuDim7.json dwave
...
Sudoku correctly solved. :)
tests: ['solver'] done
real 44m28.992s
user 1m27.995s
sys 0m7.025s
(ocean) $So the “real 44m28.992s” indicates it took 44 and a half minutes for n=7. The “Sudoku correctly solved. :)” indicates that the solution returned is checked (on a classical computer) and indeed valid.
We also had an older experiment for n=8. For that one, after (only!) 20 minutes and 50 seconds, we received an invalid sudoku solution from the hybrid method. Note that this is not a mistake in hardware or software, but the (sometimes multiple) solutions returned by D-Wave systems are solutions with low energy but it might not always be the one with the absolutely lowest energy, since the physical process of quantum sampling is a stochastic (and so a non-deterministic) one. ¹ A solution with an energy that is not the absolutely minimal one, corresponds to one where some objective function terms will not be totally minimised (to zero), corresponding to ‘constraints’ that are ‘slightly violated’.
A remedy to get around this issue of invalid solutions is to readout multiple times in one run or to rerun the whole process. However, we ran through our D-Wave quantum developer and training time budgets to do just this.
(3) Some Pre-filled Squares, Classical Gurobi:
See https://gist.github.com/PeterSels/104002759f94ea45c1dd8684d98f1d27 for the notebook gist of the test runs.
We could run all the n in [1,..7] experiments here. But the n=8 case gave no result yet even after 80 hours.
(4) Some Pre-filled Squares, Quantum D-Wave
See https://gist.github.com/PeterSels/ff076c92e8d006ae95d4605ab818a0ef for the notebook gist of the test runs.
We could run the cases n in [1..5], but for n in [6,7,8], we ran out of quantum time on our accounts.
2.4. Fitting/Embedding the QUBO on a D-Wave system
Our MIP model as well as our quantum QUBO model use n⁶ (qu)bits for a sudoku of (n*n)*(n*n) squares. This is the case because for each square, we need a number of between 1 and n*n. So we construct a full (n*n)*(n*n) sudoku-plane for each number N in 1 to n*n that holds a boolean set to true in a particular square if the number N is decided to appear in that position in the final sudoku. So, since 4⁶ = 4096, a QPU with 5000 qubits should be able to hold all 4096 qubits for a Sudoku of size (4*4)*(4*4) and certainly one of size (3*3)*(3*3), which only requires 3⁶ = 729 qubits. However the function that has to find the embedding of the problem instance graph onto the D-Wave quantum chip graph seems to fail at that, which is why we have to use the D-Wave hybrid method for n≥3.
2.5 Results from our Quantum Sudoku Experiments
Quantum calculation times for the same range of sudoku sizes as for the classical cases are shown in the below table in column ‘D-Wave total (s). Note the ‘QPU’ (for pure Quantum Processing Unit, for n in [1,2]) and the ‘Hybrid’ for n≥3 for an approach where the problem is decomposed and was then mapped to a combination of classical and quantum computers.
For each value of n, the total D-Wave time is the sum of the columns model_build, find_embedding and model_solve. Most notably is that at low values of n the quantum solver is slower and for n ≥ 6 it gets faster. The quantum calculation time reduction starts there. At n=6, we reach a time reduction of 31.2% and at n=7 one of 82.9%, a clear and large quantum advantage.
We report the files where you can verify the results we obtained in table 3 below.
3. Classical vs. Quantum Time Visualised
To get a better feel for proportions, we put our results into some charts in the Figure below. The left two charts have linear time in seconds on the vertical axis, and the right pair has a logarithmic time scale on the vertical axis. The top half represents the sudokus without preset squares and the bottom half the ones with preset squares. Blue means classical MIP Gurobi, red means quantum/hybrid QUBO D-Wave.
On the left upper corner, since time is in plain seconds, proportions are kept. Note that the number of variables when going from n=6 to n=7 goes from 6⁶=46656 to 7⁶=117649 so roughly rises by a factor 2.5 and the number of constraints rises from 4 * 6⁴ =5184 to 4 * 7⁴=9604, so by about a factor 1.8. So, loosely speaking, we could say the the complexity of the problem rises roughly by a factor 2.5 * 1.8 = 4.7.
Gurobi’s solving time goes from 1924s for=6 to 15596s for n=7, so rises by a factor 8.1. So roughly, 2 problem complexity doubling gives rise to 3 solving time doublings on a classical system.
D-Wave’s hybrid setup solving time for the same 2 doublings of problem complexity goes from 1324s for n=6 to 2669s for n=7, so almost exactly only 1 solving time doubling on a quantum system.
The charts on the right show a vertical axis on a logarithmic scale to allow for more resolution for the lower times for lower n.
Note that, independent of whether the vertical scale is in seconds or log(seconds), the horizontal axis is in fact highly compressed. Indeed consider that, if we would have put the horizontal axis in terms of fundamental problem complexity, we should have set a scale proportional to #variables*# constraints ~ n⁶ * 4 n⁴ which is equivalent to n¹⁰. That would have resulted in a much more ‘wide and flat’ curve.
Note that the reported D-Wave hybrid times still have classical computation components which may still increase harder with increasing n than their quantum sibling components. Having a pure quantum device where an n≥8 sudoku would be quickly embeddable, would be the next attainable holy grail for quantum sudoku solving!
4. Conclusions
For the super/meta-sudoku problem, the D-Wave pure quantum chips (called Quantum Processing Unit or QPU) which, today, have up to 5000 qubits can embed small sudokus up to size n=2. However, for these small problems, no quantum advantage is achieved yet.
The D-Wave hybrid approach where problems are decomposed over classical and quantum chips can embed larger problems (n≥3) and the two largest ones we could try, for sudoku dimensions 6 and 7 do show a significant quantum advantage (31.2% and 82.9% respectively). To the best of our knowledge, this is the first time this has been demonstrated on a quantum machine for meta sudokus.
Note that finding the ‘embedding’ (function find_embedding in the code) of the problem graph onto the graph representing the quantum chip (or hybrid system) is done classically and also takes time, and that this time increases with size n. With increasing n, the hope is that the quantum subsystem of the hybrid method will not rise much or not at all and that the classical time component will not dominate the total time.
If we could arrange for additional quantum time on D-Wave machines, we could try sudokus of size n≥8 again on current systems. We are also curious to see what will happen when the classical embedding-finding function gets better (and fits more on a QPU) and when the D-Wave chips with more qubits will be released in the near future.
Footnotes:
- Classical MIP solvers like Gurobi are constructed to be deterministic. This means they will always return the same answer (solution if feasible) if the input is the same and the machine is the same (or equivalently, the number of cores, number of threads and operating system is the same). This does not mean that in selecting a certain problem instance, good luck or bad luck is ruled out. For example, running Gurobi v9.1.0 on our model for a sudoku of size (8*8)*(8*8), without any pre-filled squares, we discovered it was solved in 6 seconds all-in, while the smaller case of (7*7)*(7*7) without pre-filled squares took 15596 seconds.
- O(n¹⁰) rises more slowly with increasing n than O(exp(n)) does, but it rises quickly anyway.
5. Running Quantum D-Wave Machines from your own Computer … for Free
The full Python jupyter notebooks in this article are all available on GitHub here. pip install should install all libraries needed here because it will look into the supplied file requirements.txt which holds all requirements. You can then run the notebooks yourself by adapting have_gurobi and have_dwave to what software and licenses you have available. You can setup your own free leap account at D-Wave. With pip install gurobipy (included in pip install) you can install a version of Gurobi that solves small problems up to a number of variables or constraints. For bigger ones you will need a license.
The 5 notebooks gists, containing all code but also containing all interleaved output we generated, are available at the links mentioned above in this article.
And the best thing is, you can get a minute of D-Wave quantum time (which is a lot considering that small problems are solved in milliseconds) for free, renewed each month if you register as a developer and publish your code openly, say on GitHub or so.
So now you can avoid getting seizures from trying to solve sudokus yourself.
Peter Sels and Edward Mulier, Tuesday August 31st, 2021.
All code and gists on GitHub referred to here were created by the authors of this article and all tests described where run by them.
Full Disclosure: This article is not sponsored by D-Wave, nor did we receive any payments from D-Wave ever. We admit to being enthusiastic users of their technologies, since they enable us to solve big computational challenges faster.
Copyright © 2021 Logically Yours BV.
