HOBOからQUBOへの変換

はじめに

量子アニーリングなどのアニーリングマシンを実行するには、コスト関数をQUBO(Quadratic Unconstrained Binary Optimization)という形式で定式化する必要があります。
日本語に訳すと「制約なし2値変数2次変数最適化(問題)」になります。
https://en.wikipedia.org/wiki/Quadratic_unconstrained_binary_optimization

物理の言葉で言うと2体の相互作用までを許すイジング模型として定式化する必要があります。

$$
\begin{align}
\mathcal{H} = \sum_{i,j}J_{ij}\sigma_i \sigma_j + \sum_i h_i \sigma_i
\end{align}
$$

イジング模型のもっと詳しい解説はかしゃるふぁさん (@Kashalpha)の記事を参照して下さい。

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

この記事でははスピン変数$\sigma_i \in \{-1, 1 \} $ ではなくバイナリ変数$x_i \in \{0,1\}$を使用します。
この2種類のノーテーションの間には

$$
\begin{align}
x_i = \frac{\sigma_i + 1}{2}
\end{align}
$$

の関係があります。これを利用すればスピン変数とバイナリ変数の間を行き来できます。
ただし定数倍の違いや、高次の項では新たな項が現れたりするのでそこらへんは注意が必要です。

この様な2次項までの関数をQUBOと呼びます。
3次以上の項が現れる関数のことをHOBO(Higher Order Binary Optimization)と呼びます(? 業界用語で一般的な名称かどうか不明です。)

例えば、3次以上の項としては
$$
\begin{align}
a_{ijk}x_i x_j x_k, ; b_{ijkl}x_i x_j x_k x_l, ; c_{ijklm}x_i x_j x_k x_lx_m, \cdots
\end{align}
$$

などが挙げられます。
一般的にこの様な3次以上の項を使えた方が定式化をする時の表現の自由度は高くなります。
グラフ色彩問題を例に説明します。

グラフ色彩問題を用いたHOBOの具体例

グラフ色彩問題については下記を参照してください。めっちゃ分かりやすいです。
https://qard.is.tohoku.ac.jp/T-Wave/?p=418

ハミルトニアンは

$$
\begin{align}
\mathcal{H} = A \sum_v \biggl(1 - \ \sum_i x_{v,i} \biggr)^2 + B \sum_{(u,v)\in E}\sum_i x_{u,i} x_{v,i}
\end{align}
$$

です。$u,v$がノード、$i$が色を指定するインデックスです。
$i$についての和は$K_{max}$まで走り、これが使用する色の数になります。

ここでは色の数は固定して、その色の数の範囲でグラフ上の隣り合うノードが同じ色で塗り分けられない様に塗れるかどうかを問う問題になっています。
$K_{max}=3$だと、3色で塗り分けられるかどうかを考えて、答えはYesかNoのどちらかが返ってきます。
(エネルギーが0だとYes。1以上だとNoだと思ってください。)
この様な問題をNP-Completeと呼び、答えはYesかNoになります。
$K_{max}=5$で答えがYesだとしても、「5色で塗り分けられる」ことしか分からなくて、塗り分けられる最小の色の数が分かるわけではありません。
($K_{max}$をどんどん小さくしていけば求まりますが。)

一方で塗り分けられる最小の色の数も同時に求める問題を考えます。
こちらの方が難しい問題で、NP-Hardと呼ばれています。
答えはYesかNoの二択ではなく、塗り分けられる最小の色の数が返ってきます。

私の頭で考えられるハミルトニアンは、今まで扱ってきた変数$x_{v,i}$を$x_{v}x_i$とノードと色に分離して元のハミルトニアンに代入します。すると

$$
\begin{align}
\mathcal{H} = A \sum_v \biggl(1 - x_v \sum_ix_i \biggr)^2 + B \sum_{(u,v)\in E}\sum_ix^2_i x_u x_v
\end{align}
$$

になります。ここに色を最小にするための項を追加すると、

$$
\begin{align}
\mathcal{H} = A \sum_v \biggl(1 - x_v \sum_ix_i \biggr)^2 + B \sum_{(u,v)\in E}\sum_ix_i x_u x_v +C\sum_ix_i
\end{align}
$$

になります。第2項では$x^2_i=x_i$を利用しています。
第1項は2次の2乗項なので最大4次の項が現れます。
さらに第2項に注目すると、$x$の3次項が現れるのでQUBOになっていません。
実際のアニーリングマシンを実行するにはこれを2次以下にリダクションする必要があります。

今回は下の論文から2種類の方法について解説します。
http://www.f.waseda.jp/hfs/TGBMMTTFOC.pdf
上の日本語バージョン
http://www.f.waseda.jp/hfs/miru2009.pdf

リダクション方法

HOBOからQUBOへのリダクション方法は2種類あるそうです。
今回は上記の論文の3章で開設されている、「最小値選択の変換」、「置換による変換」を説明します。
今回は3次から2次へのリダクション方法についてだけ列挙します。
高次の項については元論文を参照して下さい。

対象なる変数は$xyz$とします。

最小値選択での変換

この変換では、

$$
\begin{align}
xyz = \max_{w \in \mathbb{B}} \bigl[w(x+y+z-2)\bigr]
\end{align}
$$

となります。ここで新しいバイナリ変数$w$が導入されています。
アニーリングマシンに乗せる時に必要なビット数が1だけ増加します。
実際にこの式が正しいのかは(x,y,z)の取りうる8通りの組み合わせを調べると、

x	y	z	x + y + z - 2	xyz
0	0	0	-2	0
1	0	0	-1	0
0	1	0	-1	0
0	0	1	-1	0
1	1	0	0	0
1	0	1	0	0
0	1	1	0	0
1	1	1	1	1

$w$の最大化を考えると、$x+y+z-2$が負になった時は$w=0$、生の時は$w=1$をとる。
しかし、アニーリングマシンで最適化しようとすると$w$は逆の振る舞いをする。
(エネルギーを最小にしようとするから)

\begin{align}
\max_{x}\big[-f(x)\big] = -\min_{x} \big[f(x)\big] \\
-\max_{x}\big[f(x)\big] = \min_{x} \big[-f(x)\big]
\end{align}

を利用すると、

$$
\begin{align}
-xyz =-\max_{w \in \mathbb{B}} \bigl[w(x+y+z-2)\bigr] = \min_{w \in \mathbb{B}} \bigl[-w(x+y+z-2)\bigr]
\end{align}
$$

$\max$から$\min$へと変換できる。
論文上ではもっと一般化して、$a<0$のとき、

$$
\begin{align}
axyz = \min_{w \in \mathbb{B}} \bigl[aw(x+y+z-2)\bigr]
\end{align}
$$

としています。
アニーリングマシンを使用するときは、被最適化関数を最大化するのか最小化するのかによって、上記の$\min, \max$のどちらを選択すべきか注意が必要だと思います。

$a>0$の場合は、

$$
\begin{align}
xyz = \max_{w \in \mathbb{B}} \bigl[w(x+y+z-2)\bigr]
\end{align}
$$

に対して、$x$を$1-x$、$y$を$1-y$、$z$を$1-z$に置き換えて変形すると、

\begin{align}
&(1-x)(1-y)(1-z) = \max_{w \in \mathbb{B}} \bigl[w(1-x+1-y+1-z-2)\bigr] \\
&\Leftrightarrow -xyz =  \max_{w \in \mathbb{B}} \bigl[-w(x+y+z -1)\bigr]-\biggl( (xy+yz+zx) - (x+y+z) + 1\biggr) \\
&\Leftrightarrow xyz =  \min_{w \in \mathbb{B}} \bigl[w(x+y+z -1)\bigr] + (xy+yz+zx) - (x+y+z) + 1
\end{align}

これより、$a>0$では

$$
\begin{align}
axyz = \min_{w \in \mathbb{B}} \biggl[a \biggl\{ w(x+y+x-1) + (xy+yz+zx) - (x+y+z) + 1 \biggr\} \biggr]
\end{align}
$$

4次の項$xyzt$を変換しようとすると、$-1$が偶数回かかってしまい、$max$から$min$への変換がうまく出来ない。5次以上の場合も奇数個の積ではうまくいくが、4次以上の偶数個の積が生じるので$max$と$min$が乱立してしまいうまく変換できない。
($max$か$min$のどちらかで統一できない。)
なのでこの変換は3次までに制限されるようです。

グラフ色彩問題を最小値選択で変換する場合

$x_ix_ux_v$に対して上記の変換を施す。$a=B>0$なので、

$$
\begin{align}
Bx_ix_ux_v = \min_{w_j \in \mathbb{B}} \biggl[B \biggl\{ w(x_i+x_u+x_v-1) + (x_ix_u+x_ux_v+x_vx_i) - (x_i+x_u+x_v) + 1 \biggr\} \biggr]
\end{align}
$$

と変換して、新しい変数$w_j$を導入する。
引数$j$は$x_ix_ux_v$の引数によって新たに決まる値である。
これを元のハミルトニアンにプラスする。
新たに必要なビット数は$\sum_{(u,v)\in E}\sum_i1$の合計となる。

置換による変換

これは$xy$を1つの変数$z$で置換する方法です。
3次のバイナリ変数$xyw$を$zw$にする変換について考えます。
$z=xy$が成立するために、次の制約項をプラスします。

$$
\begin{align}
D(x, y, z) = xy -2z(x+y) + 3z
\end{align}
$$

これは$z=xy$が成立している時のみ$0$になります。
表で表すと、

x	y	z	$D(x, y, z)$
0	0	0	0
1	0	0	0
0	1	0	0
1	1	1	0

$z=xy$を満たさない場合、

x	y	z	$D(x, y, z)$
0	0	1	3
1	0	1	1
0	1	1	1
1	1	0	1

と$D(x, y, z)$が1以上になりエネルギーが増加します。

問題点

$D(x, y, z)=0$になるために非常に大きな生の数$M$を$D(x, y, z)$全体にかける必要があります。

$$
\begin{align}
MD(x, y, z) = Mxy -2Mz(x+y) + 3Mz
\end{align}
$$

日本語の解説記事から引用すると、

第一項Mxyは非常に大きな正の係数を持つ2次の項だが，多項式で表した2次擬ブール関数は，2次の項の係数が全て負であるとき，またそのときに限り劣モジュラであることが知られている [7].そのため，この項 Mxyのために変換後の関数は常に劣モジュラでなくなってしまう.また，同じ正の係数でもMは非常に大きいためか，この方法によって変換されたエネルギーはQPBOを使ってもほとんど全てのサイトにはラベルが与えられず，十分最小化できないことが実験的に確かめられる.

この辺りの意味が分かる方がいれば解説お願いします。
$M$が非常に大きい場合は$x,y,z$のビットフリップが起こりにくくなり更新があまりされず局所解に留まってしまう的なニュアンスでしょうかね。

グラフ色彩問題を置換換する場合

$x_ix_ux_v$の中の$x_ux_v$を新たに$z_j = x_ux_v$と置換すると、

\begin{align}
Bx_ix_ux_v &= Bx_iz_j + M D(x_u, x_v, z_j)\\
D(x_u, x_v, z_j) &= x_u x_v - 2z_j(x_u + x_v) + 3z_j
\end{align}

となる。
$z_j$の引数$j$は$x_ux_v$の引数によって新たに決まる値である。
新たに必要なビット数は$x_ux_v$の項の総数$\sum_{(u,v)\in E}1$である。

最後に

HOBOからQUBOへ変換する方法について説明しました。
変換時にはどうしても補助的なビットが必要になりビット数が増えてしまいます。
「最小値選択での変換」では対象となる項の総数分だけビット数が増加します。
一方で「置換による変換」では上手に置換するビットのペアを決めれば必要な補助ビット数を削減することができます。
しかし、劣モジュラ性のためうまく最適化できない場合があるようです。

さらに高階の変数項をリダクションする方法についてもこの論文では解説されています。
余力があればまた自分なりにまとめて記事を書きます。