東芝シミュレーテッド分岐マシン (SBM) による最大カット問題のベンチマーク

はじめに

こんにちは。Kawaharaです。

東芝デジタルソリューションズ株式会社が AWS Marketplace 上で公開しているシミュレーテッド分岐マシン(Simulated Bifurcation Machineｍ, 以下 SBM) を試用し、組合せ最適化問題の一つである最大カット問題 (MAX-CUT) を例題にベンチマークしてみました。

SBM

SBMは東芝デジタルソリューションズが提供している組合せ最適化ソルバーです。分岐現象・断熱過程・エルゴード過程（カオス）といった古典力学の現象を組合せ最適化に応用し、大規模な問題を高速に解くことができます。

またAWS MarketplaceでPoC版を試すことができます。

最大カット問題

SBM で実行できる問題の例として、最大カット問題 (MAX-CUT) を解いてみたいと思います。

最大カット問題とは、与えられた頂点集合において2つの部分集合に分割する際に、異なるグループの間に張られたエッジの重みを最大化する問題です。

目的関数

頂点集合 $V$ から、$S$, $T$ の２つの部分集合に分けることを考えます。

頂点$v_i \in V$が $S$, $T$ のどちらに含まれるかで、対応する2値変数(Ising変数) $s_i$ の値を次のように決めます。

$$s_i = \begin{cases}{} +1 & (v_i\in S) \\\\ -1 & (v_i\in T) \end{cases}$$

この変数を用いて目的関数を以下のように定式化できます。

$${\rm max} \ H = \frac{1}{2} \sum_{1\leq i <j \leq n} W_{ij} \left( 1 - s_i s_j \right)$$

頂点 $v_i, v_j$ について、それぞれが異なるグループに属している場合は $\left( 1 - s_i s_j \right) > 0$、同一のグループの場合は $\left( 1 - s_i s_j \right) =0$ となるので、関数 $H$ の最大化により最適な分割がなされます。

多くのアニーリングマシンでは目的関数の最大値ではなく最小値を探索するため、全体に $-1$ を掛けることで最小化問題とすることが出来ます。

$${\rm min} \ H = - \frac{1}{2} \sum_{1\leq i <j \leq n} W_{ij} \left( 1 - s_i s_j \right)$$

また、目的関数を Ising変数 $s_i \in \{-1, 1\}$ ではなく、バイナリ変数 $q_i \in \{0, 1\}$ で表す場合、下記の変数変換を用いて書き換えることが出来ます。

$$q_i = \frac{s_i+1}{2}$$

$$\begin{aligned} {\rm min} \ H &= - \frac{1}{2} \sum_{1\leq i <j \leq n} W_{ij} \left(1- (2q_i-1) (2q_j-1) \right) \\\\ &= \sum_{1\leq i <j \leq n} W_{ij} \left( 2q_iq_j -q_i - q_j \right) \end{aligned}$$

バイナリ変数 $q_i \in \{0, 1\}$ を受け取るアニーリングマシンの場合、こちらの式を用いて実行します。

SBM の実行

MAX-CUT モード

SBM には 最大カット問題専用のソルバー (MAX-CUTモード) が用意されています。このソルバーは問題を表現するテキストデータを入力することが出来ます。データは MatrixMarket フォーマットか 、最大カット問題のベンチマーク用データセットである Gset フォーマットに対応しています。

Gsetフォーマットは、以下のような空白区切りの整数列で与えられます。

    N E
    i0 j0 w_i0j0
    i1 j1 w_i1j1
    i2 j2 w_i2j2

1行目N、Eはそれぞれ問題グラフの頂点数$N$と辺の数$E$、2行目以降は頂点間$i,j$の辺の重み$w_{ij}$を表しています。

例として、下記に5頂点の最大カット問題を Gset 形式で表しました。

    5 7
    1 2 1
    1 5 1
    2 3 -1
    2 4 1
    2 5 -1
    3 4 -1
    4 5 1

これをグラフ化^[1]すると下記のようになります。

SBMのパラメータ

SBM には下記のポストパラメータを与えることが出来ます。

• steps: SBアルゴリズムの実行ステップ数。これが runs（ソルバー固有の並列数）=320の個数、並列実行されます。この実行ステップ全体を1ループと呼びます。steps を 0 に設定した場合、実行ステップ数が自動設定されます。
• loops: 上記の１ループ終了までの処理を繰り返す数。0に設定した場合、timeout時間まで繰り返します。
• timeout: 実行時間上限秒数。実行時間がこの設定値まで達すると、その時点で実行中のループは破棄され、timeout 到達時点で得られた最良値が解として返されます。
• target(オプション): 終了条件を指定することができます。マシンに与えた目的関数の値がtargetで指定した値に到達すると計算を終了します。
• stats: マシンの出力値の他に以下のような統計的情報を返すかどうか指定することができます。これは runs(並列数) * loops 個の解についての統計です。
  • none: 何も返さない
  • summary: 平均値、標準偏差、最良解の頻度
  • full: 平均値、標準偏差、全解のヒストグラム

先ほどの5頂点の問題を curl コマンドで実行してみました。接続先の 192.168.0.1 は環境に応じて書き換えてください。

 $ echo -e "5 7\n1 2 1\n1 5 1\n2 3 -1\n2 4 1\n2 5 -1\n3 4 -1\n4 5 1" | curl -i -H "Content-Type: application/octet-stream" -X POST "http://192.168.0.1:8000/solver/maxcut?steps=0&loops=0&timeout=1&stats=full" --data-binary @-

今回は下記のパラメータを使用しています。

• steps = 0
• loops = 0
• timeout = 1
• stats = full

正常に実行されると、次のようなレスポンスが得られます。

    HTTP/1.1 200 OK
    Content-Type: application/json; charset=utf-8
    Content-Length: 233
    ETag: W/"e9-6jYjTGcvinFBL6/eKnutYNzItOQ"
    Date: Mon, 19 Aug 2019 03:59:29 GMT
    Connection: keep-alive

    {"id":"r1477108799","time":1,"wait":0,"runs":396800,"steps":10,"message":"timeout","value":3,"result":[0,1,0,0,1],"stats":{"avg":2.441142,"stddev":1.372172,"histogram":[[-2,24009],[-1,4124],[0,15682],[1,12766],[2,12636],[3,327583]]}}

レスポンスの json の代表的なキーは下記のようになっています。

• result: 目的関数を最小化するQUBO変数の組み合わせ
• value: result 変数における目的関数の値
• runs: (実行されたloopsの数)*(ソルバー固有の並列数)

Gsetによるベンチマーク

次に Gset と呼ばれる最大カット問題のベンチマークセットを実行してみます。Gsetのデータセットは以下からダウンロード出来ます。

http://web.stanford.edu/~yyye/yyye/Gset/

Gsetには、G1 - G81 (飛び番あり) の全71個の問題があり、各問題の頂点数は 800 - 20000です。エッジ重み無しグラフと エッジの重み {+1, -1} の重み付きグラフの２種類があり、幾何学的な構造としては以下の３種類に分類できます。

• ランダムグラフ
• 平面グラフ
• トロイダル格子グラフ

実行コマンド

簡単のために curl コマンドにて SBM に Gset をPOSTしてみます。Gset の問題データが配置されたディレクトリに移動し、例えば問題 G1 について、steps=0(自動), loops=0(自動), timeout=10, target=11624 という設定で実行するには、以下のコマンドを入力します。

$ cat G1 | curl -i -H "Content-Type: application/octet-stream" -X POST "http://192.168.0.1:8000/solver/maxcut?steps=0&loops=0&timeout=10&target=11624" --data-binary @-

ベンチマーク結果

以下は G1 - G72 の全69問のデータに対して、steps=0(自動), loops=0(自動), timeout=10, target=(best-known) で実行した結果です。target には先行研究で知られている最良値を入力しています。

instance	node	edge	weight	type	target	time[sec]	result	diff
G1	800	19176	+1	random	11624	0.02	11624	0
G2	800	19176	+1	random	11620	0.13	11620	0
G3	800	19176	+1	random	11622	0.09	11622	0
G4	800	19176	+1	random	11646	0.02	11646	0
G5	800	19176	+1	random	11631	0.58	11631	0
G6	800	19176	+1,-1	random	2178	0.02	2178	0
G7	800	19176	+1,-1	random	2006	0.02	2006	0
G8	800	19176	+1,-1	random	2005	0.02	2005	0
G9	800	19176	+1,-1	random	2054	0.02	2054	0
G10	800	19176	+1,-1	random	2000	0.06	2000	0
G11	800	1600	+1,-1	toroidal	564	0.01	564	0
G12	800	1600	+1,-1	toroidal	556	0.02	556	0
G13	800	1600	+1,-1	toroidal	582	0.02	582	0
G14	800	4694	+1	planar	3064	10.00	3063	1
G15	800	4661	+1	planar	3050	0.06	3050	0
G16	800	4672	+1	planar	3052	0.18	3052	0
G17	800	4667	+1	planar	3047	4.64	3047	0
G18	800	4694	+1,-1	planar	992	0.16	992	0
G19	800	4661	+1,-1	planar	906	0.02	906	0
G20	800	4672	+1,-1	planar	941	0.01	941	0
G21	800	4667	+1,-1	planar	931	0.08	931	0
G22	2000	19990	+1	random	13359	0.20	13359	0
G23	2000	19990	+1	random	13344	9.99	13342	2
G24	2000	19990	+1	random	13337	0.07	13337	0
G25	2000	19990	+1	random	13340	0.45	13340	0
G26	2000	19990	+1	random	13328	0.07	13328	0
G27	2000	19990	+1,-1	random	3341	0.07	3341	0
G28	2000	19990	+1,-1	random	3298	0.20	3298	0
G29	2000	19990	+1,-1	random	3405	0.07	3405	0
G30	2000	19990	+1,-1	random	3413	0.20	3413	0
G31	2000	19990	+1,-1	random	3310	0.79	3310	0
G32	2000	4000	+1,-1	toroidal	1410	2.08	1410	0
G33	2000	4000	+1,-1	toroidal	1382	9.77	1382	0
G34	2000	4000	+1,-1	toroidal	1384	0.78	1384	0
G35	2000	11778	+1	planar	7687	10.00	7680	7
G36	2000	11766	+1	planar	7680	10.01	7675	5
G37	2000	11785	+1	planar	7691	10.00	7685	6
G38	2000	11779	+1	planar	7688	10.00	7686	2
G39	2000	11778	+1,-1	planar	2408	10.00	2407	1
G40	2000	11766	+1,-1	planar	2400	4.95	2400	0
G41	2000	11785	+1,-1	planar	2405	10.00	2404	1
G42	2000	11779	+1,-1	planar	2481	9.87	2475	6
G43	1000	9990	+1	random	6660	0.02	6660	0
G44	1000	9990	+1	random	6650	0.02	6650	0
G45	1000	9990	+1	random	6654	0.02	6654	0
G46	1000	9990	+1	random	6649	0.02	6649	0
G47	1000	9990	+1	random	6657	0.02	6657	0
G48	3000	6000	+1,-1	toroidal	6000	0.07	6000	0
G49	3000	6000	+1,-1	toroidal	6000	0.07	6000	0
G50	3000	6000	+1,-1	toroidal	5880	0.07	5880	0
G51	1000	5909	+1	planar	3848	0.94	3848	0
G52	1000	5916	+1	planar	3851	1.77	3851	0
G53	1000	5914	+1	planar	3850	10.00	3849	1
G54	1000	5916	+1	planar	3852	10.00	3851	1
G55	5000	12498	+1	random	10299	10.00	10289	10
G56	5000	12498	+1,-1	random	4017	9.80	4008	9
G57	5000	10000	+1,-1	toroidal	3494	10.00	3480	14
G58	5000	29570	+1	planar	19293	10.01	19257	36
G59	5000	29570	+1,-1	planar	6086	10.00	6067	19
G60	7000	17148	+1	random	14188	10.00	14168	20
G61	7000	17148	+1,-1	random	5796	10.01	5777	19
G62	7000	14000	+1,-1	toroidal	4870	10.01	4844	26
G63	7000	41459	+1	planar	27045	10.01	26986	59
G64	7000	41459	+1,-1	planar	8751	8.82	8728	23
G65	8000	16000	+1,-1	toroidal	5562	10.01	5532	30
G66	9000	18000	+1,-1	toroidal	6364	9.01	6324	40
G67	10000	20000	+1,-1	toroidal	6950	10.01	6906	44
G70	10000	9999	+1	random	9591	10	9522	69
G72	10000	20000	+1,-1	toroidal	7006	10	6966	40

数千程度の小さな問題については、概ね1秒以下とかなり高速にターゲットとした最適解の探索に成功していることが分かります。

上記の結果では steps を自動設定としましたが、より性能を引き出すためには steps の調整が有効と思われます。簡単に 60秒の実行時間制約においてどの程度の性能を引き出せるか、G65 を対象に簡単に調査してみました。steps は 10,000～1,000,000 の間で 変化させています。

図中の各点に付記した数字は、そのstepsのもとで実行時間60秒とした時のループ回数です。

G65では500,000以上のような大きなstepsを設定するほうが、小さなstepsのループを 多く実行するよりも効率がいいという傾向がみられました。

このように、問題によって短いループを多数回まわす方が良いか、あるいはある程度長い1ループを実行するほうが良いのか異なる場合があるため、効率よく解を得るためには最適なsteps設定を見つけることが鍵になるのかもしれません。

まとめ

SBM は最大カット問題に対し高速に近似解を得られることが体験できました。特に頂点数$N=2000$ 程度の比較的小さい問題では概ね1秒以下で最適解に到達していました。

最大カット問題は、イジング定式化された組合せ最適化問題のなかでは、制約項のない比較的解きやすい問題だと考えられています。今回はベンチマークとして最大カット問題と Gset を例に取り上げましたが、より複雑かつ実用的な組合せ最適化問題についても期待が持てるのではと思われます。また、チューニング次第でどこまで高速あるいは高精度な解が出せるのか、継続して調査していく予定です。

脚注

1. https://csacademy.com/app/graph_editor/