第３回目は、変数n個のうちmだけを１にする問題をPyQUBOを使ってアニーリングマシンで解いてみようとおもいます。
アニーリングマシンでは、イジング模型として、一番低いエネルギーの状態（ハミルトニアン）すなわち基底状態を探ることができるマシンですので、それをうまく応用することで、世の中にある問題を解決することが可能です。
これを、第１回で、量子アニーリングでは、イジング模型として、エネルギーの状態をハミルトニアン$H$で次の方程式のように表せると書き、
$$H = \sum_{i<j} J_{ij} \sigma_i \sigma_j + \sum_i h_i \sigma_i $$
そして、第２回で、QUBO形式でこの式を表す場合についてふれました。
$$H = \sum Q_{ij} q_i q_j + \sum_i b_i q_i $$
$$q_i =
{ 0, 1 } $$
アニーリングマシンで計算をするために、この形式へ変換する作業が必要になってきます。
まずはじめは、 ２つのうち１つだけを１にする問題を考えます。
$ q_0 = 1$のとき $ q_1 = 0 $ また、$ q_0 = 0$のとき $ q_1 = 1 $を
両方満たし、両方のエネルギーが最小になるような式を考えます。
$$ H = - q_0 - q_1 + 1 $$
$ q_0 = 1$, $ q_1 = 0 $ は $ H=0$
また、$ q_0 = 0$, $ q_1 = 1 $は $ H=0$
となります。
次に、３つのうち１つだけを１にする問題を考えてみます。
このときも同じように、
$ q_0 = 1$のとき $ q_1 = 0 $ また、$ q_0 = 0$のとき $ q_1 = 1 $を
両方満たし、両方のエネルギーが最小になるような式を考えます。
３個になると、Qubit間の相互作用を考えいかないとうまくいかないくなります。これを考えていき、エネルギー関数を
$$ H = 2*(q_0*q_1 + q_1 * q_2 + q_0 * q_2) - (q_0 - q_1 + q_2)+ 1 $$
とすると、
$ q_0 = 1$, $ q_1 = 0 $ $ q_2 = 0 $ は $ H=0$
$ q_0 = 0$, $ q_1 = 1 $ $ q_2 = 0 $ は $ H=0$
$ q_0 = 0$, $ q_1 = 0 $ $ q_2 = 1 $ は $ H=0$
となります。
そして、４つのうち１つだけを１にする問題を考えてみます。
こちらも同じように考えていくと、
の場合は、
$$ H = 2(q_0 * q_1 + q_1 * q_2 + q_0q_2 + q_0q_3 + q_1q_3 + q_2q_3 ) - ( q_0 + q_1 + q_2 + q_3 ) +1 $$
とすると、
$ q_0 = 1$, $ q_1 = 0 $ $ q_2 = 0 $ $ q_3 = 0$ は $ H=0$
$ q_0 = 0$, $ q_1 = 1 $ $ q_2 = 0 $ $ q_3 = 0$ は $ H=0$
$ q_0 = 0$, $ q_1 = 0 $ $ q_2 = 1 $ $ q_3 = 0$ は $ H=0$
$ q_0 = 0$, $ q_1 = 0 $ $ q_2 = 0 $ $ q_3 = 1$ は $ H=0$
となります。
実際、変数n個のうちmだけを１にするという問題はよく知られていて、
QUBOにおいては、
$$H = ( \sum_{i}^n q_{i} -m )^2 $$
とすることで表現することが可能です。
QUBOでは、 変数がq = 0 or 1 という性質を利用すると $ q^2 = s$となりますので、
$$H = \sum_{i \neq j } q_i q_j + (1-2m)\sum_{i} q_i + m^2 $$
と変形ができます。
実際２個のQubitの場合は、
$$ H = - q_0 - q_1 + 1 $$
としていましたが、
$$ H = 2q_0q_1- q_0 - q_1 + 1 $$
とすると、 ３個の場合、そして４個の場合と見ていくと、一般式に展開できるようになると思います。
これを利用して、４個のうち１つを１にする問題をPyQUBOをつかって解いてみようと思います。
$$ H = 2(q_0 * q_1 + q_1 * q_2 + q_0q_2 + q_0q_3 + q_1q_3 + q_2q_3 )
- ( q_0 + q_1 + q_2 + q3 ) +1 $$
この式を展開して整理すると
$$ H = 2 * q_0q_1 + 2 q_0q_2 + 2 * q_0q_3 - q_0 + 2 * q_1q_2 + 2q_1q_3 - q_1 + 2q_2*q_3 - q_2 - q_3 + 1 $$
となり、これで用意ができました。それでは、PyQUBOとD-Wave Leapを使って解いてみたいとおもいます。
といっても、D-Wave Leap、PyQUBO を使うために、インストールが必要です。
D-Wave Leapの実機を使うためには、D-Waveからの登録とD-WaveのOcean Toolsのインストールが必要になります。
pip install dwave-ocean-sdk
で、インストールは完了になります。
PyQUBOのインストールは、詳細はドキュメントがをみればいいのですが、pipは、 D-Wave Ocean SDKに内包されているので、さきほどのpip install dwave-ocean-sdkするだけで不要になります。
と、D-Wave LeapのWEBから、
endpoint = "https://cloud.dwavesys.com/sapi"
token = "DEV-xxxxxxxxxxxxxxxxxx"
solver_name = "DW_2000Q_2_1"
が必要です。
tokenは、API Tokenのところからコピーすることが可能です。
また、それ以外の endpoint, solver_nameはページの下の方に
と書いてありますので、こちらを利用します。
これで、用意ができたので、PyQUBO に落としていきます。
from pyqubo import Binary, solve_qubo
from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite
# Define hamiltonian
q0, q1, q2, q3 = Binary("q0"), Binary("q1"), Binary("q2"), Binary("q3")
# Hamiltonian Eq.
H = 2 * q0*q1 + 2* q0*q2 + 2 * q0*q3 - q0 + 2 * q1*q2 + 2*q1*q3 - q1 + 2*q2*q3 - q2 - q3 + 1
# Create QUBO
model = H.compile()
qubo, offset = model.to_qubo()
#print(qubo) ## QUBO形式の表示用
# Solve QUBO model by SA
solution = solve_qubo(qubo)
print(solution)
print("")
# Solve QUBO model by D-Wave
endpoint = "https://cloud.dwavesys.com/sapi"
token = "DEV-xxxxxxxxxxxxxxxxxx" ##<-- ここは自分のKeyを入れる
solver_name = "DW_2000Q_2_1"
sampler = EmbeddingComposite(DWaveSampler(endpoint=endpoint, token=token, solver=solver_name))
response = sampler.sample_qubo(qubo)
print(response)
これの計算結果は
(pyqubo) mori-mba:src morikawa$ python 4n-m-1-v1.py
{'q0': 0, 'q1': 0, 'q2': 0, 'q3': 1}
q0 q1 q2 q3 energy num_oc. chain_.
0 0 0 1 0 -1.0 1 0.0
['BINARY', 1 rows, 1 samples, 4 variables]
となります。
上の行は、 PyQUBOでおこなったSAの結果で、下方がD-Waveからの結果になっています。
４個のうち１つが１ということで、結果がちゃんとでていますが、同じ量子ビットが１でないのは、この式が特にどれを１にすると指定してないからです。
今回は、
$$H = ( \sum_{i}^n q_{i} -1 )^2 $$
式を展開してましたが、PyQUBOでは、そのままの式を
H = (q0+q1+q2+q3-1)**2
のようにしても、自動的に展開してQUBOファイルを作ってくれるとても便利なツールとなっています。
また、１回の計算だと、４個のうちどれが１になるか、わからなく、最終的には４個が同じ確率で１になるようになるはずです。 何回もこの計算をD-Waveに投げるとそういう結果がでるはずですが、 今日は長くなったので、D-Waveのことと、PyQUBOの詳細は次回にでも書きたいとおもいます。
更新履歴：
初版公開v1.0: 2019/4/15
v1.1: わかりにくい部分があったため、2019/4/16に大幅に修正をいれました。
PyQUBO開発者の棚橋さんアドバイスいただきありがとうございました。