色んなIsing模型（QUBO）を見てみよう！！

はじめに

量子アニーラーやデジタルアニーラーなどのイジングマシン（QUBO solver）が扱える問題が、QUBOと呼ばれる最適化問題です。この記事では量子アニーリングの研究に登場する様々なQUBOについて見ていくことにしましょう。

最初に、QUBOについてもうちょっと説明をします。そもそもQUBOとは何でしょう。

QUBOは**二次形式の制約なし二値変数最適化（Quadratic Unconstrained Binary Optimization）**の頭文字をとったものです。式で表すと、$$f(b):=\sum_{i,j=1}^NQ_{ij}b_ib_j$$という二次形式のコスト関数の最小値を与える$b_1,\dots,b_N$を求める問題です。ここで変数$b_i$は$0$か$1$の値をとるbinary変数です。$Q_{ij}\in\mathbb{R}$が問題を定めるパラメータになっています。

量子アニーリングの研究では、基本的に上のQUBO形式ではなく、等価な変形を施したIsing形式が使われます。Ising形式では、次の形のコスト関数を考えます。
$$H(s):=\sum_{i,j=1}^NJ_{ij}s_is_j$$
といっても全く同じ形をしていますね。違うのは変数$s_i$の取る値です。Ising形式の変数$s_i$は、$+1$か$-1$の値をとる変数でスピン変数と呼ばれます¹。 符号の反転$+\leftrightarrow -$で、変数の値を行き来できる便利な形式です。Ising形式とQUBO形式は、変換$$s_i=1-2b_i,\ J_{ij}=Q_{ij}/4$$で移り変わることができます。

Ising形式のコスト関数$H(s)$は、物理の分野ではエネルギーに対応しています。何のエネルギーかというと、「スピングラス」と呼ばれる不純物の混ざった磁性体のエネルギーです。物理で知りたいことの1つは、スピングラスのエネルギーが最も低い状態（実現される状態）がどういう状態かということです。これは$H(s)$の最適解を求める問題です。そのため$J_{ij}$の値を色々と変えた場合の振る舞いが調べられています。それを紹介していきたいと思います。

とっても簡単

では簡単な順にIsing形式のコスト関数$H(s)$（Ising模型）を見ていきましょう。

一番簡単なのは$J_{ij}$が$i,j$に依らずに一定値$J$をとる場合です。つまり、$$H_{\mathrm{MF}}(s)=J\sum_{i,j=1}^Ns_is_j=J\left(\sum_{i=1}^Ns_i\right)^2$$です。これは**平均場Ising模型（mean-field Ising model）**と呼ばれています。

右辺の形を見れば、この問題の最適解を求めるのは簡単ですね。$J<0$のときは、すべてのスピン変数が同じ値をとっている場合です。$J>0$のときは、$\sum_{i=1}^Ns_i=0$となる場合、つまり$-1$のスピン変数と$+1$のスピン変数が同数存在している場合が最適です。$J=0$のところを境にして、最適解が変化するのはちょっぴり面白いです。相転移です！

実際にあるのは？ −グラフと隣接相互作用−

平均場Ising模型ではすべてのスピン変数が、同じ$J$で結合されていました。しかし、D-waveの量子アニーラーのquantum processing unit（QPU）では、すべてのスピン変数を結合させることはできません。D-waveのQPUでは近くにあるスピン変数のペア$(s_i,s_j)$にだけ、非ゼロの$J_{ij}$を設定することができます。

近くにある、隣接しているといった関係を表すには、ノード（頂点）をエッジ（辺）で結んだグラフが使われます。ノードは、スピン変数$s_i$がある場所を表しています。以下では、ノードの集合を$V$で表すことにします。$V$は自然数でラベルづけができるような有限集合です。一方、エッジはノードとノードを結んだものです。エッジによって結ばれているノード同士は隣接している考えます。エッジは2つのノードのペア$(i,j)$により指定することができます。以下ではエッジの集合を$E$で表すことにします。つまりグラフとは、ノードの集合とエッジの集合のペアです：$G=(V,E)$。例えば、D-waveのQPUで使われているのはキメラグラフと呼ばれるグラフになっています²。

準備は以上です。

D-waveのQPUのように、$J_{ij}$を設定できるのがグラフ$G$のエッジ$E$に制限されている場合に実現できる最も簡単なコスト関数は、$$H_\text{Ising}(s)=J\sum_{(i,j)\in E}s_is_j$$です。エッジが存在するところにだけ、$J$があります。この模型が標準的なIsing模型です。

グラフがすべてのノード間にエッジがある完全グラフの場合には、前に出てきた平均場Ising模型に一致します。グラフ$G$を埋め込める空間の次元に応じて、1次元Ising、2次元Ising、、、と呼ばれたりします。

$H_\text{Ising}(s)$の最適解を求める問題は、$J<0$のときは簡単です。前と同様にすべてのスピン変数が同じ値をとっていればよいです。

しかし、$J>0$のときには一気に難しくなります。とりあえずは、隣接するスピン変数が異なる値をとっていれば、$Js_is_j=-1$となるので良さそうです。ですが、一般にすべてのスピン変数のペアについて、このような値の割り振りができるとは限りません。

例えば、三角形の頂点をノード、辺をエッジとしたグラフを考えると、この方針には問題があることがわかります。最後のスピン変数を$+1,-1$のどっちにすれば良いのか、この方針では決まりせん。このような状況は（幾何学的）フラストレーションがあると表現されます。フラストレーションが原因で問題が難しくなります。

$J>0$のときのIsing模型の最もエネルギーが低い状態を求める問題は、グラフ理論における最大カット問題（MAX-CUT）と等価になります。MAX-CUTはNP-hardであることが証明されている難しい問題です。したがって、Ising模型の最適解を求めるには、一般に指数的に長い時間が必要になります。$J<0$のときと$J>0$のときとで、随分と問題の難しさが変化しましたね。

ランダムネスを入れてみよう

より複雑で一般的な場合として、$J<0$と$J>0$が混ざっている場合を考えてみましょう。コスト関数は、$$H_\text{Rand}(s)=\sum_{(i,j)\in E}J_{ij}s_is_j$$です。今までと違って、$J_{ij}\in\mathbb{R}$はエッジごとに異なる値を取るとします。このような模型はランダムIsing模型またはEdwards–Anderson模型（EA模型）と呼ばれます。特にグラフ$G$が完全グラフの場合、つまり$$H_\text{SK}(s)=\sum_{i,j=1}^NJ_{ij}s_is_j$$の場合には、**Sherrington–Kirkpatrick模型（SK模型）**という名前がつけられています。Kirkpatrickはシミュレーテッド・アニーリング（疑似焼きなまし法）の提案者の1人としても有名です。

これらのコスト関数の最適解を見つける問題は、特殊な場合として$J>0$のIsing模型を含んでいますから、一般にNP-hardの問題です。NP-hardクラスの定義から、多くの問題を$H_\text{Rand}(s)$の最適解を見つける問題として定式化することができます。具体的なIsing形式へのマッピングについては、次の論文が参考になります。

• Andrew Lucas, Ising formulation of many NP problems. Frontiers in Phys. 2, 5 (2014).

一般にはNP-hardであるランダムIsing模型ですが、グラフ$G$が特殊な場合には比較的簡単に最適解を求めることができます。例えばグラフ$G$が平面グラフの場合には、$H_\text{Rand}(s)$の最適解を多項式時間で厳密に求めるアルゴリズムが存在します³。1つ次元を上げた3次元のグラフでは、NP-hardになります。

さて、$H_\text{Rand}(s)$をランダムIsing模型と読んでいましたが、今までの話だと、その"ランダムさ"がよくわかりません。$H_\text{Rand}(s)$がランダムIsing模型と呼ばれるのは、問題を生成するときに、通常$J_{ij}$を確率的にランダムに決めることが由来になっています。$J_{ij}$を決めるのに、よくあるのは次のような分布で$J_{ij}$をサンプリングすることです。

• ±J（bimodal）：$+1$が確率$1/2$、$-1$が確率$1/2$
• Unifrom-Range-k-Disorder：$\lbrace -k, -k+1,\dots,-1,1,\dots,k-1,k\rbrace$の一様分布
• Gaussian：$P(J_{ij})\propto e^{-J_{ij}^2/\sigma^2}$

もっと多く、より複雑に

最後にIsing模型をさらに拡張してみましょう。今までのコスト関数は二次形式（$b_i$と$b_j$の2つの積）でしたが、$b_i$と$b_j$と$b_k$の3つの積が現れるコスト関数$$H_{3\text{-spin}}(s)=\sum_{(i,j),(i,k),(j,k)\in E}J_{ijk}s_is_js_k$$を考えてみましょう。このような模型は3-spin模型と呼ばれます。一般に$p$個のスピン変数の積が現れる場合には、$p$-spin模型と呼びます。

このようにスピン変数の複数個の積も考えると問題は、さらに難しくなります。実際、平面グラフでは、ランダムIsing模型（2-spin模型）の最適解は多項式時間で求めることができましたが、3-spin模型にすると平面グラフでもNP-hard問題になります⁴。

$p$-spin模型のコスト関数をバイナリ変数で表した場合には、Higher Order Binary Optimization（HOBO）と呼ばれます。HOBOのコスト関数はQUBOのコスト関数に書き直す（リダクションする）ことができます。その詳しい方法については、Noriaki Shimada（@nori_autumn）さんの次の記事が参考になります。

• HOBOからQUBOへの変換

あとがき

この記事では、色々なQUBO（Ising模型）とその性質について説明をしました。ちょっと難しすぎる記事になってしまったかもです。

1. $\sigma_i$で書かれることも多いです。 ↩

2. 次世代のD-waveのQPUではキメラグラフよりも複雑にしたペガサスグラフが使われます。 ↩

3. 具体的には、F. Barahona, On the computational complexity of Ising spin glass models. J. Phys. A 15, 3241 (1982).で証明されました。また、Spin Glass Server（https://informatik.uni-koeln.de/spinglass/）を使うと、2次元のランダムIsing模型の最適解を求めることができます。 ↩

4. これは、Creighton K. Thomas and Helmut G. Katzgraber, Optimizing glassy p-spin models. Phys. Rev. E 83, 046709 (2011).で証明されました。 ↩