スパースモデリングを理解したい

最近、青の本*1でスパーズモデリングの勉強を始めた。

最適化の記事同様、ここではあまり式を使わずに全体を俯瞰できるように努める。

まずは種々アルゴリズムの概要を掴み、後の記事で機能する理由について数学的側面を纏めていく予定。

随時更新。

この書籍は、劣決定の連立方程式

$$ \begin{align} \boldsymbol{Ax} = \boldsymbol{b} \quad \boldsymbol{A} \in \mathbb{R}^{n \times m} \end{align} $$

についてスパースな解を得る様々な方法について理論と実践の両方を解説してくれている。 どちらかというと理論重きで、応用も画像処理のことがメイン。
「動かしたい」というよりも「理解したい」と思う人が読むべき本だと思う。

l0最小化問題（（P0）問題）

OMP, LS-OMP, WMPなどサポートを選んで解を更新することを繰り返す追跡アルゴリズム（MP）と呼ばれる手法がある。
また、最初にサポートを全て決定してから解を求めるしきい値アルゴリズムと呼ばれる手法もある。 LS-OMPなど逆行列計算が必要なアルゴリズムにおいては、行列が特異行列にならないなど注意が必要。
書籍では性能比較もしているが、検証の範囲ではOMPとLS-OMPはあまり性能差がないように思えた（個人の感想）。
計算が不安定になるくらいならLS-OMPではなく普通のOMPを使えば良いと思う。

直交マッチング追跡（Orthogonal Matching Pursuit; OMP）

1. 残差との内積が最大となるの列を選んでサポートを追加

2. 以下を解いて暫定解を更新

$$ \begin{align} ||\boldsymbol{b}-\boldsymbol{Ax}||_2^2 \quad \text{subject to Support }{\boldsymbol{x}}=S_{k} \end{align} $$

1. 得られたを用いて残差を更新

最小二乗-直交マッチング追跡（Least-squares OMP; LS-OMP）

OMPを効率的に解く方法。暫定解と追加されるサポートの候補を以下の式で一度に求める。

$$ \begin{align} \left[\begin{array}{c} \boldsymbol{x}_{S} \\ z \end{array}\right]_{\mathrm{opt}}=\left[\begin{array}{cc} \boldsymbol{A}_{S}^{\top} \boldsymbol{A}_{S} & \boldsymbol{A}_{S}^{\top} \boldsymbol{a}_{i} \\ \boldsymbol{a}_{i}^{\top} \boldsymbol{A}_{S} & ||\boldsymbol{a}_{i}||_{2}^{2} \end{array}\right]^{-1}\left[\begin{array}{c} \boldsymbol{A}_{S}^{\top} \boldsymbol{b} \\ \boldsymbol{a}_{i}^{\top} \boldsymbol{b} \end{array}\right] \end{align} $$

は最初に一度だけ計算しておけばよく、あとは最小の誤差を与えるを決定すれば良い。
との列が直交している場合はブロック対角となり、これまでに求めた解は変化せず、要素の追加だけが行われる。

マッチング追跡（MP）

上記よりもより雑な方法。サポート追加時に解の再計算を行わない。 OMPとは異なり前ステップまでに得られた解を更新せずに、新しい要素だけを追加するより単純な方法。そのため精度はOMPより低い。

弱マッチング追跡（WMP）

1. サポートに追加する次の要素の選択を準最適解で満足する方法。MPよりも雑な方法（ただしメリットもちゃんとある）
2. 最適な洗濯から倍だけ離れているいずれかの列を選択する。要は候補列と残差の最大の内積を探すのではなく、tだけ小さい閾値を超え得た最初の列を採用することにする（その分１ループあたりの計算が速くなる）。
3. 探索範囲を限定することで、探索が高速になる可能性がある（全体として高速になるという保証はない）

しきい値アルゴリズム

OMPアルゴリズムを簡単にしたもの。最初の射影で内積が大きい順にサポートを列選んでしまう。つまり、最初にの全ての列について

$$ \begin{align} \min _{z_{j}} ||\boldsymbol{a}_{j} z_{j} - \boldsymbol{b}||_{2}^{2} \end{align} $$

を計算して、ソートしておけば良い。ただ、閾値は自分で決めるので、非ゼロ要素の個数が既知であると想定していることに注意。または、許容誤差をクリアするまで選ぶサポートの数を増やしていく方法もある。

l0からl2への緩和

最小化問題は離散最適化なので扱いにくい。 そこで、 最小化問題を)]最小化問題など、滑らかな関数の問題に緩和する手法や最小化問題に置き換える方法がある。これにより一般的な最適化問題として取り扱えるようになる。ただし、大域最適解へ収束すると言う保証はなくなる。

FOCUSSアルゴリズム（Focal undetermined system solver）

反復重み付け最小二乗法（Itereative-reweighted-least-squares; IRLS）を用いて、ノルム（]）を重み付きノルムで表現する方法。 現在の近似解をと書くことにして、以下の問題を反復して解く。

$$ \begin{align} \left(M_{k}\right): \min_{\boldsymbol{x}} || \boldsymbol{X}_{k-1}^{+} \boldsymbol{x} ||_{2}^{2} \quad \text { subject to } \quad \boldsymbol{b}=\boldsymbol{Ax} \end{align} $$

ここで、

$$ \begin{align} \boldsymbol{X}_{k-1} = \mathrm{diag} \left( |\boldsymbol{x}|^q \right) \end{align} $$

とすると

$$ \begin{align} \boldsymbol{X}_{k-1}^{-1} \boldsymbol{x} &= \left( \frac{x_1}{\left| x_1^q \right|}, \ldots, \frac{x_m}{\left| x_m^q \right|} \right)^{\top} \\ || \boldsymbol{X}_{k-1}^{-1} \boldsymbol{x} ||_2^2 &= ||\boldsymbol{x}||_{2-2q}^{2-2q} \end{align} $$

であり、と考えればを模したことになる。つまり、の値によっていろいろなノルムをノルム最小化問題に置き換えることができる。特に、とすればの置き換えとなる。IRLSの「反復重み付け」は、このノルムにおけるに由来してるのだと思う（個人の感想）。実際には一般化逆行列を用いる点に注意。また、そのまま用いると対角成分にゼロが現れるので注意。計算に影響するのは選ばれたサポートに関する対角成分と暫定解の内積なので、計算上は対角ゼロ要素の部分は逆行列をとってもゼロとして仕舞えばよい。

最小化問題に置き換えたので、ラグランジュの未定乗数法で解くことができる。得られた解を使ってを更新し、再び同様の手順で計算を行うことを繰り返す。対角成分が0（つまり近似解にゼロ要素が現れたら）だと以降ずっとゼロになるので、初期解は全要素を非ゼロにする必要がある点に注意。

l0からl1への緩和

問題を問題に置き換える手法としては、基底追跡（Basis pursuit; BP）がある。ではの要素の値は影響しなかったが、では影響するのでの列を正規化しておきましょうという話である（当然列ごとの重みは変えられるが）。

計算時間はBPよりOMPが圧倒的に短くて済むらしいが、厳密解との誤差、サポート間距離などの評価尺度ではIRLS（FOCUSSアルゴリズム）やBPの方が断然良い。

個人的には次元削減したいのか、解が欲しいのかで使い分ければ良いのだろうと思った。

一般的な緩和法

を単純にに置き換えてしまう方法。気をつけるべきはではxの要素の値は影響しなかったが、以外では影響する点である。

つまり、行列の列を正規化しておくなど、スケーリングを行っておく必要がある。

制約条件の緩和（等式制約から不等式制約へ）

等式制約条件が厳しすぎるので、許容誤差を導入して問題を以下のように緩和する。

$$ \begin{align} \left(P_{0}^{\epsilon}\right): \min _{\boldsymbol{x}}||\boldsymbol{x}||_{0} \text { subject to }||\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x}||_{2} \leq \epsilon \end{align} $$

一意性と最適性の保証がなくなるので、に応じた解の安定性について考える必要がある。

解の安定性（最もスパースな場合）

以下の問題設定で、近似解を求めたいとする。

$$ \begin{align} \boldsymbol{x}_{0}^{\epsilon}=\operatorname{argmin}_{\boldsymbol{x}} || \boldsymbol{x} ||_{0} \quad \text { subject to } \quad || \boldsymbol{b}-\boldsymbol{A} \boldsymbol{x} ||_{2} \leq \epsilon \end{align} $$

この式において、はノイズと解釈できる。すなわち、ノイズの大きさだけ解が不安定になる。非ゼロ要素が１つだけしかない最もスパースな解について考えると、一定の範囲までは、選ばれるサポートは一意だが、ノイズが大きくなるとサポートの一意性すらなくなる。詳細は書籍p96, 図5.2を参照すること。

追跡アルゴリズムの拡張

先述した追跡アルゴリズムを、問題を解くことからを解く問題へと緩和する。

OMP

停止条件をとすれば

$$ \begin{align} ||\boldsymbol{b} - \boldsymbol{Ax}||_2 \le \epsilon \end{align} $$

が満たされるまで解ベクトルの非ゼロ要素を追加していけば良い。

BPDN (Basis pursuit denoising; BPDN)

OMPの場合とほとんど同様だが、ノルムをノルムへ置き換えて以下の問題を解く。
（行列の各列は正規化されているものとする。）

$$ \begin{align} \left(P_{1}^{\epsilon}\right): \min _{\boldsymbol{x}}||\boldsymbol{x}||_{1} \text { subject to }||\boldsymbol{b}-\boldsymbol{Ax}||_{2} \leq \epsilon \end{align} $$

この問題は、ラグランジュ乗数を用いて以下の制約なし最適化を解くことと同じである。

$$ \begin{align} \left(Q_{1}^{\lambda}\right): \min _{\boldsymbol{x}} \lambda||\boldsymbol{x}||_{1}+\frac{1}{2}||\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x}||_{2}^{2} \end{align} $$

IRLSによる解法

を解く単純な手法として反復再重み付け最小二乗法（IRLS）がある。FOCUSSアルゴリズムで用いている手法だが、枠組みが同じだけで計算の仕方は異なるので注意。

とすると、とと書ける。これはノルムを重み付きノルムで表現したことになる。これを利用するとは以下のように置き換えられる。

$$ \begin{align} \left(M_{k}\right): \min _{\boldsymbol{x}} \lambda \boldsymbol{x}^{\top} \boldsymbol{X}_{k-1}^{-1} \boldsymbol{x}+\frac{1}{2}||\boldsymbol{b}-\boldsymbol{Ax}||_{2}^{2} \end{align} $$

二次最適化になっており、標準的な線形代数を用いた反復計算を用いて解くことができる。具体的には、この目的関数をで微分して0とおくと

$$ \begin{align} \left(2 \lambda \boldsymbol{X}_{k-1}^{-1}+\boldsymbol{A}^{\top} \boldsymbol{A}\right) \boldsymbol{x}=\boldsymbol{A}^{\top} \boldsymbol{b} \end{align} $$

が得られるので、各ステップでについて数値最適化を行い、得た近似解で重み行列を更新することを繰り返すだけである。解の変化がしきい値以下になったら計算を終了する。は経験的に決定される。

LARSアルゴリズム

以下の目的関数を最小化することを考える。

$$ \begin{align} f(\boldsymbol{x})=\lambda||\boldsymbol{x}||_{1}+\frac{1}{2}||\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x}||_{2}^{2} \end{align} $$

この式の劣微分（劣勾配集合）

$$ \begin{align} \partial f(\boldsymbol{x})=\left\{\boldsymbol{A}^{\top}(\boldsymbol{A} \boldsymbol{x}-\boldsymbol{b})+\lambda \boldsymbol{z}\right\}, \quad \forall \boldsymbol{z}=\left\{\begin{array}{cc} +1 & x[i]>0 \\ {[-1,+1]} & x[i]=0 \\ -1 & x[i]<0 \end{array}\right. \label{lars_subderiv} \end{align} $$

が0となるとを同時に求める。なお、この時のを決める処理をソフトしきい値処理として知られている（ハードしきい値処理は0にするか非ゼロにするかの二値化処理のこと）。

以下の話は、解パス図においてを大きい方から小さい方へ動かすことを想像するとわかりやすい。

1. において、である。また、このとき

$$ \begin{align} \operatorname{argmin}_{\boldsymbol{x}} f(\boldsymbol{x}) = \boldsymbol{0} \end{align} $$

より、

$$ \begin{align} \partial f(\boldsymbol{x}) &= \boldsymbol{0} \\ -A_{\top}\boldsymbol{b} + \lambda \boldsymbol{z}_{\lambda} &= \boldsymbol{0} \\ \boldsymbol{z}_{\lambda} &= A^{\top}\boldsymbol{b} / \lambda \label{z_inf} \end{align} $$

が導かれる。すなわち、ならば、であることが保証され、式\eqref{lars_subderiv}より最適解はのままであることがわかる。

さらにを小さくしていき、

$$ \begin{align} \lambda = ||A^{\top} \boldsymbol{b}||_{\infty} \end{align} $$

となったとき、これよりさらにをさらに小さくしようとすると（すなわちの要素の最大値よりも小さくしようとすると）、のままでは式\eqref{lars_subderiv}のの要請を満たさない。そこで、

$$ \begin{align} z_{\lambda}[i] = \operatorname{sign} (x_{\lambda}[i]) \end{align} $$

を導入して、改めてを考えると

$$ \begin{align} 0=& \boldsymbol{a}_{i}^{\boldsymbol{T}}(\boldsymbol{a} x-\boldsymbol{b})+\lambda z_{\lambda}[i] \\ & \Rightarrow x_{\lambda}[i]=\frac{\boldsymbol{a}_{i}^{\top} \boldsymbol{b}-\lambda z_{\lambda}[i]}{\boldsymbol{a}_{i}^{\top} \boldsymbol{a}_{i}}=\frac{\boldsymbol{a}_{i}^{\top} \boldsymbol{b}-\lambda \operatorname{sign}\left(x_{\lambda}[i]\right)}{\boldsymbol{a}_{i}^{\top} \boldsymbol{a}_{i}} \end{align} $$

このとき、におけるの条件式\eqref{z_inf}は

$$ \begin{align} A^{\top}(A\boldsymbol{x} - \boldsymbol{b}) + \lambda \boldsymbol{z} &= 0 \\ \boldsymbol{z} &= \frac{A^{\top}(\boldsymbol{b} - A\boldsymbol{x})}{\lambda} \end{align} $$

という式に置き換えられている。

さらにの値を減らしていくと、現在の解はより、

$$ \begin{align} A^{\top}(A\boldsymbol{x} - \boldsymbol{b}) + \lambda \boldsymbol{z} &= 0 \\ \boldsymbol{x} = (A^{\top}A)^{-1} (A^{\top}\boldsymbol{b} - \lambda\boldsymbol{z}) \end{align} $$

と書け、サポートを使って表せば

$$ \begin{align} \left\{\begin{array}{c} \boldsymbol{x}^{S}_{\lambda} = (A_{S}^{\top}A_{S})^{-1} (A_{S}^{\top}\boldsymbol{b} - \lambda\boldsymbol{z}_{\lambda}^{S}) \\ \boldsymbol{z}_{\lambda}^{S} = \operatorname{sign} (\boldsymbol{x}_{\lambda}^{S}) \\ \end{array}\right. \end{align} $$

と書ける。

を小さくしていく過程で、のサポート外のが]の範囲から出ることがある。この場合はこの要素をサポートに追加し、解を更新する。逆に、解の非ゼロだった要素が0になってしまうことがある。この場合は該当するサポートを除いて解を更新する。

ここで気をつけるのは、サポートの追加、除外は１つずつ行っている事である。複数同時に更新しなければいけない状況ではアルゴリズムに修正が必要とのこと。

以上がLARS（least angle regression stagewise; 最小角回帰）の一連の流れである。なお、LARSは解の経路がの大域的最適解になることが保証されている。　

書籍ではIRLSとLARSの性能比較を行っており、低次元ではLARSの性能が圧倒的である。ただ、LARSの計算オーダーがなので、高次元では使い物にならないとも言っている。

書籍では得られる近似解について、OMPとLARSの比較も行っている。

OMPはサポートの再現率が良く、LARSは誤差が小さく精度が良くなる傾向があるとのこと。また、LARSでは正則化を行っているので、非ゼロ要素の値を小さく見積もる傾向がある。選ばれたサポートだけを使って改めて最小二乗法を用いるという方法もあるらしい。

*1:Michael Elad（2010）「スパースモデリング」共立出版

ガウス過程を理解したい２

ガウス過程と機械学習の第5章、補助変数法でよくわからなくなってしまったので原著を読んだのでまとめます。スパース近似ガウス過程回帰を統一的視点で眺めるという内容になっています。

@matsueushiさんが既に綺麗にまとめてくれているのですが、数弱の自分は式をパッと見ただけでは理解できなかったので、もうちょっと細かいところもまとめていこうと思います（ほぼ翻訳）。

ガウス過程回帰

まずガウス過程回帰についておさらいです。ガウス分布の定義は

任意次元の入力に対する出力の同時分布が多変量ガウス分布に従うもののこと。

でした（厳密な言い方では無いと思いますが）。ガウス過程をベイズ的に眺めると、これから紹介する補助変数法の理解がすっきりします。

導入の詳細は前回の記事を参照してください。

入力と、対応する出力が次のような関係だとします。

\begin{aligned} y_{i}=f\left(\mathbf{x}_{i}\right)+\varepsilon_{i}, \quad \text { where } \quad \varepsilon_{i} \sim \mathcal{N}\left(0, \sigma_{\text {noise }}^{2}\right) \end{aligned}

はノイズを表し はノイズの分散です。以降ではこの関係が前提であるとします。

では、ガウス過程回帰をベイズ的に眺めていきます。ガウス過程回帰とは結局のところ、関数について次のような事前分布を導入することと同義です。

\begin{aligned} p\left(\mathbf{f} | \mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{n}\right)=\mathcal{N}(\mathbf{0}, K) \end{aligned}

ここで、行列は成分がとなるカーネル行列です。つまり、関数自体が確率変数であると考えるのがガウス過程です。そして、データからその事後分布を推定するのがガウス過程回帰ということになります。
　では、関数の事後分布を導出していきます。ベイズの定理より、

\begin{aligned} p\left(\mathbf{f}, \mathbf{f}_{*} | \mathbf{y}\right) &= \frac{p\left(\mathbf{f}, \mathbf{f}_{*}\right) p(\mathbf{y} | \mathbf{f})}{p(\mathbf{y})} \\ p\left(\mathbf{f}_{*} | \mathbf{y}\right) &= \int p\left(\mathbf{f}, \mathbf{f}_{*} | \mathbf{y}\right) \mathrm{d} \mathbf{f} \\ &= \frac{1}{p(\mathbf{y})} \int p(\mathbf{y} | \mathbf{f}) p\left(\mathbf{f}, \mathbf{f}_{*}\right) \mathrm{d} \mathbf{f} \end{aligned}

ここで、

\begin{aligned} p\left(\mathbf{f}, \mathbf{f}_{*}\right | \mathbf{y}) &= \mathcal{N}\left(\mathbf{0},\left[\begin{array}{ll} K_{\mathrm{f}, \mathrm{f}} & K_{*, \mathrm{f}} \\ K_{\mathrm{f}, *} & K_{*, *} \end{array}\right]\right) \\ p(\mathbf{y} | \mathbf{f}) &= \mathcal{N}\left(\mathbf{f}, \sigma_{\mathrm{noise}}^{2} I\right) \end{aligned}

です。条件付き多変量ガウス分布の式から

\begin{aligned} p\left(\mathbf{f}_{*} | \mathbf{y}\right)=\mathcal{N}\left(K_{*, \mathbf{f}}\left(K_{\mathbf{f}, \mathbf{f}}+\sigma_{\text {noise }}^{2} I\right)^{-1} \mathbf{y}, K_{*, *}-K_{*, \mathbf{f}}\left(K_{\mathbf{f}, \mathbf{f}}+\sigma_{\text {noise }}^{2} I\right)^{-1} K_{\mathbf{f}, *}\right) \end{aligned}

となることがわかります。

もちろんこれを直接計算できることに越したことはないのですが、カーネル行列の計算量はです。逆行列に至ってはの計算量とメモリが必要になり、データが多いと計算量が爆発し使い物にならないという問題が発生します。そこで、同時分布を近似しようというのが補助変数法です。ガウス過程回帰だけでなく殆どのスパース近似で用いられる方法とのことです。

補助変数法（Inducing Variable Methodl; IVM）

計算量を減らす近似の手段として、補助変数（inducing variable）を同時分布に導入します。ガウス過程において補助変数を導入するということは、対応する補助入力（inducing inputs）を導入したこと同義です。事後分布では積分消去されるのがミソですが、ただ積分消去するだけでは計算量が増えるだけになってしまいます。そこで、の導入に際しての条件付き独立を仮定し、以下のような近似を行います。

\begin{aligned} p\left(\mathbf{f}_{*}, \mathbf{f}\right) \simeq q\left(\mathbf{f}_{*}, \mathbf{f}\right)=\int q\left(\mathbf{f}_{*} | \mathbf{u}\right) q(\mathbf{f} | \mathbf{u}) p(\mathbf{u}) \mathrm{d} \mathbf{u} \end{aligned}

このを決めようというのが補助変数法です。

また、登場する式のリファレンスとしても便利なので、厳密な条件付き分布についても一度整理しておきましょう。厳密な条件付き分布は

\begin{aligned} p(\mathbf{f} | \mathbf{u}) &=\mathcal{N}\left(K_{f, u} K_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \;\; K_{\mathbf{f}, \mathbf{f}}-Q_{\mathbf{f}, \mathbf{f}}\right) \\ p\left(\mathbf{f}_{*} | \mathbf{u}\right) &=\mathcal{N}\left(K_{*, \mathbf{u}} K_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \;\;K_{*, *}-Q_{*, *}\right)　\\ Q_{\mathbf{a}, \mathbf{b}} &\triangleq K_{\mathbf{a}, \mathbf{u}} K_{\mathbf{u}, \mathbf{u}}^{-1} K_{\mathbf{u}, \mathbf{b}} \end{aligned}

です（に関してはこの後導出します）。

The Subset of Data (SoD) Approximation

データの一部だけを使う方法です。データが多すぎるなら削ればよいということです。これは近似と呼ぶのかわかりませんが、計算量を減らす最も簡単な手法ですね。ただ、せっかく取れているデータを使わないのはもったいありません。以降で紹介する4つ（説明するのは3つ）の方法では、データを捨てずに事後分布を近似していきます。

The Subset of Regressors (SoR) Approximation

SoRは以下の事前重みを持つ線形モデルと解釈することができます。

\begin{aligned} f_* = K_{*, \mathbf{u}}\mathbf{w_u}, \quad \text{with} \quad p(\mathbf{w_u}) = \mathcal{N}\left(\mathbf{0}, K_{\,\mathbf{u, u}}^{-1}\right) \end{aligned}

事前分布の分散をとしているところが重要で、ガウス過程の事前分布が自然と現れます。

\begin{aligned} \mathbf{u}=K_{\mathbf{u}, \mathbf{u}} \mathbf{w}_{\mathbf{u}} \Rightarrow\left\langle\mathbf{u} \mathbf{u}^{\top}\right\rangle= K_{\mathbf{u}, \mathbf{u}}\left\langle\mathbf{w}_{\mathbf{u}} \mathbf{w}_{\mathbf{u}}^{\top}\right\rangle K_{\mathbf{u}, \mathbf{u}}=K_{\mathbf{u}, \mathbf{u}} \end{aligned}

したがって、より、

\begin{aligned} \mathbf{f}_{*}=K_{*, \mathbf{u}} K_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \quad \text { with } \quad \mathbf{u} \sim \mathcal{N}\left(\mathbf{0}, K_{\mathbf{u}, \mathbf{u}}\right) \end{aligned}

となり、事後分布と見比べると自然な形で式が得られたことがわかります。以上より、は

\begin{aligned} q_{\mathrm{SoR}}(\mathbf{f} | \mathbf{u}) &= \mathcal{N}\left(K_{\mathrm{f}, \mathbf{u}} K_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{0}\right) \\ q_{\mathrm{SoR}}\left(\mathbf{f}_{*} | \mathbf{u}\right) &= \mathcal{N}\left(K_{*, \mathbf{u}} K_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{0}\right) \end{aligned}

で与えられ。

\begin{aligned} q_{\mathrm{SoR}}\left(\mathbf{f}, \mathbf{f}_{*}\right)&=\mathcal{N}\left(\mathbf{0},\left[\begin{array}{ll} Q_{\mathrm{f}, \mathrm{f}} & Q_{\mathrm{f}, *} \\ Q_{*, f} & Q_{*, *} \end{array}\right]\right) \\ Q_{\mathbf{a}, \mathbf{b}} &\triangleq K_{\mathbf{a}, \mathbf{u}} K_{\mathbf{u}, \mathbf{u}}^{-1} K_{\mathbf{u}, \mathbf{b}} \end{aligned}

となることがわかります。このも一応導出しておきましょう。

\begin{aligned} cov\left(q_{\mathrm{SoR}}(\mathbf{f} | \mathbf{u}), q_{\mathrm{SoR}}(\mathbf{f}_* | \mathbf{u})\right) &= K_{\mathbf{f,u}}K_{\mathbf{u,u}}^{-1}\mathbf{u}\left(K_{\mathbf{f,u}}K_{\mathbf{u,u}}^{-1}\mathbf{u}\right)^{\top} \\ &= K_{\mathbf{f,u}}K_{\mathbf{u,u}}^{-1}\mathbf{u}\mathbf{u}^{\text{T}}K_{\mathbf{u,u}}^{-1}K_{\mathbf{u,*}} \\ &= K_{\mathbf{f,u}}K_{\mathbf{u,u}}^{-1}K_{\mathbf{u,*}} \end{aligned}

また、近似された事後分布は、の分散共分散行列の対角成分にノイズを足し、条件付き多変量ガウス分布の式を使えば

\begin{aligned} q_{\mathrm{SoR}}\left(\mathbf{f}_{*} | \mathbf{y}\right) &=\mathcal{N}\left(Q_{*, \mathbf{f}}\left(Q_{\mathbf{f}, \mathbf{f}}+\sigma_{\text {noise }}^{2} I\right)^{-1} \mathbf{y}, Q_{* *}-Q_{*, \mathbf{f}}\left(Q_{\mathbf{f}, \mathbf{f}}+\sigma_{\text {noise }}^{2} I\right)^{-1} Q_{\mathbf{f}, *}\right)\\ &=\mathcal{N}\left(\sigma^{-2} K_{*, u} \Sigma K_{\mathbf{u}, \mathbf{f}},\;\; K_{*, u} \Sigma K_{\mathbf{u}, *}\right) \\ \Sigma &= \left(\sigma^{-2} K_{\mathbf{u}, \mathbf{f}} K_{\mathbf{f}, \mathbf{u}}+K_{\mathbf{u}, \mathbf{u}}\right)^{-1} \end{aligned}

となります。の2つ目の式の形は実装するときにメモリに優しい形になっています（導出は最後にやります）。理解するだけならの1つ目の式がわかれば十分です。式をよく見るとわかるのですが、SoRは結局のところカーネル関数をと置いたガウス過程に対応します（ここ重要）。

さて、自然な形でが決められたわけですが、この近似によって生じるデメリットもあります。

メリット

厳密にガウス過程（これから厳密なガウス過程ではなくなっていくのでメリットとした）

デメリット

たまたま得られたデータのばらつきが小さかった時に予測分散が過剰に小さくなるなど、補助変数を介する条件付き分布を導入したことによる自由度の制約を受けます。

以下の手法はこのデメリットを緩和しようとするもので、分散が過剰に小さくならないよう、さらに前提を追加していきます（理論的には破綻するので厳密にはガウス過程ではなくなる）。

The Deterministic Training Conditional (DTC) Approximation

DTCは、以下のような近似を行います。

\begin{aligned} q_{\mathrm{DTC}}(\mathbf{f} | \mathbf{u}) &= \mathcal{N}\left(K_{\mathrm{f}, \mathbf{u}} K_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{0}\right) \\ q_{\mathrm{DTC}}\left(\mathbf{f}_{*} | \mathbf{u}\right) &= p\left(\mathbf{f}_{*} | \mathbf{u}\right) \end{aligned}

式を見ればわかる通り、SoRとの本質的な違いはとしている部分だけです。すなわち、近似された同時分布は

\begin{aligned} q_{\mathrm{DTC}}\left(\mathbf{f}, \mathbf{f}_{*}\right)=\mathcal{N}\left(\mathbf{0},\left[\begin{array}{ll} Q_{\mathrm{f}, f} & Q_{\mathrm{f}, *} \\ Q_{*, \mathrm{f}} & K_{*, *} \end{array}\right]\right) \end{aligned}

となります。の中にが混ざっていることからもわかる通り、この時点で厳密にはガウス過程ではなくなっています。あとは、SoRの時と同様に事後分布を計算すれば、

\begin{aligned} q_{\mathrm{DTC}}\left(\mathbf{f}_{*} | \mathbf{y}\right) &=\mathcal{N}\left(Q_{*, \mathbf{f}}\left(Q_{\mathbf{f}, \mathbf{f}}+\sigma_{\text {noise }}^{2} I\right)^{-1} \mathbf{y}, K_{*, *}-Q_{*, \mathbf{f}}\left(Q_{\mathbf{f}, \mathbf{f}}+\sigma_{\text {noise }}^{2} I\right)^{-1} Q_{\mathbf{f}, *}\right.\\ &=\mathcal{N}\left(\sigma^{-2} K_{*, u} \Sigma K_{\mathbf{u}, \mathbf{f}} \mathbf{y},\;\; K_{*, *}-Q_{*, *}+K_{*, \mathbf{u}} \Sigma K_{*, \mathbf{u}}^{\top}\right) \end{aligned}

となります。SoRの事後分布と見比べると、分散共分散行列の式中のがになっていることがわかります。が正定値であることから、SoRよりもDTCの法が予測分の分散が大きくなります。

メリット

SoRよりも予測分散が大きくなることが保証される。

デメリット

厳密なガウス過程ではなくなる。

The Fully Independent Training Conditional (FITC) Approximation

スパースガウス過程とも呼ばれ、以下のように近似を行います。

\begin{aligned} q_{\mathrm{FTTC}}(\mathbf{f} | \mathbf{u}) &= \prod_{i=1}^{n} p\left(f_{i} | \mathbf{u}\right) = \mathcal{N}\left(K_{\mathrm{f}, \mathbf{u}} K_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \operatorname{diag}\left[K_{\mathrm{f}, \mathrm{f}}-Q_{\mathrm{f}, \mathrm{f}}\right]\right)\\ q_{\mathrm{FITC}}\left(f_{*} | \mathbf{u}\right) &= p\left(f_{*} | \mathbf{u}\right) \end{aligned}

DTCと異なるのは、の分散共分散行列の部分です。FITCがDTCより優れている点は、と見比べるとよくわかります。FITCは対角成分以外を捨ててはいるものの、厳密なガウス過程の分散共分散行列の情報を使って近似を行っています（それでもとしてるので厳密なガウス過程ではない）。

後は今まで通り計算していくと、

\begin{aligned} q_{\mathrm{FITC}}\left(\mathbf{f}, f_{*}\right) &= \mathcal{N}\left(\mathbf{0},\left[\begin{array}{cc} Q_{\mathrm{f}, \mathbf{f}}-\operatorname{diag}\left[Q_{\mathrm{f}, \mathrm{f}}-K_{\mathrm{f}, \mathrm{f}}\right] & Q_{\mathrm{f}, *} \\ Q_{*, \mathrm{f}} & K_{*, *} \end{array}\right]\right) \\ q_{\mathrm{FITC}}\left(f_{*} | \mathbf{y}\right) &=\mathcal{N}\left(Q_{*, \mathbf{f}}\left(Q_{\mathrm{f}, \mathbf{f}}+\Lambda\right)^{-1} \mathbf{y}, K_{*, *}-Q_{*, \mathbf{f}}\left(Q_{\mathrm{f}, \mathbf{f}}+\Lambda\right)^{-1} Q_{\mathrm{f}, *}\right) \\ &=\mathcal{N}\left(K_{*, \mathbf{u}} \Sigma K_{\mathbf{u}, \mathbf{f}} \Lambda^{-1} \mathbf{y}, K_{*, *}-Q_{*, *}+K_{*, \mathbf{u}} \Sigma K_{\mathbf{u}, *}\right) \\ \Sigma &= \left(K_{\mathbf{u}, \mathbf{u}}+K_{\mathbf{u}, \mathbf{f}} \Lambda^{-1} K_{\mathbf{f}, \mathbf{u}}\right)^{-1} \\ \Lambda &=\operatorname{diag} \left[K_{\mathrm{f}, \mathrm{f}}-Q_{\mathrm{f}, \mathrm{f}}+\sigma_{\mathrm{noise}}^{2} I\right] \end{aligned}

となります。

The Partially Independent Training Conditional (PITC) Approximation

書籍「ガウス過程と機械学習」の取り扱い範囲を超えているので今回は取り扱いません。ただ、ここまで読み進めてくれた方ならわかると思いますが、さっきは相関を捨てていましたが、ある程度相関を拾ってあげることで近似精度を上げようという方法です。

最後に

メモリに優しい形の式を導出します。見やすさのために式を再掲しておきます。

\begin{aligned} q_{\mathrm{SoR}}\left(\mathbf{f}_{*} | \mathbf{y}\right) &=\mathcal{N}\left(Q_{*, \mathbf{f}}\left(Q_{\mathbf{f}, \mathbf{f}}+\sigma_{\text {noise }}^{2} I\right)^{-1} \mathbf{y}, Q_{* *}-Q_{*, \mathbf{f}}\left(Q_{\mathbf{f}, \mathbf{f}}+\sigma_{\text {noise }}^{2} I\right)^{-1} Q_{\mathbf{f}, *}\right) \\ &=\mathcal{N}\left(\sigma^{-2} K_{\mathbf{*, u}} \Sigma K_{\mathbf{u,f}},\;\; K_{\mathbf{*, u}} \Sigma K_{\mathbf{u, *}}\right) \\ \end{aligned}

[導出]

上記で与えられたを使って

\begin{aligned} \Sigma &= \left(\sigma^{-2} K_{\mathbf{u}, \mathbf{f}} K_{\mathbf{f}, \mathbf{u}}+K_{\mathbf{u}, \mathbf{u}}\right)^{-1} \\ \Sigma^{-1} &= \sigma^{-2} K_{\mathbf{u}, \mathbf{f}} K_{\mathbf{f}, \mathbf{u}}+K_{\mathbf{u}, \mathbf{u}} \\ \sigma^{2} \Sigma^{-1} &= K_{\mathbf{u}, \mathbf{f}} K_{\mathbf{f}, \mathbf{u}}+\sigma^{2} K_{\mathbf{u}, \mathbf{u}} \\ \sigma^{2} \Sigma^{-1} K_{\mathbf{u}, \mathbf{u}}^{-1} &= K_{\mathbf{u}, \mathbf{f}} K_{\mathbf{f}, \mathbf{u}} K_{\mathbf{u}, \mathbf{u}}^{-1} + \sigma^{2} I \\ \sigma^{2} \Sigma^{-1} K_{\mathbf{u}, \mathbf{u}}^{-1} K_{\mathbf{u}, \mathbf{f}} &= K_{\mathbf{u}, \mathbf{f}} K_{\mathbf{f}, \mathbf{u}} K_{\mathbf{u}, \mathbf{u}}^{-1} K_{\mathbf{u}, \mathbf{f}} + \sigma^{2} K_{\mathbf{u}, \mathbf{f}} \\ \sigma^{2} \Sigma^{-1} K_{\mathbf{u}, \mathbf{u}}^{-1} K_{\mathbf{u}, \mathbf{f}} &= K_{\mathbf{u}, \mathbf{f}} Q_{\mathbf{f}, \mathbf{f}} + \sigma^{2} K_{\mathbf{u}, \mathbf{f}} \\ \sigma^{2} \Sigma^{-1} K_{\mathbf{u}, \mathbf{u}}^{-1} K_{\mathbf{u}, \mathbf{f}} &= K_{\mathbf{u}, \mathbf{f}} \left( Q_{\mathbf{f}, \mathbf{f}} + \sigma^{2} I \right) \\ K_{\mathbf{u}, \mathbf{u}}^{-1} K_{\mathbf{u}, \mathbf{f}} &= \sigma^{-2} \Sigma K_{\mathbf{u}, \mathbf{f}} \left( Q_{\mathbf{f}, \mathbf{f}} + \sigma^{2} I \right) \\ \end{aligned}

くどすぎますが、徹底的に式変形を書いておきました。

これを使って を変形していきましょう。

\begin{aligned} Q_{*, \mathbf{f}}\left(Q_{\mathbf{f}, \mathbf{f}}+\sigma^{2} I\right)^{-1} &= K_{*, \mathbf{u}} K_{\mathbf{u, u}}^{-1} K_{\mathbf{u,f}} \left(Q_{\mathbf{f}, \mathbf{f}}+\sigma^{2} I \right)^{-1} \\ &= K_{\mathbf{*,u}} \sigma^{-2} \Sigma K_{\mathbf{u}, \mathbf{f}}\left(Q_{\mathbf{f}, \mathbf{f}}+\sigma^{2} I\right) \left(\mathbf{Q}_{\mathbf{f}, \mathbf{f}}+\sigma^{2} \mathbf{I}\right)^{-1} \\ &=\sigma^{-2} K_{*, \mathbf{u}} \Sigma K_{\mathbf{u}, \mathbf{f}} \end{aligned}

では分散共分散行列も導出していきましょう。

\begin{aligned} Q_{*,*}-Q_{*, \mathbf{f}}\left(Q_{\mathbf{f}, \mathbf{f}}+\sigma_{\text {noise }}^{2} I\right)^{-1} Q_{\mathbf{f}, *} \end{aligned}

において、は既に求めたので、

\begin{aligned} Q_{*,*}-Q_{*, \mathbf{f}}\left(Q_{\mathbf{f}, \mathbf{f}}+\sigma_{\text {noise }}^{2} I\right)^{-1} Q_{\mathbf{f}, *} &= K_{\mathbf{*,u}} K_{\mathbf{u,u}}^{-1} K_{\mathbf{u,*}} -\sigma^{-2} K_{\mathbf{*,u}} \Sigma K_{\mathbf{u,f}} K_{\mathbf{f,u}} K_{\mathbf{u,u}}^{-1} K_{\mathbf{u,*}} \\ &= K_{\mathbf{*,u}} \left\{I - \sigma^{-2} \Sigma K_{\mathbf{u,f}} K_{\mathbf{f,u}} \right\} K_{\mathbf{u,u}}^{-1} K_{\mathbf{u,*}} \\ &= K_{\mathbf{*,u}} \Sigma \left\{\Sigma^{-1} - \sigma^{-2} K_{\mathbf{u,f}} K_{\mathbf{f,u}} \right\} K_{\mathbf{u,u}}^{-1} K_{\mathbf{u,*}} \\ &= K_{\mathbf{*,u}} \Sigma \left\{\sigma^{-2} K_{\mathbf{u,f}} K_{\mathbf{f,u}} + K_{\mathbf{u,u}} - \sigma^{-2} K_{\mathbf{u,f}} K_{\mathbf{f,u}} \right\} K_{\mathbf{u,u}}^{-1} K_{\mathbf{u,*}} \\ &= K_{\mathbf{*,u}} \Sigma \left\{K_{\mathbf{u,u}} \right\} K_{\mathbf{u,u}}^{-1} K_{\mathbf{u,*}} \\ &= K_{\mathbf{*,u}} \Sigma K_{\mathbf{u,*}} \end{aligned}

と導けます。以上でメモリに優しい式の形を導出できました。

LaTeX数式からはてなブログ数式へ一発で変換する

勉強するときはlatexまたはmarkdownでまとめているのですが、これをはてなブログへ転記しようと思うと相当に面倒くさい（特に数式周りが）。

数式の記述方法については以下のブログ記事でとてもよくまとめてくれている。

はてなブログで数式を書く - 七誌の開発日記

数式に関して少し発見があったのでメモ。

上記まとめブログの「裏技」にて、

【注意】2020年2月4日現在、この方法は使えなくなったようです。 どこか一箇所で [ tex: ～ ] を使用すれば、独立した数式に $$ ～ $$ や [ ～ ] を使用できるようです。

とあるのですが、2020年2月29日現在は [ tex: ～ ] を書いておけば一応$$~$$でも数式は表示されました。ただ、$$が残ってしまうという問題がありました。

そこで、$$を削除したところ普通に数式表示されました。推奨される方法ではないとは思いますが、便利なので使っていこうかと思います。

自動化した

Markdown→はてなブログ記事へのドキュメント変換スクリプトをPythonで書きました。 使い方は

python <path_to_script> <path_to_markdown>

です。とりあえず自分が使う用に書いた（今書いている別の記事の為）ので汚い&必要最低限の変換しか考慮していないですが、数式変換には便利かと思います。

import argparse
import os
import re

parser = argparse.ArgumentParser()
parser.add_argument("path", help="text file including latex equations.")

args = parser.parse_args()
path = args.path

def convert(text):
    link_indices = [m.span() for m in re.finditer("\[[あ-んぁ-んァ-ン一-龥a-zA-Z0-9!-/:-@[-`{-~]*\]([a-zA-Z0-9!-/:-@[-`{-~]*)", text)][::-1]
    links = []
    for i,idx in enumerate(link_indices):   # 逆順になっていることに注意
        links.append(text[idx[0]:idx[1]])
        text = text[:idx[0]] + "@@@" + text[idx[1]:]

    text = text.replace("\\\\", "\\\\\\")
    text = text.replace("_", "\\_")
    text = text.replace("^", "\\^")
    text = text.replace("[", "\\[")
    text = text.replace("]", "\\]")
    text = text.replace("\\{", "\\\\{")
    text = text.replace("\\}", "\\\\}")
    text = text.replace("*", "\\*")

    # equations
    text = text.replace(r"$$", r"")
    N = len([m.span() for m in re.finditer("$$", text)][::-1])
    if N > 0:
        text = """[tex: \displaystyle]
""" + text
    
    indices = [m.span() for m in re.finditer("\$", text)][::-1]
    for i,idx in enumerate(indices):   # 逆順になっていることに注意
        if i % 2 == 0:
            # $ -> ]
            text = text[:idx[0]] + "]" + text[idx[1]:]
        elif i % 2 == 1:
            text = text[:idx[0]] + "[tex: \\displaystyle " + text[idx[1]:]


    link_indices = [m.span() for m in re.finditer("@@@", text)][::-1]
    for i,idx in enumerate(link_indices):   # 逆順になっていることに注意
        text = text[:idx[0]] + links[i] + text[idx[1]:]

    return text


with open(path) as f:
    text = f.read()

newtext = convert(text)

with open("hatena.txt", "w") as f:
    f.write(newtext)

数式表示には上記の方法を利用しています。

ガウス過程を理解したい１

時間が空きましたが、ガウス過程を理解するために多変量ガウス分布について学習してきました（以下を参照）。

hotoke-x.hatenablog.com

hotoke-x.hatenablog.com

ガウス過程と機械学習の第5章、補助変数法でよくわからなくなってしまったので原著を読みました。

原著論文の紹介をいきなりしても良いのですが、論文のモチベーションがわかるようにまずはガウス過程の導入から行っていきます。

登場する用語の整理

ガウス過程に入る前に、用語の整理をしておきます。

1. 動径基底関数

   $$ \begin{align} \exp \left(- \frac{\left(x - \mu \right)^{2}}{\sigma^{2}} \right) \end{align} $$

   のような形の関数のこと。（オイラーの公式で三角関数に分解できることを想像すると、名前に納得感があると思った。）

2. 動径基底関数（radial basis function, RBF）回帰
   線形回帰モデルの基底関数に、動径基底関数を用いたモデル。

   $$ \begin{align} y = \sum_{h=-H}^{H} w_{h} \exp \left(- \frac{\left(x - \mu_h \right)^{2}}{\sigma^{2}} \right) \label{reg:rbf} \end{align} $$

   素直にこの式に従うとデータが低次元じゃないと計算できません（次元の呪い）。たとえば、1次元なら任意の曲線を表現するためにの値を細かくとれば良いですが、2次元で任意の局面を表現するためには細かいメッシュを切る必要があります。

ガウス過程とは

任意次元の入力に対する出力の同時分布が多変量ガウス分布に従うもののこと。

ガウス過程のポイント

1. に関数を導入したことで、線形回帰モデルの任意の点（つまり無限点）の基底関数を表現可能になった
2. 式\eqref{reg:rbf}において、について全範囲]を解析的に積分することで、観測された次元に限って他を積分消去（周辺化）してしまう（扱うのは有限次元で済む）
3. 周辺化によってデータ点がある次元だけで事後分布を表現できる

無限点とは何のことか、周辺化が何をしているのかのイメージをつかみたければガウス過程と機械学習 p66の図3.5を参照。ちゃんと理解したい人は再生核ヒルベルト空間について調べましょう。

ちゃんと式で考えていきます。

ちゃんと式で考えていきます。入力に対して

$$ \begin{align} \left(\begin{array}{c}{\widehat{y}_{1}} \\ {\widehat{y}_{2}} \\ {\vdots} \\ {\widehat{y}_{N}}\end{array}\right)&=\left(\begin{array}{cccc}{\phi_{0}\left(\boldsymbol{x}_{1}\right)} & {\phi_{1}\left(\boldsymbol{x}_{1}\right)} & {\cdots} & {\phi_{H}\left(\boldsymbol{x}_{1}\right)} \\ {\phi_{0}\left(\boldsymbol{x}_{2}\right)} & {\phi_{1}\left(\boldsymbol{x}_{2}\right)} & {\cdots} & {\phi_{H}\left(\boldsymbol{x}_{2}\right)} \\ {\vdots} & {\vdots} & {\ddots} & {\vdots} \\ {\phi_{0}\left(\boldsymbol{x}_{N}\right)} & {\phi_{1}\left(\boldsymbol{x}_{N}\right)} & {\cdots} & {\phi_{H}\left(\boldsymbol{x}_{N}\right)}\end{array}\right)\left(\begin{array}{c}{w_{0}} \\ {w_{1}} \\ {\vdots} \\ {w_{H}}\end{array}\right) \\ \boldsymbol{\widehat y} &= \boldsymbol{\Phi w} \end{align} $$

とすると、

$$ \begin{align} \mathbb{E}\left[\boldsymbol{y}\right] &= \mathbb{E}\left[\boldsymbol{\Phi w}\right] = \boldsymbol{\Phi} \mathbb{E}\left[\boldsymbol{w}\right] = \boldsymbol{0} \\ \Sigma &=\mathbb{E}\left[\boldsymbol{y y}^{\mathrm T}\right]-\mathbb{E}[\boldsymbol{y}] \mathbb{E}[\boldsymbol{y}]^{\mathrm T} \\ &= \mathbb{E}\left[(\boldsymbol{\Phi} \boldsymbol{w})(\boldsymbol{\Phi} \boldsymbol{w})^{\mathrm T}\right] \\ &= \boldsymbol{\Phi} \mathbb{E}\left[\boldsymbol{w} \boldsymbol{w}^{\mathrm T}\right] \boldsymbol{\Phi}^{\mathrm T} = \lambda^{2} \boldsymbol{\Phi} \Phi^{\mathrm T} \end{align} $$

より、

$$ \begin{align} y \sim \mathcal{N} \left(\boldsymbol{0}, \lambda^{2}\boldsymbol{\Phi \Phi^{\mathrm T}} \right) \end{align} $$

となります。ここでは、としました。について知らなくてよい点で強力です。

分散共分散行列（今回のようなの形の行列をグラム行列ともいう）の成分を見てみると、

$$ \begin{align} \boldsymbol{K}_{n, n'} = \lambda^{2} \phi(\boldsymbol{x_{n}})^{\mathrm T}\phi(\boldsymbol{x_{n'}}) \end{align} $$

となっており、基底関数の内積の定数倍となっています。は分散共分散行列なので、正則な対称行列であることに注意しましょう。

ちなみに、データ点が増えすぎると計算量とメモリが爆発するので一般には計算方法を工夫する必要があります。その手法として補助変数法があります（後述）。

カーネルトリック

上記の内積を、とのカーネル関数（kernel function）と呼び、以下のように定義する。

$$ \begin{align} k(\boldsymbol{x}_{n}, \boldsymbol{x}_{n '}) = \phi(\boldsymbol{x_{n}})^{\mathrm T}\phi(\boldsymbol{x_{n'}}) \end{align} $$

は積分消去されましたが、このままだと特徴ベクトルが無限次元の場合表現できません。 そこで、を直接表現することを避け、内積の結果の式だけを利用します（カーネルトリック）。

ガウス過程回帰

カーネル関数

以下のような観測モデルを考えます。

$$ \begin{align} \left\{\begin{array}{l} y_n &= f\left(\boldsymbol{x_n}\right) + \epsilon_n \\ \epsilon &\sim \mathcal{N}(0, \sigma^2) \end{array}\right. \end{align} $$

すると、が与えられたとき、の確率分布は

$$ \begin{align} p(\mathbf{y} | \mathbf{X}) &=\int p(\mathbf{y}, \mathbf{f} | \mathbf{X}) d \mathbf{f} \\ &=\int p(\mathbf{y} | \mathbf{X}) p(\mathbf{f} | \mathbf{X}) d \mathbf{f} \\ &=\int p(\mathbf{y} | \mathbf{f}) p(\mathbf{f} | \mathbf{X}) d \mathbf{f} \\ &=\int \mathcal{N}\left(\mathbf{y} | \mathbf{f}, \sigma^{2} \mathbf{I}\right) \mathcal{N}(\mathbf{f} | \boldsymbol{\mu}, \mathbf{K}) d \mathbf{f} \end{align} $$

なお、を利用しました。ガウス分布の再生性より

$$ \begin{align} p(\mathbf{y} | \mathbf{X}) &= \int \mathcal{N}\left(\boldsymbol{\mu},\boldsymbol{K} + \sigma^2 \boldsymbol{I} \right) \end{align} $$

となります。ガウス分布の再生性は積率母関数を使えば容易に証明可能です。以上より、はガウス過程に従い、そのカーネル関数は

$$ \begin{align} k^{\prime}\left(\mathbf{x}_{n}, \mathbf{x}_{n^{\prime}}\right)=k\left(\mathbf{x}_{n}, \mathbf{x}_{n^{\prime}}\right)+\sigma^{2} \delta\left(n, n^{\prime}\right) \end{align} $$

となることがわかります。はクロネッカーのデルタです。

予測（欠損値補完など）

非観測点における出力を予測するためには、を含めたカーネル行列

$$ \begin{align} \left(\begin{array}{l}{\mathbf{y}} \\ {y^{*}}\end{array}\right) \sim \mathcal{N}\left(\mathbf{0}, \left(\begin{array}{cc} {\boldsymbol{K}} & {\boldsymbol{k}_{*}} \\ {\boldsymbol{k}_{*}^{T}} & {k_{* *}} \end{array}\right)\right) \end{align} $$

を用いてをサンプリングすればよいです。ここで、

$$ \begin{align} \left\{\begin{array}{l}{\boldsymbol{k}_{*}=\left(k\left(\boldsymbol{x}^{*}, \boldsymbol{x}_{1}\right), k\left(\boldsymbol{x}^{*}, \boldsymbol{x}_{2}\right), \ldots, k\left(\boldsymbol{x}^{*}, \boldsymbol{x}_{N}\right)\right)^{T}} \\ {k_{* *}=k\left(\boldsymbol{x}^{*}, \boldsymbol{x}^{*}\right)}\end{array}\right. \end{align} $$

です。条件付き多変量ガウス分布の公式より、

$$ \begin{align} \left\{\begin{array}{l}{p\left(y^{*} | \boldsymbol{x}^{*}, \mathcal{D}\right)=\mathcal{N}\left(\boldsymbol{k}_{*}^{T} \boldsymbol{K}^{-1} \boldsymbol{y}, k_{* *}-\boldsymbol{k}_{*}^{T} \boldsymbol{K}^{-1} \boldsymbol{k}_{*}\right)} \\ {\mathcal{D}=\left\{(\boldsymbol{x}_1, y_1), (\boldsymbol{x}_2, y_2), \ldots, (\boldsymbol{x}_N, y_N) \right\}} \end{array}\right. \end{align} $$

が得られます。なお、予測値（欠損値など）が複数点ある場合でも同様です。詳細は書籍p85を参照のしてください。

基本はここまでです。しかし、このままではデータの量が増えるのに伴って計算量が爆発してしまいます。これを避けるために計算を効率化する補助変数法と呼ばれる計算方法があります。

補助変数法については次回まとめます。

参考

1. 持橋 大地、大羽 成征、（2019）「ガウス過程と機械学習」講談社
2. Quinonero-C, ela, J. & Rasmussen, C. E. A Unifying View of Sparse Approximate Gaussian Process Regression. 6, 1939–1959 (2005).

多変量ガウス分布を理解したい２

前回の記事で、多変量ガウス分布の正規化定数を導出しました。

hotoke-x.hatenablog.com

今回は

1. 多変量ガウス分布の周辺化
2. 条件付きガウス分布

を行います。

前回同様、数学の準備をしてから本題に入ります。

数学の準備

ブロック行列

ブロック行列の積は、通常の行列と同様に計算できます。しかしブロック行列の逆行列を求めようとすると一工夫必要です。 覚えておくと便利なので、ついでに行列式を得てから逆行列を求めます。

行列式

をそれぞれブロック行列を構成する要素の行列とするとき、

$$ \begin{align} \left( \begin{array}{cc} A & B \\ C & D \\ \end{array} \right) = \left( \begin{array}{cc} I & O \\ CA^{-1} & I \\ \end{array} \right) \left( \begin{array}{cc} A & O \\ O & D-CA^{-1}B \\ \end{array} \right) \left( \begin{array}{cc} I & A^{-1}B \\ O & I \\ \end{array} \right) \end{align} $$

のように分解できます。導出は以下の通りです。

【導出】
まず、左辺からを消去します。

$$ \begin{align} \left( \begin{array}{cc} A & B \\ C & D \\ \end{array} \right) \left( \begin{array}{cc} I & A^{-1}B \\ O & I \\ \end{array} \right) &= \left( \begin{array}{cc} A & -B+B \\ C & D-CA^{-1}B \\ \end{array} \right) \\ &= \left( \begin{array}{cc} A & O \\ C & D-CA^{-1}B \\ \end{array} \right) \end{align} $$

次にを消去します。

$$ \begin{align} \left( \begin{array}{cc} I & O \\ -CA^{-1}I & I \\ \end{array} \right) \left( \begin{array}{cc} A & O \\ C & D-CA^{-1}B \\ \end{array} \right) &= \left( \begin{array}{cc} A & O \\ O & D-CA^{-1}B \\ \end{array} \right) \end{align} $$

以上より、

$$ \begin{align} \left( \begin{array}{cc} I & O \\ -CA^{-1} & I \\ \end{array} \right) &\left( \begin{array}{cc} A & B \\ C & D \\ \end{array} \right) \left( \begin{array}{cc} I & -A^{-1}B \\ O & I \\ \end{array} \right) = \left( \begin{array}{cc} A & O \\ O & D-CA^{-1}B \\ \end{array} \right) \\ \left( \begin{array}{cc} A & B \\ C & D \\ \end{array} \right) &= \left( \begin{array}{cc} I & O \\ -CA^{-1} & I \\ \end{array} \right)^{-1} \left( \begin{array}{cc} A & O \\ O & D-CA^{-1}B \\ \end{array} \right) \left( \begin{array}{cc} I & -A^{-1}B \\ O & I \\ \end{array} \right)^{-1} \\ &= \left( \begin{array}{cc} I & O \\ CA^{-1} & I \\ \end{array} \right) \left( \begin{array}{cc} A & O \\ O & D-CA^{-1}B \\ \end{array} \right) \left( \begin{array}{cc} I & A^{-1}B \\ O & I \\ \end{array} \right) \end{align} $$

ここで、より、

$$ \begin{align} \left| \begin{array}{cc} A & B \\ C & D \\ \end{array} \right| &= \left| \begin{array}{cc} I & O \\ CA^{-1} & I \\ \end{array} \right| \left| \begin{array}{cc} A & O \\ O & D-CA^{-1}B \\ \end{array} \right| \left| \begin{array}{cc} I & A^{-1}B \\ O & I \\ \end{array} \right| \\ &= \left| \begin{array}{cc} A & O \\ O & D-CA^{-1}B \\ \end{array} \right| \\ &= \left|A\left(D-CA^{-1}B \right) \right| = \det A \det \left(D-CA^{-1}B \right) \end{align} $$

とブロック行列の行列式が求まりました。

逆行列

では、少し寄り道をしたところでブロック行列を分解した式から逆行列を求めていきます。

$$ \begin{align} \left(\begin{array}{cc} A & B \\ C & D \\ \end{array} \right)^{-1} &= \left(\left( \begin{array}{cc} I & O \\ CA^{-1} & I \\ \end{array} \right) \left( \begin{array}{cc} A & O \\ O & D-CA^{-1}B \\ \end{array} \right) \left( \begin{array}{cc} I & A^{-1}B \\ O & I \\ \end{array} \right)\right)^{-1} \\ &= \left( \begin{array}{cc} I & A^{-1}B \\ O & I \\ \end{array} \right)^{-1} \left( \begin{array}{cc} A & O \\ O & D-CA^{-1}B \\ \end{array} \right)^{-1} \left( \begin{array}{cc} I & O \\ CA^{-1} & I \\ \end{array} \right)^{-1} \end{align} $$

表記が大変なので、として、

$$ \begin{align} \left(\begin{array}{cc} A & B \\ C & D \\ \end{array} \right)^{-1} &= \left( \begin{array}{cc} I & -A^{-1}B \\ O & I \\ \end{array} \right) \left( \begin{array}{cc} A^{-1} & O \\ O & S^{-1} \\ \end{array} \right) \left( \begin{array}{cc} I & O \\ -CA^{-1} & I \\ \end{array} \right) \\ &= \left( \begin{array}{cc} A^{-1} & -A^{-1}BS^{-1} \\ O & S^{-1} \\ \end{array} \right) \left( \begin{array}{cc} I & O \\ -CA^{-1} & I \\ \end{array} \right) \\ &= \left( \begin{array}{cc} A^{-1} + A^{-1}BS^{-1}CA^{-1} & -A^{-1}BS^{-1} \\ -S^{-1}CA^{-1} & S^{-1} \\ \end{array} \right) \end{align} $$

さらに、が正則なら、逆行列の補助定理

$$ \begin{align} A^{-1} + A^{-1}BS^{-1}CA^{-1} = \left(A-BD^{-1}C \right)^{-1} \end{align} $$

より、

$$ \begin{align} \left(\begin{array}{cc} A & B \\ C & D \\ \end{array} \right)^{-1} &= \left( \begin{array}{cc} A^{-1} + A^{-1}BS^{-1}CA^{-1} & -A^{-1}BS^{-1} \\ -S^{-1}CA^{-1} & S^{-1} \\ \end{array} \right) \\ &= \left( \begin{array}{cc} \left(A - BD^{-1}C \right) & -A^{-1}B\left(D-CA^{-1}B\right)^{-1} \\ -\left(D-CA^{-1}B\right)^{-1}CA^{-1} & \left(D-CA^{-1}B\right)^{-1} \\ \end{array} \right) \end{align} $$

と書けます。ちなみに多変量ガウス分布を考える場合は、このブロック行列は分散共分散行列に対応するため半正定値です。そのためランク落ちさえしなければは正則になります。

多変量ガウス分布の周辺化

$$ \begin{align} p({\boldsymbol x}_{1})= \int p \left({\boldsymbol x}_{1}, {\boldsymbol x}_{2} \right) {\mathrm d}{\boldsymbol x}_{2} = {\mathcal N} \left({\boldsymbol x}_{1}, \Sigma_{1} \right) \end{align} $$

を確かめます。ここで、は、

$$ \begin{align} \left(\begin{array}{c} {\boldsymbol x}_{1} \\ {\boldsymbol x}_{2} \\ \end{array} \right) &\sim {\mathcal N} \left( \left( \begin{array}{c} {\boldsymbol 0} \\ {\boldsymbol 0} \\ \end{array} \right), \left( \begin{array}{cc} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \\ \end{array} \right) \right) \end{align} $$

だとします。データを適当に2分割して表示しただけですね。また、

$$ \begin{align} \Lambda = \left( \begin{array}{cc} \Lambda_{11} & \Lambda_{12} \\ \Lambda_{21} & \Lambda_{22} \\ \end{array} \right) &= \left( \begin{array}{cc} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \\ \end{array} \right)^{-1} \end{align} $$

とすれば、

$$ \begin{align} p \left(\begin{array}{c} {\boldsymbol x}_{1} \\ {\boldsymbol x}_{2} \\ \end{array} \right) &\propto \exp \left\{ - \frac{1}{2} \left( \begin{array}{c} {\boldsymbol x_{1}} \\ {\boldsymbol x_{2}} \\ \end{array} \right)^{\mathrm T} \left( \begin{array}{cc} \Lambda_{11} & \Lambda_{12} \\ \Lambda_{21} & \Lambda_{22} \\ \end{array} \right) \left( \begin{array}{c} {\boldsymbol x_{1}} \\ {\boldsymbol x_{2}} \\ \end{array} \right) \right\} \end{align} $$

を確かめればよいです。2分割以上の場合でも同様に示せます。

は対称行列なので、と、2次形式の平方完成を使って

$$ \begin{align} {\mathcal L} &= - \frac{1}{2} \left( \begin{array}{c} {\boldsymbol x_{1}} \\ {\boldsymbol x_{2}} \\ \end{array} \right)^{\mathrm T} \left( \begin{array}{cc} \Lambda_{11} & \Lambda_{12} \\ \Lambda_{21} & \Lambda_{22} \\ \end{array} \right) \left( \begin{array}{c} {\boldsymbol x_{1}} \\ {\boldsymbol x_{2}} \\ \end{array} \right) \\ &= {\boldsymbol x_{1}}^{\mathrm T} \Lambda_{11} {\boldsymbol x_{1}} + {\boldsymbol x_{1}}^{\mathrm T} \Lambda_{12} {\boldsymbol x_{1}} + {\boldsymbol x_{2}}^{\mathrm T} \Lambda_{21} {\boldsymbol x_{1}} + {\boldsymbol x_{2}}^{\mathrm T} \Lambda_{22} {\boldsymbol x_{2}} \\ &= {\boldsymbol x_{1}}^{\mathrm T} \Lambda_{11} {\boldsymbol x_{1}} + 2 {\boldsymbol x_{1}}^{\mathrm T} \Lambda_{21} {\boldsymbol x_{2}} + {\boldsymbol x_{2}}^{\mathrm T} \Lambda_{22} {\boldsymbol x_{2}} \\ &= \left({\boldsymbol x}_2 + \Lambda_{22}^{-1}\Lambda_{21}{\boldsymbol x}_{1} \right)^{\mathrm T} \Lambda_{22} \left({\boldsymbol x}_2 + \Lambda_{22}^{-1}\Lambda_{21}{\boldsymbol x}_{1} \right) - {\boldsymbol x}_1^{\mathrm T} \Lambda_{21}^{\mathrm T}\Lambda_{22}^{-1}\Lambda_{21}{\boldsymbol x}_{1} + {\boldsymbol x}_{1}^{\mathrm T}\Lambda_{11}{\boldsymbol x}_{1} \end{align} $$

ここで第一項の二次形式はで積分されたとき、前回の記事の多変数のガウス積分を用いれば、に依存する定数項になります。以上より、

$$ \begin{align} p({\boldsymbol x}_{1}) &= \int p \left({\boldsymbol x}_{1}, {\boldsymbol x}_{2} \right) {\mathrm d}{\boldsymbol x}_{2} \\ &= \exp \left\{ - \frac{1}{2} \left( {\boldsymbol x}_1^{\mathrm T} \Lambda_{21}^{\mathrm T}\Lambda_{22}^{-1}\Lambda_{21}{\boldsymbol x}_{1} + {\boldsymbol x}_{1}^{\mathrm T}\Lambda_{11}{\boldsymbol x}_{1} \right) \right\} \exp \left( const. \right) \\ &\propto \exp \left\{ - \frac{1}{2} {\boldsymbol x}_1^{\mathrm T} \left( \Lambda_{11} - \Lambda_{21}^{\mathrm T}\Lambda_{22}^{-1}\Lambda_{21} \right) {\boldsymbol x}_{1} \right\} \end{align} $$

以上より、

$$ \begin{align} p({\boldsymbol x}_{1}) &= {\mathcal N} \left({\boldsymbol 0}, \left( \Lambda_{11} - \Lambda_{21}^{\mathrm T}\Lambda_{22}^{-1}\Lambda_{21} \right)^{-1} \right) \end{align} $$

である。式 (17)、(18)より、

$$ \begin{align} \left( \Lambda_{11} - \Lambda_{21}^{\mathrm T}\Lambda_{22}^{-1}\Lambda_{21} \right)^{-1} &= \Sigma_{11} \end{align} $$

となることがわかる。式 (20)と見比べると、ちょうどに関するものが消えていることがわかります。

条件付き多変量ガウス分布

$$ \begin{align} \begin{array}{l} p\left({\boldsymbol x}_{2} | {\boldsymbol x}_{1} \right) &\propto p\left({\boldsymbol x}_{1}, {\boldsymbol x}_{2} \right)\\ &\propto \exp \left(-\frac{1}{2} \left\{ \left(\begin{array}{l} {\boldsymbol x}_{1}-{\boldsymbol \mu}_{1} \\ {\boldsymbol x}_{2}-{\boldsymbol \mu}_{2} \end{array}\right)^{T} \left(\begin{array}{ll} {\boldsymbol \Lambda}_{11} & {\boldsymbol \Lambda}_{12} \\ {\boldsymbol \Lambda}_{21} & {\boldsymbol \Lambda}_{22} \end{array} \right) \left(\begin{array}{l} {\boldsymbol x}_{1}- {\boldsymbol \mu}_{1} \\ {\boldsymbol x}_{2}- {\boldsymbol \mu}_{2} \end{array}\right) \right\}\right) \\ &=\exp \left(-\frac{1}{2}\left\{\left(x_{1}-\mu_{1}\right)^{\mathrm{T}} \mathbf{\Lambda}_{11}\left(x_{1}-\boldsymbol{\mu}_{1}\right)+\left(x_{1}-\boldsymbol{\mu}_{1}\right)^{\mathrm{T}} \mathbf{\Lambda}_{12}\left(x_{2}-\boldsymbol{\mu}_{2}\right)+\left(x_{2}-\boldsymbol{\mu}_{2}\right)^{\mathrm{T}} \mathbf{\Lambda}_{21}\left(x_{1}-\boldsymbol{\mu}_{1}\right)+\left(x_{2}-\boldsymbol{\mu}_{2}\right)^{\mathrm{T}} \mathbf{\Lambda}_{22}\left(x_{2}-\boldsymbol{\mu}_{2}\right)\right\}\right)\\ &\propto \exp \left[-\frac{1}{2}\left\{\left(x_{2}-\mu_{2}\right)^{\mathrm{T}} \mathbf{\Lambda}_{22}\left(x_{2}-\mu_{2}\right)+2\left(x_{1}-\mu_{1}\right)^{\mathrm{T}} \mathbf{\Lambda}_{21}\left(x_{2}-\mu_{2}\right)\right\}\right]\\ &\propto \exp \left[-\frac{1}{2}\left\{x_{2}^{x} \mathbf{\Lambda}_{2 x} x_{2}-x_{2}^{T} \mathbf{\Lambda}_{2} \mu_{2}-\mu_{2}^{T} \mathbf{\Lambda}_{2} x_{2}+\mu_{2}^{T} \mathbf{\Lambda}_{22} \mu_{2}+2\left(x_{1}-\mu_{1}\right)^{T} \mathbf{\Lambda}_{21} x_{2}-2\left(x_{1}-\mu_{1}\right)^{T} \mathbf{\Lambda}_{21} \mu_{2}\right\}\right] \end{array} \end{align} $$

に関する項だけ残せば

$$ \begin{align} &\propto \exp \left(-\frac{1}{2}\left\{x_{2}^{\mathrm{T}} \mathbf{\Lambda}_{22} x_{2}-x_{2}^{\mathrm{T}} \mathbf{\Lambda}_{22} \boldsymbol{\mu}_{2}-\boldsymbol{\mu}_{2}^{\mathrm{T}} \boldsymbol{\Lambda}_{22} x_{2}+2\left(x_{1}-\boldsymbol{\mu}_{1}\right)^{\mathrm{T}} \boldsymbol{\Lambda}_{21} x_{2}\right\}\right)\\ &=\exp \left[-\frac{1}{2}\left\{x_{2}^{\mathrm{T}} \mathbf{\Lambda}_{22} x_{2}-2 x_{2}^{\mathrm{T}} \mathbf{\Lambda}_{22} \boldsymbol{\mu}_{2}+2 x_{2}^{\mathrm{T}} \mathbf{\Lambda}_{21}\left(x_{1}-\mu_{1}\right)\right\}\right]\\ &=\exp \left[-\frac{1}{2}\left\{x_{2}^{\mathrm{T}} \mathbf{\Lambda}_{22} x_{2}-2 x_{2}^{\mathrm{T}}\left(\mathbf{\Lambda}_{22} \boldsymbol{\mu}_{2}-\mathbf{\Lambda}_{21}\left(x_{1}-\mu_{1}\right)\right)\right\}\right] \end{align} $$

これを平方完成すれば、

$$ \begin{align} \propto \exp \left(-\frac{1}{2}\left\{\left(x_{2}-\Lambda_{22}^{-1}\left(\Lambda_{22} \mu_{2}-\Lambda_{21}\left(x_{1}-\mu_{1}\right)\right)\right)^{T} \Lambda_{22}\left(x_{2}-\cdots\right)\right\}\right) \end{align} $$

以上より、

$$ \begin{align} p\left(x_{2} | x_{1}\right) & \sim \mathcal{N}\left(\Lambda_{22}^{-1}\left(\Lambda_{22} \mu_{2}-\Lambda_{21}\left(x_{1}-\mu_{1}\right)\right), \Lambda_{22}^{-1}\right) \\ &=\mathcal{N}\left(\mu_{2}-\Lambda_{22}^{-1} \Lambda_{21}\left(x_{1}-\mu_{1}\right), \Lambda_{22}^{-1}\right) \end{align} $$

ここで、ブロック行列の逆行列を用いれば、

$$ \begin{align} p\left(x_{2} | x_{1}\right) & \sim \mathcal{N}\left(x_{2}-S\left(-S^{-1} \Sigma_{21} \Sigma_{11}^{-1}\right)\left(x_{1}-\mu\right), S\right) \\ &=\mathcal{N}\left(\mu_{2}+\Sigma_{21} \Sigma_{11}^{-1}\left(x_{1}-\mu_{1}\right), \Sigma_{22}-\Sigma_{21} \Sigma_{11}^{-1} \Sigma_{21}\right) \end{align} $$

となります。平均については、で標準化してから、で線形変換していると考えるとしっくりくる気がします（個人の感想）。

多変量ガウス分布を理解したい１

ガウス分布こと正規分布は、以下の式で定義される最も基本的かつ重要な確率分布です。

$$ \begin{align} p(x)= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left\{ -\frac{(x - \mu)^2 }{ 2 \sigma^2 } \right\} \end{align} $$

は確率変数、は平均、は分散です。

上で示した単変量の正規分布の扱いは簡単だけど、多変量になった途端扱いが大変だなぁと思ったので正規化定数の導出までを今回はまとめます。

多変量ガウス分布の定義

次元の確率変数ベクトルが平均、共分散行列の多変量ガウス分布に従うとき、その確率密度関数は

$$ \begin{align} p({\boldsymbol x} | {\boldsymbol \mu}, {\boldsymbol \Sigma})= \frac{1}{\left( \sqrt{2 \pi} \right)^n \sqrt{\det \Sigma}} \exp \left\{ -\frac{1}{2} \left({\boldsymbol x}- {\boldsymbol \mu} \right)^{\mathrm T} \Sigma^{-1} \left({\boldsymbol x}- {\boldsymbol \mu} \right) \right\} \end{align} $$

で定義される。

多変量ガウス分布を理解するための数学の準備

多変量ガウス分布をすっきり理解するために、いくつか数学の準備をしておきましょう。

単変数のガウス積分

を正の定数として、変数に関して $$ \begin{align} I = \int_{-\infty}^{\infty} e^{-ax^2} \mathrm{d}x \end{align} $$ で書かれる積分をガウス積分といいます。

これが $$ \begin{align} \int_{-\infty}^{\infty} e^{-ax^2} \mathrm{d}x = \sqrt{\frac{\pi}{a}} \end{align} $$

となることを確認します。

$$ \begin{align} I^2 &= \left( \int_{-\infty}^{\infty} e^{-ax^2} \mathrm{d}x \right)^2 \\ &=\left( \int_{-\infty}^{\infty} e^{-ax^2} \mathrm{d}x \right)\left( \int_{-\infty}^{\infty} e^{-ay^2} \mathrm{d}y \right) \\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-ax^2-ay^2} \mathrm{d}x \mathrm{d}y \end{align} $$

(6)から(7)へはフビニの定理を用いました。ここで、

$$ \begin{align} x &= r \cos \theta \\ y &= r \sin \theta \end{align} $$

と置くと、ヤコビアンは

$$ \begin{align} \left| \begin{array}{cc} \frac{\mathrm{d}x}{\mathrm{d}r} & \frac{\mathrm{d}x}{\mathrm{d}\theta} \\ \frac{\mathrm{d}y}{\mathrm{d}r} & \frac{\mathrm{d}y}{\mathrm{d}\theta} \\ \end{array} \right| = \left| \begin{array}{cc} \cos \theta & r \sin \theta \\ \sin \theta & r \cos \theta \\ \end{array} \right| = r \cos^2 \theta + r \sin^2 \theta = r \end{align} $$

となるので、

$$ \begin{align} (7) &= \int_{0}^{\infty} \int_{0}^{2\pi} e^{-ar^2} \cdot r \mathrm{d}\theta \mathrm{d}r \\ &= 2\pi \int_{0}^{\infty} re^{-ar^2} \mathrm{d}r \\ &= 2\pi \left[\frac{e^{-ar^2}}{-2a} \right]_{0}^{\infty} = 2\pi \left(0 - \frac{1}{-2a} \right) = \frac{\pi}{a} \\ I &= \sqrt{\frac{\pi}{a}} \end{align} $$

二次形式の平方完成

$$ \begin{align} -&\frac{1}{2} {\boldsymbol x}^{\mathrm T} A {\boldsymbol x} + {\boldsymbol b}^{\mathrm T}{\boldsymbol x} \\ &= -\frac{1}{2} \left(A^{1/2} {\boldsymbol x} - A^{1/2}{\boldsymbol b} \right)^{\mathrm T} \left(A^{1/2} {\boldsymbol x} - A^{1/2}{\boldsymbol b} \right) + \frac{1}{2} \left( {\boldsymbol b}^{\mathrm T} A^{-1} {\boldsymbol b} \right) \end{align} $$

多変数のガウス積分

数学的準備ができたところで多変数のガウス積分をやっていきます。 をの正則な対称行列、を次元の定数のベクトルとして、次元の変数ベクトルに関して、

$$ \begin{align} I &= \int \exp\left\{ -\frac{1}{2} {\boldsymbol x}^{\mathrm T} A {\boldsymbol x} + {\boldsymbol b}^{\mathrm T} {\boldsymbol x} \right\} \mathrm{d}{\boldsymbol x} \\ \end{align} $$

の形の積分を、多変数のガウス積分と呼びます。これが

$$ \begin{align} I &= \frac{1}{\left( \sqrt{2 \pi} \right)^n \sqrt{\det A}} \exp \left\{ -\frac{1}{2} {\boldsymbol b}^{\mathrm T} A^{-1} {\boldsymbol b} \right\} \end{align} $$

となることを確認します。

$$ \begin{align} I &= \int \exp \left\{-\frac{1}{2} {\boldsymbol x}^{\mathrm T} A {\boldsymbol x} + {\boldsymbol b}^{\mathrm T} {\boldsymbol x} \right\} {\mathrm d}{\boldsymbol x} \end{align} $$

2次形式の平方完成より

$$ \begin{align} I &= \int \exp \left\{-\frac{1}{2} \left(A^{1/2}{\boldsymbol x} - A^{-1/2}{\boldsymbol b} \right)^{\mathrm T} \left(A^{1/2}{\boldsymbol x} - A^{-1/2}{\boldsymbol b} \right) + \frac{1}{2}{\boldsymbol b}^{\mathrm T}A^{-1}{\boldsymbol b} \right\} {\mathrm d}{\boldsymbol x} \\ &= \exp \left\{\frac{1}{2}{\boldsymbol b}^{\mathrm T}A^{-1}{\boldsymbol b} \right\} \int \exp \left\{-\frac{1}{2} \left(A^{1/2}{\boldsymbol x} - A^{-1/2}{\boldsymbol b} \right)^{\mathrm T} \left(A^{1/2}{\boldsymbol x} - A^{-1/2}{\boldsymbol b} \right) \right\} {\mathrm d}{\boldsymbol x} \end{align} $$

ここで、

$$ \begin{align} {\boldsymbol y} &= A^{1/2}{\boldsymbol x} - A^{-1/2}{\boldsymbol b} \end{align} $$

とおくと、

$$ \begin{align} {\boldsymbol x} &= A^{1/2}{\boldsymbol y} + A^{-1}{\boldsymbol b} \end{align} $$

なので、この時のヤコビアンは

$$ \begin{align} \frac{\partial {\boldsymbol x}}{\partial {\boldsymbol y}} &= \det A^{-1/2} = \left(\det A \right)^{-1/2} \end{align} $$

より、

$$ \begin{align} I & = \exp \left\{\frac{1}{2}{\boldsymbol b}^{\mathrm T}A^{-1}{\boldsymbol b} \right\} \int \exp \left\{-\frac{1}{2} {\boldsymbol y}^{\mathrm T} {\boldsymbol y} \cdot \left(\det A \right)^{-1/2} \right\} {\mathrm d}{\boldsymbol y} \\ &= \left(\det A \right)^{-1/2} \cdot \exp \left\{\frac{1}{2}{\boldsymbol b}^{\mathrm T}A^{-1}{\boldsymbol b} \right\} \int \exp \left\{-\frac{1}{2} {\boldsymbol y}^{\mathrm T} {\boldsymbol y} \right\} {\mathrm d}{\boldsymbol y} \\ &= \frac{1}{\left( \sqrt{2 \pi} \right)^n \sqrt{\det A}} \exp \left\{ -\frac{1}{2} {\boldsymbol b}^{\mathrm T} A^{-1} {\boldsymbol b} \right\} \end{align} $$

最後の式変形には、1変数のガウス積分を用いました。

多変量ガウス分布では、より、

$$ \begin{align} I & = \frac{1}{\left( \sqrt{2 \pi} \right)^n \sqrt{\det \Sigma}} \end{align} $$

となり、多変量ガウス分布の正規化定数が求まりました。

あまりに記事が長くなりすぎるのでいったんここで区切ります。 次から実際に分布を弄って、周辺化、条件付き分布を導出しようと思います。

参考

www.amazon.co.jp

Wikipedia: 多変量正規分布

ガウス積分の公式の２通りの証明 | 高校数学の美しい物語

多変数のガウス積分 | 高校数学の美しい物語

ヤコビアンに関してはヨビノリたくみさんのyoutube動画がお勧め。

EMアルゴリズムを理解したい【最尤法から】

よく忘れるので一旦ちゃんとEMアルゴリズム*1についてまとめておきます。
最尤法について振り返ってからEMアルゴリズムの説明に入っていきます。
あまり計算に自信はありませんので、間違いを見つけたら教えていただけると幸いです。

最尤法

観測事象の確率密度関数から尤度を構成し、その確率密度関数のパラメータの推測を尤度の最大化によって行う方法です。
さっそく説明をしていきたいところですがその前に、尤度について振り返りましょう。

尤度についてはこちらをどうぞ

hotoke-x.hatenablog.com

それでは最尤法についてですが、実は一度正規分布を例に一度記事を書いています。

hotoke-x.hatenablog.com

ここでは特定の分布を仮定せず、一般的な最尤法について書いていきます。

---

を母数、 を確率変数の実現値、を母集団が従う確率密度関数とします。
この時、個の観測値の同時確率関数は として、

$$ \begin{align} p({\boldsymbol x} | {\boldsymbol \theta}) = \prod_{i=1}^{n} p(x_i|{\boldsymbol \theta}) \end{align} $$

で与えられます。これを、同時確率分布ではなく、母数の関数として考えると

$$ \begin{align} L_n ({\boldsymbol \theta} |{\boldsymbol x}) = \prod_{i=1}^{n} p(x_i | {\boldsymbol \theta}) \end{align} $$

と書け、これを尤度関数と呼びます。

最尤法は、式(2)の尤度を最大化するを推定する方法です。この時のを最尤推定値と呼び、と書くのが慣習的となっています。また、の関数、を最尤推定量と呼びます。忘れがちなのでちゃんと区別しましょう。

現実問題では、最尤推定量を陽に得ることが難しいことが多いため、対数尤度

$$ \begin{align} l \left({\boldsymbol \theta}, {\boldsymbol x} \right) = \log L \left({\boldsymbol \theta}, {\boldsymbol x} \right) = \sum_{i=1}^{n} \log p \left( x_i | {\boldsymbol \theta} \right) \end{align} $$

を数値的に最適化してを求めることが多く、その最適化法の一つがEMアルゴリズムがあります。最適化手法はもちろん他にも色々あります。その他についてはこちら。

hotoke-x.hatenablog.com

EMアルゴリズム

実現値が得られない観測に対する確率変数をとすると、観測データの確率密度関数は

$$ \begin{align} p ({\boldsymbol x}|{\boldsymbol \theta}) = \int p ({\boldsymbol x},{\boldsymbol y}|{\boldsymbol \theta}) \mathrm{d}{\boldsymbol y} \end{align} $$

と書けます。「観測できない点は考えられる点を全部調べる」という気持ちが積分に表れていますね。

なぜ観測できていない点をわざわざ式に埋め込むかというと、不完全データの対数尤度関数より完全データの対数尤度関数の方が扱いやすいからです。

式の上では完全データとはいえ、実際には観測できてない確率変数が混ざっているのでちょっと扱いを工夫します。それがEMアルゴリズムです。

「なぜそれで良いのか」という細かいことは置いといて、まずはEMアルゴリズムの枠組みをおさらいしましょう。

1. パラメータの初期値を適当な値に設定（スケールくらいは当たりをつけて置くと良い）。
2. 観測値とパラメータの推定値から、観測できないについて条件付き期待値を算出。これを疑似的なの観測値とする。疑似観測から疑似完全データを構成（Eステップ）
3. 疑似完全データに基づいて対数尤度を最大化するパラメータを推定（Mステップ）
4. を新たなパラメータ初期値としてEステップ、Mステップを繰り返す。

さて、いろいろ気になる点があります。

1. 疑似的な観測値を構成するのに条件付き期待値を使ってよいのか
2. 疑似的な完全データの尤度を最大化することが、観測データの尤度を最大にしたことになるのか
3. 疑似的な完全データによる推測と、完全な観測による推測とで乖離はないのか（バイアスがあるのではないかということ）

1つ目は今回の記事では解説しきれそうにないので後日追記ということにして、2つ目の疑問から確認していきます。

疑似的な完全データの尤度を最大化することが、観測データの尤度を最大にしたことになるのか

観測が与えられたときの尤度は

$$ \begin{align} l \left({\boldsymbol \theta}, {\boldsymbol x} \right) &= \log p \left({\boldsymbol x}|{\boldsymbol \theta} \right) \\ &= \log p \left({\boldsymbol x, \boldsymbol y}|{\boldsymbol \theta} \right) - \log p \left({\boldsymbol y}|{\boldsymbol x, \boldsymbol \theta} \right) \end{align} $$

となります。しかし困りました。観測できないが式に含まれています。仕方がないので適当なとして完全データの条件付き期待値を計算しましょう（つまり観測していない部分を周辺化するということ）。適当にを使っても問題ない理由もこの後わかります。

$$ \begin{align} l \left({\boldsymbol \theta}, {\boldsymbol x} \right) &= \int p \left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta^{(k)}} \right) \left\{ \log p \left({\boldsymbol x, \boldsymbol y}|{\boldsymbol \theta} \right) - \log p \left({\boldsymbol y}|{\boldsymbol x, \boldsymbol \theta} \right) \right\} \mathrm{d}{\boldsymbol y} \\ &= \mathrm{E}_{\boldsymbol \theta^{(k)} } \left[ \log p \left({\boldsymbol x, \boldsymbol y}|{\boldsymbol \theta} \right) \right] - \mathrm{E}_{\boldsymbol \theta^{(k)} } \left[ \log p \left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta} \right) \right] \end{align} $$

これはとの関数になっているので、新たに関数を導入して

\begin{align} l \left({\boldsymbol \theta}, {\boldsymbol x} \right) &= Q \left({\boldsymbol \theta}, {\boldsymbol \theta^{(k)}} \right) - H \left( {\boldsymbol \theta}, {\boldsymbol \theta^{(k)}} \right) \end{align}

と書くことができます。また、は、によって定まる関数を引数にとる汎関数であるという見方もできます。

さて、ここでの式を改めて眺めてみます。は

$$ \begin{align} Q \left({\boldsymbol \theta}, {\boldsymbol \theta^{(k)}} \right) &= \int p \left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta^{(k)}} \right) \log p \left({\boldsymbol x, \boldsymbol y}|{\boldsymbol \theta} \right) \mathrm{d}{\boldsymbol y} \end{align} $$

です。これは、事前分布に基づいて完全データの対数尤度の期待値を計算しているだけです。完全データの条件付き期待値を用いた対数尤度の最適化ができれば、実際に完全データが得られた場合の対数尤度の値へ収束することがわかります。初期値に適当なを使って問題ないのは、「完全データの条件付き期待値の対数尤度が最適化可能である」という前提からなのでした（実際は局所解へ陥る可能性があることに注意）。

いやいや関数は放っておいていいのかよと思われるかもしれませんが、この後説明します。

疑似的な完全データによる推測と、完全な観測による推測とで乖離はないのか（バイアスがあるのではないかということ）

当然バイアスはあります。実は、疑似データを用いたことによる現実との乖離度合いを表しているのがになります。式を眺めてみましょう。

$$ \begin{align} H \left( {\boldsymbol \theta}, {\boldsymbol \theta^{(k)}} \right) &= - \int p \left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta^{(k)}} \right) \log p \left({\boldsymbol y}|{\boldsymbol x, \boldsymbol \theta} \right) \mathrm{d}{\boldsymbol y} \\ \end{align} $$

ある程度勉強してきた方ならピンとくるのではないでしょうか。クロスエントロピーの式になっていますね。つまり、理想的な最適解が求まったとすればこの値は自己エントロピーに一致し最小となります。

対数尤度を最大にするのにこの項は最小になってよいのかと思われる方がいるかもしれませんが、その分、完全データの条件付き期待値の対数尤度が大きくなります。そもそも「疑似データを用いたことによる現実との乖離度合いを表している」だけなのです。

疑似的な観測値を構成するのに条件付き期待値を使ってよいのか（pending）

まず、これは個人の見解であり正確な解釈ではない可能性があることをご理解下さい。
実はEMアルゴリズムは暗に指数分布族と呼ばれる確率密度関数を使うことを仮定しているのだと思います。
指数分布で考えると、対数尤度をとるときにexponentialをきれいに打ち消せるのです。
すると全体が見えてきて、十分統計量をその条件付き期待値で置き換えて考える問題に帰着します。
今回は実用重視の記事にしたいので、数理統計の細かいところは次の機会もしくは追記という形にしたいと思います。

参考

小西 貞則、越智 義道、大森 裕浩（2008）「計算統計学の方法」朝倉書店

qiita.com

*1:Dempster, A. P., Laird, N. M. and Rubin, D. B. (1977), Maximum Likelihood from Incomplete Data Via the EM Algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39: 1-22. doi:10.1111/j.2517-6161.1977.tb01600.x

neovimでAnacondaの環境ごとにpython補完させる

環境：Ubuntu18.04 LTS

Anaconda (Miniconda) のインストール

neovim用環境の作成
conda create -n neovim python=3
pip install neovim

neovimのインストール

sudo add-apt-repository ppa:neovim-ppa/stable
sudo apt-get update
sudo apt-get install neovim

Deinの導入

mkdir -p $XDG_CONFIG_HOME/nvim/dein/toml
touch $XDG_CONFIG_HOME/nvim/dein/toml/init.viml
touch $XDG_CONFIG_HOME/nvim/dein/toml/dein.toml
touch $XDG_CONFIG_HOME/nvim/dein/toml/dein_lazy.toml

init.vimを編集。 典型的なinit.vimなのでコピペで動くと思います。

let g:python3_host_prog = $HOME.'/miniconda3/envs/neovim3/bin/python'

" initialize
if &compatible
  set nocompatible
endif
set runtimepath+=~/.cache/dein/repos/github.com/Shougo/dein.vim
if dein#load_state('~/.cache/dein')
  call dein#begin('~/.cache/dein')
  call dein#load_toml('$XDG_CONFIG_HOME/nvim/dein.toml', {'lazy': 0})
  call dein#load_toml('$XDG_CONFIG_HOME/nvim/dein_lazy.toml', {'lazy': 1})
  call dein#end()
  call dein#save_state()
endif
filetype plugin indent on
syntax enable

if dein#check_install()
  call dein#install()
endif

Deopleteの導入 dein.tomlへ以下を記入

[[plugins]] repo = 'Shougo/deoplete.nvim'
hook_add = 'let g:deoplete#enable_at_startup = 1'

Deoplete-jediの導入 activateしているAnacondaのpython_pathを自動取得してdeopleteへ反映

dein_lazy.tomlへ以下を記入（遅延読み込みさせる）

[[plugins]]
repo = 'deoplete-plugins/deoplete-jedi'
on_ft = 'python'
hook_add = """
let g:deoplete#sources#jedi#python_path = system("which python")[:-2]
"""

system("which python") でpython pathを取得し、[:-2]で改行文字を落としています。

source activateまたはconda activateで環境に入ってneovimを起動すれば、自動的にその環境での補完をしてくれます。

Julia tips #10: マクロを使って動的に構造体を定義する

Juliaでゼロから学ぶDeep Learningという記事を書きました。

hotoke-x.hatenablog.com

実はこちらの実装、ほぼ本家のPythonそっくりになるように書かれています。

というのも私がDeep Learning初心者なので、自己流で進めすぎるとどこが間違っているのかわからなくなってしまうためです。

「Juliaらしさ」も未だつかめておらず、見比べたらわかるんですがほぼPythonです。

実装に四苦八苦する中でマクロについてちょっと調べる機会があったので記事にしておこうと思った次第です。

ただし、Documentを全て読んだわけではないですし、今でも色々わかっていないということはご了承下さい。

（忘れないうちに記事にして消化しておきたかった）

構造体を動的に定義したい

とりあえず何も考えずに以下のコードを実行（ver. 1.1.0です）。

macro make_struct(struct_name, schema...)
    fields=[:($(entry.args[1])::$(entry.args[2])) for entry in schema]
    field_names=[:($(entry.args[1])) for entry in schema[1:2]]
    obj = :OBJ
    __init = [:($obj.$name=$name) for name in field_names] 
    esc(quote mutable struct $struct_name
        $(fields...)
            function (::Type{$(struct_name)})($(field_names...))
                $obj = new()
                $(__init...)
                $obj
            end
        end
    end)
end
println(@macroexpand @make_struct STRUCT_NAME (x,Integer) (y, Integer) (z, Integer))
@make_struct STRUCT_NAME (x, Integer) (y, Integer) (z, Integer)
STRUCT_NAME(1,2)

STRUCT_NAMEという構造体が定義されたはずです。

このマクロは前半でコードの部品生成、後半部分で部品の埋め込みを行っています。

部品生成部分にコメントをつけました。

macro make_struct(struct_name, schema...)
    # 部品（Expr）生成###################################
    fields=[:($(entry.args[1])::$(entry.args[2])) for entry in schema]   # [x::Integer, y::Integer, z::Integer]を生成
    field_names=[:($(entry.args[1])) for entry in schema[1:2]]   # [x::Integer, y::Integer]を生成（内部コンストラクタの引数に使う）
    obj = :OBJ   # このOBJという名前を内部コンストラクタで使う
    __init = [:($obj.$name=$name) for name in field_names]   # [OBJ.x=x, OBJ.y=y]を生成
    # 部品（Expr）生成###################################

    esc(quote mutable struct $struct_name
        $(fields...)
            function (::Type{$(struct_name)})($(field_names...))
                $obj = new()
                $(__init...)
                $obj
            end
        end
    end)
end
println(@macroexpand @make_struct STRUCT_NAME (x,Integer) (y, Integer) (z, Integer))
@make_struct STRUCT_NAME (x, Integer) (y, Integer) (z, Integer)
STRUCT_NAME(1,2)

埋め込み部分が何をしているかは

println(@macroexpand @make_struct STRUCT_NAME (x,Integer) (y, Integer) (z, Integer))

でわかります。

おわりに

ゼロから始めるDeep Learningを実装する中でJuliaの辞書が遅いということに気がついたので、「動的にstructを定義できたら速いのでは？」と思ったので今回のようなことを思いついたのでした。

ちなみに、「Juliaの辞書が遅い」というのは正確には間違っています。

Juliaの辞書は型指定無しもしくは抽象型で作った時にPythonより遅くなるようです。具象型であればPythonより高速です。

これは後ほど記事にする予定です。

Juliaでゼロから作るDeep Learning

もともとPython使いの自分には、ゼロから作るDeep LearningをPythonでやっても面白くない！

ということでJuliaでやってみた。

速度が犠牲にならないように多少は考慮したつもりだが、本家Python版の方が速いかも。

なんにも考えなくても速度がでるNumpy凄い...

Juliaは数学的に自然な形でコードが書けて嬉しいが、行列演算をちょっと工夫して書かないとNumpyより遅くなる。

また、実装していて気がついたのだがどうやらdictionaryへのアクセスがJuliaは遅いらしい？

discourse.julialang.org

dictionaryは便利だけど、できるだけstructを使って回すのが良さそうです。

最後までちゃんと読まなかったようで、Juliaの方が早いとの事。失礼しました。

#Julia言語 1秒(Python2,3)と0.18秒(Julia v0.7)を比較するとJuliaの側が明らかに圧倒的に速い。

知識が足りないせいで根本的におかしなことを言っている人がフルボッコにされているhttps://t.co/CA7lrk1Ur1

を引用しながら、「dictionaryへのアクセスがJuliaは遅いらしい？」と書くのは変。

— 黒木玄 Gen Kuroki (@genkuroki) 2019年6月17日

twitter.com

ただ、何にも工夫しないと結局遅いようなので後ほど試してまとめます。

以下にコードがあります。

商用・非商用問わず自己責任でご利用下さい。

github.com