シミュレーテッド分岐マシン(SBM)で巡回セールスマン問題を解く

はじめに

こんにちは。Kawaharaです。

前回の記事では、最大カット問題をSBMで実行してみました。

今回は、組み合わせ最適化問題の具体例として取り上げられることの多い「巡回セールスマン問題(TSP)」を、SBMで解いてみました。

巡回セールスマン問題(TSP)

巡回セールスマン問題とは、セールスマンが都市を訪問する際、全ての都市に1度ずつ訪問する最短経路を求める問題です。

都市の数を$N$とすると解の候補となる経路の数は$N!$通りあり、 $N$ が増えると組み合わせの数は爆発的に増えてしまうため、全探索で最適解を求めるのは不可能になってしまいます。

この問題をQUBO形式で表現し、SBMで解いてみようと思います。

目的関数

TSPのQUBO表現での目的関数を作成します。詳しい説明はこちらのページを参照してください。

今回用いる方法では、「何番目に」「どの都市に」訪問するかを1つのQUBO変数とします。

例えば $N=4$の場合において、「B→A→C→D」の経路となるQUBO変数の組み合わせを考えます。以下の表は、このQUBO変数の組で「何番目に訪問するか」を行に、「どの都市に訪問するか」を列に対応するように並べたものです。

	A	B	C	D
1番目	0	1	0	0
2番目	1	0	0	0
3番目	0	0	1	0
4番目	0	0	0	1

都市数$N$に対しては $N^2$個のQUBO変数が必要になります。

都市$(i,j)$ 間の距離を$d_{i,j}$として、「$\alpha$番目に都市$i$に訪問するかどうか」を表すQUBO変数$n_{\alpha, i}$を用いると、目的関数は次式のようになります。

$$H = \frac{1}{2} \sum_{\alpha, i, j} d_{i,j} n_{\alpha, i} (n_{\alpha+1, j} + n_{\alpha-1, j}) + k_1 \sum_{i} (\sum_{\alpha} n_{\alpha, i} -1)^2 + k_2 \sum_{\alpha} (\sum_{i} n_{\alpha, i} -1)^2$$

第1項が最小化したい総距離を表す項、第2項と第3項が制約条件の項で、$k_1$, $k_2$は制約条件の強さを表す重み係数です。

通常のシミュレーテッド・アニーリング（SA）において、得られる解が制約を満たせるのならば、制約条件の項にかかる重みはできるだけ小さいほうが良いと思われます。なぜなら、ある局所解から別の局所解に移動する際、その移りやすさを高められると考えられるからです。

本記事では重みの値として$\frac{1}{2} \max(d_{i,j})$を用います。

$$k_1 = k_2 = \frac{1}{2} \max(d_{i,j})$$

SBMのアルゴリズムはSAとは本質的に異なりますが、同じ定式化を用いているためこのような工夫が効いてくる可能性があります。実際に試したところ、重みの調整が解の品質に効いてくることを確認しました。

TSPLIB

TSPLIBは巡回セールスマン問題のベンチマークセットの1つです。問題データはリンク先からダウンロードできます。

今回は以下の2問題について評価を行いました。それぞれ14都市、24都市の巡回セールスマン問題です。

• burma14
• gr24

SBMでの実行

ISINGソルバー

今回TSPを解くにあたっては、SBMのISINGソルバーを用います。ISINGソルバーは一般的なQUBOモデルの問題を解くソルバーです。ISINGソルバーで与える問題データは、QUBO多項式の係数を以下のMatrixMarketフォーマットと呼ばれる形式に従うテキストデータで与えます。

%%MatrixMarket matrix coordinate real symmetric
N N E
i0 j0 w_i0j0
i1 j1 w_i1j1
i2 j2 w_i2j2

1行目は与える行列が実対称行列であることを宣言します。

2行目は行列の次元$N$と要素の個数$E$を記述します。言い換えると、イジングマシンでのQUBO変数とその相互作用係数(非ゼロ行列要素)の個数にあたります。

3行目以降は、行列要素のインデックスと値です。対象行列のうち下三角要素、つまり$i \geq j$の要素のみを記述するよう指定されています。

一例として、ランダムに生成した6都市問題のQUBO多項式をMatrixMarket化したデータはこのようになります。

'%%MatrixMarket matrix coordinate real symmetric
36 36 396
1 1 -1.526
2 1 1.526
2 2 -1.526
3 1 1.526
3 2 1.526
3 3 -1.526
4 1 1.526
4 2 1.526
4 3 1.526
4 4 -1.526
5 1 1.526
5 2 1.526
(以下略)

実行コマンド

上記のフォーマットで作成した問題データをcurlコマンドで送信し、計算結果を取得してみます。接続先の 192.168.0.1 は環境に応じて書き換えてください。

curl -i -H "Content-Type: application/octet-stream" -X POST "http://192.168.0.1:8000/solver/ising?steps=500&loops=0&timeout=10&dt=0.2" --data-binary @problem.data

SBMのパラメータ(新パラメータ補足)

SBM PoC版の新バージョンがリリースされ、新たに以下のマシンパラメータが設定できるようになりました。

• dt：微分方程式を symplectic Euler methodで解くときの時間幅。
• SBM論文における(12)式、(13)式の$\Delta t$。
• C：SBMハミルトニアンを構成する定数。
• SBM論文における(2)式の$\xi_0$。

dtのデフォルト値は1.0、Cはデフォルトでは自動調整となっているのですが、問題によってはこれを適切に設定することで、より良い解が得られやすくなることがあります。今回は、このdt、Cの最適な設定値を探索してみました。

またstepsについても数種類の設定値を試してみました。

実行結果

パラメータdtの調整

dtを0.01 ~ 1.0の範囲で変化させて実行し、得られたQUBO変数の値からTSP経路の距離を計算しました。

その他のマシンパラメータの設定値は次のとおりです。

• steps: 100, 1000, 10000
• loops: 0 (自動設定)
• target: 既知の最適値
• timeout: 10 (burma14), 20 (gr24)
• C: 0 (自動設定)

burma14

上の図は横軸dt、縦軸を距離としてプロットしたグラフです。
縦軸の下端は最適値(3323)となっています。

14都市の問題burma14では、実行した範囲では概ね最適解に到達することができました。特に0.1<dt<0.3の領域において最適解が得られやすいことを示しています。

gr24

gr24 の問題データは距離行列のみで都市の位置が与えられていないため、経路図の出力はできませんでした。

gr24では、各dtに対して10回実行してその平均値と標準偏差を調べてみました。
図はTSPの解の距離の平均値をプロットしたもので、標準偏差をエラーバーとしています。
こちらも、縦軸の下端値は最適値(1272)としています。

dt>0.5の領域では、TSPの実行可能解は得られませんでした。

0.1<dt<0.3の領域でstepsごとの結果と比較すると、steps=10000で最適解からの差分が小さくなっています。
stepsについては、ある程度大きい値を設定したほうがよい解が得られやすい傾向が見られました。

実行時間の上限を固定したとき、stepsとrunsはトレードオフの関係にあるため、最適なstepsを設定することで精度のよい解を得られやすくなる可能性があります。

この実験では、10回の平均値が最小となるパラメータの組は(steps=10000, dt=0.26)でした。

パラメータCの調整

新バージョンで追加されたもう一つのパラメータであるCについても、様々な設定値を試してみました。

まず、先のdtの実験でのCの自動設定値を確認してみます。

図は各dtに対し10回分の実行結果からCの自動設定値の平均を計算しプロットしたものです。
dt=0.3以下までは0.00175付近、それ以降では0.0020~0.0022の値がよく設定されていたことがわかりました。

そこで0.00025~0.005の範囲で0.00025刻みで値をとり、先と同じく各10回実行してみました。
C以外のパラメータの設定は、前節で平均的に良い解が得られていた設定値(steps=10000, dt=0.26)にしています。

自動調整での設定値より少し小さい値であるC=0.00150の結果が平均値・最小値ともに良く、最適解にも到達しています。
また、C=0.003以上では実行可能解は得られませんでした。

このようにCの調整も解の向上に寄与する可能性があることがわかりました。

まとめ

今回実験したTSPの問題と、得られた最良値の結果を示します。

問題	最適値	steps	dt	C	実行時間[sec]	runs	平均値	最良値	最適値差分
burma14	3323	100	0.13	auto	0.01	960		3323	0.00%
		1000	0.42	auto	0.01	320		3323	0.00%
		10000	0.26	auto	0.07	320		3323	0.00%
gr24	1272	100	0.34	auto(0.00215)	20.0	1360640	1525.6	1346	5.82%
		1000	0.44	auto(0.00224)	20.0	217280	1400.5	1314	3.30%
		10000	0.40	auto(0.00223)	20.0	22720	1362.1	1290	1.42%
		10000	0.26	auto(0.00173)	20.0	22720	1346.2	1317	3.54%
		10000	0.26	0.00150	14.53	16960	1317.1	1272	0.00%

steps, dt, Cの各種マシンパラメータを調整することで、TSPの最適解に近い解を得られることがわかりました。また、今回扱った問題においては0.2<dt<0.4の範囲で調整することで比較的良い解が得られやすい傾向が見られました。

参考

• TSPのQUBO定式化について
• TSPLIB
• SBM論文