概要
PRML第7章の関連ベクトルマシン（RVM; relevance vector machine）による回帰をpythonで実装．
コードと実験結果をまとめたJupyter notebook．
RVMのモデル
- ベイズ線形回帰の事前分布を疎な解が得られるように修正したもの
与えられた入力ベクトル$\mathbf{x}$に対する実数値の目標変数$t$の条件付き確率分布
$$
p(t|\mathbf{x}, \mathbf{w}, \beta) = \mathcal{N}(t|y(\mathbf{x}), \beta^{-1})
$$
ここで$\beta=\sigma^{-2}$はノイズの精度パラメータ（ノイズの分散の逆数），平均値は次の線形モデルで定義される．
$$
y(\mathbf{x})=\sum_{i=1}^{M}w_{i}\phi_{i}(\mathbf{x})=\mathbf{w}^{T}\mathbf{\phi}(\mathbf{x})
$$
基底関数として個々の訓練データを一方の引数としたカーネル関数を用いる．
$$
y(\mathbf{x})=\sum_{n=1}^{N}w_{n}k(\mathbf{x}, \mathbf{x}_{n})+b
$$
パラメータの数は全部で$M=N+1$である．尤度関数は次式で与えられる．
$$
p(\mathbf{t}| \mathbf{X}, \mathbf{w}, \beta) = \prod_{n=1}^{N} p(t_{n}|\mathbf{x}_{n}, \mathbf{w}, \beta)
$$
パラメータベクトル$\mathbf{w}$の事前分布として，平均0のガウス事前分布を用いる．RVMにおいては個々の重みパラメータ$w{i}$ごとに異なった超パラメータ$\alpha{i}$を用いる．__つまり，重みに対する事前分布は次のようになる．
$$
p(\mathbf{w}|\mathbf{\alpha}) = \prod_{i=1}^{M} \mathcal{N} (w_{i}|0, \alpha_{i}^{-1})
$$
ここで，$\alpha_{i}$は対応する重みパラメータ$w_{i}$の精度を表し，$\mathbf{\alpha} = (\alpha_{1}, \dots ,\alpha_{M})^{T}$である．これらの超パラメータについてエビデンスを最大化すると，大部分の超パラメータは無限大になり，対応する重みパラメータの事後分布を零一点に集中する．すると，これらのパラメータに対応する基底関数（対応するデータ点との距離を表すカーネル関数）は予測において何の役割も果たさないため，取り除くことができ，疎なモデルが得られる．
重みベクトルに対する事後確率は再びガウス分布となり，次の形で表される．
$$
p(\mathbf{w}|\mathbf{t}, \mathbf{X}, \mathbf{\alpha}, \beta) = \mathcal{N}(\mathbf{w}|\mathbf{m}, \mathbf{\Sigma})
$$
ここで平均，および共分散は次の式で与えられる．
$$
\mathbf{m} = \beta \mathbf{\Sigma} \mathbf{\Phi}^{T} \mathbf{t}
$$
$$
\mathbf{\Sigma} = \left( \mathbf{A} + \beta \mathbf{\Phi}^{T} \mathbf{\Phi} \right)^{-1}
$$
ただし，$\mathbf{\Phi}$は，$i=1, \dots ,N$について要素$\Phi_{ni}=\phi_{i}(\mathbf{x_{n}})$を，そして$n=1, \dots ,N$について要素$\Phi_{nM}=1$を持つ$N \times M$の計画行列であり，$\mathbf{A}=\rm{diag}(\alpha_{i})$である．
$\mathbf{\alpha}$と$\beta$の値は__エビデンス近似（evidence approximation）__としても知られている第二種の最尤推定によって求められる．第二種の最尤推定を行うためには，まず重みパラメータについて積分を行う．
$$
p(\mathbf{t}|\mathbf{X}, \mathbf{\alpha}, \beta) = \int p(\mathbf{t}|\mathbf{X}, \mathbf{w}, \beta)p(\mathbf{w}|\mathbf{\alpha}) d\mathbf{w}
$$
この式は2つのガウス分布のたたみ込み積分となっているため，解析的に積分を実行することができ，以下のように対数尤度が求まる．
$$
\ln p(\mathbf{t}|\mathbf{X}, \mathbf{\alpha}, \beta)
= \ln \mathcal{N} (\mathbf{t}|\mathbf{0}, \mathbf{C})
= -\frac{1}{2}{ N \ln (2 \pi) + \ln |\mathbf{C}| + \mathbf{t}^{T}\mathbf{C}^{-1}\mathbf{t}}
$$
ここで$\mathbf{t} = (t_{1}, \dots ,t_{N})^{T}$である．また$N \times N$行列$\mathbf{C}$を次のように定義した．
$$
\mathbf{C} = \beta^{-1}\mathbf{I} + \mathbf{\Phi} \mathbf{A}^{-1} \mathbf{\Phi}^{T}
$$
得られた対数尤度の微分を0とおくことで，次のように超パラメータの更新式を得る．
$$
\begin{split}
\alpha_{i}^{new} &= \frac{\gamma_{i}}{m_{i}^{2}} \
(\beta^{new})^{-1} &= \frac{|\mathbf{t} - \mathbf{\Phi m}|^{2}}{N - \Sigma_{i}\gamma_{i}}
\end{split}
$$
ここで，$\gamma_{i}$は対応する重みパラメータ$w_{i}$が，データにどれだけよく特定されたかを表す量で，次の式で定義される．
$$
\gamma_{i} = 1 - \alpha_{i} \Sigma_{ii}
$$
超パラメータの学習は，以上の結果を用いて，次のように進める．まず，適当に選んだ$\mathbf{\alpha}$と$\beta$の初期値から事後確率の平均$\mathbf{m}$と共分散$\mathbf{\Sigma}$を推定する．次に得られた値から超パラメータを推定する．この過程を適当な収束条件が満たされるまで交互に繰り返す．
実験
- PRMLの付録Aの人工データ（正弦関数）にRVM回帰を最適化
- ガウスカーネルの標準偏差$\sigma=0.2$．
訓練データ
結果
デザイン行列と関連ベクトル
出力$\mathbf{t}$に関連のないベクトルは$\alpha$の値が$\infty$となる．つまり$\alpha^{-1}$の値の大きい列が関連ベクトル．
RVM回帰結果
考察
- 逐次的疎ベイジアン学習アルゴリズムも実装しようかと思ったけど力つきた．
- SVMよりスパースな解が得られやすいので，SVMでスパースな解が得られず予測が遅くなる問題に応用したい．
- アルゴリズムが割に簡単なので，スパースモデリングにも使えそう．