NP完全問題を量子コンピュータと同じ方式で解いてみる

概要
自然数分割問題
D-WAVEの計算方式
自然数分割問題をqbsolvを使って解く
まとめ

概要

NP完全問題の一つである自然数分割問題をD-WAVEが公開しているOSSを使用してD-WAVEを用いて問題を解く時と同じ方式で解いてみます．

自然数分割問題

ある自然数の集合が与えられたときに，その自然数の集合をそれぞれの集合に含まれる自然数の和が等しい2つの集合に分割する問題です．

例えば以下のような自然数の集合が与えられたとします．

8 2 4 3 5 11 6 9 1 7 10

この場合はともに合計が33となる次の2つの集合に分割することができます．

A:2, 5, 6, 9, 1, 10
B:8, 4, 3, 11, 7

このような分割を見つけるのがこの問題です．

D-WAVEの計算方式

D-WAVEの計算のアーキテクチャは次のようになっています．
f:id:kadora:20180516203706p:plain
(出典: Software | D-Wave Systems )

この図の通り，まずは解きたい問題をQUBOに変換しそれをD-WAVEのAPIに投げることでD-WAVEがそのQUBOを解いてくれるという感じになっています．

なので問題をQUBOに変換し，そのQUBOの解から問題の解答を作成することで，実際にD-WAVEを使用する時と同じような感じで問題を解くやり方をなぞることができます．

イジングモデル

QUBOについて説明する前にまずイジングモデルについて説明します．イジングモデルとは次のようなエネルギー関数で表されるモデルのことです．

$\displaystyle H = \sum_{i \ne j}J_{ij}\sigma_i\sigma_j+\sum_ih_i\sigma_i (\sigma = \pm 1)$

ここで $\sigma$ は入力変数を， $J_{ij}$ は2つの入力変数間の相互作用のパラメータを， $h_i$ は，各入力変数のパラメータを表します．

QUBO

先程のイジングモデルは，入力変数が $\pm 1$ を取るモデルだったのでそれを計算機での計算に都合が良いように0か1を入力変数として持てるモデルに変換します．
イジングモデルのエネルギー関数で， $\sigma=2(q-1)$ を使って変数変換することで，等価なQUBOのエネルギー関数を得ることができます．

イジングモデルやQUBOについては次のページも参考にしてください．

quantum.fixstars.com

D-WAVEのマシンはこのQUBOのエネルギー関数を最小にするような入力変数の組を計算の結果として出力します．

qbsolv

qbsolvはD-WAVEが公開しているこのQUBOのソルバです．
github.com
実際にはD-WAVEが一度に扱えるビット数に限りがあるのでこのqbsolvは問題のサイズが大きいときに問題を分割することも行っています．

設定をすることで実際にD-WAVEマシンを使用することもできますが，設定しない場合にはclassicalなソルバでQUBOを解いてくれます．今回は残念ながら使用できるD-WAVEマシンがないので特に設定はしないで使用します．

自然数分割問題をqbsolvを使って解く

実際に問題を解いているコードは次のリポジトリになります．
github.com
パスが通っている場所にqbsolvがある状態でリポジトリ内のexample.shを実行することで実際に実行して問題を解くことができます．qbsolvのビルド方法などはqbsolvのリポジトリを参考にしてください．

以下で実際に問題を解く際のステップを順に追っていきます．

自然数分割問題をイジングモデルに変換

まずは自然数分割問題をqbsolvを使用して計算できる形式へ変換するためにイジングモデルへ変換します．

一発でQUBO形式に変換出来るのであればそれでもいいのですが，今回は参考にした論文がイジングモデルへの変換までなのでまずはイジングモデルへ変換します．

参考にした論文はこちらです．

Ising formulations of many NP problems

https://arxiv.org/pdf/1302.5843.pdf

与えられたi番目の自然数が分割した時にどちらのグループに属するかを $\sigma_i\in \{1,-1\}$ , 自然数の値を $s_i$ として表すとハミルトニアンが次のように表せます．

$\displaystyle H = A(\sum_{i=1}^N\sigma_i s_i)^2$

このハミルトニアンは $\sigma_i=1$ と $\sigma_i=-1$ の2つのグループに分割された自然数の合計が等しくなるときに0となり最小値となるのでこのハミルトニアンが最小値を取るときに問題の答えを表しているのでこのハミルトニアンが問題を解くために最小値を求めたいハミルトニアンとなります．

イジングモデルからQUBOへ変換

前述の通り $\sigma=2(q-1)$ を用いて変換することでイジングモデルと等価なQUBOへ変換することができます．
これを用いて変換した結果が次のとおりです．

$\displaystyle H = 2 \sum_i^{N-1}\sum_{j=i+1}^N s_i s_j q_i q_j - \sum_i^N s_i (\sum_{j=1}^N s_j - s_i) q_i$

ハミルトニアンの最小値を与えるqの値を求めるのが目的なので正の定数倍やqに関わらない項を除外しています．

サンプルのコードではこちらのコードで問題の入力からQUBOへの変換を行っています．

qbsolvの実行

先程求めた式により問題を解くために最小値を求めたい式が得られたので問題の入力からqbsolvを用いて解くQUBOを生成してqbsolvで解きます．
qbsolvの実行を行っているのはexample.sh内の以下の部分です．

qbsolv -i partition.qubo -o partition.qbout -v1 -n 8

iで解く対象のQUBOが記述されたファイルを指定し，oで結果を出力するファイルを指定します．

vは出力の詳細のレベルを，nは最小値を求めるための計算の繰り返しの回数を指定します．

その他にもオプションを指定することができますが，詳しくはqbsolv自体のページに説明があるので興味がある方はそちらを参照してください．

github.com

qbsolvの結果から問題の解答へ変換

qbsolvで出力されたファイルを見てみると次のようになっています．

11 bits,  find Min, SubMatrix= 47, -a o, timeout=2592000.0 sec
01001011101
-1089.00000 Energy of solution
0 Number of Partitioned calls, 1 output sample
 0.00182 seconds of classic cpu time

このうち2行目がハミルトニアンを最小にするqの値を，3行目がハミルトニアンの最小値を表しています．
分割可能であれば，最終的な分割の解答は2行目と入力を用いて，2行目のビットが0になっている位置の自然数の集合と，1になっている位置の自然数の集合に分割するとそれが解答となっています．
分割可能であるかどうかは分割した集合のそれぞれの合計が等しいかどうかで判定します．

ここまでの一連の流れをサンプルコードのexample.shを実行すると全て実行してくれて，結果が次のようになります．

% ./example.sh
Partitioning result

GroupA:[2, 5, 6, 9, 1, 10]

GroupB:[8, 4, 3, 11, 7]

まとめ

今回は実際にアニーリング方式の量子コンピュータで問題を解く際に必要となるイジングモデル（QUBO）への問題の変換と，その解を解釈して問題の解へ変換することをやってみました．

実際の計算部分は量子コンピュータを使うことはできませんでしたが，解かせる問題の形式などは同じものを使用しているため実際に量子コンピュータを使用する際にも同じような感じで利用することができるはずです．

D-WAVEに限らず，幾つかの企業がアニーリング方式の量子コンピュータを開発しているのでそれらの量子コンピュータを使用する際にも今回のようなやり方で使用できるはずです．

問題をイジングモデルへ変換するのが中々大変ですが，ゲート方式の量子コンピュータよりもこちらのアニーリング方式の量子コンピュータの方が実用化が早いと言われているので，このようなやり方に慣れておくといざ本当に量子コンピュータが使えるようになった時にすぐにガンガン量子コンピュータを使っていけるようになるかもしれません．