Solving the world’s most difficult SUDOKU problem using ising model on javascript
Introduction
In Japan recently to solve SUDOKU on ising model which is the basic structure of the quantum annealing machine was a little bit a fashion.
This time I will try implement this on javascript with SA algorithm
The world’s most difficult SUDOKU problem.
I was just heard about the world’s most difficult SUDOKU problem. This time we try to solve this problem on the program with SA.
http://www.independent.com.mt/articles/2012-07-29/news/introducing-the-worlds-hardest-sudoku-313820/
Interface
Let’s begin with constructing interface to make good design.
#html
<table class="table table-bordered" id="main"></div>#js
for(var i=0;i<9;i++){
var row = '<tr>';
for(var j=0;j<9;j++){
row += '<td></td>';
};
row += '</tr>';
$('#main').append(row);
}; $('td').height(10);#css
body{font-family:courier;}
.table-bordered{border:3px solid black;}
.table-bordered tr:nth-child(3n){border-bottom:2px solid black;}
.table-bordered td:nth-child(3n){border-right:2px solid black;}
Input interface
Now we have input number for better input of number.
<input type="number" pattern="\d*">
But it looks not good on smartphones.
setup some css texts…
input{-webkit-appearance:none;border:0;background:none;text-align:center;font-weight:bold;line-height:35px;}
Submit button
Next we prepare submit button to solve the SUDOKU problem…
Looks great.
Interface for the answer
The input text expressed on black text and now we have annealed answer with skyblue texts.
Let’s solve on ising model
We can solve this problem using graph coloring algorithm on ising model. The main step is simple. We have to think about constraint to solve on ising model.
- Select just one number for each cell.
- Each row and each col has array of numbers from 1 to 9. Doubled number is no allowed.
- Neighbor 9 blocks also have 1–9 array of numbers.
Select 1 from 9 on each cell
This time we are using totally 9*81=729 qubits
from blueqat import opt
a = opt.opt()
a.qubo = opt.sel(9,1)print(a.qubo)
#=>
[[-1. 2. 2. 2. 2. 2. 2. 2. 2.]
[ 0. -1. 2. 2. 2. 2. 2. 2. 2.]
[ 0. 0. -1. 2. 2. 2. 2. 2. 2.]
[ 0. 0. 0. -1. 2. 2. 2. 2. 2.]
[ 0. 0. 0. 0. -1. 2. 2. 2. 2.]
[ 0. 0. 0. 0. 0. -1. 2. 2. 2.]
[ 0. 0. 0. 0. 0. 0. -1. 2. 2.]
[ 0. 0. 0. 0. 0. 0. 0. -1. 2.]
[ 0. 0. 0. 0. 0. 0. 0. 0. -1.]]a.sa()
#=>
[0, 0, 0, 0, 0, 0, 1, 0, 0]print(a.J)
#=>
[[3.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5]
[0. 3.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5]
[0. 0. 3.5 0.5 0.5 0.5 0.5 0.5 0.5]
[0. 0. 0. 3.5 0.5 0.5 0.5 0.5 0.5]
[0. 0. 0. 0. 3.5 0.5 0.5 0.5 0.5]
[0. 0. 0. 0. 0. 3.5 0.5 0.5 0.5]
[0. 0. 0. 0. 0. 0. 3.5 0.5 0.5]
[0. 0. 0. 0. 0. 0. 0. 3.5 0.5]
[0. 0. 0. 0. 0. 0. 0. 0. 3.5]]
The upper matrix is called QUBO matrix on 01 binary number. The lower matrix is called Ising matrix on -1 and +1 of values. Now we have simple Ising matrix to realize selecting just 1 number from 9 numbers, implemented on qubits.
On javascript this is like the code below.
var L1 = 110;
for(var i=0;i<81;i++){
for(var j=0;j<9;j++){
h[i][j] = 3.5*L1;
}
}
//console.log(h);
for(var i=0;i<81;i++){
for(var j=0;j<9;j++){
for(var k=0;k<9;k++){
J[i][i][j][k] = 0.5*L1;
}
}
}
Cost function
To solve the problem we need ising model and cost function. Now we are using monte carlo simulation with metropolis method.
var anneal = setInterval(function(){
for(var kk=0;kk<50000;kk++){
var x = Math.floor(Math.random()*81);
var y = Math.floor(Math.random()*9);
//select 1 number from 9
Now we want to see the movement of the searching algorithm so using setInterval and have 50000 iteration on each temperature.
//select 1 number from 9
var dE = -2*h[x][y];
for(var i=0;i<9;i++){
if(i!=y){
dE += -2*J[x][x][y][i]*q[x][i]
}
}
Constraint on rows
Now we can think about constraint on rows using anti-ferro magnetic parameter on the qubits. With L2 hyperparameter the code is,
var L2 = 3;
//rows
for(var i=0;i<9;i++){
for(var j=0;j<9;j++){
for(var k=0;k<9;k++){
for(var l=0;l<9;l++){
J[j+i*9][k+i*9][l][l] += L2;
}
}
}
}
It’s simple. And evaluation of cost function is ,
//rows
for(var i=0;i<9;i++){
if(x<9 && x!=i){
dE += -2*J[x][i][y][y]*q[i][y]
}
}
for(var i=9;i<18;i++){
if(8<x && x<18 && x!=i){
dE += -2*J[x][i][y][y]*q[i][y]
}
}
for(var i=18;i<27;i++){
if(17<x && x<27 && x!=i){
dE += -2*J[x][i][y][y]*q[i][y]
}
}
for(var i=27;i<36;i++){
if(26<x && x<36 && x!=i){
dE += -2*J[x][i][y][y]*q[i][y]
}
}
for(var i=36;i<45;i++){
if(35<x && x<45 && x!=i){
dE += -2*J[x][i][y][y]*q[i][y]
}
}
for(var i=45;i<54;i++){
if(44<x && x<54 && x!=i){
dE += -2*J[x][i][y][y]*q[i][y]
}
}
for(var i=54;i<63;i++){
if(53<x && x<63 && x!=i){
dE += -2*J[x][i][y][y]*q[i][y]
}
}
for(var i=63;i<72;i++){
if(62<x && x<72 && x!=i){
dE += -2*J[x][i][y][y]*q[i][y]
}
}
for(var i=72;i<81;i++){
if(71<x && x<81 && x!=i){
dE += -2*J[x][i][y][y]*q[i][y]
}
}
The constraint on cols
This is almost the same as rows
//cols
for(var i=0;i<9;i++){
for(var j=0;j<9;j++){
for(var k=0;k<9;k++){
for(var l=0;l<9;l++){
J[j*9+i][k*9+i][l][l] += L2;
}
}
}
}
And cost evaluation
//cols
for(var i=0;i<9;i++){
if(x%9==0 && x!=i*9){
dE += -2*J[x][i*9][y][y]*q[i*9][y]
}
}
for(var i=0;i<9;i++){
if(x%9==1 && x!=i*9+1){
dE += -2*J[x][i*9+1][y][y]*q[i*9+1][y]
}
}
for(var i=0;i<9;i++){
if(x%9==2 && x!=i*9+2){
dE += -2*J[x][i*9+2][y][y]*q[i*9+2][y]
}
}
for(var i=0;i<9;i++){
if(x%9==3 && x!=i*9+3){
dE += -2*J[x][i*9+3][y][y]*q[i*9+3][y]
}
}
for(var i=0;i<9;i++){
if(x%9==4 && x!=i*9+4){
dE += -2*J[x][i*9+4][y][y]*q[i*9+4][y]
}
}
for(var i=0;i<9;i++){
if(x%9==5 && x!=i*9+5){
dE += -2*J[x][i*9+5][y][y]*q[i*9+5][y]
}
}
for(var i=0;i<9;i++){
if(x%9==6 && x!=i*9+6){
dE += -2*J[x][i*9+6][y][y]*q[i*9+6][y]
}
}
for(var i=0;i<9;i++){
if(x%9==7 && x!=i*9+7){
dE += -2*J[x][i*9+7][y][y]*q[i*9+7][y]
}
}
for(var i=0;i<9;i++){
if(x%9==8 && x!=i*9+8){
dE += -2*J[x][i*9+8][y][y]*q[i*9+8][y]
}
}
It’s quite easy.
By adjusting 2 hyperparameters we can solve the ising model problems.
The constraint on 9 neighbors
We think about 9 neighbors also.
//9 neighbors
for(var j=0;j<3;j++){
for(var k=0;k<3;k++){
for(var l=0;l<9;l++){
J[0+k*3+j*27][10+k*3+j*27][l][l] += L2;
J[0+k*3+j*27][11+k*3+j*27][l][l] += L2;
J[0+k*3+j*27][19+k*3+j*27][l][l] += L2;
J[0+k*3+j*27][20+k*3+j*27][l][l] += L2;
J[1+k*3+j*27][9+k*3+j*27][l][l] += L2;
J[1+k*3+j*27][11+k*3+j*27][l][l] += L2;
J[1+k*3+j*27][18+k*3+j*27][l][l] += L2;
J[1+k*3+j*27][20+k*3+j*27][l][l] += L2;
J[2+k*3+j*27][9+k*3+j*27][l][l] += L2;
J[2+k*3+j*27][10+k*3+j*27][l][l] += L2;
J[2+k*3+j*27][18+k*3+j*27][l][l] += L2;
J[2+k*3+j*27][19+k*3+j*27][l][l] += L2;
J[9+k*3+j*27][1+k*3+j*27][l][l] += L2;
J[9+k*3+j*27][2+k*3+j*27][l][l] += L2;
J[9+k*3+j*27][19+k*3+j*27][l][l] += L2;
J[9+k*3+j*27][20+k*3+j*27][l][l] += L2;
J[10+k*3+j*27][0+k*3+j*27][l][l] += L2;
J[10+k*3+j*27][2+k*3+j*27][l][l] += L2;
J[10+k*3+j*27][18+k*3+j*27][l][l] += L2;
J[10+k*3+j*27][20+k*3+j*27][l][l] += L2;
J[11+k*3+j*27][0+k*3+j*27][l][l] += L2;
J[11+k*3+j*27][1+k*3+j*27][l][l] += L2;
J[11+k*3+j*27][18+k*3+j*27][l][l] += L2;
J[11+k*3+j*27][19+k*3+j*27][l][l] += L2;
J[18+k*3+j*27][1+k*3+j*27][l][l] += L2;
J[18+k*3+j*27][2+k*3+j*27][l][l] += L2;
J[18+k*3+j*27][10+k*3+j*27][l][l] += L2;
J[18+k*3+j*27][11+k*3+j*27][l][l] += L2;
J[19+k*3+j*27][0+k*3+j*27][l][l] += L2;
J[19+k*3+j*27][2+k*3+j*27][l][l] += L2;
J[19+k*3+j*27][9+k*3+j*27][l][l] += L2;
J[19+k*3+j*27][11+k*3+j*27][l][l] += L2;
J[20+k*3+j*27][0+k*3+j*27][l][l] += L2;
J[20+k*3+j*27][1+k*3+j*27][l][l] += L2;
J[20+k*3+j*27][9+k*3+j*27][l][l] += L2;
J[20+k*3+j*27][10+k*3+j*27][l][l] += L2;
}
}
}
It is also very simple and easy to set constraint.
if(x==0 || x==3 || x==6 || x==27 || x==30 || x==33 || x==54 || x==57 || x==60){
dE += -2*J[x][x+10][y][y]*q[x+10][y];
dE += -2*J[x][x+11][y][y]*q[x+11][y];
dE += -2*J[x][x+19][y][y]*q[x+19][y];
dE += -2*J[x][x+20][y][y]*q[x+20][y];
}
if(x==1 || x==4 || x==7 || x==28 || x==31 || x==34 || x==55 || x==58 || x==61){
dE += -2*J[x][x+8][y][y]*q[x+8][y];
dE += -2*J[x][x+10][y][y]*q[x+10][y];
dE += -2*J[x][x+17][y][y]*q[x+17][y];
dE += -2*J[x][x+19][y][y]*q[x+19][y];
}
if(x==2 || x==5 || x==8 || x==29 || x==32 || x==35 || x==56 || x==59 || x==62){
dE += -2*J[x][x+7][y][y]*q[x+7][y];
dE += -2*J[x][x+8][y][y]*q[x+8][y];
dE += -2*J[x][x+16][y][y]*q[x+16][y];
dE += -2*J[x][x+17][y][y]*q[x+17][y];
}
if(x==9 || x==12 || x==15 || x==36 || x==39 || x==42 || x==63 || x==66 || x==69){
dE += -2*J[x][x-8][y][y]*q[x-8][y];
dE += -2*J[x][x-7][y][y]*q[x-7][y];
dE += -2*J[x][x+10][y][y]*q[x+10][y];
dE += -2*J[x][x+11][y][y]*q[x+11][y];
}
if(x==10 || x==13 || x==16 || x==37 || x==40 || x==43 || x==64 || x==67 || x==70){
dE += -2*J[x][x-10][y][y]*q[x-10][y];
dE += -2*J[x][x-8][y][y]*q[x-8][y];
dE += -2*J[x][x+8][y][y]*q[x+8][y];
dE += -2*J[x][x+10][y][y]*q[x+10][y];
}
if(x==11 || x==14 || x==17 || x==38 || x==41 || x==44 || x==65 || x==68 || x==71){
dE += -2*J[x][x-11][y][y]*q[x-11][y];
dE += -2*J[x][x-10][y][y]*q[x-10][y];
dE += -2*J[x][x+7][y][y]*q[x+7][y];
dE += -2*J[x][x+8][y][y]*q[x+8][y];
}
if(x==18 || x==21 || x==24 || x==45 || x==48 || x==51 || x==72 || x==75 || x==78){
dE += -2*J[x][x-17][y][y]*q[x-17][y];
dE += -2*J[x][x-16][y][y]*q[x-16][y];
dE += -2*J[x][x-8][y][y]*q[x-8][y];
dE += -2*J[x][x-7][y][y]*q[x-7][y];
}
if(x==19 || x==22 || x==25 || x==46 || x==49 || x==52 || x==73 || x==76 || x==79){
dE += -2*J[x][x-19][y][y]*q[x-19][y];
dE += -2*J[x][x-17][y][y]*q[x-17][y];
dE += -2*J[x][x-10][y][y]*q[x-10][y];
dE += -2*J[x][x-8][y][y]*q[x-8][y];
}
if(x==20 || x==23 || x==26 || x==47 || x==50 || x==53 || x==74 || x==77 || x==80){
dE += -2*J[x][x-20][y][y]*q[x-20][y];
dE += -2*J[x][x-19][y][y]*q[x-19][y];
dE += -2*J[x][x-11][y][y]*q[x-11][y];
dE += -2*J[x][x-10][y][y]*q[x-10][y];
}
Main loop
And the main setting all done. Just thinking about the metropolis and monte carlo.
dE *= q[x][y];
if(Math.exp(-dE/T) > Math.random()){
q[x][y] = -q[x][y]
}
};
//display of numbers
for(var i=0;i<81;i++){
for(var j=0;j<9;j++){
if(q[i][j]==1){
$('#answer_'+Math.floor(i/9)+'_'+i%9+'_'+(j+1)).html(j+1);
}else{
$('#answer_'+Math.floor(i/9)+'_'+i%9+'_'+(j+1)).html('');
}
}
}
T*=0.9;
//console.log(T);
At the end of calculations.
If it is solved, it will finish. If not solved it continue solving until solved.
if(T<0.7){
var check = 0;
//horizontal check
for(var m=0;m<9;m++){
var E1 = 0
for(var i=0+9*m;i<9+9*m;i++){
for(var j=0;j<9;j++){
if(q[i][j] == 1){
E1 += j;
}
}
}
console.log(E1);
if(E1 != 36 ){
check++;
}
}
//vertical check
for(var m=0;m<9;m++){
var E1 = 0
for(var i=0;i<9;i++){
for(var j=0;j<9;j++){
if(q[i*9+m][j] == 1){
E1 += j;
}
}
}
console.log(E1);
if(E1 != 36 ){
check++;
}
}
if(check !=0){
T=100;
for(var i=0;i<81;i++){
for(var j=0;j<9;j++){
q[i][j] = Math.floor(Math.random()-0.5)*2+1;
}
}
//clearInterval(anneal);
}else{
clearInterval(anneal);
};
}
},1)
})
Done
That’s all. Let’s solve actual problem on puzzle book.
Let’s just solve.
The input is like the reference, just input the number on the app and submit SUDOKU button
Solved
Finally it is solved on iPhone.
https://i.gzn.jp/img/2010/08/22/hardest_sudoku/answer.jpg
The code
I put the demo site on github. Please try to use the sudoku program from demo site below.
https://minatoyuichiro.github.io/sudoku_js_ising/
And the whole code is shared on github.