D-Waveでボルツマン機械学習する際のテクニック集

# はじめに

量子アニーリングマシンな D-Wave には、組み合わせ最適化問題への応用が期待されていますね。
ただ実際のところ理想的な断熱時間発展はできないので、最適解が得られるとは限りません。むしろ現状だと「[D-Waveによる量子アニーリングで３行４列のピクロスをノリで解いてみた](https://qiita.com/piyo7/items/5e2cf2a92ddbe016572b)」のような問題サイズでも正答率 1/3 とかそんなものみたいです。色々とチューニングの余地はあると思うのですが、最適化問題を解かせるのはなかなか厳しそうという感触です。

そこで厳密な最適解を得るのは諦めて、 D-Wave によるサンプリングがボルツマン分布に従っているとみなし、ボルツマン機械学習する際の重いサンプリングコストを担わせようという流れがあります。
といってもボルツマン分布に従う理論的な後ろ盾はないようなのですが、ディープラーニングの事前学習として制限付きボルツマンマシンを学習する際に D-Wave を利用したところうまくいったという論文「[Application of Quantum Annealing to Training of Deep Neural Networks](https://arxiv.org/abs/1510.06356)」が面白そうだったので輪読してみました。
ありがたいことに D-Wave の実機に合わせた実用的なテクニックが解説されていたので、それらを紹介しつつ D-Wave の実行パラメータとどう対応しているか調べてみたよ、というのがこの記事です。MDR社の「[D-Waveで深層学習の基礎となるRBMのボルツマン学習を実行してみた](https://blog.mdrft.com/post/258)」と合わせて読むと分かりやすいかも。

# キメラグラフへの埋めこみ

論文の図が分かりやすいので、そのまま引用します。

<img width="856" alt="Figure 3" src="https://qiita-image-store.s3.amazonaws.com/0/55208/2a5c09fb-82ce-c99d-ab3a-97428b0a1d6a.png">

(a) が D-Wave のチップ図です。緑ないし青の棒が量子ビットで、たとえば赤線でマークしたものが 1 量子ビットになります。これが横 4 本と縦 4 本で交差してセルを構成します。それを 4 セル描いたのが左上図ですね。 D-Wave のチップについて詳しく知りたい方は「[Introduction to the D-Wave Quantum Hardware](https://www.dwavesys.com/tutorials/background-reading-series/introduction-d-wave-quantum-hardware)」を参照してください。
さて量子ビットは接続点で相互作用するのですが、量子ビットを点として相互作用を線として描きなおしたのが右上図です。いわゆるキメラグラフですね。ここで 1 セルがすでに 4x4 の制限付きボルツマンマシンの形をしていますが、もっと大きい 8x8 のような制限付きボルツマンマシンをどう表現するのかというと、

<img width="740" alt="Figure 4" src="https://qiita-image-store.s3.amazonaws.com/0/55208/c2cd61dc-5dc7-21f8-592a-40e7937fe9ed.png">

というふうに、同じ方向に繋がっている量子ビットたちをチェインとして、一つのノードとして扱ってしまいます。すなわち、それらの子ビットたちが同じスピン状態になるように相互作用をかけます。よって上図は 8x8 の制限付きボルツマンマシンを 32 量子ビットに埋めこんだことになります。同様に元論文の実験では 32x32 の制限付きボルツマンマシンを 512 量子ビットに埋めこんでいます。制限付きボルツマンマシンの埋めこみはかなりシンプルなのですが、それでもわりと膨らんではしまいますね……。

## ライブラリによる埋めこみ

では具体的にどうやって D-Wave 実機に埋めこむかというと、公式ライブラリ dwave-system の [FixedEmbeddingComposite](https://docs.ocean.dwavesys.com/projects/system/en/latest/reference/composites/fixed_embedding.html) を使うのが良さそうです。引数として、イジングモデルのノードをどこの量子ビットたちにマッピングするかの情報を渡します。自動でマッピングしてほしい場合は [EmbeddingComposite](https://docs.ocean.dwavesys.com/projects/system/en/latest/reference/composites/embedding.html) を使いましょう。
なおイジングモデルの重みは、チェイン内で均等割り付けされるようです。

```python
        b = bias / len(chain)
...
        b = bias / len(available_interactions)
```

[dwave-system/dwave/embedding/transforms.py](https://github.com/dwavesystems/dwave-system/blob/f1cc9a956b53a9af03d9ac85b44fd6ff26bbd165/dwave/embedding/transforms.py)

## 壊れた量子ビット

D-Wave の量子ビットはいくつか壊れていて、元論文では 512 量子ビット中 8 量子ビットほど壊れていたそうです。実際 D-Wave の [List Solvers](https://docs.dwavesys.com/docs/latest/c_rest_api_3.html#list-remote-solvers) API を訊いてみると、いくつか欠番あるのが分かります。世知辛いですね。
元論文では、制限付きボルツマンマシンは結合いくつか抜けてもわりとうまく学習してくれるので、気にせず埋めこんでいいとあります。なるほど……。

## エネルギーのずれ

キメラグラフ埋めこみ時にチェインの制約が入るので、制限付きボルツマンマシンのエネルギーと D-Wave 上のエネルギーがずれてしまう問題があります。あるのですが制約が守られているかぎりはエネルギーが定数ずれるだけなので、分布としては変わらないので大丈夫というお話です。
……しかし残念ながらチェインは破れるので、これは近似になります。

## チェインの破れ

サンプル値のチェインが破れていた（スピン状態が揃っていない）場合、制限付きボルツマンマシンの埋めこみをどう引き戻そうかという問題があります。
ここで多数決をして多い状態の占める割合が閾値を超えていたら採用するという方法を取ると、閾値を下げれば下げるほど採用できるサンプル数は増えるけれど正確さが落ちることになります。元論文では閾値を下げた方が学習全体の精度は良かったとのことでした。

なお [FixedEmbeddingComposite](https://docs.ocean.dwavesys.com/projects/system/en/latest/reference/composites/fixed_embedding.html) の場合 [sample](https://docs.ocean.dwavesys.com/projects/system/en/latest/reference/composites/generated/dwave.system.composites.FixedEmbeddingComposite.sample.html) メソッドの引数 `chain_break_fraction` を True にするかデフォルトのままにしておくと、レスポンスに壊れたチェインの割合を付与してくれます。

```python
    if chain_break_fraction:
        data['chain_break_fraction'] = broken_chains(record.sample, chain_idxs).mean(axis=1)[idxs]
```

[dwave-system/dwave/embedding/transforms.py](https://github.com/dwavesystems/dwave-system/blob/f1cc9a956b53a9af03d9ac85b44fd6ff26bbd165/dwave/embedding/transforms.py)

# 重み誤差の補正

D-Wave の実機には重み設定に誤差があるらしく、元論文では次の 4 パターンでサンプリングして平均をとることで補正したそうです。ほむ……。

- そのまま
- 可視層の重みと相互作用の重みを反転して、可視層のスピン状態を反転して読む
- 隠れ層の重みと相互作用の重みを反転して、隠れ層のスピン状態を反転して読む
- 可視層の重みと隠れ層の重みを反転して、可視層と隠れ層のスピン状態を反転して読む

## APIによる反転

D-Wave の Solver API には [num_spin_reversal_transforms](https://docs.dwavesys.com/docs/latest/c_solver_1.html#num-spin-reversal-transforms) というパラメータがあって、これはランダムに量子ビットを反転するパターンを指定した数だけ生成してくれるみたいです。
[FixedEmbeddingComposite](https://docs.ocean.dwavesys.com/projects/system/en/latest/reference/composites/fixed_embedding.html) の [sample](https://docs.ocean.dwavesys.com/projects/system/en/latest/reference/composites/generated/dwave.system.composites.FixedEmbeddingComposite.sample.html) を使っている場合は、可変長キーワード引数 `**parameters` 経由で `num_spin_reversal_transforms=2` のように渡せばいいはず。

# 逆温度の推定

ここまで D-Wave によるサンプリングがボルツマン分布に従っているとして話を進めてきたのですが、その逆温度はいくつやねんという問題があります。ボルツマン分布における逆温度とは自然対数の肩に乗ってるやつで、

```math
\begin{align}
E(x) &= - \sum_{ij} J_{ij} x_i x_j - \sum_i h_i x_i \\
B(x) &= \frac{1}{Z(\beta)}e^{- \beta E(x)}
\end{align}
```

の $\beta$ ですね。理論的な後ろ盾がない悲しみで、この値は実験的に推定するしかないみたいです。
元論文では 5x5 8x8 の制限付きボルツマンマシンでは 4.5, 12x12 14x14 16x16 では 3, 32x32 では 2 と推定していました。問題サイズが大きくなるほど逆温度が下がる、すなわちより低いエネルギーのサンプルが得られにくくなるということで、直感に合っていますね。

## APIによる後処理

D-Wave の Solver API には [beta](https://docs.dwavesys.com/docs/latest/c_solver_1.html#beta) というパラメータがあります。お、好きな逆温度を設定できるじゃーんと飛びつきたくなりますが、これは注意が必要で [postprocess](https://docs.dwavesys.com/docs/latest/c_solver_1.html#postprocess) を `sampling` に設定した時に効くパラメータです。つまり実機からのサンプリングには影響せず、古典コンピュータによる後処理で使われます。
後処理のフローは「[Sampling Postprocessing](https://docs.dwavesys.com/docs/latest/c_post-processing_1.html#sampling-postprocessing)」で説明されていて、この `SampleSubgraph` の中身は 1 量子ビットごとの MCMC （マルコフ連鎖モンテカルロ法）みたいですね。ここ D-Wave の API で隠さず、ライブラリ化してほしいなぁ……。
この beta を動した実験結果が [Sampling Tests and Results](https://docs.dwavesys.com/docs/latest/c_post-processing_4.html) に載っていて、各指標が後処理によってどのくらい変化するかグラフ化されています。 D-Wave 実機によるサンプリングが本当にボルツマン分布に従っていて、その逆温度と同じ値を後処理でも設定できていれば、後処理しても同じボルツマン分布からのサンプリングとみなせるので、なるべく指標は変わってないほうが良い設定という見方になります。

# 感想

量子アニーリングは、解きたい問題をイジングモデルに変換してキメラグラフへ埋めこめば OK というイメージでしたが、やはり実機で行うとなると細かな問題が色々とあって、でもそのあたり D-Wave は一つ一つ API やライブラリで対応しようとしていて流石だなぁと感じました。これまで応用分野として組み合せ最適化問題や制約充足問題を期待していましたが、たしかに現状では機械学習と一番相性がいいのかもしれませんね。確率的モデルは多少のノイズを呑んでくれたりするので。
あとは本当にボルツマン分布に従うかどうかですねー。元論文でも believe, conjecture, assume とあって悩ましさが滲みでていますが、そのあたりも研究が進むといいなーと思います。