量子アニーリング利用のための定式化方法紹介

はじめに

量子コンピュータのうち、量子アニーリング方式は、ゲート方式とは異なり、整数の組み合わせ最適化問題に特化したマシンであり、商用マシンとしては、D-Waveが有名です。

この記事では、D-waveのような、量子アニーリング方式のマシン（以下、量子アニーリングマシン）を用いた、組み合わせ最適化問題の求解の流れ、特に、問題をどのようにイジングモデルやQUBOとして定式化するかについて、重点的に説明したいと思います。

量子アニーリングマシンを用いた最適化問題の求解の流れ

ある問題に対して、量子アニーリングマシンを実行し、解を得るまでの流れは以下の通りです。

1. 問題の定義
2. イジングモデル/QUBOとして定式化
3. マシンに入力・実行
4. 解の取得

手順の2番にあるように、問題を定義した後に、イジングモデル/QUBOとして定式化する必要があり。この記事では、この方法について、重点的に説明します。

それでは、最初に定式化の説明に入る前に、イジングモデルやQUBOについて、説明します。

イジングモデルとは

イジングモデルは物理学の一分野「統計力学」における最もシンプルなモデル¹であり、組み合わせ最適化問題を量子アニーリングを利用して、求解するためには、イジングモデルのハミルトニアン（エネルギー関数）として、問題を表現する必要があります。

具体的に、エネルギー関数$H(σ)$は以下に示す通りです。

\displaylines{
H(σ)=-\sum_{i<j}J_{ij}σ_iσ_j-\sum_{i=0}^{N-1} h_iσ_i \\
σ_0,...,σ_{N-1} =±1
}

$σ$はイジング変数と呼ばれ、取りうる値が-1、または+1の2値変数となっています。また、$J_{ij}$は相互作用、$h_i$は磁場と呼ばれる定数であり、最適化問題の関数の係数に相当します。

量子アニーリングマシンの実行によって、イジングモデルのエネルギーが最も低い状態、つまり、$H(σ)$を最小にする$N$個の変数$σ_0,...,σ_{N-1}$の組み合わせを求めることができます。

ところで、組み合わせ最適化問題は目的関数と制約条件を別の式として表すことが多いですが、イジングモデルのエネルギー関数として表現するためには、緩和問題として、制約条件を目的関数の罰則項として実装し、目的関数と制約条件を1つの2次多項式として、表す必要があります。

QUBOとは

以上で述べた通り、量子アニーリングで問題を求解するためには、イジングモデルのエネルギー関数、つまり、変数の値が$(-1,1)$である1つの2次多項式として定式化する必要があります。

しかし、最適化問題では$(0,1)$のバイナリ変数を用いて、定式化することが良くあるため、量子アニーリングマシンによる求解の際にも、以下のように、バイナリ変数$q_i(=0,1)$をイジング変数$σ_i(=±1)$に変換することで、
$$q_i=\frac{σ_i+1}{2}$$
問題をバイナリ変数の2次多項式として、定式化した場合も、量子アニーリングマシンでの求解が可能となっています。

この、バイナリ変数による、2次多項式問題をQUBOといいます。

OUBOはQuadratic Unconstrained Binary Optimizationの略（和訳：2次制約無し2値最適化問題）であり、イジングモデルと同じく、2値変数の2次多項式問題を示しますが、先ほど述べた通り、バイナリ変数を利用して表現することが特徴です。なお、Unconstrained（制約無し）は、こちらもイジングモデルで説明した通り、制約条件を目的関数の罰則項として実装し、目的関数と制約条件をそれぞれ別々の式として表現しないことを示しています。

具体的に、QUBOでは以下のように定式化を行います。

\displaylines{
H(q)=-\sum_{i,j}Q_{ij}q_iq_j\\
q_0,...,q_{N-1} =0,1
}

$q$はバイナリ変数であり、イジングモデルと同様に、マシンの実行によって、$H(q)$を最小にする$N$個の変数$q_0,...,q_{N-1}$の組み合わせをを求めます。また、$Q_{ij}$は係数を示します。

量子アニーリングを利用して、最適化問題を求解する際は、これまで紹介したイジングモデル、QUBOを問題に合わせて選択し定式化を行うといった流れとなります。なお、この記事ではQUBOとしての定式化を中心に紹介します。

どうやって、問題をイジングモデル/QUBOとして、定式化するか。

実は、「Ising formulations of many NP problems²」という論文で、代表的な組み合わせ最適化問題のイジングモデル、またはQUBOとしての定式化方法がまとめられています。（英語が苦手な方は、Fixstars社のホームページ³がおすすめです。）

そのため、具体的な問題の定式化方法はそちらを参照していただくとして、この記事では、定式化時に少し工夫が必要な、以下の条件の問題について、説明したいと思います。

• 制約条件が存在する時
  • 制約条件が等号の時
  • 制約条件が不等号の時
• 問題の変数が$(0,1)$、$(-1, +1)$のどちらでもない時
• 制約条件や目的関数に3次以上の式が含まれる時

制約条件が存在する時

前提

例えば、以下の制約条件の下で、目的関数を最小化する、$x_0,x_1,x_2,x_3$の組み合わせを求める問題を考えます。

目的関数：$4q_0q_1-3q_1q_3+2q_2-q_3$
制約条件①：$2q_0 + q_1+q_2=2$
制約条件②：$q_0 + q_1 + q_2+ q_3≦3$

ここで、$x_0,x_1,x_2,x_3$はバイナリ変数（取り得る値が0または1）であるとします。

量子アニーリングマシンを用いて、求解を行うために、制約条件を目的関数の項に含めた1つの2値変数の2次多項式として定式化する必要があります。この問題は、バイナリ変数で定義されているため、QUBOとして定式化します。

定式化例

結論から言うと、上で定義した問題は以下のように変形して、定式化されます。

\displaylines{
H(q)=4q_0q_1-3q_1q_3+2q_2-q_3 \\
+A(2q_0 + q_1+q_2-2)^2 \\
+A((q_0 + q_1 + q_2+ q_3 - r_1-2r_2-3r_3)^2+ (1-(r_0+r_1+r_2+r_3))^2)
}

それぞれの1行目の項が目的関数、2行目が制約条件①、3行目が制約条件②を示しています。
$A$は制約項の重み（制約の強さ）を示す定数です。($r$については後で説明します。)

前提で示した式と、$H(q)$の各項を見比べると、目的関数は最小化したい式そのものであるため、QUBOでも元の式と同様に定式化されています。

一方、制約条件は2つとも元の式から、異なる形式となっていますが、これらは共通して、**制約条件を満たす際は値が0、逆に破った際にはペナルティが発生する（値が大きくなる）**ように定式化されています。

それでは、詳細について、制約条件が等号の時（制約条件①）と、不等号の時（制約条件②）にわけて説明します。

制約条件が等号の時

等号で表される制約を、QUBOとして、定式化する際は、移項して「$=0$」の形を作り、両辺を2乗することが多いです。

具体的に、制約条件①は以下のように、等号で表されていますが、
$$2q_0 + q_1+q_2=2$$
これを、右辺の定数2を、左辺に移項して、両辺を2乗すると、以下のようになります。
$$(2q_0 + q_1+q_2-2)^2=0$$
この式の左辺は、2次多項式であり、かつ、$2q_0 + q_1+q_2=2$を満たすとき、値が$0$となり、満たさないときは正の値をとる（=ペナルティを発生させる）ため、制約条件①を表現するQUBOの項として、採用されています。

制約条件が不等号の時

制約条件が、左辺が2値変数の多項式、右辺が定数である不等式で表現される場合、「制約条件の左辺の和が$k$がどうか」を示す、2値変数(今回は$r_k$)を追加で導入し、QUBO（もしくはイジングモデル）として定式化を行います。

ここで、$r_k$の意味は以下の通りです。

r_k= \left\{
\begin{array}{ll}
1 & (左辺の値がkである。) \\
0 & (左辺の値がkではない。)
\end{array}
\right.

具体的に、例として挙げた問題の制約条件②をQUBOとして、定式化するまでの流れを説明します

まず以下の制約条件②を表す不等式に対し、
$$q_0 + q_1 + q_2+ q_3≦3$$$r_k$を導入して、上の不等式と同じ意味を持つ等式制約に変形します。

\displaylines{
q_0 + q_1 + q_2+ q_3 = 0×r_0+1×r_1+2×r_2+3×r_3\\
 r_0+r_1+r_2+r_3=1
 }

上の式は、$q_0 + q_1 + q_2+ q_3=k$の時に、$r_k=1$が成立し、$q_0 + q_1 + q_2+ q_3≠k$の時、$r_k=0$が成り立ちます。
具体的に、$r_k$と、等式 $q_0 + q_1 + q_2+ q_3= r_1+2r_2+3r_3$の対応を以下の表に示します。

$r_0$	$r_1$	$r_2$	$r_3$	$q_0 + q_1 + q_2+ q_3=r_1+2r_2+3r_3$
$1$	$0$	$0$	$0$	$q_0 + q_1 + q_2+ q_3 = 0$
$0$	$1$	$0$	$0$	$q_0 + q_1 + q_2+ q_3 = 1$
$0$	$0$	$1$	$0$	$q_0 + q_1 + q_2+ q_3 = 2$
$0$	$0$	$0$	$1$	$q_0 + q_1 + q_2+ q_3 = 3$
$0$	$1$	$0$	$1$	$q_0 + q_1 + q_2+ q_3 = 4$ → $r_0+r_1+r_2+r_3=2$のため制約を満たさない

特に、表の最後の行に見ていただければ、2つの等式制約が、元の不等式制約$q_0 + q_1 + q_2+ q_3≦3$と同じ意味であることが分かると思います。

あとは、「制約条件が等式の時」で説明した通り、各等式に対し、移項して「$=0$」の形を作り、両辺を2乗した式の左辺の和をとることで、冒頭でも説明したQUBO形式の式を得ることができます。

$$(q_0 + q_1 + q_2+ q_3 - r_1-2r_2-3r_3)^2+ (1-(r_0+r_1+r_2+r_3))^2$$

問題の変数が2値変数でない時

次に、問題の変数が2値変数でない、つまり整数変数で定義されているときを考えます。
この時、問題をQUBO（もしくはイジングモデル）として定式化するために、整数変数を2値変数で置き換える必要があります。

例えば、以下の組み合わせ最適化問題を考えます。

目的関数：$4X_0-3X_1$
制約条件：$0≦X_0≦2,0≦X_1≦7$

ここで、$X_0,X_1$はそれぞれ、$0≦X_0≦2,0≦X_1≦7$を満たす整数とします。そのため、このままでは、2値変数の2次多項式であるQUBOとして、定式化できません。

そこで、$X_0,X_1$をうまく2値変数で表現する必要があるのですが、今回は、以下の2種類の方法を紹介します。

• 2値変数の和として表現
• 2値変数を利用し、2進数で表現

2値変数の和で表現

整数変数を、バイナリ変数$q_k$の和で表現します。
具体的には、以下のように置き換えます。

\displaylines{
X_0 = q_0 + q_1+ q_2\\
X_1 = q_3 + q_4 + q_5 + q_6 + q_7 + q_8 + q_9
}

整数変数の表現に利用する2値変数の数は、制約条件にある、変数の定義域によって、決定されます。($0≦X_0≦2$なので、$X_0$は3つ、$0≦X_1≦7$、なので$X_1$は7つの2値変数を利用している。)

以上の2値変数表現を利用して、最終的に以下のQUBO形式の式を得ることができます。

4(q_0 + q_1+ q_2)-3(q_3 + q_4 + q_5 + q_6 + q_7 + q_8 + q_9)

2値変数を利用し、2進数で表現

整数変数を2値変数の和として表現という方法は、多くの2値変数を必要とします。特に、$X_1$を表現するだけ7変数を利用しており、コスパが悪いように感じてしまいます。

そこで、変数の数を節約するために、整数変数がある条件を満たすのならば、バイナリ変数を利用し、2進数で表現するという方法が有効です。

ここで、「ある条件」とは、整数変数を$X$、$m,n$ともに整数とすると、以下のいずれかを満たすことです。

\displaylines{
2^m≦X≦2^n-1　(m<n)\\
0≦X≦2^n-1　
}

つまり、整数変数が、取りうる最小値が2の累乗、または0、最大値が(2の累乗)-1を満たす、2値変数の2進数表現の利用によって、利用する変数の数を減らすことができます。

例えば、$X_1$の定義域は、$0≦X_1≦7→0≦X_1≦2^3-1$ であり、上の条件を満たすため、以下のように表現できます。

X_1 = q_3×2^2 + q_4×2 + q_5  

ここで、$q_k$は2進数の$k$桁目の値を示しています。

具体的に、$q_k$と$X_1$の対応は以下の表のとおりです。

$q_3$	$q_4$	$q_5$	$X_1$
$0$	$0$	$0$	$0$
$0$	$0$	$1$	$1$
$0$	$1$	$0$	$2$
$0$	$1$	$1$	$3$
$1$	$0$	$0$	$4$
$1$	$0$	$1$	$5$
$1$	$1$	$0$	$6$
$1$	$1$	$1$	$7$

上で説明した$X_1$のバイナリ変数の2進数での表現を反映させると、目的関数は以下のようになります。

目的関数：$4(q_0 + q_1+ q_2)+3(q_3×2^2 + q_4×2 + q_5)$

$X_1$を2値変数による2進数で表現した結果、利用する変数の数を4つ減らすことができました。

このように、整数変数を2進数で表現することで、比較的少ない2値変数で、定式化ができます。

制約条件や目的関数に3次以上の式が含まれる

QUBOやイジングモデルは2次多項式であるため、目的関数や制約条件が2次、および1次の問題をこの記事では扱ってきましたが、3次以上の問題に関しても、うまく2次多項式の形式に落とし込むことで、量子アニーリングマシンで扱うことができます。

具体的には、3次以上の項に含まれる2つの変数の積を新たな変数で置き換えることで、QUBOとして、定式化を行います。

例えば、以下の3次⁴の項を目的関数に含む問題を考えます。

目的関数：$q_0q_1q_2+q_0q_1-q_2$
制約条件：無し

ここで、$q$はバイナリ変数とします。以上の問題をQUBOとして定式化することを考えます。

まず、$q_1q_2=q_3$として、目的関数を以下のように置き換えます。

$$q_0x_3+q_0q_1-q_2$$

これにより、見かけ上は目的関数が2次となりましたが、新たに、$q_1q_2=q_3$を問題の制約条件として加える必要があります。しかし、「制約条件が等号の時」で紹介した方法で、定式化しようとすると、$(q_1q_2-q_3)^2$となり、4次式となってしまいます。

そのため、$q_1q_2=q_3$を満たす際は値が0、逆に満たさないときは値が大きくなるような、QUBO形式の定式化を、自分で独自に行う必要があります。

今回は書籍⁵で紹介されている定式化を示します。
$$3q_3+q_1q_2-2q_1q_3-2q_2q_3$$

$q_1,q_2,q_3$と上の式の値の対応を以下表に示します。

$q_1$	$q_2$	$q_3$	$$3q_3+q_1q_2-2q_1q_3-2q_2q_3$$
$0$	$0$	$0$	$0$
$0$	$1$	$0$	$0$
$1$	$0$	$0$	$0$
$1$	$1$	$1$	$0$
$1$	$1$	$0$	$1$
$0$	$0$	$1$	$3$
$0$	$1$	$1$	$1$
$1$	$0$	$1$	$1$

$q_1q_2=q_3$を満たす際（表の上4列）は、$3q_3+q_1q_2-2q_1q_3-2q_2q_3=0$で、満たさないときは、式が0より大きい数をとるということが、表を見るとわかると思います。

以上をまとめると、QUBOとして以下のように定式化できます。
$$q_0x_3+q_0q_1-q_2+A(3q_3+q_1q_2-2q_1q_3-2q_2q_3)$$

ちなみに・・・

今まで、量子アニーリングを利用して、求解を行うために、問題をQUBOとして定式化する必要がある！！ということを述べてきましたが、実はこの定式化をサポートする、SDKやライブラリがいくつか提供されています。

例えば、D-waveにも対応している、Amplify SDK(Fixstars社)を利用する場合、マシンの実行の際の入力（コード上の変数 model）を以下のように直感的に定義できます。

• 制約条件が等式、不等式で表現される時

目的関数：$4q_0q_1-3q_1q_3+2q_2-q_3$
制約条件①：$2q_0 + q_1+q_2=2$
制約条件②：$q_0 + q_1 + q_2+ q_3≦3$

制約条件が等式、不等式で表現される時

from amplify.constraint import equal_to, less_equal
from amplify import gen_symbols, BinaryPoly

q = gen_symbols(BinaryPoly, 4) #バイナリ変数を4つ生成

cost = 4*q[0]*q[1] - 3*q[1]*q[3] + 2*q[2] - q[3] # 目的関数
c1 = equal_to(2*q[0] + q[1] + q[2], 2) # 制約条件①
c2 = less_equal(q[0] + q[1] + q[2] + q[3], 3) # 制約条件②

A=4 # 制約項の重み
model = cost+A*(c1+c2)　# マシンの入力の設定

• 3次以上の式が含まれる時

目的関数：$q_0q_1q_2+q_0q_1-q_2$

3次以上の式が含まれる時

q = gen_symbols(BinaryPoly, 3) #バイナリ変数を3つ生成

cost = q[0]*q[1]*q[2] + q[0]*q[1] - q[2] # 目的関数

model = cost #マシンの入力の設定

そして、最後に以下を記述⁶することで、マシンを実行できます。

result = solver.solve(model) # マシンの実行

このように、この記事で紹介したQUBO（イジングモデル）の定式化の方法を知らなくても、Amplify等のツールの利用によって、量子アニーリングマシンの実行ができてしまいます。（楽ですね！）

となると、この記事の意味は・・・ということになりますが、QUBO（イジングモデル）としての定式化方法を知っておくことは、量子アニーリングマシンを実行する上で、悪いことではない（むしろ、知っておくべき？）ことだと思うので、是非、この記事で紹介したような定式化方法を理解したうえで、量子アニーリングを利用していただければと思います。

まとめ

この記事では、量子アニーリングマシンを用いた、組み合わせ最適化問題の求解の流れ、特に、「問題を、どのようにイジングモデルやQUBOとして定式化するか」について、重点的に説明しました。

最後に、改めて、記事で紹介した、定式化の際にに少し工夫が必要な条件と、その定式化する上でのポイントを以下に示します。

• 制約条件が存在する時
  →制約条件を満たす際は値が0、逆に破った際にはペナルティが発生するように定式化

  • 制約条件が等号の時
    →移項して「=0=0」の形を作り、両辺を2乗する

  • 制約条件が不等号の時
    →制約条件の左辺の和がkkがどうか」を示す、2値変数を追加で導入

• 問題の変数が$(0,1)$、$(-1, +1)$のどちらでもない時
  　→整数変数を2値変数で置き換える（2進数表現もアリ）

• 制約条件や目的関数に3次以上の式が含まれる時
  　→3次以上の項に含まれる2つの変数の積を新たな変数で置き換える

記事は以上になります。最後までお読みいただきありがとうございました。
もし、なにか質問、ご指摘等ございましたら、コメント欄までお願いいたします。

1. https://annealing-cloud.com/ja/tutorial/1.html ↩

2. https://www.frontiersin.org/articles/10.3389/fphy.2014.00005/full ↩

3. https://amplify.fixstars.com/ja/techresources/research/ising-model-formulation/ ↩

4. 4次以上の項を含む場合、変数の置き換えを再帰的に行うことで、定式化をできます。これについて、こちらが参考になるかと思います。 ↩

5. 書籍：量子アニーリングの基礎（https://www.kyoritsu-pub.co.jp/bookdetail/9784320035386） ↩

6. 説明の都合上省略していますが、実際はソルバー実行の前に、アクセストークンや実行するマシンの設定、また、Solverオブジェクトの定義をする必要があります。詳しくは、Amplify SDKの公式のチュートリアル などをご覧ください。 ↩