量子アニーリングマシンD-Waveで良い解を出すためのコツ（前編）

# この記事の概要

これまでエンジニアにとって身近に感じることが難しかった量子コンピュータですが、AWSが提供するAmazon Braketのローンチを皮切りにエンジニアにとっても身近なものとなりつつあると感じています。

そこで、量子アニーリングマシンとして世界的に最も有名でありAmazon Braketでも利用ができる**D-Wave**を使ってみたのですが、良い解を求めるためにチューニングする際のコツがありそうだったのでまとめてみました。

要約すると、

* チューニングポイントは「制約項の係数」「チェーン強度」「サンプリング回数」
* 制約項の係数を調整するにはPyQUBOのプレースホルダーをうまく活用する
* チェーン強度は強くするべきだが、QUBO内の値を意識しながら強くしすぎないように調整する
* サンプリング回数は多ければ多いほどいいが、金額の問題や上限が10000回であることに注意する

といった内容になります。

なお、私自身量子コンピュータについて勉強中の身ですので、間違ったことを書いている可能性もありますが、その場合は優しく指摘いただければと思います。

# D-Waveを使うにあたって
## D-Waveについて
<details><summary>知っている人は飛ばしてOKです</summary><div>
ここ最近は汎用的な量子計算ができるとされている「量子ゲートマシン」が界隈を賑わせている印象がありますが、D-Waveはそれとは異なり、組合せ最適化計算に特化した「量子アニーリングマシン」です。
（もしかしたら表現が不適切かもしれないですが、量子ゲートマシンと量子アニーリングマシンの対応関係は、汎用的な処理ができるパソコンと四則演算に特化した電卓の対応関係に似ているような気がしています。）

組合せ最適化というジャンルで代表的に上げられるものが[巡回セールスマン問題](https://ja.wikipedia.org/wiki/%E5%B7%A1%E5%9B%9E%E3%82%BB%E3%83%BC%E3%83%AB%E3%82%B9%E3%83%9E%E3%83%B3%E5%95%8F%E9%A1%8C)です。
このような爆発的に解の組合せが増えてしまうために逐次的に解を探索することが難しい問題を得意としているものが、D-Waveのようなマシンとなります。

そんなD-waveは「D-Wave Advantage」という5760Qubitを搭載したマシンをリリースしており、トポロジ（Qubit同士の接合方式のようなもの）もより効率の良い方式となってきています。
これにより、今まで解くことができなかった規模の大きな問題にも対応できるようになりつつあります。
とは言っても、10都市強の巡回セールスマン問題を解くことができる程度ですが…。

そんなD-Waveですが、PoCではあるものの企業に利用されてきていたり、今後更に搭載Qubitを増やしていく予定となっていたりと良いニュースも増えています。
</div></details>

というわけで、D-Waveを軽く使ってみたいなと思い触ってみたものの、なかなかいい解を導けないぞ？と思い、色々調べながら試行錯誤した結果なんとなく導き出したコツをここにまとめます。

## D-Waveで求解をする流れ
流れとしては以下の通りです。

* 問題を定式化する
* 定式化したものをQUBOに変換する★1
* QUBOを入力としてD-Waveに求解させる
  * Qubitに問題の埋め込みをするための前処理が行われる★2
  * 量子トンネル効果を利用して求められた変数の組合せを複数回観測する★3
* 観測結果のなかで最もエネルギーが低い変数の組合せを最適な組み合わせとして出力する

この中で今回工夫してみる部分が★のところです。

★1では、制約項をQUBOに入れ込むときの係数の調整を行います。
★2では、Qubitの埋め込みをするときに、1つの変数を複数のQubitに割り当てるような動作をするのですが、論理的に同じ変数として割り当てられたQubitがバラバラにならないようにチェーン強度というものを調整します。
★3では、量子効果は確率的な現象なので必ずしも一発で正しい解が出るとは限らないので、なるべく多く観測を行うよう設定をします。

詳細については後述していきます。

# 試しにやってみる
## この記事での前提
### 環境
プラットフォーム：AWS（AmazonBraket）
ソルバ：D-Wave Advantage_system1.1
ノートブックインスタンス：ml.t3.medium

### 求解対象
今回は以下のような6都市のTSPを求解します。

```
cities = [
    ('0', (0,0)),
    ('1', (17,8)),
    ('2', (23,0)),
    ('3', (27,15)),
    ('4', (70,20)),
    ('5', (52,0))
]
```
こんな感じの格好の配置となります。
![2021-02-26_094704.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214621/293ff3ba-5226-9f73-b0bb-a5d794a4b5bc.png)

## ★1 制約項の係数を調整する
定式化した問題をQUBOに落とし込む際に制約項もQUBOに入れる必要があります。
そもそもQUBOというものは「制約なし二次形式二値変数最適化（Quadratic Unconstrained Binary Optimization）」の略です。（詳しくは[こちら](https://qard.is.tohoku.ac.jp/T-Wave/?glossary=qubo-%E5%88%B6%E7%B4%84%E3%81%AA%E3%81%97%E4%BA%8C%E6%AC%A1%E5%BD%A2%E5%BC%8F%E4%BA%8C%E5%80%A4%E5%A4%89%E6%95%B0%E6%9C%80%E9%81%A9%E5%8C%96)）
…え？制約なしって言ってるのに制約を考える必要があるの？と思った方は勘が鋭いですね。
実はここで言っているのは、制約を別式として明示的に設定しませんという意味なのです。
つまり、もし制約がある場合は**制約を制約項として目的関数の中にうまく含める**必要があります。

定式化についての細かい解説は他に譲らせていただきますが、今回TSPを解くにあたっての制約はおなじみの以下2つです。

* 一度に二箇所以上へ訪問しない

```math
\sum_{a}^{n}x_{a,i} = 1
```

* 各都市必ず一度だけ訪問する

```math
\sum_{i}^{n}x_{a,i} = 1
```

一方、TSPにおける目的関数は距離を最小化するような式となるので以下のような式となります。
この値は小さければ距離が短いことを示すので、小さければ小さいほどよいということになります。
（これを距離項と言ったりしますね。）

```math
\frac{1}{2}\sum_{i,j,a}^{n}d_{i,j}x_{a,i}(x_{a+1,j}+x_{a-1,j})
```

この式に対して先程の2つの制約を制約項として加えるとすると、それぞれの制約項は次の通りとなります。
この値がそれぞれ**0となっていること**が**制約を満たしていること**を示します。

```math
\sum_{i}^{n}(\sum_{a}^{n}x_{a,i} - 1)^2\\
\sum_{a}^{n}(\sum_{i}^{n}x_{a,i} - 1)^2
```

あとは距離項と制約項を足したハミルトニアン$H$をQUBO行列化します。
ハミルトニアン$H$はこのように書くことができます。

```math
H=\frac{1}{2}\sum_{i,j,a}^{n}d_{i,j}x_{a,i}(x_{a+1,j}+x_{a-1,j})+\lambda(\sum_{i}^{n}(\sum_{a}^{n}x_{a,i} - 1)^2)+\sum_{a}^{n}(\sum_{i}^{n}x_{a,i} - 1)^2)
```

前置きが長くなりましたが、ここで重要となってくるのが$\lambda$のチューニングです。
制約項の係数である$\lambda$を調整する際に便利なのが**[PyQUBO](https://github.com/recruit-communications/pyqubo)**です。

通常、QUBO行列を生成するときにfor文を何個も書いてその中に制約項の係数を入れたり入れなかったりしながらQUBO行列を生成するのですが、可読性が悪く、どこに制約項の係数がかかってくるのかとか数式で言うとどこの部分を指しているのかとかがわかりにくいソースコードとなります。
これにより、制約項の係数をチューニングしたくても簡単にできない点が課題となります。
しかし、PyQUBOを使うことで数式をほぼそのままの形でQUBO化してくれる上に、制約項の係数のチューニングをしやすい仕掛けも備えています。
その仕掛けが**プレースホルダー機能**です。
プレースホルダー機能を使ってQUBO化するまでのソースはこんな感じです。

```python
# 一度に二箇所以上へ訪問しない
time_const = 0.0
for i in range(n_city):
    time_const += Constraint((sum(x[i, j] for j in range(n_city)) - 1)**2, label="time{}".format(i))

# 各都市必ず一度だけ訪問する
city_const = 0.0
for j in range(n_city):
    city_const += Constraint((sum(x[i, j] for i in range(n_city)) - 1)**2, label="city{}".format(j))

# 距離コスト
distance = 0.0
for i in range(n_city):
    for j in range(n_city):
        for k in range(n_city):
            d_ij = dist(i, j, cities)
            distance += 1/2 * d_ij * x[k, i] * (x[(k+1)%n_city, j]+x[(k-1)%n_city, j])

# ハミルトニアン生成
A = Placeholder("A")
H = distance + A * (time_const + city_const)

# モデルコンパイル
model = H.compile()

# QUBO生成
feed_dict = {'A': 100.0}
qubo, offset  = model.to_qubo(feed_dict=feed_dict)
```

このように`A`というプレースホルダーを使って制約項の係数$\lambda$を仮で置くことができます。
こうすることで後から微妙に調整してもう一回QUBOを生成したいとか、何パターンか$\lambda$を設定してQUBOを生成して求解結果がいいものを採用したいとかいう場合に、ここ↓のブロックだけを実行すれば良くなります。

```
# QUBO生成
feed_dict = {'A': 100.0}
qubo, offset  = model.to_qubo(feed_dict=feed_dict)
```

で、私自身D-Waveを実際に動かしてみるとなかなかいい解にはたどり着かず、ケースバイケースで制約項の係数をちょっとずつ調整をしたり、制約項の係数を変えた場合の解の比較をするような場面も結構あり、PyQUBOの機能に度々助けられました。
こんな感じで、**PyQUBO**はとても便利なのです。

# 前編はここまで
後編では★2、★3の部分についての解説をしていきます。
特に個人的に激アツなのは★2の部分なのでぜひご覧になっていただければと思います:bow:
（[後編](https://qiita.com/izey0306/items/83cc3b67abba79007a95)に続く）

なお、量子コンピュータ関係の他の記事は、[こちらのページで一覧にしています](https://qiita.com/shnchr/items/41fa7fd1284d74007a7c)。
よろしければこちらの記事も参照いただければと思います！