TSP DP: 巡回セールスマン問題 / 動的計画法 / メモ化再帰

前回の再帰関数と動的計画法の続きとして、今回は巡回セールスマン問題（TSP/Travelling Salesman Problem）の厳密解を求める動的計画法（DP/Dynamic Programming）を試してみました。

再帰関数の例題によく出てくるフィボナッチ数などは単純なので分かりやすいのですが、TSPになると途端に難しくなり、中身を理解するのに予想以上に時間がかかってしまいました。

今回もまた理解した内容をできるだけ自前でプログラミングしてみたので、正式な動的計画法になっているかはわかりませんが一応機能します（DPといってもいろいろやり方があるようなので、どれがいいのかよく分からない）。

調べてみると、ビッド演算子を使うbit DPという方法が便利そうですが、今回はビッド演算子を使わず通常のforループとset()関数で試してみました。

驚きの数学 巡回セールスマン問題

posted with amazlet at 19.05.19

ウィリアム・J・クック
青土社
売り上げランキング: 100,713

Amazon.co.jpで詳細を見る

下準備：（Python 3.6、Jupyter Notebook使用）

まずノードの集合（簡単に4ノード）として、

S = {0, 1, 2, 3}

をPythonのセット関数（集合関数）で用意します。

0からスタートして残りの{1, 2, 3}を巡って最後にまた0に戻ってくることにします。同時に各ノード間の距離テーブルdist[a][b]も用意しておきます。

dist = [[0, 1, 5, 4]
　　　　[1, 0, 3, 6]
　　　　[5, 3, 0, 2]
　　　　[4, 6, 2, 0]]

そして距離を求める再帰関数として、

TSP(a, S-{a}, b)

を用意します。aからスタートして、ノード集合S-{a}を全て巡ってbに戻るという設定です。

実際の数値を代入すれば、

TSP(0, {1,2,3}, 0)

になります。

S-{a}というのは集合関数で、

S-{0} = {0, 1, 2, 3} - {0}　= {1, 2, 3}

となり、集合Sから要素{0}を除去するという意味です。ちなみに要素を付け足すにはプラス（+）ではなくS.add(0)になります。

TSP DPの考え方：

DPの特長として、全体を部分に切り分けて答えを求めていく「分割統治法」というのがあります。その際、部分的な答えを記憶保存しておく「メモ化」というテクニックを使います。
その際求めたい答え（TSPなら全体的な道順）を後回しにして、とりあえず目先の部分的な処理（次の行き先）だけを考えますが、この再帰的な仕組みが理解しにくい。必ずしも再帰ループを使わなくてもいいらしいですが、部分を解いて結果的に全体を解くという流れにすると再帰的になってしまいます。再帰的な途中処理の結果をメモ化で記録して再利用することで効率があがるようです。
まずは、最短経路の合計距離だけを求める式を考えてみます（道順についてはあとで求めます）。

先ほど用意した

TSP(0, {1,2,3}, 0)

というのは、左の0から出発し、次に順番はどうであれ中継地点{1, 2, 3}をそれぞれ一回ずつ通って目的地0に到着するということを表現しています。

TSP(出発地, {中継地点}, 目的地)

というわけです。最終的にはこの式から距離が求まりますが再帰的な処理をしていくので、いきなりは答えがでてきません。
{1, 2, 3}の道順の組み合わせは6通りありますが、とりあえず次に向かう先だけを考えます。そうすると1か2か3の3通りだけを考えることになります。残り2つは気になるけど後回しにしておきます。

次の手続きとして、もし中継地点のうち1に移動すると0から1までの距離dist[0][1]と、次の出発点1として残り{2, 3}が中継地点となります。

それを式に表すと、

TSP(0, {1, 2, 3}, 0) = dist[0][1] + TSP(1, {2, 3}, 0)

ということになります。全体の道順をいきなり解くというより、次の行き先だけを解いて残りは保留（後で解く）という感じです。次の行き先までの距離dist[0][1]は距離テーブルから求められますが、TSP(1, {2, 3}, 0)の距離はこのままでは計算できません。ちなみに、この段階では出発点0から最も近い中継地点に向かうことがベストとは限りません。出発点0から一番遠い地点であっても残りの経路TSP(1, {2, 3}, 0)で挽回することがありうるので、とりあえず式だけ組んでおくという感じです。
そして次の可能性を考えますが、もし1の次に2に向かうなら、

TSP(1, {2, 3}, 0) = dist[1][2]　+ TSP(2, {3}, 0)

と表すことができます。これは1から2へ向かって（距離は計算可能で確定）、2から3へ向かうということになります（予定）。このように確定と予定の二つで式をつくっていきます。あとは中継地点が{3}しかないので、3を次の出発点にして0へ向かうということになりますが、同様の手順を踏んでTSP(2, {3}, 0)の部分を、

TSP(2, {3}, 0) = dist[2][3] + TSP(3, {}, 0)

と表すことができ、さらに

TSP(3, {}, 0) = dist[3][0]

となり、中継地点はないので3から0へ直接移動するということになります。
この一連の流れでは、{}内の要素が一つずつ減っていき空集合になるまで進めていくと、最終的にはdist[a][b]として距離を求められます。また、中継地点のうちどれか一つが次のステップでの出発点となり、左側の変数aは変化し続けます。右側の変数bは最終目的地（出発地でもある）なのでずっと固定値です。
こうやって未定の部分を可能性として掘り下げていくと距離テーブルをつかって実際の値を求めていくことができます。{}内の要素がゼロになったとき、この再帰ループが終わるので、それを終了条件としておけばいいことになります。
あとで気づきましたが、TSP(2, {3}, 0)を求める段階でこの距離はdist[2][3] + dist[3][0]とわかるので、{}内の要素が一つのとき再帰ループの終了条件とし、それ以降のTSP(3, {}, 0)までループさせる必要がないかもしれません。

しかし、上記で行ったのは一つのパターンであって、実際はもっと多くの組み合わせがあります。最初のTSP(0, {1, 2, 3}, 0)のときは、0から次に向かう先は1か2か3の3通りあり、そのうちの最短経路を求めたいのでmin()関数を使って書き換えると、

TSP(0, {1, 2, 3}, 0) = min(dist[0][1] + TSP(1, {2, 3}, 0), 
                           dist[0][2] + TSP(2, {1, 3}, 0), 
                           dist[0][3] + TSP(3, {1, 2}, 0))

になります。これらのうちの最小のものを採用することになります。

しかし、このままではTSP()の部分が計算できないので、先ほどやったように以下のパターン全てを掘り下げて計算していきます。

TSP(1, {2, 3}, 0) = min(dist[1][2] + TSP(2, {3}, 0),
                        dist[1][3] + TSP(3, {2}, 0))

TSP(2, {1, 3}, 0) = min(dist[2][1] + TSP(1, {3}, 0),
                        dist[2][3] + TSP(3, {1}, 0))

TSP(3, {1, 2}, 0) = min(dist[3][1] + TSP(1, {2}, 0),
                        dist[3][2] + TSP(2, {1}, 0))

同様に右辺式後半の各TSP()部分が未確定なので、

TSP(2, {3}, 0) = dist[2][3] + TSP(3, {}, 0)
TSP(3, {2}, 0) = dist[3][2] + TSP(2, {}, 0)
TSP(1, {3}, 0) = dist[1][3] + TSP(3, {}, 0)
TSP(3, {1}, 0) = dist[3][1] + TSP(1, {}, 0)
TSP(1, {2}, 0) = dist[1][2] + TSP(2, {}, 0)
TSP(2, {1}, 0) = dist[2][1] + TSP(1, {}, 0)

となり、さらに

TSP(3, {}, 0) = dist[3][0]
TSP(2, {}, 0) = dist[2][0]
TSP(1, {}, 0) = dist[1][0]

となり、{}が空集合になるまで繰り返すと最後はdist[a][b]だけで値を求められます。あとはここから元に戻っていけば最初の答えがわかるというわけです。このときメモ化をつかって、各段階でTSP(a, S-{a}, b)を記録しておきます。そうすることで、一度計算したパターンについては、また掘り下げて計算しなくて済むようになります。その分計算量が減ってより高速になるというわけです。今回の例ではノード数が少ないため重複する部分も少ないですが（今回の場合は最後の3パターンが一つ前の段階で重複しているだけ）、ノード数が増えるほど重複計算も増えるので効果がでてきます。

実際は以下のようにfor文や再帰処理を使うので、ここまで書き出す必要はなくなります。

TSP（再帰）：
まずはメモ化なしで再帰的な処理だけで組んでみます。

def TSP(a, S, b):
    if len(S) == 0:
        return dist[a][b]

    d_min = float('inf')
    for s in S - {a}:
        d = dist[a][s] + TSP(s, S - {a, s}, b)
        if d < d_min:
            d_min = d

    return d_min

%time print('Total Length:', TSP(0, S, 0))

関数内の最初の処理としては、len(S)で中継地点の数をカウントし、もしゼロ（空集合）であれば、dist[a][b]の値を返すようにしてあります（再帰ループの終了条件）。ここで一旦リターンされますが、内部的にまだ再帰処理が残っている場合は、残りの処理を続けるためにまたループし始めます。
Sが空集合でない場合は、続くforループで複数のパターンを処理させます。
関数TSP(a, S, b)の部分は、最初TSP(0, {0,1,2,3}, 0)であり、関数内でS-{a}={0,1,2,3}-{0}={1,2,3}にしてからforループさせています。次の処理となるforループ内では、S-{a, s}={0,1,2,3}-{0,s}で前回の出発点0と次の目的地sを差し引いて処理し最小値を取り出しています。
もっと完結に書けば、

def TSP(a, S, b):
    if len(S) == 0:
        return dist[a][b]
    return min([dist[a][s] + TSP(s, S - {a, s}, b) for s in S - {a}])

%time TSP(0, S, 0)

このくらいシンプルになりますが、手順が分かりにくいので最初のほうが理解しやすいかと。
再帰処理の場合、終了条件（今回の場合は中継地点Sが空集合になるまで）がなければ無限ループになってしまうので、まずは終了条件を何にするのか、そして内部処理では、終了条件に向かうために何かを変化させていけなければいけないのですが、その部分を何にするかというのが最初は分かりにくいという感じです。終了条件さえ決められれば、while文でループさせてもいいかもしれませんが、自己代入で再帰的に処理するほうが余計な変数など減っていいのかもしれません。

ちなみに、このコードでは最短経路の合計距離（厳密解）だけを求めています。再帰的とはいえ結局はすべてのパターンを計算しています。

このままだと計算量は(n - 1)!なので、10ノードだと1秒くらいかかってしまいます（遅い）。ということから、前回参考にしたメモ化（Memoization）を加えて計算が重複しないように工夫してみました。

TSP DP（メモ化）：

先ほどのコードにメモ化を追加して、一度計算した内容は遡って計算しないように記録しておき、その都度直接呼び出せるようにします。これで一応動的計画法になるようです。
例えば、TSP(1, {2, 3}, 0)の値は、TSP(2, {3}, 0)とTSP(3, {2}, 0)、TSP(3, {}, 0)やTSP(2, {}, 0)、そしてdist[3][0]とdist[2][0]まで遡りつつ比較最小値をとらなければいけませんが、一度最小値の計算結果を求めてあれば、その値をmemoから直接一回で呼び出せます。

memo = {}

def TSP_DP(a, S, b):
    if len(S) == 0:
        memo[(a, tuple(S - {a}), b)] = dist[a][b]
        return dist[a][b]

    d_min = float('inf')
    for s in S - {a}:
        if (s, tuple(S - {a, s}), b) not in memo:
            memo[(s, tuple(S - {a, s}), b)] = TSP_DP(s, S - {a, s}, b)
        d = dist[a][s] + memo[(s, tuple(S - {a, s}), b)]  
        if d < d_min:
            d_min = d

    return d_min

%time print('Total Length:', TSP_DP(0, S, 0))

今回はPythonの辞書機能を使用して、キーをそのままTSP()の引数にしています。辞書dict()では、任意の数や文字を直接キーに割り当てられますが、リスト[]やセット{}は登録できないようです。タプル()なら登録できるので、セット{}をタプル()に変えて登録しています。

結果的には、10ノードで1秒もかかっていた処理速度が34.5 msまで短縮できました。理論上の計算量は2ˆn * nˆ2になるそうです。(n - 1)!よりは高速ですが、ノード数15を超えたあたりからきつくなってくるので、DPとはいえ数十、数百ノードのTSPを処理させることは不可能です。

追記：
その後少しメモ化の部分を変更してみました（以下）。

memo = {}

def TSP_DP(a, S, b):
    S = S - {a}
    if len(S) == 0:
        memo[(a, tuple(S), b)] = dist[a][b]
        return dist[a][b]
    
    if (a, tuple(S), b) in memo:
        return memo[(a, tuple(S), b)]
    
    d_min = min([dist[a][s] +  TSP_DP(s, S - {s}, b) for s in S])
    memo[(a, tuple(S), b)] = d_min

    return d_min

%time print('Total Distance:', TSP_DP(0, S, 0))

辞書機能を使っているのは同じですが、最小値を得たあとに記憶させておくことにしました。この方が効率がいいかもしれません。今回はmin()関数を使い一行forループにし、それとS = S - {a}を最初の方で書いておきもう少しすっきりさせました。

メモ化するためのライブラリfunctools.lru_cacheもありますが、集合関数set()や{}あるいはlist[]を含めてしまうと適用できないようでbit DPにして整数で集合を表現した方がよさそうです（その方がより高速になる）。

経路の出力：

先ほどのコードでは最短距離（厳密解）しか求められないので、経路も求められるコードを追加してみました。経路順はメモ化した内容から割り出せます。最短候補のノードを一つずつ拾い上げてつなげていく感じです。

P = [0]  # starting from to "0"

for i in range(num - 2, -1, -1):
    d_min = np.inf
    p_min = None
    for m in memo:
        if len(m[1]) == i and  set(P) | {m[0]} | set(m[1]) == S:
            d = dist[P[-1]][m[0]] + memo[m]
            if d < d_min:
                d_min = d
                p_min = m[0]
    P.append(p_min)

print('Tour:', P + [0])
plot_path(P + [0])

S={0,1,2,3}の場合のmemoの中身は、

{(1, (), 0): 0.5761741880674179,
 (1, (2,), 0): 0.6102626177520184,
 (1, (2, 3), 0): 1.8033372646495835,
 (1, (3,), 0): 1.5373375235172968,
 (2, (), 0): 0.3706271783270686,
 (2, (1,), 0): 0.8158096274923677,
 (2, (1, 3), 0): 1.7769729629422466,
 (2, (3,), 0): 1.5637018252246337,
 (3, (), 0): 0.7446829688310912,
 (3, (1,), 0): 1.3688287427536237,
 (3, (1, 2), 0): 1.4029171724382241,
 (3, (2,), 0): 1.189646034720611}

となっており、最初の状態(0, {1,2,3}, 0)以降の状態が記録されています。実際for m in memoで呼び出されるのは、(1, (2, 3), 0)というキーの部分だけですが、そのキーを入力すれば対応する距離（計算結果）を一発で呼び出せます。{}がタプルに変換されて()になっています。ちなみに各キーに対応する値（dist[a][b]による距離の計算結果）はnp.random.random()で生成した値を使っているので（以下のGistに記載してあります）このページ最初に設定した値とは異なっています。
最初の出発地点は0、中継地点は{1,2,3}になるので、次の移動先は{1,2,3}のどれかになります。
つまり次のメモ化内容は、出発地が{1,2,3}のどれかで、中継地点の数が2個のものになります。
中継地点{1,2,3}の3通り分だけforループ処理して、さらに検索条件として中継地点の数が2個len(S)==2のものを抜き出します。そうすると、

(1, (2, 3), 0): 1.8033372646495835
(2, (1, 3), 0): 1.7769729629422466
(3, (1, 2), 0): 1.4029171724382241

のうちどれかになります。
TSPの経路として最小値を選択するには、これらの距離だけでは判断できません。これらに一つ前の距離（開始点0からの距離）を足した合計距離で最短のものを選択します。よって式は、

dist[0][1] + memo[(1, (2, 3), 0)]: 1.8033372646495835
dist[0][2] + memo[(2, (1, 3), 0)]: 1.7769729629422466
dist[0][3] + memo[(3, (1, 2), 0)]: 1.4029171724382241

となり、dist[i][j]は距離テーブルから実際の値が求まるので、それにmemo[m]から呼び出した距離を合計して最短経路を選びます。
この結果、最小値のものが、

dist[0][3] + memo[(3, (1, 2), 0)]: 1.4029171724382241

となったとします。
memo[(出発点, (中継地点), 最終到着点)]であるので、出発点3は前回のステップから見れば0からの移動先です。この値3を配列Pにappendしておき最終的な経路を記録していきます。
ということから、さらに次の移動先は中継地点の数が1個で、もうすでに訪れた3が含まれていないということが条件になります。今回の検索条件は

if len(m[1]) == i and  set(P) | {m[0]} | set(m[1]) == S:

としています。mはforループで抜き出したmemo内のキーです。(m[0], m[1], 0)は(出発地, 中継地点, 到着地)に対応しているので、len(m[1])は中継地点の数となります。set(P) | {m[0]} | set(m[1]) == Sの部分は、格納した経路地点と出発地と中継地点を合わせれば常に{0,1,2,3} == Sになるという条件です。 そうすると、次のステップは中継地点数が1個のもので、最終経路の一部としてPに格納されている3が含まれていないもので、出発地が1か2のもの、中継地点も1か2のものとなります。そうすると、

 (1, (2,), 0): 0.6102626177520184
 (2, (1,), 0): 0.8158096274923677

次の移動先はこれら二つに絞られ、さらにこれらにdist[i][j]を加えてから最小値の方を選択します。よって距離の合計式は、

 dist[3][1] + memo[(1, (2,), 0)]: 0.6102626177520184
 dist[3][2] + memo[(2, (1,), 0)]: 0.8158096274923677

となって、

dist[3][1] + memo[(1, (2,), 0)]: 0.6102626177520184

が選ばれたとします。 そうすると1が次の移動先となり、1を配列Pにappendしておきます。いまのところ0→3→1となります。つぎは消去法で2が移動先になりますが、一応先ほどと同じ条件で検索してみます。 次の条件は中継地点数が0個のもので、Pに格納されている[0,3,1]以外のものとなります。そうすると

(2, (), 0): 0.8158096274923677

が抜き出され、Pには2がappendされて[0,3,1,2]になります。 最後に0を加えて[0,3,1,2] + [0] = [0,3,1,2,0]にして最終的な巡回経路にします。

以下はn=10で試したものです。一応厳密解の経路が形成されました。

画像上部のTour:が経路順です。厳密解を見つけることができるのでいいのですが、TSPの規模が大きくなると無理なので、今回はTSPを材料に動的計画法やメモ化について学んでみたという感じです。他のことに応用できるか分かりませんが、それなりの収穫があったと思います。

TSPに関してはこれまでいくつかのヒューリスティックな方法を試しましたが、次は整数計画法を試してみたいです。

プログラミングコンテストチャレンジブック [第2版]　～問題解決のアルゴリズム活用力とコーディングテクニックを鍛える～

posted with amazlet at 19.05.19

秋葉拓哉 岩田陽一 北川宜稔
マイナビ
売り上げランキング: 10,015

Amazon.co.jpで詳細を見る

この↑「プログラミングコンテストチャレンジブック（通称：蟻本）」には様々なアルゴリズムが書いてあります。TSP　DPについても書いてありますが、やはりbit DPを使うと便利らしい。
以下のGitの最後にはmemoizationとbit DPを使ったコードも書いておきました。

続き：TSP DP（その２）　bit DP / 巡回セールスマン問題 / 動的計画法
関連：Traveling Salesman Problem:巡回セールスマン問題について（まとめ）