@yuichirominato 2017.12.03更新

量子アニーリング、イジングモデルとフレームワーク

QUBO イジング シミュレーテッドアニーリング 量子アニーリング

はじめに

最近は従来型のユニバーサルゲートモデルと呼ばれる量子コンピュータの他に、カナダのD-wave社のような量子アニーリングモデル、またNTTの量子ニューラルネットワークのようなイジングマシンもできてました。また、富士通のような従来型の半導体技術を使用してイジングモデルを解くようなマシンも出てきています。これらイジングモデルは同一のフレームワークで動く事実上のデファクトスタンダードのような規格があり、それを理解することでイジングモデルのアプリケーションの開発やハードウェアの開発などを理解することができます。

量子コンピュータ、特に量子アニーリングを学び始める際にいきなり量子アニーリングを学ぶよりも、その元となっているシミュレーテッドアニーリング(以下SA)を学ぶのがいいと思います。またイジングモデルという物理モデルを学ぶ必要もあるので、分けて紹介したいと思います。

SAとは?

「SAアルゴリズムは、解を繰り返し求め直すにあたって、現在の解のランダムな近傍の解を求めるのだが、その際に与えられた関数の値とグローバルなパラメータ T(温度を意味する)が影響する。そして、前述の物理過程との相似によって、T(温度)の値は徐々に小さくなっていく。このため、最初はTが大きいので、解は大胆に変化するが、Tがゼロに近づくにつれて収束していく。最初は簡単に勾配を上がっていけるので、山登り法で問題となるようなローカルな極小に陥ったときの対処を考える必要がない。」
引用元:wikipedia

SAは温度Tを使いグラフの大域的最適解を求めるのに適しているようです。

簡単なグラフでの確認

適当なグラフで温度T=100を書記として,$x$をランダムからスタートして動作をみてみます。
温度Tが大きい時は赤い丸の動きは激しいですが、温度が低くなるにつれグラフの低いところに移動をして動きが少なくなって来ます。
ここでは$x$をランダムで選択し、それに対応する$y$の値が小さくなれば無条件で更新、大きくなる場合でも差分$dy$を求めて$\mathrm{e}^{-dy/T}$の確率で更新するメトロポリス法で更新判定を行い、それを温度Tが十分低くなるまで繰り返しています。

イジングモデル

実際の量子アニーリングマシンにはイジングモデルという物理モデルがハードウェアとして実装されており、シミュレーテッドアニーリング(もしくはシミュレーテッド量子アニーリング)と組み合わせることで私たちのコンピュータ上でも計算をシミュレートすることができます。カナダのD-wave社のマシンはキメライジングという少し変わった並び方をしていますが、最後に紹介します。

統計力学において、イジング模型(英: Ising model、イジングモデルとも言う)とは二つの配位状態をとる格子点から構成され、最隣接する格子点のみの相互作用を考慮する格子模型
引用元:wikipedia

格子状のモデルのようです。

1次元イジング

1次元イジングモデルは量子ビットが一列に並んでいます。量子ビットはスピンを考え、上向き+1(ここでは赤)、下向き-1(同様にグレー)の値を取ります。
7.png

2次元イジング

2次元イジングは考え方は1次元イジングと同様で、隣接する量子ビットが増えます。
8.png

キメライジング(キメラグラフ)

縦横4量子ビットの井型のユニットセルがたくさん並んだ形状を構成するキメライジングです。1量子ビットに最大6の隣接(交差)があります。
9.png

SA+イジング

シミュレーテッドアニーリング(以下SA)とイジングを組み合わせて計算します。計算が速いのでモンテカルロ法+メトロポリス法という乱数を使用する方法を採用します。

イジングモデルのエネルギー関数

SAでイジングモデルの大域最適解を探索するには、$J_{ij}$と呼ばれる量子ビット間の結合の強さと、量子ビットにかかる磁場項$h_i$から構成されるエネルギー関数(コスト関数、ハミルトニアン)にSAを適用します。

$
E=-\sum_{i,j}J_{ij}q_iq_j-\sum_ih_iq_i
$

簡単な例で初期エネルギーを計算

1次元のイジングモデルで5量子ビット$q_0$から$q_4$が並んだ例を用意します。今回は簡単のため量子ビット間の結合の強さ$J_{ij}=1$,また、量子ビットに個々にかかる磁場項を$h_i=0$とします。また、初期状態の量子ビット$q_i$の値は+1(ここでは赤)もしくは-1(同様にグレー)のランダムでスタートします。

1-1-768x102.png

この例では初期エネルギーはエネルギー関数の式を使って以下のように計算できます。$J_{ij}$はすべて+1、$h_i$はすべて0なので式を簡単にします。

$
E_0 = -(J_{01}q_0q_1+J_{12}q_1q_2+J_{23}q_2q_3+J_{34}q_3q_4)-(h_0q_0+h_1q_1+h_2q_2+h_3q_3+h_4q_4) \\
= -(q_0q_1+q_1q_2+q_2q_3+q_3q_4)
$

今回は$q_0=1$,$q_1=1$,$q_2=1$,$q_3=-1$,$q_4=1$なので、エネルギーは上記を計算して、
$E_0 = 0$となります。

SAで量子ビットを反転させていく

次にSAを始めます。今回は初期温度を $T_0=5$と設定してみます。
まず量子ビットをランダムで一つ選びます。今回は$q_0$を選びました。下記のとおり青枠をつけてみました。
次に、この量子ビットを反転させた場合のエネルギーを計算します。

2.png

現在$q_0=1$なので反転させた場合には$q_0′ = -1$

3.png

その場合のエネルギー$E’$は反転させた場合の$q’$の値を考えて、
$E’ = 2$

次に反転させる前、反転させた後のエネルギー差を考慮して、本当に反転させるかどうかを判断します。今回、エネルギー差$dE = E’-E_0 = 2$なので、メトロポリス法で反転確率 $e^{-dE/T}$ を使用して、反転させるかどうかの確率を計算します。計算すると$e^{-dE/T} = 0.6703200460356393$なので、この確率で量子ビットを反転させます。プログラム的には適当に乱数を0-1の間で発生させ、この確率と比べて確率が大きければ反転、小さければ反転させないでそのままになります。量子ビットを反転する、しないをきめて実行した後、また次の量子ビットをランダムで選択し、同じことを繰り返します。

温度を下げてまた繰り返し

現在初期温度を5に設定しました。ある程度上記の操作を繰り返した後に、今度は温度を下げます。そしてまた同じように繰り返していきます。そうすると温度が高い時からだんだん低くなるにつれて動きが少なくなり、安定した状態に近づき解を求めることができます。

デモ

$T=5$から$T=0.02$までだんだんと温度を下げながら計算を行ってみます。

4.gif

キメラグラフ(キメライジング)とは?

キメラグラフは通常のイジングモデルの形状ではなく、4×4の格子状の量子ビット回路がたくさん並んでいる形状をしています。量子ビットは細長い形をしており、最大で隣接するもしくは交差する量子ビットの数が「6」となっています。この4×4の格子状の量子ビット回路を1ユニットセルと呼び、そのユニットセルが縦横にたくさん並ぶことで大きな回路を構成しています。2048量子ビットのマシンでは、この4×4のユニットセルが縦横16×16並んでおり、8×16×16=の2048量子ビットの回路を構成しています。

長い結合とキメラグラフ

キメラグラフは量子ビット同士の最大接続数は6ですが、キメラグラフで完全グラフを作るなど、接続数を作るためにはキメラグラフの量子ビット同士を接続するテクニックを使用します。実際の量子コンピュータでは結合する部分で多少のエラーが蓄積しますが、理論的には完全グラフを構成することができます。下記の赤い部分の量子ビット間にカプラーという結合がありますので、その結合の値を量子ビットが同じ値をとるように設定するとキメラグラフで長い結合を作ることができます。

1.png

完全グラフとキメラグラフ

キメラグラフで完全グラフを作成するには、下記のように回路をくの字に折りながら長距離の結合を三角形状に作り、重ねていきます。

2.png

これをさらにキメラで表現すると三角形のエリアを使って完全グラフを計算します。

3.png

QA離散化

イジングモデルにSAを適用してきましたが、こちらを量子アニーリング(シミュレーテッド量子アニーリング、以下SQA)に展開します。量子ビットの重ね合わせを数学的な処理を用いて近似し、並列処理に落とし込みます。

量子項を考慮したイジングモデルのエネルギー関数

量子コンピュータや量子アニーリングマシンでは量子ビットは重ね合わせの状態と呼ばれる状態に落ちています。既存の計算機で計算する場合にはこの重ね合わせ状態を並列処理に落とし込みます。業界標準の方法となっている経路積分と鈴木トロッター分解という操作を行うと、下記のエネルギー関数の式を得ることができます。

$E=-\sum_{i,j}\sum^m_{k=1}\frac{J_{ij}}{m}q_{i,k}q_{j,k}-\sum_i\sum^m_{k=1}\frac{h_i}{m}q_{i,k}-\frac{T}{2}\sum_{i,j}\sum^m_{k=1}lncoth(\frac{\Gamma}{mT})q_{i,k}q_{i,k+1}$

$m$はトロッタと呼ばれる並列数。$k$は何番目のトロッタか。$\Gamma$は磁場の強さです。SAのエネルギー関数から難しくなった感じがありますが、

$E=-\sum_{i,j}J_{ij}
q_i
q_j
-\sum_ih_iq_i
$

最初の2項が$m$で割られて、最後に何か長い式がついてます。
量子アニーリングでは最後の式に入っている$\Gamma$を最初に強く設定して、それをSAの温度のように段々弱めて最後は0にすることで計算を行います。

虚時間軸とトロッタ、並列離散近似の可視化

式だけだと理解しづらいので、絵にして書いて見ます。1次元イジングを考えた場合、SAの場合一列だったイジングが縦にたくさんレプリカが並びます。このレプリカを作る方向を虚時間軸とよび、レプリカの通し番号をトロッタと呼んでいます。問題に設定されたパラメータは$J_{ij}$や$h_i$はすべてのトロッタで同じです。1次元の問題は量子系では2次元に展開されます。下記ではトロッタ数m=4の場合を想定しています。

3.png

量子モンテカルロ法と計算

ここまでいくと、あとはSAの時と同じように量子ビットを一つ選んで反転させるかどうかの試行を繰り返し、温度$T$の代わりに磁場の強さ$\Gamma$を初期値から段々弱めていきます。エネルギー関数からエネルギーを計算する際に量子項の計算だけ面倒ですが、そこだけ頑張れば大丈夫です。

デモ

$\Gamma$=20から$\Gamma$=0.01まで段階的に下げていった場合の量子アニーリングの様子を可視化して見ます。量子アニーリングでは基本的に縦方向の量子ビットが揃えば計算が成功です。

2.gif

そして実際の問題を解くために

ここまでで一通り学ぶと量子アニーリングを使うことができます。実際の問題ではNのサイズを大きくして、問題に合わせて$J_{ij}$や$h_i$を自分で作って見て設定します。また、並列数を大きくして量子シミュレートをすることでより量子アニーリングの効果を実感することができると思います。

量子コンピュータは繰り返し計算につよい

量子コンピュータは超電導で冷却をしたり、値の設定やキャリブレーションなどがとても難しい。アナログコンピュータのイジングマシンでは特にそうである。アナログコンピュータなので動きは予測不能で、解がもとまる時もあれば、もとまらない時もある。なのでデジタル計算機のように一度計算すれば解がもとまるものではなくて、一度に何十回、何百回と計算を回して、その結果からいい値を選ぶようになる。解の精度が低いが、繰り返し計算に強いので、解を取り出して、自分の設定した問題の式に当てはめていい答えが出たものを選ぶようになる。

検算や解のチェック方法

量子コンピュータのイジングマシンではまずイジングモデルというものにアプリケーションをマッピングする。その設定したイジングのパラメータと出てきた答えを組み合わせることで、自分の計算したもののハミルトニアン・エネルギーを計算できる。イジングマシンは最小値問題を解くための機械なので、自分の設定したパラメータと出てきた量子ビットの値からエネルギーを古典計算機で計算し、エネルギーの低いものを答えとして選ぶ。実際、答えは本当に合っているのかどうかは実は古典計算機では計算できないので、解精度がある程度保証されているアルゴリズムSDPとかと比較したり、SAと比較したりして確からしさをチェックする。

量子コンピュータは消費電力が低い

D-wave社の動作する消費電力は25kWと言われているが、このほとんどが冷凍機の動作に使用されているという。実際アナログコンピュータなのでチップ自体にかかる消費電力はほぼ0であるので、繰り返し計算をしても冷凍機の消費電力以上のものは基本的にかからないので、将来的にも消費電力の低減が期待できる。

現状

上記イジングモデルに対して、開発を行いたいと思った場合にどうすればいいのでしょうか。現状はイジングモデルを開発するための開発者キットやフレームワークの公開は見つけるのが困難です。これまで量子アニーリング界隈で先行していたカナダ・バンクーバーのベンチャー企業である1Qbit社( https://1qbit.com )が公開していたQDKというイジングフレームワークと開発者キットがあったのですが、今年それが更新されず非公開になってしまいました。さらに、現在全世界でイジングモデルのマシンを開発しているのはカナダと日本のみなので、この両国で開発キットが見つからない状況では開発は少し厳しい状況ですので、開発キットの登場を待つ必要があります。

QUBOフレームワーク

イジングモデルはQUBOと呼ばれるフレームワーク(ミドルウェア)に落とされます。

これは、

1、問題をグラフ構造で2量子ビット間の相互作用を非対角項に、磁場項と呼ばれるバイアス項を対角項に書いた行列の形で提供する。
2、{0,1}のビットで考えられた問題を{-1,1}に変換する。

という二つの操作をします。私たちがイジングモデルのアプリケーションを作成する場合には現状の古典計算機で考えるように{0,1}で考えてもらいますが、量子コンピュータの内部では{-1,1}で計算を行うために変換が必要です。

QUBOから実機への実装

イジングモデルの量子コンピュータやシミュレータはそれぞれ異なった回路形状を持っています。QUBOmatrixへ一度落とされた問題は、これら実機の回路形状に最適化された形で変換されます。例えばD-wave社はキメラグラフと呼ばれる回路形状を持っています。これは44の井形のユニットセルと呼ばれる単位を2次元に展開しており、1量子ビットで最大6の結合数を持つ特殊な形をした回路です。通常QUBOmatrixはすべての量子ビットからすべての量子ビットへの完全結合の形で書かれているので、これらQUBOmatrixの結合をD-waveのキメラグラフに落とし込むためにいくつかのステップを通じて変換されます。D-wave社の2017年現在での量子ビット数は2000量子ビットなので、この量子ビットを超える問題であったり、量子ビット数は超えないけれど問題の形によって変換をかけると超えることがあるので、その場合にはqbsolvと呼ばれるフレームワークでさらに問題分割を行なって複数回の計算を通じて解を求めるなどの工夫をしています。NTTの量子ニューラルネットワークは完全結合ですので、QUBOmatrixはそのまま使用できますが、現状のマシンはJijのみのサポートなので、QUBOmatrixの対角成分は計算ができません。~~富士通のデジタルアニーラは古典計算機&FPGAでの実装ですので、D-waveマシンやNTTマシンに比べて自由度が高いと推測されますが、市場に出回っている情報があまりないため評価は現状ではあまりできません。~~(2017.12.3追記)富士通のデジタルアニーラの論文を見つけました。基本的に1024bitの実装で完全結合、磁場項hもあるので、10241024のQUBOmatrixなら分解なしで実装ができます。

アプリケーションからQUBOへの実装

では、QUBOに落とし込めればマシンを問わず問題を解くことができることが何と無く分かりました。QUBOに落とし込むまではどのようにするのでしょうか。QUBOに落とし込むには自分で直接考えてダイレクトに行列を作る方法もありますが、一般的には「ハミルトニアン」とか「コスト関数」とか「エネルギー関数」と呼ばれる関数をグラフ問題として作ります。このハミルトニアンは相互作用係数(結合荷重)のJijと磁場項(バイアス)h、そしてそれぞれの量子ビットqiをつかって式を作ります。問題によって色々制約条件というものがあり、制約条件を満たすようにJij,h,qi,qjなどを駆使して式を立てます。複雑な式でもいいですが、式を立てた後は展開をして解いていきます。ハミルトニアンを作れば必ずしもQUBOmatrixに落とし込めるわけではないので、例えばq1q2q3のように3量子ビットの掛け算などが出てくることがあります。これらは多体問題とよばれ、QUBOでは直接は計算できないので、2体の問題に落とし込む必要があります。これらは多少数学的な知識が必要になります。数学変換などを活用してハミルトニアンを2体問題に落とし込めれば、QUBOmatrixにいれることができ、各種量子コンピュータやシミュレータで解くことができます。

pythonとSymPyを使ってQUBOの実装を見て見たいと思います。

量子アニーリングは最小値を求めるアルゴリズムなので、通常二次計画問題に落とし込まれた二次式を式変換したりします。また、ビットの表現は問題を作成する際には{0,1}ビットですが、物理モデルに落とす部分は{-1,1}なので、その変換が毎回発生します。

$(q_0+2q_1-q_3-q_4)^2 = q_0^2+4q_1^2+q_3^2+q_4^2+4q_0q_1-2q_0q_3-2q_0q_4-4q_1q_3-4q_1q_4+2q_3q_4$

SymPyの使い方(展開)

簡単に確認して見ました。式展開をしたい場合と因数分解をしたい場合は下記の通りでした。expandを使うと展開できるようです。


from sympy import *

x = Symbol("x")
f = (x+1)**2

print(expand(f))

x**2 + 2*x + 1

すごい、、、今までの苦労はなんだったんだ、、、

SymPyの使い方(因数分解)


from sympy import *

x = Symbol("x")
f = (x+1)**2

print(factor(expand(f)))

(x + 1)**2

便利じゃないですか、、、

早速量子アニーリング関連の二次計画問題の展開をば。。。

このブログでやった量子アニーリングでの1+1の展開をSymPyでやって見たいと思います。
https://qiita.com/mdrft/items/0277d6f695031dd1f186

$(q_0+2q_1−q_3−q_4)^2$

この式を展開したいだけなんです。。。


from sympy import *

q0,q1,q3,q4 = symbols("q0 q1 q3 q4")
f = (q0+2*q1-q3-q4)**2

print(expand(f))

q0**2 + 4*q0*q1 - 2*q0*q3 - 2*q0*q4 + 4*q1**2 - 4*q1*q3 - 4*q1*q4 + q3**2 + 2*q3*q4 + q4**2

すべてが一瞬です、、、

$q_0^2+4q_1^2+q_3^2+q_4^2+4q_0q_1−2q_0q_3−2q_0q_4−4q_1q_3−4q_1q_4+2q_3q_4$

では、QUBO変換は

QUBO変換はビットを{0,1}から{-1,1}に変換します。その際には$q=\frac{q’+1}{2}$を実行します。

変数の代入は下記のコマンドsubsで、、、

f1 = q0 + 4*q0*q1 - 2*q0*q3 - 2*q0*q4 + 4*q1 - 4*q1*q3 - 4*q1*q4 + q3 + 2*q3*q4 + q4
f2 = f1.subs([(q0,(q00+1)/2),(q1,(q11+1)/2),(q3,(q33+1)/2),(q4,(q44+1)/2)])
print(simplify(f2))

実行したところ、、、下記が得られました。

q00*q11 - q00*q33/2 - q00*q44/2 + q00/2 - q11*q33 - q11*q44 + q11 + q33*q44/2 - q33/2 - q44/2 + 2

バッチリです。これの係数を確認することで簡単に量子アニーリングに取りかかれます。

繰り返し計算、正答率、サンプリング

各種イジングマシンは元となっている考え方がメタヒューリスティックスと呼ばれる確率・統計力学的な手法となっており、組み合わせ最適化問題といえども現状のマシンでは正解がバシッと一度で決まることは稀です。ですが、この類のマシンは一度問題を設定すると繰り返しの計算に強いので、一度で100回や1000回の答えを戻すことができます。それらの解答群の分布を使用したり、その中の最適な解答を選んだりと柔軟に対応して解答を処理します。

まとめ

現在量子コンピュータは黎明期とありまだ構造がシンプルです。また毎年文言やシステムの仕様が少しずつ決まっていきます。今回紹介したアニーリング型、イジング型を量子コンピュータとして考慮しないこともありますが、解ける問題が比較的スパコンが苦手なグラフ構造の組み合わせ最適化問題なので、フレームワークやアプリケーションが現在どんどん考えられている状況です。ミドルウェアやフレームワークの構造も比較的単純で、ユーザー側のインターフェイスから順番に、

・SDK、各種インターフェイス
・ハミルトニアン
・QUBO
・各種変換ミドルウェア(ハードウェア依存)
・量子コンピュータ、イジングシミュレータ

のような状況です。アプリケーションを開発したい場合にはまずハミルトニアンからQUBOへの落とし込みをかんがえ、SDKやAPIなどでインターフェイスを起こすことができますので、ハミルトニアンの式を立てる部分とQUBOを理解すれば進めることができます。ミドルウェアを学ぶのは結構大変で、QUBO自体を調べることになり、それぞれのハードウェアの特性を理解した上で、実際のハードへの実装を理解する必要があるので、時間がかかります。

各種マシンはそれぞれの得意分野や解いているアルゴリズムが違うので、問題によっては疎密やサイズ、問題の傾向によって正答率が変わってきます。ハードウェアの特性に合わせた組み合わせ最適化問題までをチューニングできるようになればかなりアプリケーション開発も有利に進めることができます。

あわせて読みたい