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