Notes of LMDP (1)  - 線形ベルマン方程式

※大人の事情で以前削除した、個人的なメモ書きの再掲です。参考になれば……

# 概要
MDPのベルマン方程式

```math
V\left ( s \right ) = \min_{a} \left \{ \ell \left ( s, a \right ) + \sum_{s^{\prime}} p \left ( s^{\prime} | s, a \right ) V \left (  s^{\prime} \right ) \right \} 
```

から、線形ベルマン方程式

```math
z(s) =  \exp \Bigl( - q\left ( s \right ) \Bigr)  \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) z(s^{\prime})
```

を導出



# マルコフ決定過程(MDP)
### MDPにおけるベルマン方程式
通常のMDPでは

- 状態集合$S$に属する状態$s$における状態価値関数$V(s)$
- 状態$s$で行動集合Aに属する行動$a$を実行することで得るコスト$\ell (s,a)$
- 状態$s$で行動$a$を実行して次状態$s^{\prime}$に遷移する状態遷移確率$p \left ( s^{\prime} | s, a \right ) $

を用いた、次のベルマン方程式を（価値反復法などによって）解きます。

```math
V\left ( s \right ) = \min_{a} \left \{ \ell \left ( s, a \right ) + \sum_{s^{\prime}} p \left ( s^{\prime} | s, a \right ) V \left (  s^{\prime} \right ) \right \} 
```
### min邪魔問題
このベルマン方程式でV(s)を計算するためには、$\min_a$、すなわち、｛｝の中身が最小となるような行動$a$を探して計算することを繰り返し行うしかない。
普通見かけるベルマン方程式は$\max$が使われていると思いますが、今回は報酬$r$の代わりにコスト$\ell$で表現していて、$\ell$は$r$の符号を反転させたもの、すなわち$\ell(s,a) = - r(s,a)$です。コストはなるべく小さくしたいので、$\min$。

## 線形可解マルコフ決定過程(LMDP)
**線形可解マルコフ決定過程(** Linearly-solvable MDP; **LMDP)** とは、エマヌエル・トドロフ（アメリカ）によって提案されたマルコフ決定過程の亜種です。別名を**カルバック・ライブラー制御(KL制御)**ともいいます。以下、主に[トドロフの論文(PDF)](https://homes.cs.washington.edu/~todorov/papers/TodorovNIPS06.pdf)に従う形で解説します。

### 行動aから制御uへ
まず、行動$a$を、制御$\mathbf{u}$というベクトルで置き換えます。行動$a$が「上・下・左・右」のような抽象的な概念だったのに対し、このベクトルは、$\mathbf{u} \in \mathbb{R}^{|S|}$、状態sの数の大きさを持つ、実数値のベクトルです。状態$s$における$\mathbf{u}$である$\mathbf{u}_s$を

```math
\mathbf{u}_s = u(s^{\prime}|s)
```

として定義します。$\mathbf{u}_s$は全状態$s$の数だけあるベクトルです。例えば、$3 \times 3$のグリッドワールドなら9個の要素を持ちます。

### 2つの仮定
LMDPでは、2つの大きな仮定を置きます。

#### 1.状態遷移確率の定義

MDPで、

- 状態遷移確率$p \left ( s^{\prime} | s, a \right ) $

を定義しましたが、LMDPでは「行動$a$」$\rightarrow$「状態$s$についての制御$\mathbf{u}_s$」に変えた、

- 状態遷移確率$p ( s^{\prime} | s, \mathbf{u}_s  ) $

を、次の式で定義します。

```math
p\left ( s^{\prime} | s, \mathbf{u}_s \right ) = \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right )
```

ここで、$\bar{p}$は、「uncontrolled probability」（非制御遷移確率）とか「passive dynamics」（受動ダイナミクス）とか呼ばれる**「制御をしないときの、もともとの状態遷移確率」**です。LMDPでは、状態遷移確率を、「もともとの遷移確率が、制御を実行することで変化したもの」と捉えるのです。

もちろん、もともと遷移しない組み合わせ、すなわち$ \bar{p}\left ( s^{\prime} | s \right )=0$ならば$p\left ( s^{\prime} | s, \mathbf{u}_s \right ) =0$です。
#### 2. コスト関数の定義
LMDPでは、コスト関数$\ell(s,\mathbf{u}_s)$を、次の式で定義します。

```math
\ell \left ( s, \mathbf{u}_s \right ) = q\left ( s \right ) + D_{\rm{KL}} \Bigl (  p\left ( s^{\prime} | s, \mathbf{u}_s \right ) \Bigl|\Bigr| \bar{p}\left ( s^{\prime} | s \right ) \Bigr )
```

第一項の$q(s)$は、制御によらず、ある状態に到達した時点で受け取る状態そのもののコストで、「**状態依存のコスト**」といいます。
第二項は、制御の入った遷移確率$p$と、制御が入っていないときの遷移確率$\bar{p}$の[カルバック・ライブラーダイバージェンス](https://ja.wikipedia.org/wiki/%E3%82%AB%E3%83%AB%E3%83%90%E3%83%83%E3%82%AF%E3%83%BB%E3%83%A9%E3%82%A4%E3%83%96%E3%83%A9%E3%83%BC%E6%83%85%E5%A0%B1%E9%87%8F)です。つまり、$p$と$\bar{p}$の確率分布の差を表しており、制御によって大きく変えられるほど、この項の数値は大きくなります。無茶な遷移をしようとするほど大きくなるので、そういう遷移を選択しないようにする効果があるようです。
カルバック・ライブラーダイバージェンスの定義

```math
D _ { \mathrm { KL } } ( P \| Q ) = \sum _ { i } P ( i ) \log \frac { P ( i ) } { Q ( i ) }
```

から

```math
\ell \left ( s, \mathbf{u}_s \right )  = q\left ( s \right ) + \sum_{s^{\prime}: \bar{p}(s^{\prime}|s) \neq 0} p \left (  s^{\prime} | s, \mathbf{u}_s \right ) \log  \frac{p \left (  s^{\prime} | s, \mathbf{u}_s \right )}{\bar{p}(s^{\prime}|s)}
```
となります。

ここで、状態遷移確率の定義

```math
p\left ( s^{\prime} | s, \mathbf{u}_s \right ) = \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right )
```

から、式変形して両辺のlogをとると

```math
\mathbf{u}_s = \log \frac{ p\left ( s^{\prime} | s, \mathbf{u}_s \right ) }{ \bar{p}\left ( s^{\prime} | s \right ) }
```
となるので、これをコスト関数の定義の式に代入すると、

```math
\ell \left ( s, \mathbf{u}_s \right )  = q\left ( s \right ) + \sum_{s^{\prime}} p \left (  s^{\prime} | s, \mathbf{u}_s \right ) \mathbf{u}_s
```
と変形できます。

## 打倒min　～minが邪魔なら中身を予めminしてしまえ～

さて、冒頭のベルマン方程式

```math
V\left ( s \right ) = \min_{a} \left \{ \ell \left ( s, a \right ) + \sum_{s^{\prime}} p \left ( s^{\prime} | s, a \right ) V \left (  s^{\prime} \right ) \right \} 
```

を、行動$a$を制御$\mathbf{u}_s$にした上で、先程の定義

```math
\begin{align}
p\left ( s^{\prime} | s, \mathbf{u}_s \right ) &= \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right ) \\
\ell \left ( s, \mathbf{u}_s \right )  &= q\left ( s \right ) + \sum_{s^{\prime}} p \left (  s^{\prime} | s, \mathbf{u}_s \right ) \mathbf{u}_s
\end{align}
```

を代入すると、

```math
\begin{align}
V\left ( s \right ) &=& \min_{\mathbf{u}_s} \left \{ q\left ( s \right ) + \sum_{s^{\prime}} p \left (  s^{\prime} | s, \mathbf{u}_s \right ) \mathbf{u}_s + \sum_{s^{\prime}} p \left ( s^{\prime} | s, \mathbf{u}_s \right ) V \left (  s^{\prime} \right ) \right \} \nonumber \\
&=& \min_{\mathbf{u}_s} \left \{ q\left ( s \right ) + \sum_{s^{\prime}} p \left (  s^{\prime} | s, \mathbf{u}_s \right ) \Bigl( \mathbf{u}_s + V \left (  s^{\prime} \right ) \Bigr) \right \} \nonumber \\
&=& \min_{\mathbf{u}_s} \left \{ q\left ( s \right ) + \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right ) \Bigl( \mathbf{u}_s + V \left (  s^{\prime} \right ) \Bigr) \right \} 
\end{align}
```


のように変形できます。
ここで、例の邪魔な$\min$の中身を見てみると、$q(s)$は制御に関係ないのでスルーでき、後半の

```math
\sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right ) \Bigl( \mathbf{u}_s + V \left (  s^{\prime} \right ) \Bigr)
```

を、$u$について最小化していることになります。これを数学的に最小化してしまえばいいじゃないか！

### ラグランジュの未定乗数法による攻略
ここで、極値の候補を求めることができる「[ラグランジュの未定乗数法](https://ja.wikipedia.org/wiki/%E3%83%A9%E3%82%B0%E3%83%A9%E3%83%B3%E3%82%B8%E3%83%A5%E3%81%AE%E6%9C%AA%E5%AE%9A%E4%B9%97%E6%95%B0%E6%B3%95)」を用います。

まず、極値をみつけたい関数$f(x)$と制約条件$g(x)$で、ラグランジュ関数$L$をつくります。

```math
L(x, \lambda) = f(x) - \lambda g(x)
```

このとき、Lのxについての偏微分と、Lのλについての偏微分がともに「0」になるような、つまり

```math
\frac{\partial L}{\partial x}   = 0 \\
\frac{\partial L}{\partial \lambda} = 0
```

となる$x$が極値の候補なのです。

今回、$f(x)$にあたるのが

```math
\sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right ) \Bigl( \mathbf{u}_s + V \left (  s^{\prime} \right ) \Bigr)
```
です。

また、制約条件$g(x)$には、状態遷移確率$p$は確率なので、全部足すと1になること、つまり

```math
\sum_{s^{\prime}} p\left ( s^{\prime} | s, \mathbf{u}_s \right ) =1
```

を、状態遷移確率の定義から変形して

```math
\sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right )  = 1 
```

として用います（式に代入するときは=0の形にします）。

こうして、次のラグランジュ関数を作り出すことができます。

```math
\mathcal{L} \Bigl( \mathbf{u}_s, \lambda_s \Bigr)= \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right ) \Bigl( \mathbf{u}_s + V \left (  s^{\prime} \right ) \Bigr) + \lambda_s \Bigl( \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right ) - 1 \Bigr)
```

これを解きましょう！

## ラグランジュ関数の偏微分を求めよう

#### uについての偏微分

```math
\frac{\partial \mathcal{L}}{\partial \mathbf{u}_s} = \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right ) \Bigl( \mathbf{u}_s + V \left (  s^{\prime} \right ) \Bigr)
```

$\bar{p} (s^{\prime}|s) \neq 0$のとき，$\frac{\partial \mathcal{L}}{\partial \mathbf{u}_s}=0$を満たす唯一の解は

```math
\mathbf{u}_s + V(s^{\prime})+\lambda_s + 1 = 0
```

であり，ゆえに

```math
\mathbf{u}_s = - V(s^{\prime}) - \lambda_s - 1
```

です。
このときの$\mathbf{u}_s=\mathbf{u}^{*}_s$とします。

$\bar{p} (s^{\prime}|s) > 0$かつ

```math
\exp \left ( \mathbf{u}^{*}_s \right ) > 0
```

なので、二階偏導関数は

```math
\left. \frac{ {\partial}^2 \mathcal{L} } {\partial \mathbf{u}_s \partial \mathbf{u}_s} \right|_{\mathbf{u}_s = \mathbf{u}^{*}_s} = \bar{p}(s^{\prime}|s) \exp \left ( \mathbf{u}^{*}_s \right ) > 0
```

となります。つまり、$\mathbf{u}^{*}_s$において極小なのです！

#### λについての偏微分
次に、$\lambda_s$についての偏微分を求めましょう。

```math
\frac{\partial \mathcal{L}}{\partial \lambda_s}=\sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right ) - 1=0
```

なので、uについての偏微分の項目で求めた式を変形して

```math
\begin{align}
&\mathbf{u}_s + V(s^{\prime})+\lambda_s + 1 = 0 \\
&\mathbf{u}_s = - V(s^{\prime})-\lambda_s - 1
\end{align}
```

として、先程のλについての式に代入してあげると、

```math
\begin{align}
\sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left \{  - V(s^{\prime}) - \lambda_s - 1 \right \} - 1 &= 0 \nonumber \\
\sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left \{  - V(s^{\prime}) \right \} \exp (  - \lambda_s - 1  ) &= 1 \nonumber \\
\sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left \{  - V(s^{\prime}) \right \}  &= \exp ( \lambda_s + 1 ) \\
\log \Bigl( \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left \{  - V(s^{\prime}) \right \} \Bigr)  &= \lambda_s + 1 \\
\lambda_s  &= \log \Bigl( \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left \{  - V(s^{\prime}) \right \} \Bigr)  - 1
\end{align}
```

と、$\lambda_s$が求まります。

これを、$\mathbf{u}^{*}_s$の式

```math
\mathbf{u}^{*}_s = - V(s^{\prime})-\lambda_s - 1
```

に代入してあげると、

```math
\mathbf{u}^{*}_s = - V(s^{\prime}) -  \log \Bigl( \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left \{  - V(s^{\prime}) \right \} \Bigr)
```

を得ます。

## 線形ベルマン方程式の導出

### min邪魔問題の突破

最小化したかったのは、

```math
\sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right ) \Bigl( \mathbf{u}_s + V \left (  s^{\prime} \right ) \Bigr)
```

でしたね。$\mathbf{u}^{*}_s$のとき、この項は

```math
\sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}^{*}_s \right ) \Bigl( \mathbf{u}^{*}_s + V \left (  s^{\prime} \right ) \Bigr) \nonumber
```

となります。
このとき、$\mathbf{u}^{*}_s$で極小なので、例のベルマン方程式

```math
V\left ( s \right ) =\min_{\mathbf{u}^{*}_s} \left \{ q\left ( s \right ) + \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}^{*}_s \right ) \Bigl( \mathbf{u}^{*}_s + V \left (  s^{\prime} \right ) \Bigr) \right \} 
```

のうち、$\min$の中身を考えると、

- 第1項はそもそも「状態依存のコスト」なので$\mathbf{u}^{*}_s$に関係ない
- 第2項はラグランジュの未定乗数法で極小値であるとわかった

ので、$\mathbf{u}^{*}_s$のとき最小になります。

**つまりminを外せます！！**

```math
V(s) =  q\left ( s \right ) + \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}^{*}_s \right ) \Bigl( \mathbf{u}^{*}_s + V \left (  s^{\prime} \right ) \Bigr)
```

min邪魔問題を突破しました。

### 線形ベルマン方程式への変形

さて、後少しです。ラグランジュの未定乗数法の制約条件に使った

```math
\begin{align}
\sum_{s^{\prime}} p\left ( s^{\prime} | s, \mathbf{u}_s \right ) &=1 \\
\sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left ( \mathbf{u}_s \right )  &= 1 
\end{align}
```

を先程の式に入れると

```math
V(s) =  q\left ( s \right ) +  \Bigl( \mathbf{u}^{*}_s + V \left (  s^{\prime} \right ) \Bigr)
```

ラグランジュの未定乗数法で求めた

```math
\mathbf{u}^{*}_s = - V(s^{\prime}) -  \log \Bigl( \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left \{  - V(s^{\prime}) \right \} \Bigr)
```

を代入して

```math
\begin{align}
V(s) &=  q\left ( s \right ) +  \Bigl( \mathbf{u}^{*}_s + V \left (  s^{\prime} \right ) \Bigr) \\
&=  q\left ( s \right ) +  \Bigl( - V(s^{\prime}) -  \log \Bigl( \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left \{  - V(s^{\prime}) \right \} \Bigr) + V \left (  s^{\prime} \right ) \Bigr) \\
&=  q\left ( s \right )  -  \log \Bigl( \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left \{  - V(s^{\prime}) \right \} \Bigr)
\end{align}
```

最後に、両辺の指数をとって

```math
\begin{align}
\exp \left ( - V(s) \right ) &=  \exp \Bigl( - q\left ( s \right ) +  \log \Bigl( \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left \{  - V(s^{\prime}) \right \} \Bigr) \Bigr) \\
\exp \left ( - V(s) \right ) &=  \exp \Bigl( - q\left ( s \right ) \Bigr)  \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) \exp \left \{  - V(s^{\prime}) \right \}
\end{align}
```

ここで、$\exp \left ( - V(s) \right )$を、$z(s)$に置き換えた


```math
z(s) =  \exp \Bigl( - q\left ( s \right ) \Bigr)  \sum_{s^{\prime}} \bar{p}\left ( s^{\prime} | s \right ) z(s^{\prime})
```


を、**線形ベルマン方程式**と呼びます。

$z(s)$を、**desirability function**といい、日本語だと「**好適度関数**」などと訳されます。これは、今回$V(s)$が**コスト（なるべく減らしたい）の期待値**になっているので、この符号を反転させている$z(s)$は、**増えれば増えるほどよい**値となるためです。

さて、この線形ベルマン方程式を行列の形に書き換えてみましょう。

- $\mathbf{G}$ : 状態数$\times$状態数の対角行列、対角成分$G(i,i)=\exp ( -q(i) )$
- $\mathbf{\bar{P}}$ : 状態数$\times$状態数の行列、$\bar{P}(i,j)=\bar{p}(j|i)$
- $\mathbf{Z}$ : 状態数の大きさのベクトル、$Z(i)=z(i)$

を定義したとき、線形ベルマン方程式は

```math
\mathbf{Z} = \mathbf{G} \mathbf{\bar{P}} \mathbf{Z}
```

という式になります。

つまり、状態依存のコスト$q(s)$と、もともとの遷移確率$\bar{p}(s^{\prime}|s)$が予めわかっているなら、**上の式で固有値$\mathbf{Z}$を求めれば良いだけです！**

ちなみに、状態価値が知りたければ、$z(s)$を求めれば、$V(s)=-\log(z(s))$で計算できます。

なお、自分がよく理解できていないだけかもしれませんが、トドロフの論文では、固有値計算に[べき乗法](https://ja.wikipedia.org/wiki/%E3%81%B9%E3%81%8D%E4%B9%97%E6%B3%95)を紹介しています。どうやら固有値zは絶対値が最も大きな固有値のようです。
べき乗法を使うやり方では、k=0,1,2,......のとき$\mathbf{Z}_0 = 1$（全要素が1）として、

```math
\mathbf{Z}_{k+1} = \mathbf{G} \mathbf{\bar{P}} \mathbf{Z}_{k}
```

を、$\mathbf{Z}_k$が収束するまで計算すればよいです。

### 最適状態遷移確率
$z(s)$が計算できれば、最適な状態遷移確率$p^{*}$は次の式で求められます。

```math
p^{*} \left(s^{\prime} | s \right )= \frac{ \bar{p} \left ( s^{\prime}|s \right ) z(s^{\prime}) }{ \sum_{s^{\prime}} \bar{p} \left ( s^{\prime}|s \right ) z(s^{\prime}) }
```

MDPでいうところの最適方策$\pi^{*}$ですね。これで状態$s$における遷移先$s^{\prime}$を（確率的に）選択してやればよいわけです。結局$\mathbf{u}_s$は式変形以外に使わなくてよいのですね。

## 参考文献

原論文はこちらです。

・[Todorov, E. (2007). "Linearly-solvable Markov decision problems." In Advances in neural information processing systems (pp. 1369-1376).](http://papers.nips.cc/paper/3002-linearly-solvable-markov-decision-problems.pdf)

勉強する際に特に参考にしたもの。

・[松原崇充(2015)「確率最適制御の最近の動向一確率推論による解法」システム／制御／情報・Vol．59　No．10　pp．369−374](https://www.jstage.jst.go.jp/article/isciesci/59/10/59_KJ00010110285/_pdf)
・[内部英治(2016)「線形可解マルコフ決定過程を用いた順・逆強化学習」日本神経回路学会誌・23巻1号](https://www.jstage.jst.go.jp/article/jnns/23/1/23_2/_article/-char/ja)
・[鹿島久嗣(2007)「Linearly-solvable Markov decision problems / Emanuel Todorov (UCSD)」(ppt) 機械学習研究グループ T-PRIMAL](http://www.ms.k.u-tokyo.ac.jp/T-PRIMAL/OpenSeminar/NIPS2006/Kashima.ppt)