Notes of LMDP (3)  - LMDPの逆強化学習(1) OptV, OptQ

論文 "[Inverse Optimal Control with Linearly-Solvable MDPs](https://homes.cs.washington.edu/~todorov/papers/DvijothamICML10.pdf)"のメモです。

##線形可解マルコフ決定過程（LMDP）
詳しくは[Notes of LMDP (1) - 線形ベルマン方程式](https://qiita.com/deepirl_learner/items/e05f9a73d87df34389e2)を参照。

この論文で出てくるところをおさらい。

- **状態依存のコスト関数** $q(s)$ : 強化学習でおなじみ報酬$R(s)$の符号反転版で、小さいほどよい。状態にのみ依存する(state-dependent)。
- **制御遷移確率** $\pi (s^{\prime}|s)$ : 制御器やエージェントは、従来のMDPの$s^{\prime} \sim \pi(s,a)$の代わりに次状態への遷移確率を直接指定できる（ $s^{\prime} \sim \pi (\cdot |s)$ ）。
- **非制御遷移確率** $\bar{p}(\cdot|s)$ : システム固有の状態遷移確率。言うなればマルコフ過程の遷移確率みたいなもの。論文では受動ダイナミクス（passive dynamics）と呼ばれる。

これらにより、コスト関数 $\ell(s,\pi(\cdot|s) )$ は

```math

\ell(s,\pi(\cdot|s)) = q(s) + D_{\operatorname{KL}} \left \{ \pi(\cdot|s) || \bar{p}(\cdot | s)\right \}

```

となる。$D_{\operatorname{KL}}$はKLダイバージェンス。

また、状態価値関数$V(s)$、**好適度関数** $z(s) = \exp( -V(s) ) $を定義するとき、最適方策$\pi^{*}$は

```math

\pi^{*}(s^{\prime}|s) = \frac{ \bar{p}(s^{\prime}|s) z(s^{\prime}) }{\sum_k \bar{p}(k|s) z(k) }

```

となる。また、$z$についての線形ベルマン方程式

```math

\lambda z(s) = \exp (-q(s)) \sum_k \bar{p}(k|s) z(k)

```

が成り立つ。$\lambda$は"first-exit"（ゴールに着いたら終了=episodic task）の問題では$\lambda=1$、無限区間平均コスト問題では$\lambda=\exp(-平均コスト)$。

# LMDPでの逆強化学習(IRL)
LMDPでのIRLの初の論文（多分）が
K. Dvijotham, E. Todorov (2010) "[Inverse Optimal Control with Linearly-Solvable MDPs](https://homes.cs.washington.edu/~todorov/papers/DvijothamICML10.pdf)"
である。論文の名前としては「逆最適制御」(IOC)だが、これはIRLに同じ。

## 使うデータセット
だいたいのIRLアルゴリズムはエキスパートの軌跡一本づつを所与とするが、紹介するアルゴリズム群ではエキスパート、すなわち最適方策$\pi^{*}$からサンプリングされた$N$個の状態遷移のペア

```math

\{ s_{(n)}, s^{\prime}_{(n)} \}_{n=1,\cdots N} \\
s^{\prime}_{(n)} \sim \pi^{*}(\cdot | s_{(n)})

```

を利用する。つまり、軌跡の長さなどに影響されないし、断片的な軌跡でも使える。

## OptV

### zの最尤推定
zをエキスパートデータから最尤推定することを考える。なお、$\bar{p}$は所与とする。$z$についての損失関数は、次の負の対数尤度で与えられる。

```math

L [ z(\cdot) ]=-\sum_n \ln z(s^{\prime}_{(n)}) + \sum_n \ln \sum_{s^{\prime}}\bar{p}(s^{\prime}|s)z(s^{\prime})

```

これを$z$について最小化すればよいのだが、この$L$は非凸であるために、一応実験してみたものの「局所解に陥りがちな上に遅かった」という。

### Vについての書き直し＋訪問カウントの導入
そこで、この$L$を$V$についての関数に書き直してやると、log-sum-expと線形関数の和になるために凸になる。さらに、「データポイントの数が状態数を上回った時」（状態の次元数よりデータ数が多いという意味？）により高速に$L$を計算するテクニックとして、「訪問カウント」を式に導入する。これはデータセット中の特定の状態の数をカウントしたもので、現状態$s$についてと次状態$s^{\prime}$についての2つを定義する。


```math
a(s^{\prime}) \cdots s^{\prime}_{(n)} = s^{\prime} の回数 \\
b(s) \cdots s_{(n)} = s の回数
```

これらにより、$L$を次のように書き直すことができる。

```math
L[V(\cdot)] = \sum_s^{\prime} a(s^{\prime})V(s^{\prime}) + \sum_s b(s) \ln \sum_{s^{\prime}}\bar{p}(s^{\prime}|s)\exp\{ 
 -V(s^{\prime})\}
```

この損失を最小化するような$V$を推定すればよい。このアルゴリズムを**OptV**と呼ぶ。著者らはこの式のヘッセ行列を解析的に求めたうえで、バックトラック法の直線探索によるニュートン法で実装したそうである。

ここで疑問に思ったのが、「$V$は求まるけどどうやって$q$求めるんだ」というところ。論文を見てみると、

> Once we have an estimate $\hat{V}$, we can obtain $\hat{z} = \exp(−\hat{V})$ , $\hat{\pi}^{*}$ from (2), and $\hat{q}$ from (4).
※表記を一部変更

とある。式(4)とは上述した線形ベルマン方程式である。どうやら、線形ベルマン方程式を$q$について変形すると

```math

z(s) = \exp (-q(s)) \sum_k \bar{p}(k|s) z(k)　 \\
\exp (-q(s)) = \frac{ z(s) }{ \sum_k \bar{p}(k|s) z(k) } \\
q(s) = -\ln \frac{ z(s) }{ \sum_k \bar{p}(k|s) z(k) } \\
q(s) = -\ln \frac{ \exp \{ -V(s) \} }{ \sum_k \bar{p}(k|s) \exp \{ -V(k) \} } 

```

が導出できるので、この式の$V$に推定した$\hat{V}$を入れてやれば$\hat{q}$が求まるっていう算段らしい。たぶん。

## OptQ
OptVは状態価値関数$V$を推定してやることで状態依存のコスト関数$q$を推定するって話だったけど、そんなじれったいことせずに、直接$q$を推定しようっていうアルゴリズムがOptQである。詳しい説明は論文にあるので、ここでは損失関数を挙げる。
状態遷移対が長さ$T(k)$の軌跡群$\varsigma^{(k)}$からサンプルされたと仮定する。軌跡$k$の時刻$t$における状態を$s_{t}^{(k)}$とするとき、$q$についての損失関数$L$は次の式で表される。

```math
L[q(\cdot)] = \sum_k \left ( \ln z(s_{0}^{(k)}) + \sum_{t=0}^{T(k)} q(s_{t}^{(k)})\right )
```
この式で$q$を推定するアルゴリズムが**OptQ**である。推定のたびに、$z$を計算しなければならず、OptVよりは遅いという。でも、個人的には、例えばモデルフリーRLであるQ学習のLMDP版、Z学習(Z-Learning)で$z$を求めることもできると考えられるので、OptVのように、明示的に$\bar{p}$を所与とする必要がない点では優れているかもしれない、と思ったり。

# MaxEntIRL（最大エントロピー逆強化学習）との関係

論文2.4節にちょっと興味深い話が載っているので注釈をつけながら概訳しておく。

> **MaxEntIRL**アルゴリズム (Ziebart+ 2008) は特徴から報酬を推定する。ここでは状態ごとにデルタ関数による「特徴」をもつルックアップテーブルの場合を考える。

この「ルックアップテーブル」なるものは、おそらくone-hotベクトル化することを言っている。

>MaxEntIRLは密度推定アルゴリズムであり、観測された状態訪問頻度に一致する最大エントロピー分布を探す。モーメントマッチング制約をもつ最大エントロピー分布は、指数型分布族に属することが知られる。

ここは[情報幾何学による最大エントロピー原理の再現](https://shuhoyo.hatenablog.com/entry/information-science-principle-of-maximum-entropy)という労作のページが参考になるかもしれない。

>したがって、MaxEntIRLは[次の分布]族内で観測された軌跡の確率を最大化する$q(\cdot)$を見つける問題に帰着する:

式(12)

```math
p_{\operatorname{MaxEnt}}(\varsigma|s_0) \propto \exp \left ( -\sum_{t=0}^T q(s_t) \right )

```
> ボトルネックは、再帰的手順によって行われる最適化の各ステップで、分配関数を計算する点である。

これは

初期状態$s_0$から軌跡$\varsigma$をたどる確率$p(\varsigma|s_0)$　は　$\varsigma$をたどったときの累積報酬$ -\sum_{t=0}^T q(s_t)$に指数をかけたもの　に比例する

という意味だと思う。要するにエキスパートは報酬が多くもらえる軌跡をたどるよねって話。


> 直感的には、MaxEntIRLはIRLの手法と似ている。しかしながら、現在に至るまで、MaxEntIRLはどのような順方向の最適制御問題を反転しているのか、そもそもそのような問題が存在するのかすら不透明であった。

そもそも、逆強化学習は、強化学習の「逆」である。MaxEntIRLは確かにIRLっぽいことはしているのだが、じゃあMaxEntIRLの逆の問題って具体的には何なのさ？っていう疑問は残ったままだった。

> 今、我々はこの疑問に答えることができる。式(12)を式(6)と比較するとき、軌跡の確率は受動ダイナミクスが一様の場合に等しいことがわかる。
つまり、MaxEntIRLは、一様な受動ダイナミクスをもつLMDPの逆にあたる手法である。

式(6)というのは、LMDPにおいて、最適方策からある軌跡$\varsigma$が生成する確率を示す次の式である。

```math
p^{*}(\varsigma|s_0) = \frac{\bar{p}(\varsigma|s_0)\exp \left ( -\sum_{t=0}^{T} q(s_t) \right ) }{z(s_0)}
```

ここで、$z(s_0)$は定数とみなせるので、

```math
p^{*}(\varsigma|s_0) \propto \bar{p}(\varsigma|s_0)\exp \left ( -\sum_{t=0}^{T} q(s_t) \right )
```

が成り立つ。これと、式(12)

```math
p(\varsigma|s_0) \propto \exp \left ( -\sum_{t=0}^T q(s_t) \right )

```

を比べると、たしかに非制御遷移確率$\bar{p}(\varsigma|s_0)$を一様分布としたバージョンになっている。MaxEntIRLは「LMDPの一様分布版」の逆だったのだ！（ΩΩΩ<な、なんだってー!?）

> 実際、Ziebart+ (2008) で分配関数を計算するために使用される再帰[計算]は、Todorov (2007; 2009b)の好適度関数を計算するための反復法に非常に似ている。
それぞれの再帰は、順問題を解くことと計算上同等である。その結果として、MaxEntIRLとOptQの両方がOptVよりも遅く、その上MaxEntIRLはOptQの特殊なケースである。 
MaxEntIRLの一様な受動ダイナミクスに対する制限は、物理システムのモデリングにおいて特に問題となる。

[Ziebart+ (2008)、すなわちMaxEntIRLの論文](http://new.aaai.org/Papers/AAAI/2008/AAAI08-227.pdf)のAlgorithm 1に載っている分配関数$Z_{s_i}$の計算式

1. $Z_{s_i,0} = 1$とする

2. 以下をN回再帰的に計算

```math
Z_{a_{i, j}} = \sum_k P(s_k|s_i,a_{i,j})\exp \{ R(s_i|\theta) \} Z_{s_k} \\
Z_{s_i} = \sum_{a_{i,j}} Z_{a_{i,j}}
```

は、まさに[Todorov (2007)](https://homes.cs.washington.edu/~todorov/papers/TodorovNIPS06.pdf)で好適度関数$z$をかんたんに計算する一例の「繰り返し手法」(iterative method)として示されている

```math
\mathbf{z}_0 = 1 \\
\mathbf{z}_{k+1} = \mathbf{G} \mathbf{\bar{p}} \mathbf{z}_{k}
```
に他ならない。似ているどころか基本的には同じ。
こうなると、MaxEntIRLは各軌跡の生成確率が一様だと仮定していることが問題になりそうなことがわかる。一様分布を仮定するより、何らかの分布が事前に存在することを知ることができるなら、それを利用できたほうが正確な推定ができるのは当然とも言える。