TSP（巡回セールスマン問題）について調べていくと、派生的に色々なアルゴリズムに出くわします。それと、ここしばらくはPythonを使っていますが、Python自体の書き方においてもいまだに発見があります（便利というか奥が深いというか）。
TSPにおいては、初期経路生成のアルゴリズムとその後の経路最適化アルゴリズムの二段階に分かれます。まずは、初期経路生成について色々試しています。
前回は凸包（convex hull）による初期経路を試してみましたが、単純な貪欲法よりはまともな初期経路ができあがりました。結論的には、凸包の挿入法による初期経路形成のほうが今回のMST（minimum spanning tree）による方法よりもいいような気がしますが、実際どうなのかまだ検証していないので分かりません（追記：その後簡単な比較検証してみた結果、MSTを用いたほうが優れていました）。いずれにせよMSTは様々な分野において応用が効きそうなので勉強してみようかと。
今回は、最小全域木（MST: minimum spanning tree）を使って初期経路を生成しようと思います。MSTはそれ自体では枝分かれしたグラフであり、そのままではTSPの巡回路としては使えないのですが、MSTを生成することで様々なことに応用できるようです。枝分かれしたそれぞれの辺は最小値になるように選択しているため、TSPにおける経路の上限値と下限値などがわかるようです。MSTの生成アルゴリズムとしては、主にクラスカル法とプリム法の二つがあるようです。
クラスカル法：
全ての辺の組み合わせを計算して最短の辺から順番にリストに加えていきます。
統合：
リストに加えて行く際に、すでに選択した辺に繋がる場合は統合します。辺は二つのノード（点）なので、新たに加える辺のどちらか一方のノードが、すでに選択した辺のノードと一致する場合は同じグループとみなします。今回の場合はリスト内にグループごとの子リストをつくって仕分けしていきます。
グループ追加：
加えようとする辺の二つのノードが、既に選択されている辺のどれにも一致しない場合は統合されないため、新たな子グループ（リスト）として親リスト内に追加します。
却下：
グループに新たな辺を加える際には閉路ができないようにします。加えようとする辺のノードが二つともすでに選択した辺と一致してしまう場合は閉路とみなすことができるので選択を却下します。
このように統合、グループ追加、却下を繰り返していくと、全てのノードが繋がって一つのグループになります。
set()関数：
ある数値（ノードの番号）が、あるリスト内の数値に一致するかどうかは、集合関数であるset()を使うといいようです。
A = [1, 2, 3, 4, 5]
B = [0, 2]
set(A) & set(B):
とすれば、
{2}
という出力があり（2がどちらにも共通という）、リストAに対してリストBの部分一致を確かめることができます。つまりBのどちらか一方がAに所属するということがわかるというわけです。
set(A) >= set(B)
なら、Bの全てがAに含まれるかどうかがわかります。
ということで、この一連の流れをコーディングしてみました。実際は計算量を軽減するための工夫も必要ですが、今回は結果だけを求めています。
ノード50個の場合。
このように最短の辺で構成された枝分かれのグラフができあがります。
すべてのノードへ訪問可能にするために必要な経路の合計が最短ということです。
この段階では、
[[4,15],[23,45]...]
というような二つのノードを組にした辺のリストとして順序は無関係に並んでいます。
プリム法：
もう一つのMST生成アルゴリズムであるプリム法です。
プリム法では、まず任意のノードを選び、それに最も近いノードを次に選びます。選んだノードはリストに加えていきますが、一組の数値（辺：[a, b]）として加えていきます。次は、選択した二点のノードのそれぞれに近いノードを残りのノードから選びだし、より小さいほうだけをリストに加えます。このようにリストに加えたノードに最も近いノードを残りのノードから一つずつ選んで辺として加えていきます。
クラスカル法に比べアルゴリズムは単純ですが、毎回選択したノード数分の距離を総当たりで計算させているため計算量がかさんでしまいます。あらかじめ距離の総当たりテーブルを用意しておけばいいのかもしれません。
TSP初期経路生成：
MSTを元にTSPのための初期経路を生成します。主には２種類あり、今回はDouble TreeアルゴリズムとChristofidesアルゴリズムをコーディングしてみました。大まかな内容を理解してから勝手にアルゴリズムを組み立てているので、厳密なアルゴリズムにはなっていません。
オイラー閉路：すべての辺を通り開始点に戻ることができる。
ハミルトン閉路：すべてのノードを通り開始点に戻ることができる。
Double Treeアルゴリズム：
このアルゴリズムではMSTの上にもう一つ同じ経路を重ね合わせて二重経路にします。MSTを経路というよりは壁と見立てて、その周囲を壁沿いに移動していくイメージです。そうすることで、どの経路も往復路になり、行き止まりになる枝の先端部分であっても戻ることができるため全てのノードに到達可能となります。TSPでは同じノードを二回通ることはできませんが、まずはグラフ理論におけるオイラー路（すべての辺を通るグラフ、同じノードを複数回訪れてもよい）として考え、さらに循環する場合はオイラー閉路とみなしてすべてのノードへの経路を考えます。そのため単純計算すれば、最大でMSTの二倍の距離になります。最終的にはTSPの経路に変換するため同じノードをスキップします。
結果的にはあまりきれいな経路にはなりませんが、このあとの最適化によってまともな経路になることに期待します。
Christofidesアルゴリズム：
このアルゴリズムはやや複雑で、Double Treeでは2倍長の経路であったのに対し、グラフ理論における最適マッチングを適用することで加える経路を半分に減らすことができ、合計1.5倍の経路で巡回可能になるようです。
次数が奇数のノードを選択する：
まず、MST上で接続されている辺の数が奇数のノードをすべてピックアップします。例えば三つに分岐しているノードなどです。枝の先端も接続されている辺の数は一本だけなので奇数ノードとみなします（ノードにおける辺の接続数のことを「次数」という）。
ノードにおける奇数と偶数の違いは、そこにたどり着く経路とそこから出ていく経路の数として見るといいようです。同じ経路を二回通らないという設定であれば、三分岐点の場合は、その地点から出発して戻ってきたあと再度出発すれば二度と戻ってこれないということになり、往復経路のどちらか一方が足りないということになります。そのためこれら奇数のノードに対してもう一本経路を追加してあげれば往復可能、通過可能、巡回可能なノードになります。
最適マッチング（minimum weight perfect matching）：
奇数ノードだけを抜き出し、それらだけで最適マッチングによる辺（サブ経路：M）の追加を行います。サブ経路Mが追加されることで、行き止まりとなる枝の先端ノードや分岐点が巡回可能になります。
左図（上図）はノード数20個のMST。この中から次数が奇数のノードを選択します。この場合、三又の分岐点（3箇所）以外に枝の先端（5箇所）も含め8箇所あります。そして最適マッチングによりそれらの8点を4本（8ノードの半数）の辺としてつなぎます（右図/下図）。最適マッチングによる辺は互いに接続せず離れ離れに配置されます。連結した辺というよりは破線のように一つ飛ばしで辺を構成していく感じです。この4辺をMとして、MSTの辺と重ね合わせます。これで全てのノードにおいて偶数の経路と接続されたこととなります。
このときminimum-weightになるように、合計の辺の長さが短くなるような組み合わせを選びます。8点あるので組み合わせの数は7*5*3*1=105パターンになります。もし10点ある場合は9*7*5*3*1=945というように数が増えるほど計算量が膨大になっていきます。今回の場合そのまま計算させるとかなり時間がかかってしまったので（ベタに計算させたコードはこちら）、グラフ理論ライブラリであるNetworkXのmax_weight_matching()関数を利用しました（短時間で計算してくれます）。ただし、NetworkXにはminimum_weight_matchingの関数はないので、パラメータをマイナス反転して結果だけを利用しています。
尚、最適マッチングの計算方法としては、Blossom algorithmやHungarian algorithmを参考にするといいようです。
最適マッチングによって辺Mが追加された結果、枝の先端は行き止まりになることなく、また三又の分岐点も巡回可能なノードになります。この状態から新たにTSP用に経路を構築していきます。
初期経路生成：
MSTを元に、Double TreeやChristofidesアルゴリズムを経てTSPに合わせた経路に変換します。前段階まではオイラー閉路（全ての辺を通る、同じノードを数回通ってもよい）であったので、同じノードを一回しか通ることができない経路にするため、複数回通るノードをスキップします。
MSTや最適マッチングによって生成された辺の集合は、
[[12,20],[4,8],[6,22]...]
などと順序は無関係に並んでいるので、これらを
[[12,20],[20,23],[23,5]...]
というように前後の番号がつながるように並び替えます。
そして最終的には、
[0,12,20,23...18,15,0]
というような重複しないノードのリストに変換してTSP経路をつくります。
基本的には、各辺は[a,b]というような一組のノード番号によってリストに入っており、次につながる辺は[b,c]もしくは[c,b]となるため前後を逆にした二通りを検索し、つながる辺が見つかれば追加していきます。
しかし、そのまま各辺を追加してくと、開始点によっては途中で道が途絶えてしまうことがあります。これは、開始点と現在地とで局所的ループが発生して、それ以外の経路へ移動できない状態になっているということです。局所的ループによって行き止まりになったら、現在加えた辺を経路リストの最後から先頭に移動して現在地を一つ手前に戻します。それでも接続可能な辺を見つけられなければ、接続可能な辺が見つかるまでこの手順を繰り返します。
こうすることで、どの地点から開始しても現在地と開始点をずらしながら通り過ぎてしまった分岐点に遡ることが可能となり、見逃していた辺（経路）を追加していくことができます。
そして順路ができあがったらノード0を起点に並べ替えて、最後に再度ノード0を付け足して巡回路とします。
ノード数20の場合。左（上）がDouble Tree。右（下）がChristofides。
最終的に重複するノードをスキップさせるためか経路が部分的に交差してしまうことがあります。二つを何回か比較すると、Christofidesのほうがより安定した経路を形成します。
そして最適化させると以下のようになります。交差している部分は解消されほぼ厳密解に近い経路になります。
TSPにおいては、初期経路生成のアルゴリズムとその後の経路最適化アルゴリズムの二段階に分かれます。まずは、初期経路生成について色々試しています。
前回は凸包（convex hull）による初期経路を試してみましたが、単純な貪欲法よりはまともな初期経路ができあがりました。結論的には、凸包の挿入法による初期経路形成のほうが今回のMST（minimum spanning tree）による方法よりもいいような気がしますが、実際どうなのかまだ検証していないので分かりません（追記：その後簡単な比較検証してみた結果、MSTを用いたほうが優れていました）。いずれにせよMSTは様々な分野において応用が効きそうなので勉強してみようかと。
今回は、最小全域木（MST: minimum spanning tree）を使って初期経路を生成しようと思います。MSTはそれ自体では枝分かれしたグラフであり、そのままではTSPの巡回路としては使えないのですが、MSTを生成することで様々なことに応用できるようです。枝分かれしたそれぞれの辺は最小値になるように選択しているため、TSPにおける経路の上限値と下限値などがわかるようです。MSTの生成アルゴリズムとしては、主にクラスカル法とプリム法の二つがあるようです。
クラスカル法：
全ての辺の組み合わせを計算して最短の辺から順番にリストに加えていきます。
統合：
リストに加えて行く際に、すでに選択した辺に繋がる場合は統合します。辺は二つのノード（点）なので、新たに加える辺のどちらか一方のノードが、すでに選択した辺のノードと一致する場合は同じグループとみなします。今回の場合はリスト内にグループごとの子リストをつくって仕分けしていきます。
グループ追加：
加えようとする辺の二つのノードが、既に選択されている辺のどれにも一致しない場合は統合されないため、新たな子グループ（リスト）として親リスト内に追加します。
却下：
グループに新たな辺を加える際には閉路ができないようにします。加えようとする辺のノードが二つともすでに選択した辺と一致してしまう場合は閉路とみなすことができるので選択を却下します。
このように統合、グループ追加、却下を繰り返していくと、全てのノードが繋がって一つのグループになります。
set()関数：
ある数値（ノードの番号）が、あるリスト内の数値に一致するかどうかは、集合関数であるset()を使うといいようです。
A = [1, 2, 3, 4, 5]
B = [0, 2]
set(A) & set(B):
とすれば、
{2}
という出力があり（2がどちらにも共通という）、リストAに対してリストBの部分一致を確かめることができます。つまりBのどちらか一方がAに所属するということがわかるというわけです。
set(A) >= set(B)
なら、Bの全てがAに含まれるかどうかがわかります。
ということで、この一連の流れをコーディングしてみました。実際は計算量を軽減するための工夫も必要ですが、今回は結果だけを求めています。
def kruskal(path): dist = [] idx = [] for i in range(len(path)): for j in range(i+1, len(path)): dist.append(np.abs(XY[i] - XY[j])) idx.append([i,j]) arg = np.argsort(dist) edge_list = [idx[a] for a in arg] plist = [] group = [] for edge in edge_list: if set(edge) & set(sum(group, [])): # either in all groups merge_list = [] for i in range(len(group)): if set(edge) <= set(group[i]): # both in one group continue elif set(edge) & set(group[i]): # either in one group merge_list.append(i) if len(merge_list) > 0: group[merge_list[0]] = list(set(sum([group[m] for m in merge_list] + [edge], []))) [group.pop(i) for i in merge_list[1:]] plist.append(edge) else: group.append(edge) plist.append(edge) return plist最初に全ての辺の組み合わせを計算しなければいけないので、ノード数が増えるとかなりきつくなりそうです。この部分に工夫が必要そうです。今回はベタに計算しています。
このように最短の辺で構成された枝分かれのグラフができあがります。
すべてのノードへ訪問可能にするために必要な経路の合計が最短ということです。
この段階では、
[[4,15],[23,45]...]
というような二つのノードを組にした辺のリストとして順序は無関係に並んでいます。
プリム法：
もう一つのMST生成アルゴリズムであるプリム法です。
プリム法では、まず任意のノードを選び、それに最も近いノードを次に選びます。選んだノードはリストに加えていきますが、一組の数値（辺：[a, b]）として加えていきます。次は、選択した二点のノードのそれぞれに近いノードを残りのノードから選びだし、より小さいほうだけをリストに加えます。このようにリストに加えたノードに最も近いノードを残りのノードから一つずつ選んで辺として加えていきます。
クラスカル法に比べアルゴリズムは単純ですが、毎回選択したノード数分の距離を総当たりで計算させているため計算量がかさんでしまいます。あらかじめ距離の総当たりテーブルを用意しておけばいいのかもしれません。
def nearest_node_dist(path, index): best_dist = np.inf best_node = None for p in path: if p != index: d = abs(XY[p] - XY[index]) if d < best_dist: best_dist = d best_node = p return best_node, best_dist def prim(path): copy = list(path) if path[0] == path[-1]: copy = list(path[:-1]) node, _ = nearest_node_dist(copy, path[0]) plist = [[path[0], node]] copy.remove(path[0]) copy.remove(node) while len(copy): best_dist = np.inf best_edge = None for p in set(np.ravel(plist)): node, dist = nearest_node_dist(copy, p) if dist < best_dist: best_dist = dist best_edge = [p, node] plist.append(best_edge) copy.remove(best_edge[1]) return plist今回は毎回各ノードに対して最も近いノードとの距離計算をさせているためnearest_node_dist()という関数を別途用意してあります。
浅野 孝夫
近代科学社
売り上げランキング: 444,789
近代科学社
売り上げランキング: 444,789
TSP初期経路生成：
MSTを元にTSPのための初期経路を生成します。主には２種類あり、今回はDouble TreeアルゴリズムとChristofidesアルゴリズムをコーディングしてみました。大まかな内容を理解してから勝手にアルゴリズムを組み立てているので、厳密なアルゴリズムにはなっていません。
オイラー閉路：すべての辺を通り開始点に戻ることができる。
ハミルトン閉路：すべてのノードを通り開始点に戻ることができる。
Double Treeアルゴリズム：
このアルゴリズムではMSTの上にもう一つ同じ経路を重ね合わせて二重経路にします。MSTを経路というよりは壁と見立てて、その周囲を壁沿いに移動していくイメージです。そうすることで、どの経路も往復路になり、行き止まりになる枝の先端部分であっても戻ることができるため全てのノードに到達可能となります。TSPでは同じノードを二回通ることはできませんが、まずはグラフ理論におけるオイラー路（すべての辺を通るグラフ、同じノードを複数回訪れてもよい）として考え、さらに循環する場合はオイラー閉路とみなしてすべてのノードへの経路を考えます。そのため単純計算すれば、最大でMSTの二倍の距離になります。最終的にはTSPの経路に変換するため同じノードをスキップします。
結果的にはあまりきれいな経路にはなりませんが、このあとの最適化によってまともな経路になることに期待します。
Christofidesアルゴリズム：
このアルゴリズムはやや複雑で、Double Treeでは2倍長の経路であったのに対し、グラフ理論における最適マッチングを適用することで加える経路を半分に減らすことができ、合計1.5倍の経路で巡回可能になるようです。
次数が奇数のノードを選択する：
まず、MST上で接続されている辺の数が奇数のノードをすべてピックアップします。例えば三つに分岐しているノードなどです。枝の先端も接続されている辺の数は一本だけなので奇数ノードとみなします（ノードにおける辺の接続数のことを「次数」という）。
ノードにおける奇数と偶数の違いは、そこにたどり着く経路とそこから出ていく経路の数として見るといいようです。同じ経路を二回通らないという設定であれば、三分岐点の場合は、その地点から出発して戻ってきたあと再度出発すれば二度と戻ってこれないということになり、往復経路のどちらか一方が足りないということになります。そのためこれら奇数のノードに対してもう一本経路を追加してあげれば往復可能、通過可能、巡回可能なノードになります。
最適マッチング（minimum weight perfect matching）：
奇数ノードだけを抜き出し、それらだけで最適マッチングによる辺（サブ経路：M）の追加を行います。サブ経路Mが追加されることで、行き止まりとなる枝の先端ノードや分岐点が巡回可能になります。
このときminimum-weightになるように、合計の辺の長さが短くなるような組み合わせを選びます。8点あるので組み合わせの数は7*5*3*1=105パターンになります。もし10点ある場合は9*7*5*3*1=945というように数が増えるほど計算量が膨大になっていきます。今回の場合そのまま計算させるとかなり時間がかかってしまったので（ベタに計算させたコードはこちら）、グラフ理論ライブラリであるNetworkXのmax_weight_matching()関数を利用しました（短時間で計算してくれます）。ただし、NetworkXにはminimum_weight_matchingの関数はないので、パラメータをマイナス反転して結果だけを利用しています。
尚、最適マッチングの計算方法としては、Blossom algorithmやHungarian algorithmを参考にするといいようです。
def min_weight_perfect_matching(edges): # find odd-degree nodes from MST(edges) odd = [] for i in range(num): cnt = 0 for e in np.ravel(edges): if i == e: cnt += 1 if cnt % 2 == 1: odd.append(i) print('Odd-Degree Edges:', odd) # get combinations of edges from the odd-degree nodes combs = list(combinations(odd, 2)) # calculate weights on the edges weight = [abs(XY[c[0]]-XY[c[1]]) for c in combs] # negating the weights to get minimum-weight from the max_weight_matching function edge_weight = [(combs[i][0], combs[i][1], {'weight': -weight[i]}) for i in range(len(combs))] # find minimum-weight perfect matching(negative max_weight_matching) by NetworkX library G = nx.Graph(edge_weight) M = nx.max_weight_matching(G, True, 'weight') M = [list(m) for m in M] G.clear() return M
最適マッチングによって辺Mが追加された結果、枝の先端は行き止まりになることなく、また三又の分岐点も巡回可能なノードになります。この状態から新たにTSP用に経路を構築していきます。
初期経路生成：
MSTを元に、Double TreeやChristofidesアルゴリズムを経てTSPに合わせた経路に変換します。前段階まではオイラー閉路（全ての辺を通る、同じノードを数回通ってもよい）であったので、同じノードを一回しか通ることができない経路にするため、複数回通るノードをスキップします。
MSTや最適マッチングによって生成された辺の集合は、
[[12,20],[4,8],[6,22]...]
などと順序は無関係に並んでいるので、これらを
[[12,20],[20,23],[23,5]...]
というように前後の番号がつながるように並び替えます。
そして最終的には、
[0,12,20,23...18,15,0]
というような重複しないノードのリストに変換してTSP経路をつくります。
# T, M=False: Double Tree Algorithm # T, M=M : Cristofides Algorithm def tsp_tour(T, M=False): if M: TM = (T + M).copy() else: TM = (T + T).copy() tour = [TM[0]] TM.remove(TM[0]) while len(TM): for tm in TM: if tour[-1][1] == tm[0]: # tour(n,0) == TM(0,m) tour.append(tm) TM.remove(tm) elif tour[-1][1] == tm[1]: # tour(n,0) == TM(m,0) tour.append(tm[::-1]) TM.remove(tm) if tour[0][0] == tour[-1][1]: # when meets deadend tour = [tour[-1]] + tour[:-1] newtour = [] for t in tour: if not t[0] in newtour: newtour.append(t[0]) zero = newtour.index(0) newtour = newtour[zero:] + newtour[:zero] + [0] return newtourMSTで生成した経路をT、最適マッチングによって追加された経路をMとしています。なお引数にMを渡さなければ、Double Treeアルゴリズムで経路生成します。その場合はTを2倍にした経路が初期状態となります。
基本的には、各辺は[a,b]というような一組のノード番号によってリストに入っており、次につながる辺は[b,c]もしくは[c,b]となるため前後を逆にした二通りを検索し、つながる辺が見つかれば追加していきます。
しかし、そのまま各辺を追加してくと、開始点によっては途中で道が途絶えてしまうことがあります。これは、開始点と現在地とで局所的ループが発生して、それ以外の経路へ移動できない状態になっているということです。局所的ループによって行き止まりになったら、現在加えた辺を経路リストの最後から先頭に移動して現在地を一つ手前に戻します。それでも接続可能な辺を見つけられなければ、接続可能な辺が見つかるまでこの手順を繰り返します。
こうすることで、どの地点から開始しても現在地と開始点をずらしながら通り過ぎてしまった分岐点に遡ることが可能となり、見逃していた辺（経路）を追加していくことができます。
そして順路ができあがったらノード0を起点に並べ替えて、最後に再度ノード0を付け足して巡回路とします。
最終的に重複するノードをスキップさせるためか経路が部分的に交差してしまうことがあります。二つを何回か比較すると、Christofidesのほうがより安定した経路を形成します。
そして最適化させると以下のようになります。交差している部分は解消されほぼ厳密解に近い経路になります。
比較結果（追記）：
その後、前回の凸包・挿入法と今回のMSTによる初期経路生成について、TSPライブラリのPyconcordeの結果を元に100回分のスコアを比較しました。
Pyconcorde ：100.0000 %
凸包・挿入法：103.3175 %
Double Tree ：105.5888 %
Christofides：102.7978 %
それぞれのアルゴリズムで経路生成したあとに最適化アルゴリズムを通して得られたスコアです。Pyconcordeのスコアを100%とした場合の経路長の率で、%が低いほど短距離でいい結果ということです。今回は簡単な最適化しか行っていないためもう少し改善することはできましたが、Pyconcordeより短い距離を出すことはほぼ無理。何回かに一回はPyconcordeと同じスコアになるときはあります。
凸包・挿入法は一見よさそうだったのですが、結果的にはChristofidesのほうがわずかながら優れています。凸包とはいえ仕組み的には貪欲法なので、やはりMSTを構築したほうが結果がよくなるという感じでしょうか。MSTは最短の辺で全体をつなぎとめるので、全体的な構成や骨格が最小で済み、あとは局所的な最適化さえすればいいという感じで理にかなっているのかもしれません。
関連：
ウィリアム・J・クック
青土社
売り上げランキング: 255,242
青土社
売り上げランキング: 255,242
