@yuichirominato 2017.12.03更新 261views

GPUを用いた高速リアルタイムシミュレーテッド量子アニーリングシミュレータ

GPU SA SQA シミュレータ 量子アニーリング

はじめに

先日NVIDIA社のAIスタートアップパートナーのピッチ登壇がGTC JAPAN2017であり、量子コンピュータとGPUの可能性について話をしてきたところ反響が結構ありましたので、少し書いてみようと思います。

https://www.nikkei.com/article/DGXMZO24542480S7A211C1XY0000/

通常シミュレーテッド量子アニーリング(以下SQA)を実行していると量子モンテカルロ法やその他の方法で内部で何が起きているのかわかりません。
量子計算の状態を把握するためにシミュレータで内部の状況を確認することは実務上でアニーリングを使用するなどの場合には有用です。
特にハイパーパラメータなど人間が設定しないといけないものを設定する必要があるような問題では内部状態を確認しながらハイパーパラメータを設定することで、アプリケーションの開発速度を大幅にあげることができます。

SQAに関するハイパーパラメータ

量子アニーリングの近似式をプログラミングして実問題を計算します。E=−∑i,j∑k=1mJijmqi,kqj,k−∑i∑k=1mhimqi,k−T2∑i,j∑k=1mlncoth(ΓmT)qi,kqi,k+1E=−∑i,j∑k=1mJijmqi,kqj,k−∑i∑k=1mhimqi,k−T2∑i,j∑k=1mlncoth(ΓmT)qi,kqi,k+1

普段使用するハイパーパラメータはトロッタ数と呼ばれる並列数m、磁場の強さΓΓ、温度T(もしくは逆数のββです。問題依存は量子ビット間の結合荷重であるJijJij、量子ビット自体にかかる磁場項hihi。そして問題によっては問題の中で制約条件を調整する調整パラメータλλがでてきます。

また、温度や磁場の強さを減衰させるための係数ττや温度や磁場の初期値など多数のハイパーパラメータで構成されています。モンテカルロ法を使用する場合には、またモンテカルロループを回すためのループの回数などを設定します。

これら多数のハイパーパラメータを使用して問題を解いていきます。研究用途など十分な時間が取れる場合には十分なアニーリング時間をとって計算できますが、実務ではこの時間は極めて有限で精度と速度がトレードオフな現状ではパラメータ調整を行い、問題のサイズや性質に合わせて解く必要があります。

イジングの問題を解くためのノウハウを貯めるために

とにかくイジングにマッピングをして実機で解くのが手っ取り早いですが、現状(昔から)D-wave社との契約は時間あたり2000-3000ドル(うちらのチームは有用な研究のため課金されませんでしたが、、、)、年間数千万円かかりますし、無料で使いたい場合には米国のNASA Ames研究所とやりとりをすればできますが、ちょっと承認されるのに時間がかかります。
その際にシミュレータを開発し、独自に進めれば良いですが、できるだけ色々な問題を解くためには手元に優秀なシミュレータが欲しくなります。実際D-wave社もGPUやスパコンを使用した量子シミュレートを使用して開発を進めたりクライアントワークを進めているようです(直接聞きました、、、)。

シミュレータの活用

シミュレータを活用すると近似的なアルゴリズムですが同じような問題を解くことができます。ただ、量子モンテカルロ法は多大な並列数と十分なアニーリングタイムのためのモンテカルロ計算を必要とするため、大きなマシンパワーと時間が必要になります。どの分野でもモンテカルロの収束には時間がかかります。シミュレータの開発と活用には知識と人材と時間がかかります。そこで現在スパコンに加えてこれまでCPUベースで進めていたものをGPU拡張し始めました。

まずは大きな問題を解くためにシミュレート結果の可視化

量子モンテカルロで使用されるメトロポリスの更新はGPUのような並列処理は問題依存でスパースな問題では並列化しやすいですが、効率性を考えるとCPUできちんとモンテカルロを回してあげた方が最終的に得られる精度としては手っ取り早い気がします。CPUでも十分に早いので小さい問題ではCPUで解いてみるというので十分かと思います。GPUが必要になるのは一つはシミュレート結果の可視化です。可視化を行う際にトロッタ並列数は大きい場合にはより高速な画像処理が必要になりますし、画像処理用のコードを書く必要が出ます。

例えば

大きな問題を解こうとすると温度を大幅に下げる必要が出ます。また、問題の性質によっては逆に有限温度を上げて熱混合と量子アニーリングを並列に行った方が解きやすいという時もあります。実際D-wave社は20mKという動作温度は有限温度ですし、googleが2015年に発表した論文では量子アニーリングの最小値の探索は熱による励起と量子トンネル効果の両方でおこなわれるという説明も紹介されました。また、フォルクスワーゲンの渋滞緩和の問題では問題の制約条件の調整パラメータが使用されたり、NASAのカリフォルニアの緑の分類機の問題ではここでも特徴選択時の調整パラメータが採用されているように、たくさんの調整パラメータを日々調べる必要があります。

今回開発したシミュレータ

とても単純ですが、量子モンテカルロループをCPUで回しながら、フロントエンドの可視化表現をGPUで書いたところ、それまでCPU計算+CPU描画、もしくはGPU計算+GPU描画よりも効率的に同時にシミュレートをおこなうことができました。CPU量子モンテカルロ+GPU可視化描画だけでも、とても大きな問題を描画することができます。こちらは1次元のイジングモデルの量子アニーリングで縦方向は並列トロッタなので同じ問題のレプリカです。トロッタサイズが1000くらいでも問題なく高速にリアルタイムに可視化ができます。すでに現在D-waveを使用する前の前段階で論文の確認をするなどシミュレータは海外の研究者に活用してもらっています。

改めて

最近では次期D-wave社の量子コンピュータはキメラグラフをやめるかもという話もあります。キメラグラフとは何ぞやというところから、実際にキメラグラフをつくってインタラクティブに動くコードとしてシミュレータを実装してみます。

キメラグラフとは?

キメラグラフは量子ビットを4*4の井形に組んだ形状をしているグラフ形状のことで、1量子ビットは最大で交差する4つの量子ビットと隣接接続する2つの量子ビットの最大6量子ビットと接続する形状を言います。

前段の知識SAとSQA

D-wave型のシミュレータを作る際にはSAと呼ばれる古典的なシミュレーテッドアニーリングアルゴリズムでも、量子的な振る舞いを導入したSQAシミュレーテッド量子アニーリングのどちらでもできます。せっかくなので両方使ってやって見たいと思います。

SAはイジングモデルという物理モデルを使用して、磁性体スピンの向く方向でビットを表現します。取りうる値は{-1,1}で、通常のビットの{0,1}と少し違うのに注意が必要です。今回は隣接格子同士で接続のあるシンプルな形状をやって見ます。

実装の元となる式は下記の通りです。これをjavascript実装します。E=−∑i,jJijqiqj−∑ihiqiE=−∑i,jJijqiqj−∑ihiqisa.js

//量子ビット数をN=150に設定
var N = 150;

//初期温度を3に設定
var kT = 3;

//ここでは相互作用項のJijはすべて1として試してみる。磁場項hはめんどいので0でやってます。
var jij = 1;

//量子ビットを準備して初期値として{-1,1}のどちらかを入れる
var q = [];
for(var i=0;i<N;i++){
 q[i] = Math.floor(Math.random()-0.5)*2+1;
};

//準備は終わったので早速計算。ループを回してメトロポリス法で判定しながら量子ビットを反転シミュレート。
for(var l=0;l<550;l++){
 for(var k=0;k<N;k++){

  //モンテカルロを使用して一つ量子ビットを選ぶ
  var x = Math.floor(Math.random()*N);

  //隣接格子間のエネルギーの差をとる
  var dE = (jij*2*q[x]*[(N+x-1)%N] + jij*2*q[x]*q[(x+1)%N]);

  //メトロポリス法でビットフリップ判定
  if(dE<0 || Math.exp(-dE/kT)>Math.random()){
   q[x] = -q[x];
  };
 };

 //温度を下げる。今回はループの終点がないので適度なところで判断する。
 T = T*0.99;
};

次に量子アニーリングです。量子アニーリングはE=−∑i,j∑k=1mJijmqi,kqj,k−∑i∑k=1mhimqi,k−T2∑i,j∑k=1mlncoth(ΓmT)qi,kqi,k+1E=−∑i,j∑k=1mJijmqi,kqj,k−∑i∑k=1mhimqi,k−T2∑i,j∑k=1mlncoth(ΓmT)qi,kqi,k+1

を実装します。sqa.js

//量子ビット数は同じ150で実装
var N = 150;

//初期温度は0.02で固定
var kT = 0.02;

//ガンマの初期値10
var G = 10;

//トロッタ数50
var m = 50 ;

//Jijは固定で1
var jij = 1;

//量子ビットを初期化{-1,1}で。トロッタ数50(並列数50)なので、トロッタ数分だけ異なるqの配列が初期化されます。
var q = [];
for(var j=0;j<m;j++){
 q[j] = [];
 for(var i=0;i<N;i++){
  q[j][i] = Math.floor(Math.random()-0.5)*2+1;
 };
};

//早速計算
for(var l=0;l<2000;l++){
 for(var k=0;k<50000;k++){
  var y = Math.floor(Math.random()*m); //トロッタを選択
  var x = Math.floor(Math.random()*N); //量子ビットを選択

  //エネルギー差の式の一部を計算
  var dE = (jij*2*q[y][x]*q[y][(N+x-1)%N] + jij*2*q[y][x]*q[y][(x+1)%N])/m;

  var kk = G/kT/m;//よく使うので計算しておいて、

  var kk1 = Math.exp(kk);//ここら辺はcothがないので地道に作ってます
  var kk2 = Math.exp(-kk);//ここら辺もcothの準備

  //先ほどの式にJijとhi以外の磁場項の計算。cothは地道に計算。。。
  dE += q[y][x]*(q[(m+y-1)%m][x]+q[(y+1)%m][x])*Math.log((kk1+kk2)/(kk1-kk2))/kT;
  if(dE<0 || Math.exp(-dE/kT)>Math.random()){ //メトロポリス法で判定
   q[y][x] = -q[y][x]; //ビット反転
  };
 };
 G = G*0.999;//磁場を弱める
};

こんな感じで実装できました。

早速キメラグラフの実装

キメラグラフは4*4の井形をしています。基本的に上記で見てきたSAとSQAの式を変形することで実装ができます。

キメラグラフSA

キメラグラフのSA実装は並列数がないので、2048量子ビットの大きめでやって見ましょう。Jijを1とか-1に変えると模様が変わるので楽しいです。sa.js

var N = 2048;//量子ビット数は2048
var kT = 5;//初期温度5
var q = [];
var jij = -1;//相互作用項-1にしてみました。
var hcell = 16;//横方向のユニットセル数
var hq = hcell*8;

//量子ビットの初期化
for(var i=0;i<N;i++){
 q[i] = Math.floor(Math.random()-0.5)*2+1;
};

//早速計算。基本的には上記のSAと変わらない。
for(var i=0;i<550;i++){
 for(var k=0;k<N;k++){
  var x = Math.floor(Math.random()*N);//量子ビット選ぶ
  var xcell = Math.floor(x/8)*8;//量子ビットが含まれるセルの位置特定
  //こっからはコツコツキメラグラフ形状にそってエネルギー計算
  if(x%8<4){
   var dE = jij*2*q[x]*(q[xcell+4]+q[xcell+5]+q[xcell+6]+q[xcell+7]);
   if(x<N-hq){
    dE += jij*2*q[x]*(q[x+hq]);
   }
   if(x>hq-1){
    dE += jij*2*q[x]*(q[x-hq]);
   }
  }else{
   var dE = jij*2*q[x]*(q[xcell+0]+q[xcell+1]+q[xcell+2]+q[xcell+3]);
   if(x%hq>7){
    dE += jij*2*q[x]*q[x-8];
   }
   if(x%hq<hq-8){
    dE += jij*2*q[x]*q[x+8];
   }
  }
  //メトロポリス法判定
  if(dE<0 || Math.exp(-dE/kT)>Math.random()){
   q[x] = -q[x];//ビット反転
  };
 }
 kT= kT*0.99;//温度下げる
};

キメラグラフSQA

sqa.js

var N = 512;//量子ビット数は512。並列数が多いので増やしすぎると大変。
var kT = 0.02;//初期温度固定の0.02
var G = 10;//初期のガンマ10
var q = [];
var m = 256;//トロッタ数256
var jij = 1;//Jijは1に設定。-1とかでも楽しい。
var hcell = 8;//横方向のセル数。自動計算でもいいので今回は自明なので設定。
var hq = hcell*8;

//量子ビットを全トロッタ、全量子ビットで初期化
for(var j=0;j<m;j++){
 q[j] = [];
 for(var i=0;i<N;i++){
  q[j][i] = Math.floor(Math.random()-0.5)*2+1;
  };
};

//早速計算。基本的には上記SAとほぼ同じです。
for(var l=0;l<600;l++){
 for(var k=0;k<N*m;k++){
  var x = Math.floor(Math.random()*N);
  var y = Math.floor(Math.random()*m);
  var xcell = Math.floor(x/8)*8;
  if(x%8<4){
   var dE = jij*2*q[y][x]*(q[y][xcell+4]+q[y][xcell+5]+q[y][xcell+6]+q[y][xcell+7]);
   if(x<N-hq){
    dE += jij*2*q[y][x]*(q[y][x+hq]);
   }
   if(x>hq-1){
    dE += jij*2*q[y][x]*(q[y][x-hq]);
   }
  }else{
   var dE = jij*2*q[y][x]*(q[y][xcell+0]+q[y][xcell+1]+q[y][xcell+2]+q[y][xcell+3]);
   if(x%hq>7){
    dE += jij*2*q[y][x]*q[y][x-8];
   }
   if(x%hq<hq-8){
    dE += jij*2*q[y][x]*q[y][x+8];
   }
  }
  dE = dE/m;
  var kk = G/kT/m;
  var kk1 = Math.exp(kk);
  var kk2 = Math.exp(-kk);
  dE += q[y][x]*(q[(m+y-1)%m][x]+q[(y+1)%m][x])*Math.log((kk1+kk2)/(kk1-kk2))/kT;
  if(dE<0 || Math.exp(-dE/kT)>Math.random()){
   q[y][x] = -q[y][x];
  };
 };
 G = G*0.99;
};

これらファイルの活用の仕方

今回紹介したのはコアのアルゴリズムの部分のjsです。これを実装したい場合には、htmlファイルなどに書き込む必要があります。せっかくなのでhtmlファイルの実装例も入れておきたいと思います。sa-canvas.js

<!doctype html>
<html>
<head>
<style>body{font-size:10px;}</style>
</head>
<body>
kT=<span id="TT"></span><br>
<canvas id="main" width="2000" height="2000"></canvas>

<script>
 var N = 2048;
 var kT = 5;
 var q = [];
 var jij = -1;
 var hcell = 16;
 var hq = hcell*8;

 for(var i=0;i<N;i++){
  q[i] = Math.floor(Math.random()-0.5)*2+1;
 };

 var canvas = document.getElementById('main');
 var ctx = canvas.getContext('2d');

 var anneal = setInterval(function(){

  for(var k=0;k<N;k++){
   var x = Math.floor(Math.random()*N);
   var xcell = Math.floor(x/8)*8;
   if(x%8<4){
    var dE = jij*2*q[x]*(q[xcell+4]+q[xcell+5]+q[xcell+6]+q[xcell+7]);
    if(x<N-hq){
     dE += jij*2*q[x]*(q[x+hq]);
    }
    if(x>hq-1){
     dE += jij*2*q[x]*(q[x-hq]);
    }
   }else{
    var dE = jij*2*q[x]*(q[xcell+0]+q[xcell+1]+q[xcell+2]+q[xcell+3]);
    if(x%hq>7){
     dE += jij*2*q[x]*q[x-8];
    }
    if(x%hq<hq-8){
     dE += jij*2*q[x]*q[x+8];
    }
   }
   if(dE<0 || Math.exp(-dE/kT)>Math.random()){
    q[x] = -q[x];
   };
  }

  ctx.clearRect(0,0,2000,900);
  for(var i=0;i<N;i++){
   if(q[i]==1){
    ctx.fillStyle = 'rgb(255,0,0)';
   }else{
    ctx.fillStyle = 'rgba(255,255,255,0)';
   }
   if(i%8<4){
    ctx.fillRect(3+i%8*7+Math.floor(i%hq/8)*40,Math.floor(i/hq)*40,3,30);
   }else{
    ctx.fillRect(Math.floor(i%hq/8)*40,i%8*7-25+Math.floor(i/hq)*40,30,3);
   }
  };
  kT= kT*0.999;

  document.getElementById('TT').innerHTML = kT;
  if(kT<0.02){
   clearInterval(anneal);
  }
 },1);
</script>
</body>
</html>

webGL

流石にjavascriptだけでブラウザでマシンパワーを使って量子シミュレーションをするのは無理があります。モンテカルロ法はなかなか収束するのにマシンパワーも必要なので。並列処理をブラウザで行うためのwebGLを使うとあまりチューニングしてないcudaくらいの速度はでるので、興味ある方はそちらもやって見てください。かなりでかいグラフでも計算ができます。

Recommended


Wikiへ移動