notebooks/RL_Puterman_memo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Puterman本のメモ\n",
    "\n",
    "参考：\n",
    "* [Markov Decision Processes: Discrete Stochastic Dynamic Programming](https://onlinelibrary.wiley.com/doi/book/10.1002/9780470316887)\n",
    "\n",
    "Puterman本を読んだときの面白かった話のメモです．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ２章\n",
    "\n",
    "### 2.3章：一般的なMDPについて\n",
    "\n",
    "* $S$をユークリッド空間上のボレル可測な部分集合とします．すなわち，ユークリッド空間の開集合全体から作るσ代数の元です．\n",
    "* $\\pi=\\left(d_1\\right) \\in \\Pi^{\\mathrm{MD}}$とします．$d_1(s)=a$とします．\n",
    "\n",
    "このとき，ベルマン作用素は，**もし密度$p_1(u \\mid s, a)$**が存在するなら，\n",
    "$$\n",
    "r_1(s, a)+E_s^\\pi\\left\\{v\\left(X_2\\right)\\right\\}=r_1(s, a)+\\int_S v(u) p_1(u \\mid s, a) d u\n",
    "$$\n",
    "として表現されます．この密度はDiscreteだと存在しません．その場合には，次のLebesgue-Stieltjes積分が便利です：\n",
    "\n",
    "$$\n",
    "r_1(s, a)+E_s^\\pi\\left\\{v\\left(X_2\\right)\\right\\}=r_1(s, a)+\\int_S v(u) p_1(d u \\mid s, a)\n",
    "$$\n",
    "\n",
    "さて，これが意味を持つのは，\n",
    "* $v(\\cdot) p_1(\\cdot \\mid s, a)$がルベーグ可積分なとき\n",
    "* $v(\\cdot)$が$p_1(d u \\mid s, a)$に対してLebesgue-Stieltjes可積分なとき\n",
    "  \n",
    "の２つです．さて，さっきの方策をこれに代入しましょう．\n",
    "\n",
    "$$\n",
    "r_1\\left(s, d_1(s)\\right)+\\int_S v(u) p_1\\left(u \\mid s, d_1(s)\\right) d u\n",
    "$$\n",
    "\n",
    "これは$s$についての関数ですね．特にこのベルマン作用素をマルチステップに拡張するには，\n",
    "* $r_1(s, \\cdot), p_1(u \\mid s, \\cdot)$そして$d_1(\\cdot)$\n",
    "  \n",
    "に適切な仮定がないといけません．最低でも，方策$d_1: S\\to A$は可測な関数でないとだめです（そうでないと，２ステップ目以降の積分が取れません）．\n",
    "そして，可測かどうかは$A$上の位相が重要になります．\n",
    "\n",
    "ここで，$d^*$をベルマン方程式を満たす方策とします．\n",
    "これが可測であることを示すためには，\n",
    "$$\n",
    "A_s^* \\equiv \\underset{a \\in A_s}{\\arg \\max }\\left\\{r_1(s, a)+\\int_S v(u) p_1(u \\mid s, a) d u\\right\\}\n",
    "$$\n",
    "が，$\\times_{s \\in S} A_s^*$に$S\\to A$の可測関数を選択定理によって含むことを保証しなければいけません．\n",
    "\n",
    "しかし，これには$r_1, p_1, A_s$に何らかの仮定がないといけません．実際，何の仮定もないと可測な最適方策が存在しない場合があります．\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ５章：無限ホライゾンの基本\n",
    "\n",
    "---\n",
    "\n",
    "![HD_neq_MD](figs/puterman_HD_neq_MD.png)\n",
    "\n",
    "**Problem 5.11: 方策が履歴依存＆決定的でも，対応するOccupancy measureを与えるマルコフな方策が確率的でないといけない場合があります．**\n",
    "\n",
    "上の図のMDPを考えてみましょう．次の履歴依存な方策$\\kappa=\\left(d_1^\\kappa, d_2^\\kappa, \\ldots\\right)$を考えてみます：\n",
    "\n",
    "* $s_1$で選択する行動$a_{1, 1}$と$a_{1, 3}$が次のように履歴に依存します：\n",
    "1. $d_3^\\kappa\\left(s_1, a_{1,1}, s_1, a_{1,3}, s_1\\right)=a_1$：つまり，１→３なら１を選択\n",
    "2. $d_4^\\kappa\\left(s_1, a_{1,1}, s_1, a_{1,3}, s_1, a_{1,1}, s_1\\right)=a_{1,3}$：つまり，１→３→１なら３を選択\n",
    "3. 後は同じ流れ．\n",
    "\n",
    "このとき，$\\kappa$の（割引なし）occupancy measureを考えると\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "x(s_1, a_{1, 1}) &= 1 + \\gamma^2 + \\gamma^4 + \\cdots = \\frac{1}{1 - \\gamma^2}\\\\\n",
    "x(s_1, a_{1, 3}) &= \\gamma + \\gamma^3 + \\cdots = \\frac{\\gamma}{1 - \\gamma^2}\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "このOccupancy measureは明らかにマルコフかつ決定的な方策では出せません（どちらかが０になります）．\n",
    "\n",
    "---\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ６章：割引無限ホライゾン\n",
    "\n",
    "---\n",
    "\n",
    "**定理：6.1.1**\n",
    "これはよくある結果ですが，証明方法が「良い」のでメモっときます．\n",
    "\n",
    "* $0 \\leq \\lambda < 1$\n",
    "* 任意の定常方策$d^\\infty$（すべてのステップで$d \\in D^{\\mathrm{MR}}$に従う）\n",
    "\n",
    "のとき，$v^{d^\\infty}_\\lambda$は次の式の唯一の解であり，\n",
    "$$\n",
    "v=r_d+\\lambda P_d v\n",
    "$$\n",
    "さらに次で表せます：\n",
    "\n",
    "$$\n",
    "v_\\lambda^{d^x}=\\left(I-\\lambda P_d\\right)^{-1} r_d\n",
    "$$\n",
    "\n",
    "**証明**\n",
    "\n",
    "まず，次の行列ノルムを定義します：\n",
    "$$\n",
    "\\|H\\|\\triangleq \\sup_{s \\in S} \\sum_{j \\in S} |H(j\\mid s)|\n",
    "$$\n",
    "\n",
    "このとき，明らかに\n",
    "$$\n",
    "\\left\\|P_d\\right\\|=1 \\text { and } \\lambda=\\left\\|\\lambda P_d\\right\\|\n",
    "$$\n",
    "です．後は[LA_normed_linear_space.ipynb](LA_normed_linear_space.ipynb)を見てください．\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**定理6.2.9**: 任意の状態空間での最適方策の議論をするために便利です．\n",
    "\n",
    "次を満たす方策$d^*$を\"conserving\"と呼ぶことにします\n",
    "\n",
    "$$\n",
    "L_{d^*} v_\\lambda^* \\equiv r_{d^*}+\\lambda P_{d^*} v_\\lambda^*=v_\\lambda^*\n",
    "$$\n",
    "\n",
    "すると，\n",
    "1. conservingな方策\n",
    "2. 最適方策\n",
    "\n",
    "のどちらかが存在する場合，決定的かつ定常な最適方策が存在します．\n",
    "\n",
    "**証明**\n",
    "\n",
    "証明は省略します（そんなに難しくないので）\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**定理 6.1.11**：$S$が有限もしくは可算だとします．このとき，任意の$\\varepsilon > 0$について，$\\varepsilon$-最適な決定的かつ定常方策が存在します．\n",
    "\n",
    "ベルマン方程式から，$\\mathscr{L} v_\\lambda^*=v_\\lambda^*$です（このベルマン作用素は$\\sup$バージョンです）．$\\varepsilon > 0$について，次を満たす$d_{\\varepsilon} \\in D^{\\mathrm{MD}}$を選びます：\n",
    "\n",
    "$$\n",
    "r_{d_{\\varepsilon}}+\\lambda P_{d_{\\varepsilon}} v_\\lambda^* \\geq \\sup _{d \\in D^{\\mathrm{MD}}}\\left\\{r_d+\\lambda P_d v_\\lambda^*\\right\\}-(1-\\lambda) \\varepsilon e=v_\\lambda^*-(1-\\lambda) \\varepsilon e\n",
    "$$\n",
    "\n",
    "変形して，\n",
    "\n",
    "$$\n",
    "r_{d_{\\varepsilon}} \\geq (I - \\lambda P_{d_{\\varepsilon}} )v_\\lambda^*-(1-\\lambda) \\varepsilon e\n",
    "$$\n",
    "\n",
    "$v_\\lambda^{\\left(d_\\varepsilon\\right)^{\\infty}}=\\left(I-\\lambda P_{d_\\varepsilon}\\right)^{-1} r_{d_\\varepsilon}$なので，両辺に$\\left(I-\\lambda P_{d_\\varepsilon}\\right)^{-1}$をかけると\n",
    "$$\n",
    "v_\\lambda^{\\left(d_\\varepsilon\\right)^\\infty} \\geq v_\\lambda^*-\\varepsilon e\n",
    "$$\n",
    "であり，これは$\\varepsilon$-最適方策です．\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "**例6.2.3**：状態空間を任意にとると，ベルマン作用素で$\\max$が取れたとしても，最適方策が存在しない場合があります．\n",
    "\n",
    "* 状態空間：$S=[0, 1]$\n",
    "* 行動空間：$A=[0, 1]$\n",
    "* $W\\subseteq S$：非ボレルな$S$上の部分集合（本ではprojection on $S$って言ってるけど，多分同じな気がする）\n",
    "* $B$：$[0, 1] \\times [0, 1]$の部分集合で，$S$の方は$W\\subset S$に射影した集合？\n",
    "\n",
    "を考えます．\n",
    "\n",
    "ここで，\n",
    "* 報酬：$r(s, a)=I_{\\{B\\}}(s, a)$とします．これは$(s, a) \\in B$なら$1$, そうでなければ$0$です．\n",
    "* 遷移：任意の$a \\in A$で$P(s\\mid s, a)=1$とします．つまり，どの行動を取っても同じ状態から動きません．\n",
    "* $\\pi=\\left(d_1, d_2, \\ldots\\right) \\in \\Pi^{H R}$とします．$d_1(\\cdot)$は$s$のみに依存する関数です．また，積分取りたいので可測関数だけ考えます．\n",
    "\n",
    "このとき，\n",
    "* 最初の行動は$s_1$に依存するので，$d_1$は可測かつ$r_{d_1}(s)=r(s, d_1(s))$も可測です．\n",
    "* よって，$Q=\\{s: r_{d_1}(s)=1\\}$は$S$上のボレル部分集合です．定義より，$Q$は$W$に含まれます．\n",
    "* すると，次を満たす$s_0 \\in W-Q$が存在します：\n",
    "\n",
    "$$\n",
    "v_\\lambda^{\\pi \\pi}\\left(s_0\\right)<\\lambda+\\lambda^2+\\lambda^3+\\cdots=\\lambda(1-\\lambda)^{-1}\n",
    "$$\n",
    "\n",
    "また，$v_\\lambda^{\\pi^{\\prime}}\\left(s_0\\right)=(1-\\lambda)^{-1}$を満たす（非可測な？）$\\pi^{\\prime} \\in \\Pi^{\\mathrm{HR}}$も存在します．\n",
    "よって，任意の$\\varepsilon < 1$について，\n",
    "\n",
    "$$\n",
    "v_\\lambda^\\pi\\left(s_0\\right)+\\varepsilon<v_\\lambda^*\\left(s_0\\right)\n",
    "$$\n",
    "です．後はこれを任意の$\\pi$について繰り返せば，$0 <\\varepsilon < 1$最適な方策が存在しないことがわかります．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "**6.3章**：収束レートについて\n",
    "\n",
    "* $\\{y_n\\} \\subset V$を$y^*$に収束する系列とします．つまり，$\\lim_{n\\to \\infty} \\|y_n - y^*\\|=0$です．\n",
    "* 次を満たすとき，収束の「オーダーが$\\alpha$」である，といいます．\n",
    "\n",
    "$$\\|y_{n+1}-y^*\\| \\leq K \\|y_n - y^*\\|^\\alpha$$\n",
    "\n",
    "を満たす定数$K > 0$が存在する．\n",
    "\n",
    "* 最低でも$\\alpha=1$のとき，線形収束する，といいます\n",
    "* 最低でも$\\alpha=2$のとき，二次収束する，といいます\n",
    "\n",
    "$$\n",
    "\\lim\\sup_{n\\to \\infty} \\frac{\\|y_{n+1}-y^*\\|}{\\|y_n - y^*\\|}=0\n",
    "$$\n",
    "\n",
    "ならば，これは優線形収束する，といいます．\n",
    "優線形収束は線形収束よりも一般に早いです．\n",
    "\n",
    "* オーダー$\\alpha$で収束する場合，$\\|y_{n+1}-y^*\\| \\leq K \\|y_n - y^*\\|^\\alpha$がすべての$n$で成立するような最低の$K$を，収束レートといいます．\n",
    "    * 一般に収束レートは線形収束する系列を比較するときに使います．\n",
    "\n",
    "* 「すべての$n$で成立する」，というのはちょっと強い仮定なので，次のasymptotic average rate of convergence (AARC)が便利です．\n",
    "\n",
    "$$\n",
    "\\lim \\sup_{n \\to \\infty} \\left[\\frac{\\|y_n - y^*\\|}{\\|y_0 - y^*\\|}\\right]^{1/n}\n",
    "$$\n",
    "\n",
    "* 非負実数関数$f(n)$に対して，$\\{y_n\\}$は次を満たすときに$O(f_n)$であるといいます：\n",
    "\n",
    "$$\n",
    "\\lim\\sup_{n\\to \\infty}\\frac{y_n - y^*}{f(n)}\n",
    "$$\n",
    "が有限である．**例えば価値反復法は$O(\\lambda^n)$で収束します（証明は定理6.3.3.）**．\n",
    "\n",
    "このとき，特に$y_n = y^* + O(f(n))$とかけます．もし$0 < \\beta < 1$について，$f(n)=\\beta^n$ならば，この収束はレート$\\beta$でgeometricに収束する，といいます．\n",
    "\n",
    "* 一般に（少なくともこの本では），初期値に関わらず特定のオーダーで収束するならば，それはglobalな収束である，といいます．\n",
    "* 一般に（少なくともこの本では），特定の初期値で収束するならば，それはlocalな収束である，といいます．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "**6.3.3：ガウス-ザイデル価値反復法**\n",
    "\n",
    "[RL_Gauss_Seidel_VI.ipynb](RL_Gauss_Seidel_VI.ipynb)を参照．\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A.4 行列の極限\n",
    "\n",
    "ここでの結果は基本的に有限と可算な状態空間のマルコフ連鎖に発展できます．\n",
    "\n",
    "* $\\left\\{A_n: n \\geq 0\\right\\}\n",
    "* $\\lim _{n \\rightarrow \\infty} A_n=A$はすべての要素が収束することを表します．\n",
    "\n",
    "periodicなマルコフ連鎖を考えると，この極限が存在しない場合があります．変わりに，次のCesaro limitを考えましょう．\n",
    "\n",
    "---\n",
    "\n",
    "**Cesaro limit**\n",
    "\n",
    "$$\n",
    "\\lim _{N \\rightarrow \\infty} \\frac{1}{N} \\sum_{n=0}^{N-1} A_n=A\n",
    "$$\n",
    "$$\n",
    "C-\\lim _{N \\rightarrow \\infty} A_N=A\n",
    "$$\n",
    "\n",
    "普通の極限は「オーダー０のCesaro limit」と呼ばれることがあります．\n",
    "\n",
    "極限行列を$P^*=C-\\lim P^N$と書くことにします．\n",
    "これは\n",
    "$$\n",
    "p^*(j \\mid s)=\\lim _{N \\rightarrow \\infty} \\frac{1}{N} \\sum_{n=1}^N p^{n-1}(j \\mid s)\n",
    "$$\n",
    "を意味します．ここで，$p^{n-1}(j\\mid s)$は$P^{n-1}$の要素であり，「状態$s$から$n-1$ステップで状態$i$に到達する確率」を表します．\n",
    "$E_s\\left[\\nu_j\\right]$を$s$からスタートして$j$を踏む期待回数とすると，\n",
    "$$\n",
    "E_s\\left\\{\\nu_j\\right\\}=\\sum_{n=1}^{\\infty} p^{n-1}(j \\mid s)\n",
    "$$\n",
    "が成り立ちます．よって，$p^*(j \\mid s)$は「$s$からスタートして，$j$をどれくらい占有するか？」を表します．\n",
    "\n",
    "Cesaro limitは極限が取れない行列に対して便利です．例えば，\n",
    "\n",
    "$$\n",
    "P=\\left[\\begin{array}{ll}\n",
    "0 & 1 \\\\\n",
    "1 & 0\n",
    "\\end{array}\\right]\n",
    "$$\n",
    "\n",
    "を考えると，これは$P^{2 n}=I$ and $P^{2 n+1}=P$なので，明らかに極限が存在しません．一方で，\n",
    "$$\n",
    "C-\\lim _{N \\rightarrow \\infty} P^N=P^*=\\left[\\begin{array}{ll}\n",
    "0.5 & 0.5 \\\\\n",
    "0.5 & 0.5\n",
    "\\end{array}\\right]\n",
    "$$\n",
    "になります．実際，Cesaro limitの極限は（状態空間が有限ならば）存在します：\n",
    "\n",
    "**証明**\n",
    "\n",
    "次の定理を使います（Theorem A.5）：\n",
    "（多分[Generalized Inverses in Discrete Time Markov Decision Processes](https://epubs.siam.org/doi/10.1137/0610009)に証明があるっぽい．後で読んでおこう．\n",
    "\n",
    "$S$が有限で，$P$が$m$個のrecurrentなクラスを持つとします．このとき，\n",
    "1. $1$は$P$の固有値である．また，その多重度が$m$であり，$m$個の線形独立な固有ベクトル．\n",
    "2. 次を満たす正則行列$W$が存在する（ここで，$I$は$m\\times m$の行列）：\n",
    "$$\n",
    "P=W^{-1}\\left[\\begin{array}{ll}\n",
    "Q & 0 \\\\\n",
    "0 & I\n",
    "\\end{array}\\right] W\n",
    "$$\n",
    "ここで，\n",
    "* $1$は$Q$の固有値ではない\n",
    "* $\\sigma(Q) \\leq 1$．もし$P$のrecurrentなsubchainがaperiodicなら，さらに$\\sigma(Q) < 1$が成り立つ\n",
    "* $(I - Q)^{-1}$が存在する（これは$1$が$Q$の固有値じゃないので，$(I-Q)$の固有値に$0$がないことを意味します．よって，逆行列が存在します）\n",
    "* $\\sigma(I-Q) = \\sigma(I-P)$が成り立つ\n",
    "3. 次を満たす行列$P^*$が唯一存在する：$P P^*=P^* P=P^* P^*=P^*$を満たし，$P_i^* > 0$が$\\forall i=1, \\dots, m$が成立する．さらに，\n",
    "$$\n",
    "P^*=W^{-1}\\left[\\begin{array}{ll}\n",
    "0 & 0 \\\\\n",
    "0 & I\n",
    "\\end{array}\\right] W\n",
    "$$\n",
    "\n",
    "\n",
    "これを使うと，$P^*$が存在することが言えます．上の補題から，\n",
    "$$\n",
    "P^n=W^{-1}\\left[\\begin{array}{cc}\n",
    "Q^n & 0 \\\\\n",
    "0 & I\n",
    "\\end{array}\\right] W\n",
    "$$\n",
    "です．よって，\n",
    "$$\n",
    "\\frac{1}{N} \\sum_{n=0}^{N-1} P^n=W^{-1}\\left[\\begin{array}{cc}\n",
    "\\frac{1}{N} \\sum_{n=0}^{N-1} Q^n & 0 \\\\\n",
    "0 & I\n",
    "\\end{array}\\right] W\n",
    "$$\n",
    "が成り立ちます．$I-Q$が正則なので，\n",
    "\n",
    "$$\n",
    "\\sum_{n=0}^{N-1} Q^n=\\left(I-Q^N\\right)(I-Q)^{-1}\n",
    "$$\n",
    "が言えます（ノイマン展開？）．ここで，$\\sigma(Q) \\leq 1$なので，$Q^N$は有界であり，\n",
    "$$\n",
    "\\lim _{N \\rightarrow \\infty} \\frac{1}{N} \\sum_{n=0}^{N-1} Q^n=0\n",
    "$$\n",
    "が成り立ちます．よって，上の補題から，\n",
    "\n",
    "$$\n",
    "\\lim _{N \\rightarrow \\infty} \\frac{1}{N} \\sum_{n=0}^{N-1} P^n=W^{-1}\\left[\\begin{array}{ll}\n",
    "0 & 0 \\\\\n",
    "0 & I\n",
    "\\end{array}\\right] W=P^*\n",
    "$$\n",
    "\n",
    "であり，極限行列が存在します．\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "極限行列$P^*$は次を満たします：\n",
    "\n",
    "$$\n",
    "P P^*=P^* P=P^* P^*=P^*\n",
    "$$\n",
    "\n",
    "ここで，\n",
    "\n",
    "$$\n",
    "\\left(P^*\\right)^2=P^*,\\left(I-P^*\\right)^2=\\left(I-P^*\\right) \\text {, and } P^*\\left(I-P^*\\right)=0\n",
    "$$\n",
    "\n",
    "なので，$P^*$と$(I-P^*)$は直行するProjection matrixです．\n",
    "\n",
    "---\n",
    "\n",
    "極限行列はcountableな状態空間の場合，要素が$0$になり得ることに注意しましょう．\n",
    "例えば，\n",
    "$p(s+1\\mid s)=1$が$s=0, 1, \\dots$であり，そうでないならば$p(j\\mid s)=0$のとき，極限行列$P^*$はゼロ行列になります．\n",
    "\n",
    "行の和がちゃんと１になる$P$のことを確率的と呼びます．\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "**$P^*$の計算方法**\n",
    "\n",
    "次の定理が便利です．\n",
    "\n",
    "$P$がpositive recurrent irreducible chainだとします．\n",
    "このとき，\n",
    "\n",
    "$$\n",
    "q^T = q^T P\n",
    "$$\n",
    "かつ$\\sum_{j \\in S} q(j)=1$は唯一の正の解を持ちます．\n",
    "\n",
    "この問題の解を**定常分布**と呼びます．\n",
    "\n",
    "ここで，$P^*P=P^*$であることを思い出しましょう．\n",
    "$P$がrecurrentかつirreducibleのとき，$P^*$は同じ行を持ち，\n",
    "$$\n",
    "P^* = eq^T\n",
    "$$\n",
    "と書けます．\n",
    "\n",
    "**TODO: 続き書く**\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "**Drazin逆行列**\n",
    "\n",
    "マルコフ決定過程では，次の方程式を頻繁に解きます：\n",
    "\n",
    "$$\n",
    "(I-P) v=r\n",
    "$$\n",
    "\n",
    "ここで，$(I-P)$は正則ではありません（１が$P$の固有値なので）．\n",
    "よって，この方程式は唯一の解を持ちません．\n",
    "\n",
    "上でやった補題から，\n",
    "$$\n",
    "I-P=W^{-1}\\left[\\begin{array}{cc}\n",
    "I-Q & 0 \\\\\n",
    "0 & 0\n",
    "\\end{array}\\right] W\n",
    "$$\n",
    "なので，$(I-P)^{-1}$は存在しません．そこで，一般逆行列を考えましょう．\n",
    "\n",
    "\n",
    "次の形式を持つ行列$B$を考えましょう\n",
    "$$\n",
    "B=W^{-1}\\left[\\begin{array}{ll}\n",
    "C & 0 \\\\\n",
    "0 & 0\n",
    "\\end{array}\\right] W\n",
    "$$\n",
    "ここで，$C$は正則行列とします．\n",
    "また，\n",
    "$$\n",
    "B^{\\#}=W^{-1}\\left[\\begin{array}{cc}\n",
    "C^{-1} & 0 \\\\\n",
    "0 & 0\n",
    "\\end{array}\\right] W\n",
    "$$\n",
    "としましょう．すると，\n",
    "\n",
    "$$\n",
    "B^\\# B B^{\\#}=B^\\#, \\quad B B^\\#=B^\\# B, \\quad \\text { and } \\quad B B^\\# B=B\n",
    "$$\n",
    "\n",
    "が成り立ってることがすぐにわかります．この条件を満たす行列$B^\\#$をDrazin逆行列といいます．\n",
    "\n",
    "Drazin逆行列を使うと，次が成り立ちます：\n",
    "\n",
    "1. 行列$(I-P+P^*)$は正則です．その逆行列を$Z_P$とします．これは次の式からすぐにわかります：\n",
    "$$\n",
    "I-P+P^*=W^{-1}\\left[\\begin{array}{cc}\n",
    "I-Q & 0 \\\\\n",
    "0 & I\n",
    "\\end{array}\\right] W\n",
    "$$\n",
    "2. Drazin逆行列$H_P$は次を満たします：\n",
    "$$\n",
    "(I-P)^\\#=\\left(I-P+P^*\\right)^{-1}\\left(I-P^*\\right) \\equiv H_P\n",
    "$$\n",
    "これは次の式から明らかです：\n",
    "$$\n",
    "\\begin{aligned}\n",
    "(I-P)^\\# & =W^{-1}\\left[\\begin{array}{cc}\n",
    "(I-Q)^{-1} & 0 \\\\\n",
    "0 & 0\n",
    "\\end{array}\\right] W \\\\\n",
    "& =W^{-1}\\left[\\begin{array}{cc}\n",
    "(I-Q)^{-1} & 0 \\\\\n",
    "0 & I\n",
    "\\end{array}\\right] W-W^{-1}\\left[\\begin{array}{ll}\n",
    "0 & 0 \\\\\n",
    "0 & I\n",
    "\\end{array}\\right] W \\\\\n",
    "& =\\left(I-P+P^*\\right)^{-1}-P^*=Z_P-P^*\n",
    "\\end{aligned}\n",
    "$$\n",
    "3. $H_P=C-\\lim _{N \\rightarrow \\infty} \\sum_{n=0}^{N-1}\\left(P^N-P^*\\right)$であり，$P$がaperiodicなら，Cesaro-limitはlimitになります．\n",
    "証明はちょっと厄介ですが，\n",
    "$$\n",
    "\\left(P-P^*\\right)^n=P^n-P^*\n",
    "$$\n",
    "を使えばいけます．\n",
    "\n",
    "---\n",
    "\n",
    "次の式が成立します：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "(I-P) & H_P=H_P(I-P)=I-P^* \\\\\n",
    "H_P P^* & =P^* H_P=0 \\\\\n",
    "H_P & =Z_P-P^* \\\\\n",
    "Z_P P^* & =P^* \\\\\n",
    "P^* & =I-(I-P)(I-P)^\\#\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "便利なので覚えておきましょう．\n",
    "\n",
    "---\n",
    "\n",
    "**$H_P$の解釈**\n",
    "\n",
    "aperiodicな連鎖の場合は，\n",
    "$$\n",
    "H_P=\\lim _{N \\rightarrow \\infty}\\left[\\sum_{k=0}^{N-1} P^k-N P^*\\right]\n",
    "$$\n",
    "が成り立ちます．つまり，$H_P(j\\mid s)$は，\n",
    "1. $s$からスタートして$j$を訪れる期待訪問回数\n",
    "2. 遷移確率が$P^*$のときに，$s$からスタートして$j$を訪れる期待訪問回数\n",
    "\n",
    "の差を表します．\n",
    "\n",
    "この$H_P$はバイアスの定義に使います．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}