# 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);``
``#cssbody{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.

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.

1. Select just one number for each cell.
2. Each row and each col has array of numbers from 1 to 9. Doubled number is no allowed.
3. 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 opta = 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.

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

https://github.com/minatoyuichiro/sudoku_js_ising