動的計画法と強化学習

# 強化学習の核心
[速習 強化学習 --基礎理論とアルゴリズム](https://medium.com/@zanker/finding-a-ruby-bug-with-gdb-56d6b321bc86)のまえがきに以下のような文章がある。
>本書では、我々が制御対象とするシステムが確率的であると仮定する。さらに、ここでは制御器がシステムの状態を推論する必要がないぐらい十分、詳細にシステムの状態を観測できることを仮定する。このような特徴をもった問題は、マルコフ決定過程（Markov decision process; MDP）の枠組みでうまく記述できる。このMDPを"解く"ための標準的なアプローチは動的計画法である。動的計画法は、良い制御器を見つける問題を良い価値関数を見つける問題に落とし込む手法である。しかしながら、MDPが極めて少ない状態数と行動数しかもたない単純な問題から離れると、動的計画法は実行不可能となる。ここで我々が議論する強化学習アルゴリズムは、大規模な問題では実行不可能となる動的計画法を実践的なアルゴリズムに落とし込み、大規模な問題に対して適用可能にするための方法といえる。

つまり、動的計画法で実質的に解けない問題があるわけだが、それをなんとかして解こうとするのが諸々の強化学習のアルゴリズムである。この記事では、
>動的計画法は、良い制御器を見つける問題を良い価値関数を見つける問題に落とし込む手法である。

これがどういうことなのか調べていく。

# 動的計画法とは
DP(Dynamic Programing)と略されることも多く、結構皆好きである。詳細に説明された記事も多く、以下の記事が分かりやすいと思う。

   + [プログラミングコンテストでの動的計画法](https://www.slideshare.net/iwiwi/ss-3578511)
   + [典型的な DP (動的計画法) のパターンを整理 Part 1 ～ ナップサック DP 編 ～](https://qiita.com/drken/items/a5e6fe22863b7992efdb)

一応説明しておくと、
>ある計算式に対して、一度計算した結果をメモリに記録しておき、同じ計算を繰り返し行うという無駄を避けつつ、それらを再利用することにより効率化を図ることは、プログラミングやアルゴリズムの設計を行う上で有効なアプローチとなります。その手法の一つが動的計画方です。[^1]

本質的なアイディアは、**同じ計算を二度しない**、だ。具体例を見ると分かりやすいと思う。

[^1]: プログラミングコンテスト攻略のためのアルゴリズムとデータ構造― p248

### 例1)　フィボナッチ数列 ###
1, 1, 2, 3, 5, 8, 13, 21, 34, ....
と続く数列である。f(n)をフィボナッチ数列とすると、n = 0,1の時、f(n) = 1　で、それ以外の時、f(n) = f(n - 1) + f(n - 2)である。
超ナイーブな実装は以下のようになるだろう。

```C++
int fib(int n) {
    if(n == 0 | n == 1) return 1;
    return fib(n - 1) + fib(n - 2);
}
```

この関数の計算量は O(((1 + sqrt(5)) / 2)^(n-1) ) となる（補足参照）。これではnが大きくなると、途端に計算できなくなる。上記の関数では、例えばfib(5)を計算するときに、fib(4)とfib(3)を計算する必要があるが、fib(4)では、別にfib(3)を計算することになる。つまり、一度計算した値をもう一度計算しているのである。そこで、一度計算した値を配列に保持することにする。そして、その値が必要なときには、計算ずみの値を返す。これをメモ化と呼ぶ。

```C++
/* メモ化で用いる配列 */
int dp[100] = {};

int fib(int n) {
    // 値をメモしつつ値を返す 
    if(n == 0 | n == 1) return dp[n] = 1;
    // 値がメモされていれば、その値を返す 
    if(dp[n] != 0) return dp[n];
    // 値をメモしつつ値を返す 
    return dp[n] = fib(n - 1) + fib(n - 2);
}
```

また、上の関数は再帰的に計算しているが、ボトムアップのアプローチもある。

```C++
int fib2(int n) {
    int prev = 1, now = 1;
    for(int i = 0; i < n - 1; ++i) {
        int next;
        next = now + prev;
        prev = now;
        now = next;
    }
    return now;
}
```

こちらのアプローチだと配列をメモリに持つ必要がない。ローカル変数をうまく使い回して計算している。


### 例２)  最長共通部分列###
最長共通部分列問題（Longest Common Subsequence problem: LCS）は、2つの与えられた列X = {$x _ {1}$, $x _ {2}$, ..., $x _ {m}$}とY = {$y _ {1}$, $y _ {2}$, ..., $y _ {n}$}の最長共通部分列を求める問題。ある列ZがXとYの両方の部分列であるとき、ZをXとYの共通部分列と言う。そして、共通部分列の中で最長のものを最長共通部分列と呼ぶ。例えば、X = {a, b, c, b, d, a, b}, Y = {b, d, c, a, b, a}の最長共通部分列は{b, c, a, b}となる。
以下では、与えられた2つの文字列X, Yに対して、最長共通部分列Zの長さを出力するプログラムを考える。
まず、{$x _ {1}$, $x _ {2}$, ..., $x _ {i}$}を$X _ {i}$、{$y _ {1}$, $y _ {2}$, ..., $y _ {j}$}を$Y _ {j}$と表す。
そして、動的計画法で使う配列を以下のように定義する。
`c[m + 1][n + 1]`:`c[i][j]`を$X _ {i}$と$Y _ {j}$のLSCの長さとする2次元配列
そうすると、`c[i][j]`の値はiとjの値に基づいて、以下のように計算可能。

```math
c[i][j] = \left\{
\begin{array}{ll}
0 & (i = 0 \hspace{1pt}　or \hspace{3pt}  j = 0) \\
c[i - 1][j - 1] + 1  & (i,j \gt 0 \hspace{3pt} and \hspace{3pt} x_{i} = y_{j}) \\
max(c[i][j - 1], c[i - 1][j]) & (i,j \gt 0 \hspace{3pt} and \hspace{3pt} x_{i} ≠ y_{j})
\end{array}
\right.
```

これを実装すると、以下のようになる。

```C++

static const int N = 100;

int lcs(string X, string Y) {
    int c[N + 1][N + 1];
    int m = X.size();
    int n = Y.size();
    int maxl = 0;
    // c[i][j] = 0 if i = 0 or j = 0
    for (int i = 0; i < m + 1; ++i) c[i][0] = 0;
    for (int j = 1; j < n + 1; ++j) c[0][j] = 0;

    for (int i = 1; i < m + 1; ++i) {
      for (int j = 1; j < n + 1; ++j) {
        if (X[i] == Y[j]) {
          c[i][j] = c[i - 1][j - 1] + 1;
        } else {
          c[i][j] = max(c[i - 1][j], c[i][j - 1]);
        }
          maxl = max(maxl, c[i][j]);
      }
    }
    return maxl;
}
```

計算量はO(MN)。

### 例3)  連鎖行列積###
n個の行列の連鎖$M _ {1}$, $M _ {2}$, $M _ {3}$, ..., $M _ {n}$が与えられたとき、スカラー乗数の回数が最小になるように積$M _ {1}$$M _ {2}$$M _ {3}$...$M _ {n}$の計算順序を決定する問題を連鎖行列積問題　(Matrix-Chain Multiplication problem)と言う。以下では、n個の行列について、行列$M _ {i}$の次元が与えられたとき、積$M _ {1}$$M _ {2}$$M _ {3}$...$M _ {n}$の計算に必要なスカラー乗数の最小の回数を求めるプログラムについて考える。

l×mの行列Aとm×nの行列Bの積はl×nの行列Cとなり、行列の積の定義よりl×m×n回の乗算が必要になる。今$M _ {i}$が$p _ {i-1}$行$p _ {i}$列の行列であるような、n個の行列のかけ算を考える。
例えばn = 4とすると。
4つの行列$M _ {1}$, $M _ {2}$, $M _ {3}$,$M _ {4}$があって、$M _ {1}$は$p _ {0}$行$p _ {1}$列、$M _ {2}$は$p _ {1}$行$p _ {2}$列、$M _ {3}$は$p _ {2}$行$p _ {3}$列となる。これらの積は結合法則により、どの順序からでも計算可能である。つまり
$$ (((M _ {1}M _ {2})M _ {3})M _ {4}) = (M _ {1}(M _ {2}(M _ {3}M _ {4})))$$
みたいな感じになる。そして、この掛け算の順序で積の回数が変わってくる。全ての可能な順番を調べてみる方法はO(n!)となる。そして、この問題はより小さな問題に分割可能なので、動的計画法を適用可能である。
まず、($M _ {i}$$M _ {i+1}$)を計算するための方法は1通りであり、$p _ {i-1}$×$p _ {i}$×$p _ {i+1}$の掛け算が必要になる。そして、これをコストとしてテーブルに記録していく。
次に、3つの行列の積($M _ {1}$$M _ {2}$$M _ {3}$)、...、($M _ {n-2}$$M _ {n-1}$$M _ {n}$)を計算するための最小の掛け算の回数を求める。($M _ {1}$$M _ {2}$$M _ {3}$)のコストは$min(((M _ {1}M _ {2})M _ {3}), (M _ {1}(M _ {2}M _ {3})))$とする。今これらの積($((M _ {1}M _ {2})M _ {3})$と$(M _ {1}(M _ {2}M _ {3}))$)のコストは以下のように計算可能。

$((M _ {1}M _ {2})M _ {3})$のコスト = ($M _ {1}M _ {2}$)のコスト + ($M _ {3}$)のコスト + $p _ {0}$×$p _ {1}$×$p _ {3}$

$(M _ {1}(M _ {2}M _ {3}))$のコスト = ($M _ {1}$)のコスト + ($M _ {2}M _ {3}$)のコスト + $p _ {0}$×$p _ {2}$×$p _ {3}$

ここで注目すべきなのは、($M _ {1}M _ {2}$)のコストや($M _ {2}M _ {3}$)のコストは計算する必要はなく、テーブルから参照可能だということ。ただし、($M _ {i}$)のコストは0とする。

同様に$M _ {1}$$M _ {2}$$M _ {3}$$M _ {4}$の最適解は、

$M _ {1}$のコスト + ($M _ {2}$$M _ {3}$$M _ {4}$)のコスト + $p _ {0}$×$p _ {1}$×$p _ {5}$

($M _ {1}$$M _ {2}$)のコスト + ($M _ {3}$$M _ {4}$)のコスト + $p _ {0}$×$p _ {2}$×$p _ {5}$

($M _ {1}$$M _ {2}$$M _ {3}$)のコスト + ($M _ {4}$)のコスト + $p _ {0}$×$p _ {3}$×$p _ {5}$

のうち、最小のコストになる。

具体的には次のような変数を準備する。

`m[n + 1][m + 1]`:`m[i][j]`を($M _ {i}$$M _ {i+1}$...$M _ {j}$)を計算するための最小の掛け算の回数とする2次元配列
`p[n + 1]`:$M _ {i}$が`p[i - 1]`×`p[i]`の行列となるような配列

これらの変数を用いて、`m[i][j]`は以下の式で得られる。

```math
m[i][j] = \left\{
\begin{array}{ll}
0 & (i = j ) \\
min _ {i≦k＜j}\hspace{3pt}\{m[i][k] + m[k+1][j] + (p[i-1]　p[k]p[j]) \} & (j \gt i )
\end{array}
\right.
```

実装の一例は以下のようになる。

```C++

static const int N = 100;

int mcm(int matrixNumber, int matrixInfo[]) {
    int m[N + 1][N + 1];
    
    for (int i = 1; i < matrixNumber + 1; ++i) m[i][i] = 0;
    
    for(int l = 2; l < matrixNumber + 1; ++l) {
        for (int i = 1; i < matrixNumber - l + 2; ++i) {
            int j = i + l - 1;
            // 適当な値を設定
            m[i][j] = (1 << 21);
            
            for (int k = i; k < j; ++k) {
                m[i][j] = min(m[i][j], m[i][k] + m[k + 1][j] + matrixInfo[i - 1] * matrixInfo[k] * matrixInfo[j]);
            }
        }
    }
    return m[1][matrixNumber];
}
```

計算量はO($n^3$)。


# 強化学習におけるDP
以下の文章は、Suttonの[強化学習](https://www.amazon.co.jp/%E5%BC%B7%E5%8C%96%E5%AD%A6%E7%BF%92-Richard-S-Sutton/dp/4627826613)の4章と[速習 強化学習 --基礎理論とアルゴリズム](https://medium.com/@zanker/finding-a-ruby-bug-with-gdb-56d6b321bc86)の第1章を参考にしている。
Suttonの[強化学習](https://www.amazon.co.jp/%E5%BC%B7%E5%8C%96%E5%AD%A6%E7%BF%92-Richard-S-Sutton/dp/4627826613)の4章を読むと、
>動的計画法は、良い制御器を見つける問題を良い価値関数を見つける問題に落とし込む手法である。


の意味がよく分かった。

良い制御器（以下ではエージェントと呼びます。その方がしっくりくるので）を見つけるための方法として、方策反復(policy iteration)と価値反復(value iteration)の2つがある。そして、それらをプログラムで実行するときにはいずれの手法でも、メモ化を行って価値関数を計算していく。すなわち、動的計画法によって、価値関数を効率的に計算していくのである。そして、それは良い方策へと繋がっていく。



### 方策反復 ###
ある方策$π _ {k}$がわかっている時、その方策$π _ {k}$に従う価値関数$V^{ π _ {k} }$を求める（方策評価と呼ばれるステップ）。そして、$V^{ π _ {k} }$に関して、グリーディーな方策を$π _ {k + 1}$として求めていく（方策改善と呼ばれるステップ）。これを繰り返していくことで最適な方策を得ることができる。そして、方策評価のステップでDPは使われている。価値関数を求めるときに、次のBellman方程式を更新規則として用いる。
$$ V_{k+1}\left( y\right) = \sum _{a}\pi \left( x,a\right) \sum _{y} P(x, a, y)[R(x, a, y) + \gamma V _ {k}(y) ]$$

つまり、この式によると$V_{l+1}(x)$を計算するためには$V_{l}(x)$が必要になり、$V_{l}(x)$を求めるには$V_{l - 1}(x)$が必要になる（以下ループ）。最新の価値関数を配列に記録しておくことで、巻き戻りをしなくて済む。つまりメモ化をしているのである。ということだと僕は理解した。

### 価値反復 ###
価値反復はもっとシンプルである。Bellman最適方程式を更新規則として用いて、価値関数の系列を得る。
$$ V_{k+1}\left( y\right) = max _{a}\sum _{y} P(x, a, y)[R(x, a, y) + \gamma V _ {k}(y) ]$$
これも方策反復と同様に動的計画法を用いて効率的に計算できる。



注意点としてはこれらの手法は大規模な問題（つまり、状態空間や行動空間が大きい）では、適用できない。$\sum _{y}$や$max _{a}$がもはや計算できなくなる。


# まとめ
動的計画法を使うと、良いagentを見つける問題を方策反復や価値反復を用いて、価値関数を見つける問題へとおとしこめる。

# 補足
フィボナッチ数列の計算量について軽くまとめる。
[この記事](http://www.aoni.waseda.jp/ichiji/2012/java-01/java-14-3.html)では、最終的な計算のオーダーは示されているが、計算過程について説明した記事は僕がネットを調べた限りだとない。そんなに難しくないので、ここでちょっとだけ説明する。fib(n)の計算量をC(n)とすると、[この記事](http://www.aoni.waseda.jp/ichiji/2012/java-01/java-14-3.html)を参考にして、
$$C(n) = C(n - 1) + C(n - 2) + 5$$ 
という漸化式が成り立つ。この式を特性方程式を解いて変形すると、
$$C\left( n\right) -\dfrac {1+\sqrt {5}}{2}C\left( n-1\right) =\dfrac {1-\sqrt {5}}{2}( C\left( n-1\right) -\dfrac {1+\sqrt {5}}{2}C\left( n-2\right)) + 5$$
これを変形して、
$$C\left( n\right) -\dfrac {1+\sqrt {5}}{2}C\left( n-1\right) - \dfrac {5}{2}(\sqrt {5} - 1) =\dfrac {1-\sqrt {5}}{2}[( C\left( n-1\right) -\dfrac {1+\sqrt {5}}{2}C\left( n-2\right)) - \dfrac {5}{2}(\sqrt {5} - 1) ]$$
よって、
$$C\left( n\right) -\dfrac {1+\sqrt {5}}{2}C\left( n-1\right) - \dfrac {5}{2}(\sqrt {5} - 1) = \left( \dfrac {1-\sqrt {5}}{2}\right) ^{n-1} × Const$$
nが十分大きいと、右辺は0に収束する。さらに、左辺の定数を無視すると、
$$C\left( n\right) =\dfrac {1+\sqrt {5}}{2}C\left( n-1\right)$$
これより、C(2) = 2　を使うと、
$$C\left( n\right) =2\left( \dfrac {1+\sqrt {5}}{2}\right) ^{n-1}$$