東芝の最適化マシンSBMを使ってみる

# はじめに

2019年4月、東芝が**シミュレーテッド分岐マシン**（以下、SBM）という最適化マシンを発表しました[^sbmpr][^sbmpaper][^sbmintro]。
このSBMの概念実証（PoC）版が、2019年10月末までの期間限定でAWS Marketplace上で公開されています[^sbmawspr][^sbmaws]（インスタンス代を除いて無料）。

AWS上から手軽に試せるということに惹かれて、実際にSBMのPoC版を利用してみました。

# SBMについて

SBMは、**QUBO**という形式で表現された最適化問題の近似解を求めることができる最適化マシンです。

QUBOとは、以下の形式で表現された関数を最適化する問題です:

```math
E = \sum_{i \ge j} Q_{ij} x_i x_j
```

各$x_i$は、0もしくは1をとる変数です。各$Q_{ij}$は、各$x_i$間の結合の強さを表します。
QUBOはイジングモデルと呼ばれるモデルで表現できることから、QUBOの最適化を行うマシンを**イジングマシン**と呼称することもあります。


イジングマシンとして有名なのは、D-Waveが実用化した**量子アニーリング**です。
一方で、SBMは量子技術を用いたものではなく、通常の（非量子的な）マシン上で動作します。紹介サイト[^sbmintro]では「量子分岐マシン研究から生まれた」と記されており、量子アニーリングの影響を受けて開発されたようです。

同様に、非量子的なアーキテクチャで動作する最適化マシンとして、富士通の**デジタルアニーラ**[^da]、日立の**CMOSアニーリング**[^cmos]や**モメンタム・アニーリング**[^ma]などがあります。

東芝はSBMの特徴として、従来アルゴリズムと異なり並列化による高速化が容易なことを挙げています[^sbmpr]。実際、GPUによる並列化を行い、10万変数・全結合の問題に対して数秒で良解を求められたとしています。

SBMの原理などについては、論文[^sbmpaper]を参照してください。

# SBM PoC版について

SBMの概念実証（PoC）版のAMIが、Amazon Marketplaceで2019年10月末までの期間限定で無償公開されています。

https://aws.amazon.com/marketplace/pp/B07T8STKZR/

AMI自体は無償ですが、インスタンス代はかかるので注意してください。利用にはp3.2xlargeのインスタンスが必要です。
ユーザマニュアルは上記ページに記載のリンクからダウンロード可能です。

SBMは、作成したEC2インスタンス上で動作するWeb APIとして提供されます。利用する際には作成したインスタンスにHTTPリクエストを投げる形になります。
SBMの設定パラメータなどを変更したい場合には、SSHでログインして設定ファイルを書き換える必要があります。

扱えるQUBOのサイズについては、以下の通りです（ユーザマニュアルに記載）:

* 変数の数: 最大10,000
* 結合行列の非零要素の数: 1,000,000

# SBM PoC版の基本的な使い方

SBMのPoC版には、以下の三つのAPIが用意されています。

* ISING
* Max-CUT
* Max-SAT

ここではISINGのAPIに最適化クエリを投げてみたいと思います。

例として、以下の関数を最適化する問題を考えます。

```math
E = -3 x_1^2 + 4 x_1 x_2 -4 x_2^2
```

行列で書くと、以下のようになります:

```math
E = \sum_{i \ge j} Q_{ij} x_i x_j \\
Q_{11} = -3 \\
Q_{12} = Q_{21} = 4 \\
Q_{22} = -4
```

SBMにリクエストする際には、この行列$Q$を、MatrixMarket形式[^matrixmarket]と呼ばれる形で表現する必要があります。
これは、行列の非ゼロ要素のみを列挙した以下のような形です:

```
%%MatrixMarket matrix coordinate real symmetric
[行数] [列数] [非ゼロ要素数]
[行番号1] [列番号1] [値1]
[行番号2] [列番号2] [値2]
...
```

先頭行は固定です。また、列挙する要素は、行列の下三角部分のみになります（つまり、$行番号\ge列番号$の部分）。

今回の場合、以下のようになります:

```
%%MatrixMarket matrix coordinate real symmetric
2 2 3
1 1 -3
2 1 4
2 2 -4
```

これをリクエストボディとして、以下のURLにPOSTメソッドでHTTPリクエストを行います:

```
http://[IPアドレス]:8000/solver/ising?steps=[ステップ数]&loops=[ループ数]
```

`steps`や`loops`は、SBMの動作パラメータです。`steps`は一回の実行で何ステップ回すかを表し、`loops`は全体の実行回数に関係します。`steps=0`とするとステップ数が自動調整されます。ユーザマニュアルにこれらの関係を示した図があるので、詳しくはそちらを参照してください。

`curl`を使う場合は、以下のようになります。

```
$ curl -H "Content-Type: application/octet-stream" -X POST  "http://[IPアドレス]:8000/solver/ising?steps=0&loops=100" --data-binary @test.data
```

`--data-binary`はcurlのオプションで、HTTPリクエストのリクエストボディを指定することができます。`@[ファイル名]`の形で指定することで、指定されたファイルの中身をリクエストボディに渡すことができます（ここでは、上記の行列表現が`test.data`というファイルに保存されているとします）。詳細は、curlのヘルプを参照してください。

上記のリクエストを投げると、以下のようなjson形式のレスポンスが返ってきます:

```
{"id":"r3990952542","time":0.19,"wait":0,"runs":32000,"steps":4,"message":"finished","value":-4,"result":[0,1]}
```

`result`がSBMが求めた解の変数の値で、`value`がそのときの目的関数の値です。
今回の場合、$x_1=0, x_2=1$のときに$E=-4$というのがSBMが求めた解になります。問題サイズが小さいこともあり、目的関数の厳密な最適解になっています。

### stepsについての注意
`steps=0`とするとステップ数が自動調整されますが、場合によっては自動調整では最適な精度が得られないこともあるようです。
この場合は、`steps`の値を手動で設定すると精度を調整できます。十分な精度が得られない場合には試してみてください。


# SBMで巡回セールスマン問題を解いてみる

SBMを用いて、巡回セールスマン問題を解いてみます。
巡回セールスマン問題とは、複数の地点をちょうど一度ずつ訪れる経路のうち、最も総距離の短い経路を求める問題です。

まず、巡回セールスマン問題をQUBOとして表現する必要があります。
ここでは簡単な説明だけ行います。既存の文献が多数ありますので、詳細はそちらを参照してください[^fixstars][^isingtsp][^isingnp][^studyqa]。

地点の数を$N$として、以下のような$N^2$個の変数を用意します:

```math
x_{u,i} = \begin{cases} 1 \quad\text{(地点$u$を$i$番目に訪れる)} \\ 0 \quad\text{(それ以外)} \end{cases}
```

これらの変数を使った以下の関数を最適化することで、巡回セールスマン問題を解くことができます:

```math
E = \sum_{u,v} \sum_{i} D_{uv} x_{u,i}x_{v,i+1} + B \sum_{i} \left( \sum_{u} x_{u,i} - 1 \right) ^2 + B \sum_{u} \left( \sum_{i} x_{u, i} - 1 \right)^2
```

一つ目の項は、経路の総距離を表しています。$D_{uv}$は地点間の距離です。
また、後ろの二項は、「同じ地点に二回訪れる」などの意味のない解にペナルティをかけています。$B$はペナルティの大きさを表します。

ペナルティの大きさ$B$、またSBMの動作パラメータである`steps`などはチューニングする必要があります。

## 結果

ランダム生成した巡回セールスマン問題をSBMに解かせてみました。

試行錯誤した結果、ペナルティを大体以下の値に設定したときに、SBMが良い解を出しやすかったです:

```math
B = 1.2 \cdot \max_{u, v} D_{uv}
```

また、SBMの動作パラメータについては、`steps=20`に設定し、`loops`は地点数ごとに調整しました。

以下、いくつかの地点数でSBMが出した経路です。

| 地点数 | 実行時間 | 経路 |
|:-----:|:--------:|:--:|
| 8 | 0.2秒 | ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/315364/246c4b42-cc75-58ee-8c4b-96f77ce0d2b4.png) |
| 10 | 20秒 | ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/315364/25cae21a-111b-987d-647a-092b8b652f80.png) |
| 13 | 5分 | ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/315364/80a5df92-7f80-39ad-8e72-4bc463ec75f4.png) |

このパラメータでは13地点でかなり多めの時間がかかってしまい、また出てきた経路もあまり最適ではありませんでした。
SBMで多地点の巡回セールス問題を解くためにはもう少し調査・研究が必要だと感じました。

# まとめ

10月末までの期間限定ではありますが、AWS上から手軽に試せるのは非常に魅力的だと感じています。11月以降どうなるかは分かりませんが、イジングマシンに興味がある方はぜひ試してみてください。

# 関連記事

D-Wave Leapで量子アニーリングマシンの実機に触れる
https://qiita.com/navitime_tech/items/a9ec5784901c7f8bfb57


[^sbmpr]: SBMについてのプレスリリース https://www.toshiba.co.jp/rdc/detail/1904_01.htm
[^sbmpaper]: SBMの論文 https://advances.sciencemag.org/content/5/4/eaav2372
[^sbmintro]: SBM紹介サイト https://www.toshiba-sol.co.jp/pro/sbm/index_j.htm
[^sbmawspr]: SBMのAmazon Marketplace公開についてのプレスリリース https://www.toshiba-sol.co.jp/news/detail/20190717.htm
[^sbmaws]: SBMのAmazon Marketplaceページ https://aws.amazon.com/marketplace/pp/B07T8STKZR/
[^da]: デジタルアニーラ紹介サイト https://www.fujitsu.com/jp/digitalannealer/
[^cmos]: CMOSアニーリング プレスリリース http://www.hitachi.co.jp/New/cnews/month/2015/02/0223b.html
[^ma]: モメンタム・アニーリング プレスリリース https://www.hitachi.co.jp/rd/news/topics/2019/0830.html
[^matrixmarket]: MatrixMarket format https://math.nist.gov/MatrixMarket/formats.html
[^fixstars]: 巡回セールスマン問題 – 量子コンピュータ情報サイト https://quantum.fixstars.com/introduction_to_quantum_computer/applications/tsp/
[^isingtsp]: イジング+4spin-flipの巡回セールスマンをアニーリングで解くには。 https://qiita.com/YuichiroMinato/items/fd2937cd3e593fa35412
[^isingnp]: "Ising formulations of many NP problems" https://www.frontiersin.org/articles/10.3389/fphy.2014.00005/full
[^studyqa]: "Study of Optimization Problems by Quantum Annealing" https://arxiv.org/abs/quant-ph/0205020