test

はじめに

イジングモデルで数独を解くのが少し前に流行りましたが、僕も解いてみます。ただ、インターフェイスをつけてやりたいです。やり方は多分graph coloringを応用して、各マスの相互作用をつければいいだけと思いますが、SAでやってみます。イジングモデルで解ければ実質的に量子コンピュータでも解けますので、今後への道筋をつけたいと思います。

参考

https://quantum.fixstars.com/introduction_to_quantum_computer/applications/sudoku/

https://qiita.com/gyu-don/items/0c59423b12256821c38c

さまざまな参考があります。今回はなんとなく上記をの話を聞いていたので同様適当に作ってみます。

世界一難しい数独

下記のような記事を教えてもらいました。とりあえずこれが解ければ大丈夫かなと頑張ってみます。

数学のエキスパートが3ヶ月かけて作成した「世界一難しい数独」

https://gigazine.net/news/20100822_hardest_sudoku/

見込み

もちろんD-Waveを使ってもできると思いますが、まずは初めてなので手元のSAというアルゴリズムでやってみます。応用でD-Waveや他のイジングマシンへは応用できますので、コードを参照してやってみてください。

シミュレーションはモンテカルロを使いますが、接続の節約のために少しモンテカルロ内のエネルギーコスト関数の評価部分には分岐や条件がつく方向性ですが、それは今後リファクタリングします。

インターフェイス

とにかく見た目をシンプルにしたい＆誰でも簡単に使えるようにしたいとhtml/jsでテーブル作ってやります。サクッとできました。table要素を普通に使いました。概要はこんな感じです。

#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でやればいいと思います。上記のjsでinputをappendします。キーボード入力がめんどいのですが、このようにすれば数字が入るようです。

できたと思ったらスタイルがダメでした。

ちょっとがんがったらだいたい良くなりましたとりあえず良さそうです。iOSです。

input{-webkit-appearance:none;border:0;background:none;text-align:center;font-weight:bold;line-height:35px;}

まだ解けるかどうかもわからないのに、解くボタンだけ先に用意しておきます。ボタンを増やすのが嫌なので、オシャレにSUDOKUの文字を解くボタンにしちゃいます。

正直大満足です。

答えを乗せるレイヤーをinputの後ろに入れておきます。

文字はオシャレにスカイブルーにしました。やはりみためがわかりやすいとテンション上がります。

いよいよイジング解く

解き方を考えてみます。graphcoloring問題というのを考えると、

１、９つの数字から１ますは１つ選ぶ制約条件

２、行と列はそれぞれ異なるものが入るという条件

３、９つの正方形ブロックは異なるものが入るという条件

あたりかなと思います。（違ってたら指摘してください。。。）

正直イジングのソフト制約を満たすのが大変なので、検算のシステムを入れながら量でランダム法に近く無理やり解いてしまいます。。。

やってみます。セルにidを振って管理してみます。まずは１マスごとに９量子ビットを使っています。まずは1から9までの数字のうち１つだけが選ばれる条件を考えてみます。こちらは良くやっている手法で、すべての量子ビットを足して-1をして二乗をして展開します。今回はQUBOはツールを使って、、、

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

QUBOは対角に-1、非対角に+2をいれて、Jijとしてのイジング変換後は対角に+3.5、非対角に+0.5でした。この通り実装してみます。その結果まずは簡単に可視化したモンテカルロを見てみます。とりあえずできました。

これで9C1で１マス分のハミルトニアンができました。これを81個用意します。ハミルトニアンはイジングmatrix上で上三角行列で大きな階段状のものができると思います。コードはjsでかきましたが、下記のようになりました。

ちょっと僕の実力不足でJijの対角項に本当はhを対応させないといけないのですが、頭が追いつかなくて別実装でお願いします。。。

ちなみに今回はきちんとイジングmatrixに直して計算をします。

 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;↵   }↵  } ↵ }↵

コスト関数の評価

コスト関数の評価部分も並行して進めました。構造はちょっと頭があとで混乱するといけないので、81マスをきめて、その中に連想配列で1から9までの数字を持つことにしました。なので、量子ビットの構造はq[i][j]でiがマス目、jが数字に対応します。今回は全部で81*9=729量子ビットを使用します。

近傍９量子ビットは完全結合ですが、横方向、縦方向、近傍９ますの評価は接続を制限できるので全体としては完全結合ではありません。

 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);↵   //9こから1↵

今回は動きが面白そうなので、forループではなく、setIntervalで計算をしながら表示も操作します。iterationは50000回で少し多めにして精度を確保しています。xはマス目、yは数字を選んでます。先ほどのイジングmatrixは下記のように表現されました。あとで全体にq[x][y]をかけますので、そこは省略されています。

   //9こから1↵   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]↵    }↵   }↵

次に横方向の制約条件

横方向の制約条件をかけます。これはシンプルです。行ごとに量子ビットを取り出して同じ数字同士に対応する量子ビットの結合をanti-ferromagneticという反磁性のパラメータを与えます。

調整パラメータを導入しながら、

 var L2 = 3;↵ //行ごとに同じ数字が出ないように調整↵ 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;↵    }↵   }↵  } ↵ }↵

簡単にできます。これで行ごとに同じ数字が出ません。1行目から9行目まで順番に設定していきます。

計算リソース（能力の限界も）は限られているので、コスト関数の評価部分は下記のように場合分けをしました。ルールは決まっているので、あとでリファクタリングして綺麗にできそうです。

//行は違う数字↵   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]↵    }↵   }↵

同様に縦方向の制約

同様に縦方向の制約もつけました。これも同じようにさっくりできました。

 //列ごとに同じ数字が出ないように調整↵ 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;↵    }↵   }↵  }↵ }↵

また、コスト評価部分は、

//列は違う数字↵   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]↵    }↵   }↵

特に難しいところはありません。

ここまで縦横の制約条件をつけられましたので、動かしてみるときちんと制約を満たして動かせましたが、まだ問題は解けませんでした。結構頑張ったのですが、やはりもうちょい制約条件が足りないようです。

ちなみにパラメータ調整は比較的うまくいき、コスト関数を評価しながらパラメータを２種類調整することでうまく問題にフィットするパラメータを探すことができました。

近傍９ますの制約条件

近傍９マスの中で数字がダブってはいけないという制約条件を最後につけます。これをつけることで結構綺麗に数独が解けました。条件は素直に近傍を評価します。

 //周辺9マス↵ 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;↵ }↵ }↵ }↵

三個飛ばしずつで、対応する量子ビット同士にanti-ferroを設定します。

コスト関数評価部分は、

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];↵  }↵

うまい方法を考えるのが面倒だったので力技でやりました。近傍９マスのうちのどの位置に量子ビットがあるかで、近傍９マス内部での相対位置で評価しました。これなら特に苦労せずできます。

あとはループ

あとはメトロポリス法を使ってサクサク評価してモンテカルロするだけです。ここはシンプルです。以外ですが、温度の減衰は結構ゆるい設定で平気でした。設定値はどちらかというと一回の並行状態を保つiterationの回数を増やした方が精度が上がりました。精度は最終の解のエネルギーコスト関数で評価しました。

   dE *= q[x][y];↵   if(Math.exp(-dE/T) > Math.random()){↵    q[x][y] = -q[x][y]↵   }↵  };↵ //数字の表示系↵ 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);↵

終了時の設定

終了しには幾らかの評価を入れながら、解けなかったら継続して再度解くというフローを入れました。これにより解けるまで頑張ってくれます。

if(T<0.7){↵   var check = 0;↵   //横チェック↵   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++;↵    }↵   }↵   //縦チェック↵   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)↵ })↵

できた

結構素直に制約条件入れてパラメータ調整しただけでうまく動くようになりました。試しに百円ショップでかったナンプレの初級編を３問連続で解かせてみましたが簡単に解けました。

もはや初級編は簡単にいきます。スマホで読み込んで、数字をポチポチ推していくことで誰でも簡単にできます。

早速世界で一番難しい問題を。。。

早速解いてみました。

こちらがその「世界一難しい数独」。ω-3脂肪酸のサプリメントを販売するEfamol社の依頼で、科学と応用数学の博士号を持つフィンランド人の環境科学者Arto Inkala博士が手がけた数独作成プログラムにより作成されました。

引用：https://gigazine.net/news/20100822_hardest_sudoku/

とのことですので、入力してみます。。。もうプログラムを作るのに疲れてパラメータ調整とかがめんどくて嫌だったので、初級編を解いた時のパラメータのままやってみました。。。zzz

iPhoneにいれて解いてみました。iPhoneでデモサイトを読み込みぽちぽち。。。解けました！時間が夜中の3:50です。。。zzz

答えはこちらにあるようですが、合っていると思います。

https://i.gzn.jp/img/2010/08/22/hardest_sudoku/answer.jpg

結構繰り返しをやっていたのでもしかしたらほぼランダムかも、、、と思いつつ組み合わせの数は結構な数なので、そう簡単に解けないと思いますので、解けているのかも。

コード

ふつうにGithubでコード公開してます＆デモ作りました。スマホでも動きます。下記が今のsudokuのコードがのっています。ぜひスマホで試してみてください。PCでもできます。

https://minatoyuichiro.github.io/sudoku_js_ising/

コードはgithubを作りました。拙いコードですが、もしかしたらリファクタリングするかもしれません。表示系を含んでいて最適化はしてないので、結構早くなると思います。

https://github.com/minatoyuichiro/sudoku_js_ising

世界で一番難しい問題も解けてしまったので、問題入力の手間がめんどいので、画像認識かなんかを導入したいなと思っています。もちろんイジングモデルで。その際に少しコードを綺麗にするかもしれません。

次回はD-Waveで解いてみたいと思います。

以上です！