雑記: 分散共分散行列のベイズ更新の話

2021-02-02 絵を追加しました。
いろいろな場面(カルマンフィルタ、ガウス過程回帰など、直接観測できない何かの分布をその線形変換の観測からベイズ更新する場面)で以下の問が出てくると思います。
 x \in \mathbb{R}^n は確率ベクトルでその事前分布は平均 0 で分散共分散行列が V \in \mathbb{R}^{n \times n} の多変量正規分布です。
 y \in \mathbb{R}^m も確率ベクトルで y = Hx + w です。H \in \mathbb{R}^{m \times n} は行列です。w \in \mathbb{R}^m は平均 0 で分散共分散行列が R \in \mathbb{R}^{m \times m} の多変量正規分布にしたがうノイズで x とは独立です。
いまある y が観測されました。この y の下での(=事後分布の)x の分散共分散行列はどうなりますか?
普通はベイズの定理で x確率密度関数を考えると思います。
方法1.ベイズの定理で x確率密度関数をかく

 \begin{split} \displaystyle p(x|y) \propto p(y|x)p(x) &\propto \exp \left( -\frac{1}{2}  (Hx)^\top R^{-1} Hx\right) \exp \left( -\frac{1}{2} x^\top V^{-1} x\right) \\ &= \exp \left( -\frac{1}{2} x^\top (V^{-1} + H^\top R^{-1} H) x\right) \end{split}

より x の事後分散共分散行列は  (V^{-1} + H^\top R^{-1} H)^{-1} である。
こうすると観測前は精度行列(分散共分散行列の逆行列で、多変量正規分布の密度関数の式で確率ベクトルに挟まれているやつ)が V^{-1} だったのが観測後は V^{-1} + H^\top R^{-1} H になっていて、精度が大きくなった(=分散が小さくなった)感があります。観測の精度 R^{-1} が大きいならより精度はより大きくなるし、H=0 なら(観測が x を反映しないなら)精度は全く大きくならないこともわかります。

別に以上でいいんですが、上の方法は終始 x の精度の世界で考えていて、x の分散に何があったのかわかりにくいものになっています(?)。そもそも y が観測されたら x の分布がどうなるのかというのは (x,y) の同時分布を考えるのが正攻法であるはずです(?)。なのでそれでやってみます。また、ここでは(※)の箇所で後述の補題を利用します(条件付き共分散をとる方法は色々あると思いますが、この補題は直接的に条件付き共分散はこれですというものなので利用しました。これは行列をブロック単位で UDL 分解するのと等価です)。

方法2.(x,y) の同時分布を出してから x の分散共分散行列を取り出す
多変量正規分布のモーメント母関数は M_X(t) = E(e^{t^\top X}) = \exp \Bigl( \mu^\top t + \frac{1}{2} t^\top \Sigma t \Bigr) なので、(x,y) の同時分布のモーメント母関数は以下のようになる。

\displaystyle \begin{split} E \left[ e^{ \left( \begin{array}{cc} t_1^{\top} & t_2^{\top} \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right)} \right] &= E \left[ e^{  t_1^{\top} x + t_2^\top H x + t_2^\top w} \right]= E \left[ e^{  t_1^{\top} x + t_2^\top H x} \right] E \left[ e^{ t_2^\top w} \right] \\ &= \exp \Bigl( \frac{1}{2} t_2^\top R t_2 \Bigr) \exp \Bigl( \frac{1}{2} (t_1^\top + t_2^\top H) V (t_1 + H^\top  t_2) \Bigr)  \\ &= \exp \Bigl( \frac{1}{2} \left( t_1^\top V t_1 + t_1^\top V H^\top t_2 + t_2^\top H V t_1 + t_2^\top HVH^\top t_2 + t_2^\top R t_2 \right) \Bigr) \\ &= \exp \left( \frac{1}{2} \left( \begin{array}{cc} t_1^{\top} & t_2^{\top} \end{array} \right)  \left( \begin{array}{cc} V & VH^\top \\ HV & HVH^\top + R \end{array} \right)  \left( \begin{array}{c} t_1 \\ t_2 \end{array} \right) \right)  \end{split}

よって (x,y) の同時分布は分散共分散行列が  \displaystyle  \tilde{V}_{x,y} = \left( \begin{array}{cc} V & VH^\top \\ HV & HVH^\top + R \end{array} \right) の多変量正規分布である。このとき x の分散共分散行列  \tilde{V}_{x} \tilde{V}_{x,y} の右下 m \times m ブロックに対するシューア補行列(※)に他ならず、したがって  V - VH^\top (HVH^\top + R)^{-1} HV である。
こうすると観測前は分散が V だったのが観測後は  V - VH^\top (HVH^\top + R)^{-1} HV になっていて分散が小さくなった感があります。(※)では以下の補題を利用しました(参考: Sherman-Morrison-Woodburyの公式 (Schur補行列) - いんふらけいようじょのえにっき)。
補題(シューア補行列)
正方行列 P \in \mathbb{R}^{(n+m) \times (n+m)}A \in \mathbb{R}^{n \times n}, \, D \in \mathbb{R}^{m \times m} が正方行列になるように  \displaystyle P = \left( \begin{array}{cc} A & B \\ C & D \end{array} \right) と区分けする。このとき PD も正則ならば、 S = A - BD^{-1}C も正則でP^{-1} の左上 n \times n ブロックは S^{-1} になる。この SPD に関するシューア補行列とよぶ。
証明は元の記事にありますが、記事中で示されているのは左上ブロックに対するシューア補行列の方なので、LDU 分解を UDL 分解に変更して右下ブロックに対するシューア補行列の方を示してください。同時共分散行列の LDU 分解には y の共分散が陽にあらわれ、UDL 分解には x の共分散が陽にあらわれます(以下)。なので x の共分散がほしいときは UDL 分解するのが自然なはずです。上の補題は UDL 分解の \tilde{X}A - BD^{-1}C であるといっているのと同じです。

  • LDU 分解すると以下の Yy の共分散: \displaystyle \left( \begin{array}{cc} t_1^\top & t_2^\top \end{array} \right) \left( \begin{array}{cc} I & O \\ W & I \end{array} \right) \left( \begin{array}{cc} X & O \\ O & Y \end{array} \right) \left( \begin{array}{cc} I & Z \\ O & I \end{array} \right) \left( \begin{array}{c} t_1^\top \\ t_2^\top \end{array} \right) = \left( \begin{array}{cc} t_1^\top + t_2^\top W & t_2^\top \end{array} \right) \left( \begin{array}{cc} X & O \\ O & Y \end{array} \right) \left( \begin{array}{c} t_1 + Z t_2 \\ t_2 \end{array} \right)
  • UDL 分解すると以下の \tilde{X}x の共分散: \displaystyle \left( \begin{array}{cc} t_1^\top & t_2^\top \end{array} \right) \left( \begin{array}{cc} I & \tilde{Z} \\ O & I \end{array} \right) \left( \begin{array}{cc} \tilde{X} & O \\ O & \tilde{Y} \end{array} \right) \left( \begin{array}{cc} I & O \\ \tilde{W} & I \end{array} \right) \left( \begin{array}{c} t_1 \\ t_2 \end{array} \right)= \left( \begin{array}{cc} t_1^\top & t_1^\top \tilde{Z} + t_2^\top \end{array} \right) \left( \begin{array}{cc} \tilde{X} & O \\ O & \tilde{Y} \end{array} \right) \left( \begin{array}{c} t_1 \\ \tilde{W}t_1 + t_2 \end{array} \right)


なお、方法1.の結論は  (V^{-1} + H^\top R^{-1} H)^{-1} で、方法2.の結論は  V - VH^\top (HVH^\top + R)^{-1} HV で、なんか違ってみえますが、以下の Sherman-Morrison-Woodbury の公式よりこれらは等しいです。というかむしろ方法1.と方法2.で事後分散共分散行列を求めると Sherman-Morrison-Woodbury の公式が証明できることになります(まともな証明は先の記事にありますが、LDU 分解と UDL 分解の逆行列の成分比較によります)。

Sherman-Morrison-Woodbury の公式
A \in \mathbb{R}^{n \times n}, \, B \in \mathbb{R}^{n \times m}, \, C \in \mathbb{R}^{m \times n}, \, D \in \mathbb{R}^{m \times m} について、AD - CA^{-1}BD A - BD^{-1}C も正則なとき、以下の恒等式が成り立つ。

 (A - BD^{-1}C)^{-1} = A^{-1} + A^{-1} B (D - CA^{-1}B)^{-1} C A^{-1}


では  (V^{-1} + H^\top R^{-1} H)^{-1} V - VH^\top (HVH^\top + R)^{-1} HV のどちらの表現が便利なのかというと、逆行列が取られている行列が前者は n \times n、後者は m \times m なので、m が大きければ前者、n が大きければ後者が逆行列の計算コストが低いはずです(なので推定対象の状態空間が高次元の場合は方法1.で前者を求めてから Woodbury の公式で後者にしましょうとかいわれると思います)。

それで、どうして方法1.と方法2.で逆行列をとる行列のサイズに違いが出てくるのか気になるわけですが、トートロジカルですがどこの行列が正則であることを利用しているかが違うはずです。

  • 方法1.: ベイズの定理を使う場合だと、終始「密度の積(精度の和)」で考えるので逆行列を取るのは最後の仕上げになり、x の精度行列 = n \times n 行列の逆行列をとることになります。
  • 方法2.: 同時分布の分散共分散行列を出してから補題(シューア補行列)を利用する場合だと、この操作で直接利用するのは同時共分散行列の右下 m \times m ブロックが正則であることです。
    • シューア補行列 S が存在するならば  S = A - BD^{-1}C であることを示すときに利用するのは D が正則であることだけです。
    • なお、同時共分散行列の LDU 分解の P^{-1} と成分比較することもできます。この場合は左上 n \times n ブロックの A が正則であることを利用することになり、また、出てくる答えが方法1.と同じになります(方法2.の \tilde{V}_{x,y} の各ブロックを  A^{-1}+A^{-1}B(D - CA^{-1}B)^{-1} CA^{-1} に代入)。なので、同時共分散行列を経由することではなく、そこから如何に x の共分散行列をとるかが逆行列のサイズを決めています。

方法1.

x の事前共分散行列」
 ↓ ベイズの定理
x の事後精度行列」
 ↓ 逆行列 ― 精度行列全体 n \times n が正則であることを利用
x の事後共分散行列」

方法2.

x の事前共分散行列」
 ↓ (x,y) の同時分布のモーメント母関数
(x,y) の同時共分散行列」
 ↓ 補題(シューア補行列)― 同時共分散行列の右下 m \times m ブロックが正則であることを利用
x の事後共分散行列」


まとめると(まとめではない)、(x, y) の同時共分散行列を P とおくとベイズの定理は精度行列 P^{-1} の UDL 分解(P の LDU 分解)の左上ブロックをとることに他ならず、それに Woodbury の公式を適用すると P^{-1} の LDU 分解(P の UDL 分解)の左上ブロックになるので正則を仮定するブロックが左上から右下に切り替わる、という感じがします(無論ベイズ更新の文脈ではどちらも正則でないと困るのですが正則であることを明示的に利用するのがどちらなのかが切り替わる)。

2021-02-02 追記
絵を描きました。スタートから出発してゴールを目指します。★は逆行列の操作なので大きな逆行列はなるべく避けます。ただし図は n > m のイメージで描いています。この大小関係によって適切なルートを選びます。

f:id:cookie-box:20210202184504p:plain