Fixstars Amplifyを使って最大カット問題を解いてみる

:::note info
この記事は[東大推薦Advent Calendar](https://adventar.org/calendars/8917)1日目の記事です。
:::

# はじめに
こんにちは、東京大学工学部推薦7期のGraphium(X: [@graphium142857](https://x.com/graphium142857))です。現在学部2年生で、工学部計数工学科に内定しています。
この記事では、私が前から気になっていたFixstars Amplifyを使ってみます。

# Fixstars Amplifyとは
[公式ホームページ](https://amplify.fixstars.com/ja/)によると、

> Fixstars Amplifyは、イジングマシン向けSDKと実行環境からなるクラウド基盤です。

このサービスとそれを提供する会社の名前がともにFixstars Amplifyです。Amplify SDKを使うことで、様々な最適化問題を定式化し、組合せ最適化問題に特化したイジングマシンで解くことができます。これを使うことで、Fixstars Amplify AEというアニーリングマシンや、量子コンピュータなどで、無料で簡単に最適化問題を解けます。

# Amplify SDKの使い方
## 扱える問題
以下のような形式で表せる目的関数の最小化をできます。
```math
\sum_{i<j}J_{i j}s_i s_j + \sum_{i}J_{i}s_i
```
$s_i$は二値変数で、$s_i\in\lbrace -1,+1\rbrace$のときイジングモデル、$s_i\in\lbrace 0,+1\rbrace$のときQUBOモデルと呼ばれます。

## 手元で使ってみる
### インストール
インストールは簡単で、
```bash
$ pip install amplify
```
をするだけです。

### アクセストークンの取得
公式サイトから簡単にアクセストークンを取得できます。これを使うことで、Fixstars Amplify AE等のマシンを無料で使用できます。以下のリンクから登録してください。
https://amplify.fixstars.com/ja/user/token

量子コンピュータを使ってみたい場合は、別にアクセストークンの取得が必要になります。
書類を提出することでD-Wave用のアクセストークンをFixstars Amplify経由で取得できますが、[公式サイト](https://www.dwavesys.com/build/getting-started/)でアカウント登録することでもすぐに取得できます。

### 書き方
[公式チュートリアル](https://amplify.fixstars.com/ja/docs/tutorial.html)のサンプルコードが非常にわかりやすいです。また、[公式のオンラインデモ & チュートリアル](https://amplify.fixstars.com/ja/demo)には、チュートリアルや、様々な最適化問題に対するデモコードがあり、オンラインでコードを見て環境構築なしに実行することができます。

# 最大カット問題
## 最大カット問題とは
重みつき無向グラフ$G=(V,E)$の頂点を2つの集合$V_1,V_2(V_1\oplus V_2 = v)$に分割することを考えます。各辺の重みを、$w_{i j}((i,j)\in E)$としたとき、2つの集合にまたがる辺の重みの和をカット容量と言い、これは下の式で表せます。
```math
\sum_{i\in V_1,j\in V_2}w_{ij}
```
これを最大化する問題を最大カット問題と呼びます。
カット容量を最小化する最小カット問題は多項式時間アルゴリズムが知られていますが、最大カット問題はNP困難であることが知られており、簡単に言うと高速なアルゴリズムが知られていません。

## イジングモデルでの定式化
最大カット問題をイジングマシンで解くために、イジングマシンでの定式化をします。
頂点$i$がカット$(V_1,V_2)$のどちらに属するかを、イジング変数$s_i \in \lbrace-1,+1\rbrace$を用いて、以下のように表現します。
```math
s_i = \left\{
\begin{array}{ll}
+1 & (i \in V_1) \\
-1 & (i \in V_2)
\end{array}
\right.
```
ここで、$\frac{1}{2}(1-s_i s_j)$は、$i,j$がそれぞれ$V_1,V_2$の異なる頂点集合に属しているときに1,同じ頂点集合に属しているときは-1をとるので、カット容量は以下のように表せます。
```math
\frac{1}{2}\sum_{(i,j)\in E} w_{i j}(1-s_i s_j)
```
今回はこれを最大化しますが、実装の際には-1をかけて最小化問題として解きます。

# 実際に解いてみた
## 実装
Fixstars Amplify AEで解くための実装は以下のようになります。
```python
import sys
from amplify import IsingPoly, IsingSymbolGenerator
from amplify import Solver
from amplify.client import FixstarsClient
from dotenv import load_dotenv
load_dotenv()

import os

def main():
    args = sys.argv
    input_path = args[1]

    with open(input_path, 'r') as f:
        lines = f.readlines()

    # v: 頂点数, e: 辺数
    v,e = map(int, lines[0].split())

    # イジング変数の設定
    gen = IsingSymbolGenerator()
    s = gen.array(v)

    # 目的関数の設定
    f = IsingPoly()
    for i in range(1,e+1):
        u,v,c= map(int, lines[i].split())
        f += c*(1-s[u-1]*s[v-1])
    f /= 2
    f *= -1

    # clientの設定
    client = FixstarsClient()
    client.parameters.timeout = 10000 # タイムアウト10秒
    client.token = os.getenv('FIXSTARS_AE_TOKEN')

    # 実行
    solver = Solver(client)
    result = solver.solve(f)

    for sol in result:
        solution = s.decode(sol.values)

    print(-sol.energy)

if __name__ == '__main__':
    main()
```
上のコードを使用する際は、.envという名前で以下のような内容のファイルを作り、自分で取得したAPIキーを書いてくさい。
```
FIXSTARS_AE_TOKEN = "xxxxxxxxxxxxxxxxxxxxxxx"
```
今回はD-WaveのLeap Hybrid Solverでの実行も行ってみました。実装が気になる方は、[githubのリポジトリ](https://github.com/Graphium314/amplify_playground)を見てください。

## テスト
今回は[Gset](https://web.stanford.edu/~yyye/yyye/Gset/)というデータセットを使いました。こちらは参考にしたブログ等[^1][^2]で使われていたので使用しました。これらのうち、一部のデータを適当に選んで使っています。

## 結果
Fixstars Amplify AEとD-Wave Leap Hybrid Solverでの実行結果を以下の表にまとめました。best knownはFixstarsのブログ[^1]にあったものです。
| problem    | best known | Fixstars Amplify AE | D-Wave Leap Hybrid Solver |
|:-----------|------------:|------------:|------------:|
| G1         | 11624       | 11624       | 11624       |
| G14        | 3064        | 3063        | 3061        |
| G35        | 7687        | 7674        | 7676        |
| G53        | 3850        | 3848        | 3848        |
| G72        | 7006        | 6924        | 6938        |

# まとめ
前々から気になってはいましたが、公式のチュートリアルやデモが充実していて、意外と簡単に使えてよかったです。また、計算資源を無料で使えるのも非常に良いですね。今回は定式化しやすい最大カット問題に取り組みましたが、他の問題も試してみたいと思います。

[^1]: [東芝シミュレーテッド分岐マシン (SBM) による最大カット問題のベンチマーク](https://proc-cpuinfo.fixstars.com/2019/09/sbm/)
[^2]: [最大カットのベンチマークを最適化ソルバーで解く](https://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1941-09.pdf)