n個のうちmだけを１にする問題をPyQUBOを使ってアニーリングマシンで解く

第３回目は、変数n個のうちmだけを１にする問題をPyQUBOを使ってアニーリングマシンで解いてみようとおもいます。

アニーリングマシンでは、イジング模型として、一番低いエネルギーの状態（ハミルトニアン）すなわち基底状態を探ることができるマシンですので、それをうまく応用することで、世の中にある問題を解決することが可能です。

これを、[第１回](https://qiita.com/morimorijap/items/d7afda9452036c53caed)で、量子アニーリングでは、イジング模型として、エネルギーの状態をハミルトニアン$H$で次の方程式のように表せると書き、

$$H = \sum_{i<j} J_{ij} \sigma_i \sigma_j + \sum_i h_i \sigma_i $$

そして、[第２回](https://qiita.com/morimorijap/items/9d99abed3f6284710337)で、QUBO形式でこの式を表す場合についてふれました。

$$H = \sum  Q_{ij} q_i q_j + \sum_i b_i q_i $$
$$q_i = 
\{ 0, 1 \} $$

アニーリングマシンで計算をするために、この形式へ変換する作業が必要になってきます。

まずはじめは、　２つのうち１つだけを１にする問題を考えます。

<img width="321" alt="Screen Shot 2019-04-16 at 13.42.28.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/260261/450e4e1e-a2a5-8199-f00d-322e97a6a49e.png">

$ 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 $を
両方満たし、両方のエネルギーが最小になるような式を考えます。

<img width="331" alt="Screen Shot 2019-04-16 at 13.42.34.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/260261/8543b234-2839-a99e-7a74-2ac24ab913fe.png">

３個になると、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$

となります。

そして、４つのうち１つだけを１にする問題を考えてみます。
こちらも同じように考えていくと、

<img width="314" alt="Screen Shot 2019-04-16 at 13.58.30.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/260261/7e762e38-21bd-3904-4bd7-f5d981e07168.png">

の場合は、

$$ H = 2(q_0 * q_1 + q_1 * q_2 + q_0*q_2 + q_0*q_3 + q_1*q_3 + q_2*q_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 = 2*q_0*q_1- q_0 - q_1 + 1 $$
とすると、 ３個の場合、そして４個の場合と見ていくと、一般式に展開できるようになると思います。

これを利用して、４個のうち１つを１にする問題をPyQUBOをつかって解いてみようと思います。

$$ H = 2(q_0 * q_1 + q_1 * q_2 + q_0*q_2 + q_0*q_3 + q_1*q_3 + q_2*q_3 )  
- ( q_0 + q_1 + q_2 + q3 ) +1 $$

この式を展開して整理すると

$$ H = 2 * q_0*q_1 + 2* q_0*q_2 +  2 * q_0*q_3 - q_0 +  2 * q_1*q_2 + 2*q_1*q_3 - q_1 + 2*q_2*q_3 - q_2 - q_3 + 1 $$

となり、これで用意ができました。それでは、PyQUBOとD-Wave Leapを使って解いてみたいとおもいます。

といっても、D-Wave Leap、PyQUBO を使うために、インストールが必要です。

D-Wave Leapの実機を使うためには、[D-Waveからの登録](https://www.dwavesys.com/)とD-Waveの[Ocean Tools](https://docs.ocean.dwavesys.com/en/latest/overview/install.html)のインストールが必要になります。

pip install dwave-ocean-sdk

で、インストールは完了になります。
PyQUBOのインストールは、詳細は[ドキュメント](https://pyqubo.readthedocs.io/en/latest/)がをみればいいのですが、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のところからコピーすることが可能です。

<img width="1010" alt="Screen Shot 2019-04-14 at 21.38.31.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/260261/0569ca77-e44f-b495-14af-c712a79aedfb.png">

また、それ以外の　endpoint, solver_nameはページの下の方に

<img width="561" alt="Screen Shot 2019-04-14 at 21.37.58.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/260261/571b4f23-e0c0-54d3-8e23-45224a5f6691.png">

と書いてありますので、こちらを利用します。

これで、用意ができたので、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開発者の棚橋さんアドバイスいただきありがとうございました。