動的計画法：DP（Dynamic Programming）/ メモ化再帰

今回は動的計画法（DP）と再帰的アルゴリズムについてです。これもまたTSP（巡回セールスマン問題）から派生したアルゴリズムなので勉強用に記録しておきます。要は演算を効率化して、ある程度の規模の組み合わせ問題においても計算可能にしていく工夫という感じです。

これまで、

・貪欲法
・2-OPT法
・挿入法
・凸包（ギフト包装法/グラハムスキャン法）
・最小全域木（MST：クラスカル法/プリム法）
・クリスフィードアルゴリズム（グラフ理論：Minimum Weight Perfect Matching）

など様々なアルゴリズムがでてきましたが、さらには、

・動的計画法（DP）
・メモ化再帰
・線形計画問題（LP）
・ラグランジュ緩和
・整数線形計画問題（ILP）

というアルゴリズムにも関係していくようで、TSP以外にも様々な組み合わせ問題や最適化問題に応用できるためもう少し勉強してみようかと。前回登場したMinimum Weight Perfect Matchingを見つける演算の際にも使えるのかもしれません。
特に線形計画問題（LP）はラグランジュ緩和に発展できるようなので、以前勉強したサポートベクターマシン（SVM）にも通じるし、機械学習における最適化を数理的に理解する上でも役立ちそうです。
TSPをきっかけにいろんなアルゴリズムに派生しましたが、今回はこの中でも動的計画法（DP）とメモ化再帰について試してみることにしました（まだまだ覚えることはたくさんありそう）。


再帰関数：(Python 3.6、 Jupyter Notebook使用)
簡単な例として階乗のアルゴリズムから。

def fact(n):
    if n == 0:
        return 1
    return n * fact(n-1)

関数fact(n)の内部で自身の関数を呼び出す方法。n=5なら5!=5*4*3*2*1=120。

フィボナッチ数（0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55...）の場合は、
fib_n = fib_n-1 + fib_n-2なので、

def fib(n):
    if n < 2:
        return n
    else:
        return fib(n-1) + fib(n-2)

これも関数fib(n)内で自身の関数fib(n-1)とfib(n-2)を呼び出す方法。しかしながらこの方法の場合はn=30あたりから計算に時間がかかってきて、n=35だと2.34s、n=40で26.7sもかかってしまいます。というのは、fib(n)を求めるにはfib(n-1)を求めなければならず、fib(n-1)を求めるにはfib(n-2)を求めなければいけなくなり、毎回n=1まで遡って計算しなければいけなくなり、その都度計算が重複してしまうためのようです。
nが大きい値になると現実的な時間内では計算不可能となるため、高速化するには以下の方法で。


末尾再帰（高速化）：
フィボナッチ数をサンプルとすれば、一度計算した値（一つ前と二つ前の値）を記憶させておく末尾再帰という方法。

def fib(n, n1=1, n2=0):
    if (n < 2): 
        return n1
    else:
        return fib(n-1, n1+n2, n1)

こうすることで、n=40で26.7sかかっていたのが14.1usになり、n=1000でも982usで済みます。nが大きい値であっても計算可能。
ループを使うなら以下。

def fib(n):
    n1, n2 = 1, 0
    for _ in range(2, n):
        n1, n2 = n1 + n2, n1
    return n1 + n2

試してみるとこちらのほうが速い。n=1000で57us。


メモ化（Memoization）/ 動的計画法：
もうひとつメモ化：Memoization（Memorizationではなく）という方法で一度計算した内容を重複しないように記録しておく方法。動的計画法（ボトムアップ）。再帰的ではない。

def fib(n):
    memo = [0] * (n+1)
    for i in range(n+1):
        if i < 2:
            memo[i] = i
        else:
            memo[i] = memo[i-1] + memo[i-2]
    return memo[n]

この場合、n=1000で162usなので処理は高速な方。しかし事前にn個数だけリストに変数を格納する分メモリを消費。
似たような方法で、リストの代わりにディクショナリを使った方法は以下。トップダウンで再帰的な動的計画法。

memo = {}
def fib(n):
    if n in memo:
        return memo[n]
    if n < 2:
        val = n
    else:
        val = fib(n-1) + fib(n-2)
    memo[n] = val
    return val

ディクショナリの場合はその都度必要に応じてキーと値を格納できるので、事前にn個分の値を格納するリストを用意しなくてもよい。n=1000で594us程度。

それなら普通にループを使って処理させたら？（以下）

def fib(n):
    memo = []
    for i in range(n+1):
        if i < 2:
            memo.append(i)
        else:
            memo.append(memo[i-1] + memo[i-2])
    return memo[-1]

リストを使って普通にappendしていくだけなのでわかりやすい。n=1000で190usなので結構高速。ただし、これもリストにn個分の値を格納するので、それだけメモリを使う。


ライブラリ/デコレータの使用：
またライブラリを使うことでメモ化することも可能。

from functools import lru_cache

@lru_cache()
def fib(n):
    if n < 2:
        return n
    else:
        return fib(n-1) + fib(n-2)

%time [fib(i) for i in range(987, 1001)]

functoolsライブラリをインポートし、最初に書いた演算が遅くなってしまう再帰関数にデコレータである@lru_cacheを追加すればメモ化が可能。しかしこの場合n>988だとエラー（スタックオーバーフロー）になってしまうため、forループで987から1000まで計算させれば大丈夫。
尚、%timeはJupyter Notebookの時間計測用マジックコマンド。

自前でデコレータを用意する場合は以下。

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

@memoize
def fib_mem(n):
    if n <= 2:
        return 1
    else:
        return fib_mem(n-2) + fib_mem(n-1)

この場合、n=1000で6.44us。これはかなり高速だし汎用的なので便利かもしれない。

なんとなく動的計画法や再帰関数の基本を理解してみましたが、次回はこれをTSP（巡回セールスマン問題）に適用して厳密解を求めてみたいと考えています。

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

最小全域木（MST: minimum spanning tree)、Christofides algorithmによるTSP初期経路生成

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に含まれるかどうかがわかります。
ということで、この一連の流れをコーディングしてみました。実際は計算量を軽減するための工夫も必要ですが、今回は結果だけを求めています。

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

最初に全ての辺の組み合わせを計算しなければいけないので、ノード数が増えるとかなりきつくなりそうです。この部分に工夫が必要そうです。今回はベタに計算しています。

ノード50個の場合。
このように最短の辺で構成された枝分かれのグラフができあがります。
すべてのノードへ訪問可能にするために必要な経路の合計が最短ということです。
この段階では、
[[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()という関数を別途用意してあります。

グラフ・ネットワークアルゴリズムの基礎: 数理とCプログラム

posted with amazlet at 19.02.08

浅野 孝夫
近代科学社
売り上げランキング: 444,789

Amazon.co.jpで詳細を見る

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を参考にするといいようです。

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 newtour

MSTで生成した経路をT、最適マッチングによって追加された経路をMとしています。なお引数にMを渡さなければ、Double Treeアルゴリズムで経路生成します。その場合はTを2倍にした経路が初期状態となります。
基本的には、各辺は[a,b]というような一組のノード番号によってリストに入っており、次につながる辺は[b,c]もしくは[c,b]となるため前後を逆にした二通りを検索し、つながる辺が見つかれば追加していきます。
しかし、そのまま各辺を追加してくと、開始点によっては途中で道が途絶えてしまうことがあります。これは、開始点と現在地とで局所的ループが発生して、それ以外の経路へ移動できない状態になっているということです。局所的ループによって行き止まりになったら、現在加えた辺を経路リストの最後から先頭に移動して現在地を一つ手前に戻します。それでも接続可能な辺を見つけられなければ、接続可能な辺が見つかるまでこの手順を繰り返します。
こうすることで、どの地点から開始しても現在地と開始点をずらしながら通り過ぎてしまった分岐点に遡ることが可能となり、見逃していた辺（経路）を追加していくことができます。
そして順路ができあがったらノード0を起点に並べ替えて、最後に再度ノード0を付け足して巡回路とします。

 

ノード数20の場合。左（上）がDouble Tree。右（下）がChristofides。
最終的に重複するノードをスキップさせるためか経路が部分的に交差してしまうことがあります。二つを何回か比較すると、Christofidesのほうがより安定した経路を形成します。
そして最適化させると以下のようになります。交差している部分は解消されほぼ厳密解に近い経路になります。

比較結果（追記）：

その後、前回の凸包・挿入法と今回のMSTによる初期経路生成について、TSPライブラリのPyconcordeの結果を元に100回分のスコアを比較しました。

Pyconcorde  ：100.0000 %

凸包・挿入法：103.3175 %

Double Tree ：105.5888 %

Christofides：102.7978 %

それぞれのアルゴリズムで経路生成したあとに最適化アルゴリズムを通して得られたスコアです。Pyconcordeのスコアを100%とした場合の経路長の率で、%が低いほど短距離でいい結果ということです。今回は簡単な最適化しか行っていないためもう少し改善することはできましたが、Pyconcordeより短い距離を出すことはほぼ無理。何回かに一回はPyconcordeと同じスコアになるときはあります。

凸包・挿入法は一見よさそうだったのですが、結果的にはChristofidesのほうがわずかながら優れています。凸包とはいえ仕組み的には貪欲法なので、やはりMSTを構築したほうが結果がよくなるという感じでしょうか。MSTは最短の辺で全体をつなぎとめるので、全体的な構成や骨格が最小で済み、あとは局所的な最適化さえすればいいという感じで理にかなっているのかもしれません。

関連：

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

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

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

posted with amazlet at 19.02.08

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

Amazon.co.jpで詳細を見る