はじめに
これはeeic(東京大学工学部電気電子・電子情報工学科) Advent Calendar 2020 18日目の記事です。
アドベントカレンダー、とりあえず登録はして何書こうか迷ったのですが、今年一番時間をかけた問題はなんだろうって考えたときに巡回セールスマン問題が頭に浮かんだのでそれについて書きます。
私は6〜8月にGoogleとメルカリのエンジニア育成プログラム、STEPプログラムとBuild@mercariに参加していました。プログラムについてはいつか(出来れば年内に)体験記書こうと思っています。
この2つのプログラムで出たのが**"巡回セールスマン問題のアルゴリズムを極めてください"**という課題。相当時間をかけて色々と試行錯誤したのでそこで学んだことを記事にします。
巡回セールスマン問題とは
循環セールスマン問題、英語ではtraveling salesman problemというのでよく**"TSP"**と略されます。
この問題を一文で表現すると、
「複数の都市を全て必ず一回だけ通るように巡回し、スタート地点に戻ってくるとき、最短の経路は何か？」
というものです。
郵便屋さんの例がよく使われますね。郵便局をスタートし、全ての郵便物を届け、また戻ってくる場合の最短ルートを知りたい という訳です。
どうして解くのが難しいのか
TSPはいわゆるNP-Hard問題です。
NP困難（エヌピーこんなん、英: NP-hard）とは計算量理論において、問題が「NPに属する任意の問題と比べて、少なくとも同等以上に難しい」ことである
(Wikipediaより)
都市数が大きくなると計算量が爆発して今現在の技術だとコンピュータを使っても最小値を求めることができません。そこでこのタイトルにある「極める」という言葉に繋がってくるのです。
手法
アルゴリズムについての話を始める前に、どのようにして初期ノードを生成し、プログラムを書き、実行・結果確認をしたかを説明します。自分でプログラム作成と実行まではしなくていい！アルゴリズムのみに興味がある！という方はアルゴリズムへ飛んでください。
初期ノード生成と結果のvisualizerまで自分でプログラムするのは大変です。アルゴリズムを極めるのが課題だったので、この部分についてはGoogleさんが用意してくださいました。こちらのGitHubをcloneしてきます。適当なフォルダ内で
$ git clone https://github.com/hayatoito/google-step-tsp-2018
$ cd google-step-tsp-2018/week5
とりあえず現在のランダムに設定してある初期経路を確認するには、
# 結果確認
$ ./output_verifier.py
Challenge 0
output : 3862.20
sample/random : 3862.20
sample/greedy : 3418.10
sample/sa : 3291.62
...
と結果が見られます。数字は総移動距離を表しています。
- output: 自分のコーディング結果
- sample/random: randomな経路の結果
- sample/greedy: greedyアルゴリズムの結果
- sample/sa: 前年度STEPプログラム参加者のベストスコア
です。最初は何もコーディングしていないのでoutputとsample/randomが等しくなっています。このoutputをできるだけ小さくすることが目標です。数値だけ見ていても難しいのでvisualizeしてみます。
# visualizerを起動
$ python3 -m http.server
http://localhost:8000/visualizer/build/default/index.html を開くと、下のようなvisualizerが起動するはずです。
これはノード数 N = 8 の場合です。じっと見ているとたくさんの改善ポイントが思いつくと思います。
solver_greedy.pyなどを参考にsolver.py等のファイルを自分で作成し、コーディングしていきます。その結果をoutputに入れるには、
# solverでinput_6の問題を解き, output_6.csvに保存
$ ./solver.py input_6.csv > output_6.csv
# greedy solverで全ての問題を解き、output_{i}.csvに保存
$ for i in `seq 0 6`;do ./solver.py input_${i}.csv > output_${i}.csv; done
です。これを行った後先ほどの結果確認とvisualizer起動を行えば自分のコーディング結果を確認できます。
アルゴリズム
本題です。ネット検索すると本当に色々な解法が載っており、流派も人それぞれです。ここでは私の解いた方法を語っていきます。これから説明する全てのコードはこのGitHubに載っているので変数の確認などに適宜参照してください。
アプローチ方法
最初から最短経路をズバッと求めようなんてことは出来ません。まずは1つルートを生成し、改善を重ねて最適化していきます。
- 初期経路生成
- 最適化アルゴリズムを複数試す
- 相性の良い最適化アルゴリズムを複数組み合わせて更に最適化
という手法を取ります。
初期経路生成
とりあえず巡回ルートを1つ考えてみましょう。適当にランダムに回ってもいいですし、少しだけ頭を使って「一番近い地点へ移動」を繰り返してもいいです。一番近い地点へ移動を繰り返す手法は**貪欲法(greedy)**と呼ばれます。その他にも
- 挿入法
- 最小全域木(MST: minimum spanning tree)
- オイラー路
などがあります。詳しくはこちらのまとめスライド P10〜をご覧下さい。今回は貪欲法と最小全域木の1つである**クラスカル法**を採用しました。
貪欲法
先に述べたとおり、"一番近い地点へ移動"を繰り返す手法です。
- cities: 都市(ノード)のリスト
- current_city: 現在地
- unvisited_cities: 未訪問の都市のリスト
- dist: 2都市間の距離を格納しているリスト
- tour: 回る順番に都市番号を格納しているリスト
として、
def greedy(cities,current_city,dist):
N = len(cities)
unvisited_cities = set(range(0, N))
tour = [current_city]
unvisited_cities.remove(current_city)
while unvisited_cities:
next_city = min(unvisited_cities,
key=lambda city: dist[current_city][city])
unvisited_cities.remove(next_city)
tour.append(next_city)
current_city = next_city
return tour,dist
です。next_cityを選ぶ際にdist[current_city][city]の最小値、つまりcurrent_cityからの距離が最短のcityを選んでいます。
クラスカル法
閉路ができないように重みが小さい順番から辺を選び、追加していく方法です。ここにその説明まで書くと長すぎるので、詳しくはまとめスライドP21〜参照です。
コードは下のようになります。UnionFindを用いることになるので少し複雑です。
from unionfind import UnionFind
def kruskal(dist, N):
tree = UnionFind(N)
relation = [[] for i in range(N)]
one_array_dist = list(itertools.chain.from_iterable(dist))
for k in np.argsort(one_array_dist):
i = k % N
j = int(k / N)
if not tree.issame(i, j) and len(relation[i]) < 2\
and len(relation[j]) < 2:
relation[i].append(j)
relation[j].append(i)
tree.unite(i, j)
def go_next():
for i in range(len(relation[next_city])):
if relation[next_city][i] in unvisited_cities:
return relation[next_city][i]
return None
edge = []
for i in range(N):
if len(relation[i]) == 1:
edge.append(i)
relation[edge[0]].append(edge[1])
relation[edge[1]].append(edge[0])
current_city = 0
unvisited_cities = set(range(1, N))
tour = [current_city]
next_city = relation[current_city][0]
while unvisited_cities:
unvisited_cities.remove(next_city)
tour.append(next_city)
next_city = go_next()
return tour
最適化
次に、初期経路を改善していく方法を考えます。総移動距離をとにかく短くしたいです。
- startノードを変える
- 2-opt / 3-opt (局所最適解に陥る可能性あり)
- or-1-opt / or-2-opt
- 焼きなまし法(SA: Simulated Annealing)
- 距離に重みづけ(minus mean)
ざっと思いつくのはこの辺ですかね、、これらについて見ていきましょう。もちろん他にもたくさんあります。
startノードを変える
これはわかりやすいです。最初のcurrent_cityを変えて色々試します。スタートの都市を全とおり試すなら
for _ in range(N): # スタート地点を全ノードで試す
current_city = _
ただN(都市数)が大きくなってくると処理に時間がかかり待っていられなくなるのでランダムにスタート都市を決定して実行を何回か繰り返し、その最小値を最終的な答えにするといいと思います。
current_city = random.randint(0, len(cities)-1) #最初のノードをランダムに決定
2-opt / 3-opt
2optは「もしcrossを見つけたら、swapして解消」です。
- A B - - A - B -
× ==>
- C D - - C - D -
AD + BC > AB + CDなら上図のように順番を入れ替えます。難しくはないです。
def two_opt(tour, dist):
N = len(tour)
while True:
count = 0
for i in range(N-2):
for j in range(i+2, N):
l1 = dist[tour[i]][tour[i + 1]]
l2 = dist[tour[j]][tour[(j + 1) % N]]
l3 = dist[tour[i]][tour[j]]
l4 = dist[tour[i + 1]][tour[(j + 1) % N]]
if l1 + l2 > l3 + l4:
new_tour = tour[i+1 : j+1]
tour[i+1 : j+1] = new_tour[::-1]
count += 1
if count == 0:
break
return tour
j+1はN以上になる場合があるので%Nをつけることに注意してください。
3-optの場合はここで出てくるノードが6つになります。考え方は同様ですが、少し複雑です。私はこれは実装していないのでWikipediaさんを参照してください。
or-1-opt / or-2-opt
or-1-optは「もしある点が別の2点間に入った方がよければ移動」
- A C - - A - C -
\ /
B ==> B
/ \
- D - E - - D E -
or-2-optは「もしある2点が別の2点間に入った方がよければ移動」
- A D - - A - D -
\ /
B - C ==> B - C
/ \
- D - E - - D E -
図を見ればやることはわかると思います。実装も想像通りです。
def one_or_opt(tour, dist):
N = len(tour)
is_improved = False
while True:
count = 0
for i in range(N):
i0 = i
i1 = (i + 1) % N
i2 = (i + 2) % N
for j in range(N):
j0 = j
j1 = (j + 1) % N
if j0 not in {i0, i1}:
l1 = dist[tour[i0]][tour[i1]]
l2 = dist[tour[i1]][tour[i2]]
l3 = dist[tour[j0]][tour[j1]]
l4 = dist[tour[j0]][tour[i1]]
l5 = dist[tour[j1]][tour[i1]]
l6 = dist[tour[i0]][tour[i2]]
if l1 + l2 + l3 > l4 + l5 + l6:
city = tour.pop(i1)
if i1 < j1:
tour.insert(j0, city)
else:
tour.insert(j1, city)
count += 1
if count != 0:
is_improved = True
return is_improved
insertの部分が少し複雑になります。tourには訪れる都市が順番に格納されているので、その順番に注意する必要があります。
def two_or_opt(tour, dist):
N = len(tour)
is_improved = False
while True:
count = 0
for i in range(N):
i0 = i
i1 = (i + 1) % N
i2 = (i + 2) % N
i3 = (i + 3) % N
for j in range(N):
j0 = j
j1 = (j + 1) % N
if j0 not in {i0, i1, i2}:
l1 = dist[tour[i0]][tour[i1]]
l2 = dist[tour[i2]][tour[i3]]
l3 = dist[tour[j0]][tour[j1]]
l4 = dist[tour[j0]][tour[i1]]
l5 = dist[tour[j1]][tour[i2]]
l6 = dist[tour[i0]][tour[i3]]
if l1 + l2 + l3 > l4 + l5 + l6:
if i2 == 0:
city1 = tour.pop(i1)
city2 = tour.pop(i2)
tour.insert(j0, city2)
tour.insert(j0, city1)
else:
city2 = tour.pop(i2)
city1 = tour.pop(i1)
if i1 < j1:
tour.insert(j0 - 1, city2)
tour.insert(j0 - 1, city1)
else:
tour.insert(j1, city2)
tour.insert(j1, city1)
count += 1
if count != 0:
is_improved = True
return is_improved
焼きなまし法
英語でSimulated Annealing, SAと呼ばれます。これは2-optで陥る可能性のある局所最適解を抜け出すために導入されることが多いです。
単にエネルギー状態(今回であれば経路の長さ)が低いほうに進むだけでは大域的最適解に至らず、局所的最適解に留まってしまうことがあるのでその対策として使われる手法です。
物理世界での熱による**"振動"**でよく喩えられます。物理世界は温度が高いときはよく振動し、温度が下がるにつれて安定な状態へと落ち着きます。
これと同じイメージで、
- 計算の最初のうちは経路長が長くなる場合でも確率的に経路の交換をするようにする。
- イテレーションの終りの方では経路長が長くなる時に経路の交換を行う確率が低くなるようにする。
ということです。焼きなまし法のめちゃくちゃわかりやすい動画があったのでこれをご覧下さい。実装はしてないのでコードないです、すみません。
距離に重みづけ(minus mean)
各エッジに重みづけをします。具体的には、2点間の距離を両端のノードの距離平均が大きい場合は、局所最適解に陥らないように、優先的にエッジを選択します。
def minus_mean_distance(dist, N):
ave_dist = np.mean(dist, axis=1)
dist = [[dist[i][j] - (ave_dist[i] + ave_dist[j])
for i in range(N)] for j in range(N)]
return dist
何度も出てきていますがdistは2都市間の距離を格納しているリスト、Nは都市数です。
ただ、これだとNが大きくなってきた場合のave_distを求める際の計算量が膨大になってしまいます。そこで登場するのがk-meansです。
def minus_k_mean_distance(dist, N, k):
ave_dist = np.mean(np.sort(dist, axis=1)[:, 1:min(N, k + 1)], axis=1)
dist = [[dist[i][j] - (ave_dist[i] + ave_dist[j])
for i in range(N)] for j in range(N)]
return dist
min(N, k + 1)で、Nとk+１のうち小さい方の数だけ計算考慮に入れることにしてあります。実際に計算実行する場合はこちらを使うことになります。
最適化アルゴリズムを組み合わせる
上にあげた以外にも最適化アルゴリズムは山ほどあります。気になった人は調べてください！
この記事ではここで終わりにして、ここからはこれらを組み合わせで出来る最善解を調べます。この過程で気になることは
- 基本的には最適化をたくさん行えば行うほど良い結果となる
- ただしノード数が大きくなった場合実行に時間がかかる
- 局所最適解に陥る可能性がある最適化は最初に行う必要あり
- 初期経路、スタートノードを変えて最適化を行って結果を比較してみる
などです。
これに関しては様々な組み合わせをひたすら自力で調べるしかありません。ここまでで使える手段をまとめると、
- 初期経路
- 貪欲法
- クラスカル法
- 最適化
- startノードを変える
- 2-opt
- or-1-opt
- or-2-opt
- minus k mean
です。
ひたすらコード変更&実行を繰り返しましょう。
結果
試行錯誤の様子をここに書いてもしょうがないのでずばり最善結果です。
Best Scoreは
- 初期経路: kruskal
- 距離の重み付け: minus-10-mean (実行時間の関係で10としました)
- 最適化処理: 2-opt + or-1-opt + or-2-opt
- スタートノード: min(N, 30) (これも実行時間の関係でmax30とした)
のとき。
| | N = 4 | N = 8 | N = 16 | N = 64 | N = 128 | N = 512 | N = 2048 |
| ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- |
| random | 3862.20 | 6101.57 | 13479.25 | 47521.08 | 92719.14 | 347392.97 | 1374393.14 |
| greedy | 3418.10 | 3832.29 | 5449.44 | 10519.16 | 12684.06 | 25331.84 | 49892.05 |
| Best score | 3418.10 | 3832.29 | 4994.89 | 8970.05 | 11489.79 | 21363.60 | 42712.37 |
どうでしょうか？Nが大きくなるにつれて最適化の効果がよく見えるようになったと思います。
初期経路greedyのとき
初期経路greedyに各最適化アルゴリズムを適用した結果
| | N = 5 | N = 8 | N = 16 | N = 64 | N = 128 | N = 512 | N = 2048 |
| ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- |
| greedy | 3418.10 | 3832.29 | 5449.44 | 10519.16 | 12684.06 | 25331.84 | 49892.05 |
| greedy + 2-opt | 3418.10 | 3832.29 | 4994.89 | 8970.05 | 11489.79 | 21363.60 | 42712.37 |
| greedy + 2-opt + or-1-opt | 3291.62 | 3778.72 | 4994.42 | 8656.07 | 11225.87 | 20902.75 | 41638.84 |
| greedy + 2-opt + or-1-opt + or-2-opt | 3291.62 | 3778.72 | 4994.42 | 8656.07 | 11225.87 | 20902.75 | 41638.84 |
最適化アルゴリズムは2-opt + or-1-opt + or-2-optで固定し、初期経路を変えたとき
| | N = 5 | N = 8 | N = 16 | N = 64 | N = 128 | N = 512 | N = 2048 |
| ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- |
| greedy | 3291.62 | 3778.72 | 4494.42 | 8656.07 | 11225.87 | 20902.75 | 41638.84 |
| minus_mean | 3291.62 | 3778.72 | 4494.42 | 8338.24 | 11047.07 | 21178.11 | 42104.09 |
| minus_10_mean | 3291.62 | 3778.72 | 4994.42 | 8497.71 | 11338.02 | 21219.00 | 41247.10 |
| kruskal | 3291.62 | 3778.72 | 4994.42 | 8451.45 | 10925.97 | 20894.69 | 41110.98 |
| kruskal + start node 30通り | 3291.62 | 3778.72 | 4994.42 | 8118.40 | 10539.04 | 20263.87 | 40537.49 |
以上になります。初期経路生成や、最適化アルゴリズムの組み合わせで随分と結果が変わってくることが分かったと思います。
他にも
- もっと多い種類のスタートノードで試す
- 距離の重み付け計算を工夫する
- 焼きなまし法を導入する
- プリム法でやってみる
などなど、出来る事はたくさんあります。余裕のある方は試してみてください。(良い結果が出たらぜひ教えてください笑)
終わりに
Advent Calendar登録したはいいものの、時間無さすぎて1日でかけるテーマ選んでしまった、、数ヶ月前にやった事だから昔の自分のコード読んで「読みにくっ！」てところがたくさんありました。TSPは本当に奥が深い問題で、まだまだ改良の余地があるので。私はここまでで力尽きたので誰か続編でも書いてください...
でもそうこうしているうちにいつかこれを一瞬で溶けるスーパーコンピュータが出来るんでしょうね。