ワープのある双六でゴールする時間の期待値をガチで求める -吸収マルコフ連鎖-

# 概要

「コイン投げで表が出るまで何回かかる？」「ガチャで全種類のカードをそろえるまでおいくら？」「振出しに戻るマスがあるすごろくでゴールするまでにサイコロを振る回数の期待値」といった問題を解くとき、動的計画法がよく使われます（ちょっと難しい問題だと連立方程式を解きます）。
そうです。`dp[i] = i種類のカードを持っている状態からカードをコンプするまでに支払う金額の期待値`とおいて、`dp[i] = (p*dp[i] + (1-p)*dp[i+1])+1`みたいな式を立てて`dp[n]`から逆順にdpテーブルを埋めたりする、アレです。

でも、期待値ってそもそも 回数*確率 の和ってのが定義ですよね。定義通り計算してないけど、これ大丈夫なの？そもそも期待値は有限なの？無限だったらこの議論ヤバくね？

この記事では、誰しも一度は考えたことのあるこうした疑問を「吸収マルコフ連鎖」により数学的に解決します。
（注意：これを読んでも競プロの問題が解けるようにはなりません）

# 数学的設定

- スタートがマス $0$, ゴールがマス $N$ となる合計 $N+1$ マスのすごろくがある
- $1$ 秒ごとにサイコロをふるとき、ゴールにたどり着くまでの秒数の期待値を求めたい
- マスごとに振るサイコロが異なる（「次にいけるマス」の場所や確率が異なる）
- でも過去のサイコロの目にはよらない。あくまで「今いるマス」だけで決まる
- どのマスから出発しても、運が良ければいつかはゴールにたどり着ける

## マルコフ連鎖
過去のサイコロの目にはよらず、「今いるマス」だけで「次にいけるマス」が決まるという状況は「マルコフ連鎖」と呼ばれます。マス $i$ からマス $j$ に移動する確率を $p(i,j)$ とおきましょう。この $p(i,j)$ を行列の形に並べた

```math
P = (p(i,j))_{i,j} =  
\left(
\begin{array}{ccccc}
p(0,0) & \cdots & p(0,N)\\
\vdots & \ddots &  \vdots \\
p(N,0) & \cdots  & p(N,N)
\end{array}
\right)
```

は遷移行列といいます。

大事な事実として、$P^t$ の $(i,j)$ 成分は 「$i$ を出発して $t$ 秒後に $j$ にいる確率」となります。これは、「$i$ を出発して $t$ 秒後に $j$ にいる確率 = （$i$ を出発して $t-1$ 秒後に $k$ にいる and $t$ 秒目に $k$ から $j$ に移動する）確率の $k$ についての和」という式を行列に直せばわかります。つまり、$P^t$ の $i$ 行目は、「$i$をスタートして$t$ 秒後に $j$ にいる確率」を $j$ について横に並べたものになっています（ので横に和をとると $1$ になります）。

さらに今回はマス $N$ が「ゴール」です。つまり一度ゴールに到達したらずっとそのまま、と思うと、 $j \neq N$ なら $p(N,j) = 0$, $p(N,N) = 1$ です。このような元を吸収元といい、吸収元があるマルコフ連鎖を吸収マルコフ連鎖といいます。
吸収マルコフ連鎖の場合、遷移行列 $P$ は行、列をサイズ $N$ とサイズ $1$ に区切って、

```math
P = \begin{pmatrix} Q & R \\ 0 & 1 \end{pmatrix}
```

と書けます。ここで $Q$ は $N \times N$ 行列、$R$ は $N \times 1$ の行列（縦ベクトル）、$0$ は $1 \times N$ のゼロ行列（横ベクトル）です。


## ゴールにたどり着く時間の期待値
ちょうど $t$ 秒後にゴールにたどり着くとは、$t-1$ 秒目はゴールでなく、$t$ 秒目にゴールにいることです。よって、$e_0 = (1,0,0,\dots,0)$ とおくと、$ e_0 Q^{t-1} R$ が、ちょうど $t$ 秒後にゴールにたどり着く確率です。
ここから期待値を定義通り計算できます。こちらがその値です。
$$ \sum_{t=0}^\infty te_0Q^{t-1}R = e_0\left(\sum_{t=0}^\infty tQ^{t-1}\right)R $$
さて、行列の無限和が出てきてしまいました。はたしてこれは収束するのでしょうか？$\sum_{t=0}^\infty tx^{t-1} = 1/(1-x)^2$という式は知っているのですが、似た計算ができないものでしょうか......

## 級数の収束

突然ですが、等比級数の和の公式 $ 1 + r + r^2 + \cdots = \frac{1}{1-r} $ というものがあります。これは $|r|<1$ のときのみ成り立つのでした。では、この行列バージョン $ I + Q + Q^2 + \cdots = (I-Q)^{-1} $ はいつ成り立つのでしょうか？実は、$Q$の「スペクトル半径」 $\rho(Q)$ が $\rho(Q)<1$ をみたすとき、この行列版の等比級数の和の公式が成り立つことが知られています。
スペクトル半径とは何ぞや？それは $Q$ の固有値の絶対値の最大値です。固有値は「固有ベクトル方向への拡大倍率」を表すので、どの固有値の絶対値もみんな $1$ より小さければ、各固有ベクトルは変換後に縮むから $Q^t$ の各成分もガンガン小さくなる、っていう気分ですね（正確には、ジョルダン標準形のベキの各成分が指数的に $0$ に近づくことから示せます）。

そして実は、「どのマスからもいつかはゴールにたどり着ける可能性がある」吸収マルコフ連鎖は、 $\rho(Q) < 1$ を満たします！
<details><summary>証明</summary>
$Q$の固有ベクトル $v$、固有値 $\lambda$ をとる（つまり $Qv = \lambda v$）。$v = (v_0,\dots,v_N)^\top$ として、絶対値最大の成分を $v_i$ とする（とくに $v_i \neq 0$）。「どのマスからもいつかはゴールにたどり着ける可能性がある」ので、マス $i$ を出発して $C$秒後にゴール $N$ に到達する確率が正となるような $C$ が存在する。つまり　$p^C(i,N) > 0$.これは $\sum_j Q^C(i,j) < 1$ を意味する。よって
$$
\lambda^C |v_i| = |\sum_j Q^C(i,j) v_j| \le \sum_j |Q^C(i,j)||v_j|
\le |v_i| \sum_j |Q^C(i,j)| < |v_i|
$$なので、$|\lambda| < 1$ を得る。
</details>

よって $Q$ は「等比級数の和の公式」が適用できる行列です。同様に、式 $\sum_{t=0}^\infty tx^{t-1} = 1/(1-x)^2$ の行列版についても、スペクトル半径が $1$ より小さいときに成りたつことが証明できます。よってそれを $Q$ に適用して、 $$\sum_{t=0}^\infty t Q^{t-1} = (I-Q)^{-2}$$とわかります。

## 計算をつづけると
ここまでの議論により、期待値は $ e_0(I-Q)^{-2}R $ となることが数学的に証明できました！でもまだ計算は続けられます。

いま、$\vec{1}$ を全成分が $1$ となるサイズ $N$ の縦ベクトルとします（つまり $\vec{1} = (1,1,\dots,1)^\top$、ここでは転置行列の記号として$\top$を使います ）。このとき $ (I-Q)^{-1}R = \vec{1} $ が成り立ちます。（ $Q\vec{1} + R = P \vec{1} = \vec{1}$ より。ここで$P$と書いたのは、正確には $P$ の上 $N$ 行。）

よって期待値は$ e_0(I-Q)^{-1}\vec{1} $と表せます。

## いつもの計算を導く

ここで「初期位置が i のときの吸収される期待値」を`p[i]`とおきましょう。いままでの議論から $p[i] = e_i(I-Q)^{-1} $ となります（$e_i$ は $i$ 番目だけ $1$, ほかは $0$ となる横ベクトル）。これを縦に並べて、


```math
\begin{pmatrix} p[0] \\ \vdots \\ p[N-1] \end{pmatrix}
= (I-Q)^{-1} \begin{pmatrix} 1 \\ \vdots \\ 1 \end{pmatrix},
\qquad
```

すなわち

```math
\begin{pmatrix} p[0] \\ \vdots \\ p[N-1] \end{pmatrix}
= Q\begin{pmatrix} p[0] \\ \vdots \\ p[N-1] \end{pmatrix}
+ \begin{pmatrix} 1 \\ \vdots \\ 1 \end{pmatrix}
```
がわかります。このように、普段通りの連立方程式が定義通りの計算から導けました。

# まとめ
期待値を定義通りに計算すると、その値は $ e_0(I-Q)^{-1} \vec{1} $ という有限値になります。これを変形すると、普段通りの`dp[i] = i スタートでゴールにつく秒数の期待値`に関する連立方程式が立ちます。安心して自己ループを消してDPしましょう。