TSP DP（その２） bit DP / 巡回セールスマン問題 / 動的計画法

前回のTSP DP（動的計画法による巡回セールスマン問題）では、訪問先を表現するためにset()関数を用いましたが、今回はbit DPで試してみました。
bit演算させたほうが高速になるのと、メモ化ライブラリであるfunctools.lru_cacheも使えるようになります。

ビット演算でどこを訪れていて、どこを訪れていないか（訪問フラグ）は以下のように表すようです。
まず訪問先を0, 1, 2, 3とすれば4桁の二進数になります。Pythonでは左側に0bをつければ二進数となります。

S = 0b0001

0ビット目（一番右側）を1にしておいて、ここが出発点となります。
つぎに2に訪れたら（右から0番目、1番目、2番目）、

S = 0b0101

となり、すでに訪れた場所の状態がわかります。
次に1に訪れて、

S = 0b0111

つぎに3に訪れて、

S = 0b1111

最終的にすべてが1になれば完了というわけです。
ビット演算用の記号を使いますが、実際は整数としても認識するので、

S = 1 | 1 << 2 = 0b0101

こんな感じになります。

ということから二進数を用いて前回のTSPを書き直してみました。


functools.lru_cacheを使った場合：（Python 3.6、Jupyter Notebook使用）
メモ化ライブラリであるfunctools.lru_cacheを使うと以下のようになります。
ちなみにこれは合計最短経路のみの出力です。尚、TSP/DPのアルゴリズムについては前回の方が詳しく説明してあるかと。

from functools import lru_cache

@lru_cache(maxsize=None)
def TSP_lru(S, pos):        
    if S == all_visited:
        return dist[pos][0]
    return min([dist[pos][i] + TSP_lru(S | (1 << i), i) for i in range(num) if S & (1 << i) == 0])

# @lru_cache(maxsize=None)
# def TSP_lru(S, pos):  # same as the above
#     if S == all_visited:
#         return dist[pos][0]
#
#     d_min = float('inf')
#     for i in range(num):
#         if S & (1 << i) == 0:     # if not visited 
#             d = dist[pos][i] + TSP_lru(S | (1 << i), i)
#             if d < d_min:
#                 d_min = d
#   return d_min

# Solve by TSP with functools.lru_cache
%time print('Total Distance:', TSP_lru(1, 0))
print(TSP_lru.cache_info())  # show cache info
TSP_lru.cache_clear()        # clear cache otherwise same outputs as before

メモ化されていないファンクションの上に@マークのデコレータをつけるだけでメモ化されます。
maxsize=Noneでメモリ制限なしにしています。デコレータなしだとメモ化しないためかなり遅くなります。
ただし、前回のようなset関数やリストを引数にしている場合はこのlru_cacheが使えません。タプルを引数にするなら大丈夫だと思います。今回は訪問先の集合をビットで表したので使うことができるというわけです。おそらくビット演算にするだけでも、少しは高速になっているのではないでしょうか。
それから、一度runするとcacheに内容が保持されたままになるので、再度TSR_lru(1, 0)をrunすると、前回と同じ答えを出してしまいます。ノード数を変えて試してみたい場合は、毎回cache_clear()する必要があります。何度も同じ答えを呼び出すには高速ですが、キャッシュクリアしないかぎりは同じ答えを出し続けるので、他のプログラムで使うときには注意したほうがいいと思います。


memoize関数を使った場合：
自前で書いたmemoize関数をデコレータにすることもできるようです（こちらを参考）。
合計最短経路のみの出力。

def memoize(f):
    cache = {}
    def func(*args):
        if not args in cache:
            cache[args] = f(*args)
        return cache[args]
    return func

@memoize
def TSP_memo(S, pos):  # re-run this function to clear cache otherwise same outputs as before
    if S == all_visited:
        return dist[pos][0]
    return min([dist[pos][i] + TSP_memo(S | (1 << i), i) for i in range(num) if S & (1 << i) == 0])

# Solve by TSP with memoization decorator
%time print('Total Distance:', TSP_memo(1, 0))

これも先ほどのlru_cacheと同じようなものです。
この場合もキャッシュクリアしない限りは同じ答えを出すので、TSP_memo(1, 0)をrunさせる場合は、def TSP_memo():も一緒に（jupyter上のこのセルごと）runさせたほうがいいです（このキャッシュクリアの方法がわからない）。


リストを使ったDPメモ化：
ビット化した以外はほぼ前回と同じです。
最後にメモ化した内容から経路順を割り出します。

def TSP_BDP(S, pos):
    # if visited before
    if dp[S][pos] >= 0:
        return dp[S][pos]
    
    # if all visited  
    if S == all_visited:         
        dp[S][pos] = dist[pos][0]
        return dp[S][pos]
    
    # if not visited yet
    dp[S][pos] = min([dist[pos][i] + TSP_BDP(S | (1 << i), i) for i in range(num) if S & (1 << i) == 0])
    
    return dp[S][pos]

# Solve TSP by bit DP (15 nodes: 0.75 sec, 18 nodes: 8.42 sec, 20 nodes: 44 sec)
dp = [[-1] * num for _ in range(1 << num)]    # initialize memoization
%time print('Total Distance:', TSP_BDP(1, 0))

基本的に、

dp[現在訪問状態][現在地] = dist[現在地][移動先] + TSP_BDP(次の訪問状態, 移動先)

を毎回記録/計算しています。forループのi（移動先）は、未訪問の地点があるかどうか確かめてから、あればそこへ移動する場合の距離を計算しmin()によって最短経路を選択します。
具体的な値を入れてみると、

dp[0b1011][3] = dist[3][2] + TSP_BDP(0b1011 | (1 << 2), 2)

Sについて2ビット目が0となっているので（右から数えて0,1,2番目）、まだ訪れていない地点は2ビット目だけということになり、forループでi=2のときこの処理がなされます。
そうすると、

dp[0b1011][3] = dist[3][2] + TSP_BDP(0b1111, 2)

右辺においてS=0b1111（全て訪問済み）となり、あとは現在地3から2までの距離と残りTSP_BDP(0b1111, 2)による距離を足して左辺のdp[0b1011][3] へ記録することになります。
しかし、dp[0b1011][3] へ記録する前に、TSP_BDP(0b1111, 2)を計算しなければいけないので、TSP_BDP(0b1111, 2)は次の再帰ループに突入してif S == all_visited:の条件分岐により、

dp[0b1111][2] = dist[2][0]

になり距離テーブルdist[2][0]から値を得てdp[0b1111][2]に記録します。このとき再帰ループで一つ先に移動しているので現在地pos=2、移動先は0になっています。距離が記録されたのちdp[0b1111][2]（= dist[2][0]）がリターンされたのでTSP_BDP(0b1111, 2)に戻ります。
この結果、

dp[0b1011][3] = dist[3][2] + TSP_BDP(0b1111, 2)

この式の右辺側にあるTSP_BDP(0b1111, 2)は、

TSP_BDP(0b1111, 2) = dist[2][0]

とわかったので、それを代入すると

dp[0b1011][3] = dist[3][2] + dist[2][0]

つまり、

未訪問2で現在地3 = 3から2までの距離 + 2から0（ゴール）までの距離

となってゴールから二歩手前までの距離が求まります。

訪問フラグと距離の関係（4ビットの場合、posはその時の現在地）：
dp[ob1111][pos]（4ビット1の場合）ゴールから一歩手前までの距離
dp[0b1011][pos]（3ビット1の場合）ゴールから二歩手前までの距離
dp[0b0011][pos]（2ビット1の場合）ゴールから三歩手前までの距離
dp[0b0001][pos]（1ビット1の場合）ゴールからスタートまでの距離

ということから、この仕組みを理解して経路の生成を以下で行いました。


経路の生成：
上で記録した二次元配列dpを元に経路の順番を生成します。

def gen_dp_tour():
    tour = [0]
    S = 1
    while len(tour) < num:
        d_min = np.inf
        pos = None
        for i in range(num):
            if S & (1 << i) == 0:
                d = dist[tour[-1]][i] + dp[S | (1 << i)][i]
                if d < d_min:
                    d_min = d
                    pos = i
        S |= (1 << pos)
        tour.append(pos)
    return tour + [0]

# Generate Tour
Tour = gen_dp_tour()
print('Tour:', Tour)
plot_path(Tour)

基本的には最短経路を求めるコードと似ています。前のコードのforループ内のTSP_BDP()の部分をdp[][]に置き換えて距離を計算させるかわりに直接メモリから呼び出します。再帰ループの代わりにwhileループにしています。最初出発点にいる状態から全てを訪問するまで、min関数で得た最短距離と同時に現在地も記録しておき、それをtour配列に入れていきます。
先ほどの訪問フラグと距離の関係から、dp[0b0001]にはスタートからゴールまでの合計最短距離が記録されています。これを一つずつ分解していく感じです。

dp[0b0001][pos] = dist[pos][i] + dp[S | (s << i)][i]

というのは、

dp[現在訪問状態][現在地] = dist[現在地][移動先] + dp[次の訪問状態][移動先]

なので、初期状態は0ビット目にいるため0b0001になり現在地pos=0になります。次の移動先iは複数の可能性があるのでforループで各パターンを確かめます。min()関数を使ってdist[][] + dp[][]の合計距離で最短の地点iを次の移動先にしますが、この時最短のiを経路用配列tourにappendしていきます。次のループでは再度訪問フラグをチェックして、複数ある未訪問の地点iについてまたdist[][] + dp[][]の合計距離を調べます。あとはそれの繰り返しとなります。尚、現在地posはその都度追加していった経路用配列tourの最後の値tour[-1]を参照すればいいことになります。初期状態では0から出発することになるので、tour=[0]にしています。

次へ：TSP LP (Linear Programming)：線形計画法による巡回セールスマン問題 / Subtour Elimination 部分巡回路除去

関連：Traveling Salesman Problem:巡回セールスマン問題について（まとめ）

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:巡回セールスマン問題について（まとめ）