変分近似アルゴリズム

変分近似アルゴリズムは，ベイズモデルの事後分布を近似的に求める手段の一つです．今回の記事では，変分近似アルゴリズムの導出についてまとめ，また実際に線形回帰モデルに適用してみたいと思います．

${\bf Z} = (z_1,z_2,\ldots,z_n)$ をハイパーパラメーターにもつ確率モデル $x | {\bf Z}$ を考える．[tex :x]に従うサンプル $D = \{x_1, \ldots, x_d \}$ に対する同時事後分布 ${\mathbb Z} | D$ を求めたい．同時事後分布の確率密度関数 $p({\bf Z} | D)$ は「Kullback-Leibler情報量

$\displaystyle{{\sf KL}[ q({\bf Z}) || p({\bf Z}, D)] = \int d {\bf Z} ~ q({\bf Z}) \ln \Bigl( \frac{q({\bf Z})}{p({\bf Z}, D)} \Bigr)}$

を最小化する分布 $q({\bf Z})$ 」と特徴づけられます．そこで，Kullback-Leibler情報量をなるべく小さくする $q{}^\ast({\bf Z})$ を見つける方法の一つに変分近似アルゴリズムがあります．以下その導出を行います．

まず， $q{}^\ast({\bf Z})$ として次の形で書けるものを考えます．

$\displaystyle{q{}^\ast ({\bf Z}) = \prod_{i = 1}^n q_i(z_i)\quad \cdots {\rm eq1}.}$

このテクニックは，平均場近似と呼ばれるものです．このような形で書ける $q{}^\ast ({\bf Z})$ のうちで $\displaystyle{{\sf KL}[ q({\bf Z}) || p({\bf Z}, D)]}$ を最小化する $q_i{}^\ast (z_i)$ を見つけるため変分法を利用します．

汎関数

$\displaystyle{S[ q_1, \ldots, q_n ]:= {\sf KL}[ q_1(z_1)\cdots q_n(z_n) || p({\bf Z}, D)] - \sum_{i = 1}^n \lambda _i \Bigl( \int d z_i ~ q(z_i) -1 \Bigr)}$

を考えます( $\lambda_i, \ldots, \lambda_n$ はLagrangeの未定係数)． $q_i$ に関する変分をとるため，次のように書き直す;

$\displaystyle{S[ q_1, \ldots, q_n] }$

$\displaystyle{ = \int dz_i {\mathcal L} (q_1, \ldots, q_n)} + (q_i\text{に依存しない部分})$ ，

ただし

$\displaystyle{{\mathcal L} (q_1, \ldots, q_n) = \int dz_{j \neq i}~ \prod_{j = 1}^n q_j(z_j)\ln \Bigl( \frac{ \prod_{j = 1}^n q_j(z_j)}{p(D, {\bf Z})}\Bigr) - \lambda_i q_i}.$

Euler-Lagrange方程式より，

$0 = \displaystyle{\frac{\partial {\mathcal L}}{\partial q_i}}$

$0 = \displaystyle{\frac{\partial}{\partial q_i} \int dz_{j \neq i}~ \prod_{j = 1}^n q_j(z_j)\ln \Bigl( \frac{ \prod_{j = 1}^n q_j(z_j)}{p(D, {\bf Z})} \Bigr) - \lambda_i}$

$0 = \displaystyle{\int dz_{j \neq i}~ \Bigl\{ \prod_{j \neq i} q_j(z_j)\ln \Bigl( \frac{ \prod_{j = 1}^n q_j(z_j)}{p(D, {\bf Z})}\Bigr) + \prod_{j \neq i}^n q_j(z_j) \Bigr\}-\lambda_i}$

整理して，

$\displaystyle{\ln ~q_i (z_i) = \int dz_{j \neq i}~ \prod_{j \neq i}^n q_j(z_j)\ln ~ p(D, {\bf Z}) + {\sf const}},$

となる( ${\sf const}$ は $z_i$ に無関係な定数)．

よって次のことがわかります．

( 平均場近似の停留解)

平均場近似eq1の下で汎関数 $S[ q_1, \ldots, q_n]$ の停留値問題の解 $q_1^{\ast}(z_1), \ldots, q_n^{\ast}(z_n)$ は次をみたす

$\displaystyle{\ln ~q_i ^{\ast}(z_i) = \int dz_{j \neq i}~ \prod_{j \neq i}^n q_j^{\ast}(z_j)\ln ~ p(D, {\bf Z}) + {\sf const}_i} \quad \cdots {\rm eq 2}$

,ここで ${\sf const}_i$ は $z_i$ に無関係な定数．

方程式eq2が解析的に解けない場合，反復アルゴリズムにより，近似的に $q_1^{\ast}(z_1), \ldots, q_n^{\ast}(z_n)$ を求める．例えば， $n = 2$ のときでは次のようになる

(変分近似アルゴリズム)

Step1.分布 $q_1(z_1)$ を適当に設定する．

以下Step2,Step3を十分に繰り返して( ${\sf KL}[p(D, z_1, z_2) || q_1(z_1)q_2(z_2)]$ の推移などで判断)，適当なところで打ち切る．

Step2.以下の式により $q_2(z_2)$ を更新する．

$\displaystyle{\ln q_2(z_2) = \int dz_1 q_1(z_1)\ln p(D, z_1, z_2) + {\sf const}}$

(ここで ${\sf const}$ は正規化定数とする．)

Step3.以下の式により $q_1(z_1)$ を更新する．

$\displaystyle{\ln q_1(z_1) = \int dz_2 q_2(z_2)\ln p(D, z_1, z_2) + {\sf const}}$

(ここで ${\sf const}$ は正規化定数とする．)

2.線形回帰モデルの場合

このセクションでは，線形回帰モデルを例に変分近似アルゴリズムを導出し，実際にその威力を見てみましょう．考えたいモデルは次のモデルです

$\displaystyle{y | x, \alpha , \beta , \tau \sim {\rm Normal} (\alpha + \beta x^{\sf T}, \tau ^{-1}) ,\quad \alpha \in {\mathbb R}, \beta \in {\mathbb R} {}^l , \tau \in (0, \infty ) \cdots {\rm eq7}}$

ここで， $y, x$ をそれぞれ応答変数，および説明変数，また ${\rm Normal}(\mu, \sigma ^{2})$ は平均 $\mu$ ，分散 $\sigma ^{2}$ の正規分布とします．さらに，パラメータ $\alpha, \beta, \tau$ について次の事前分布を仮定します

$\alpha \sim {\rm Normal}(\mu_0, \theta_0 ^{-1}),$

$\beta \sim {\rm Normal}(\lambda_0, \Phi_0 {}^{-1}),$

$\tau \sim {\rm Gamma}(a_0, b_0).$

ただしここで， ${\rm Gamma}(a, b)$ はガンマ分布，つまり次で確率密度が与えられる $(0, \infty)$ 上の分布です

$\displaystyle{ f_{\rm Gamma}(\tau | a, b) \propto \tau ^{a-1} e ^{- b \tau}.}$

$D = \{(y_1,x_1), \ldots, (y_d, x_d)\}$ をモデルeq7に従う独立なサンプルとして，変分近似アルゴリズムを導出します．

アルゴリズムの導出

記号:

${\mathbb X} := [ x_1, \ldots, x_d]$

とおく．また，確率変数 $Z$ に対して

${\mathbb E}_Z[ -]$

は，期待値を表す．

$\underline{\alpha {\text の更新}\ }$

更新後の $\alpha$ の分布を， $q {}^\ast(\alpha | D)$ とすれば，

$\ln q {}^\ast(\alpha | D)$

$= \displaystyle{ {\mathbb E}_{\tau, \beta} \Bigl[ \ln p(D, \alpha , \beta, \tau) \Bigr]}$

$= \displaystyle{ {\mathbb E}_{\tau, \beta} \Bigl[ \ln p(D | \alpha , \beta, \tau) + \ln p(\alpha) \Bigr]}$

$= \displaystyle{ {\mathbb E}_{\tau, \beta} \Bigl[ -\frac{1}{2} \tau \sum_{i = 1}^d (z_i -\alpha - \beta x_i^{\sf T})^2 - \frac{1}{2} \theta_0 (\alpha - \mu_0)^2 \Bigr] }$

となる．ただし $\alpha$ と無関係な定数を無視して計算している．ここで， $\tau | D$ と $\beta | D$ の分布が次で与えられるとしよう．

$\tau | D \sim {\rm Gamma} (a, b),$

$\beta | D \sim {\rm Normal}(\hat{\beta}, \Phi^{-1})$ .

すると $\alpha | D$ の更新後の分布は次のように表すことができる

$\alpha | D \sim {\rm Normal}(\mu,\theta ^{-1}).$

$\displaystyle{ \theta = d \hat{\tau} + \theta_0}$

$\displaystyle{\mu = \frac{\theta_0 \mu_0 + \hat{\tau} \sum_{i = 1}^d (z_i - \hat{\beta} x_i^{\sf T})}{d \hat{\tau} + \theta_0}}$

,ただし $\hat{\tau} = {\mathbb E}_{\tau | D}[ \tau] = \frac{a}{b}$ .

$\underline{\beta {\text の更新}\ }$

更新後の $\beta$ の分布を $q {}^\ast (\tau | D)$ とすれば，

$\ln q{}^\ast(\beta | D)$

$= \displaystyle{ {\mathbb E}_{\alpha, \tau} \Bigl[ \ln p(D, \alpha , \beta, \tau) \Bigl] }$

$= \displaystyle{ {\mathbb E}_{\alpha, \beta} \Bigl[ \ln p(D | \alpha , \beta, \tau) + \ln p(\beta) \Bigl] }$

$= \displaystyle{ {\mathbb E}_{\alpha, \beta} \Bigl[ - \frac{\tau}{2} \sum_{i = 1}^d (z_i - \alpha -\beta x_i^{\sf T})^2 - \frac{1}{2} (\beta - \mu_0) \Phi (\beta - \mu_0)^{\sf T} \Bigl]}$

となる．ここで， $\beta$ と無関係な定数を無視して計算している．ここで， $\alpha | D$ と $\tau | D$ の分布が次で与えられるとしよう．

$\alpha | D \sim {\rm Normal}(\hat{\alpha}, \theta ^{-1})$

$\tau | D \sim {\rm Gamma} (a, b)$ .

すると $\beta | D$ の更新後の分布は次のように表すことができる

$\beta | D \sim {\rm Normal}(\hat{\beta}, \Phi^{-1})$

$\displaystyle{\Phi = \Phi_0 + \hat{\tau} {\mathbb X}^{\sf T} {\mathbb X}}$

$\displaystyle{\hat{\beta} = \bigl\{ \mu_0 \Phi_0 + \hat{\tau}\sum_{i = 1}^d (z_i - \hat{\alpha} )x_i^{\sf T} \bigr\} (\Phi_0 + \hat{\tau} {\mathbb X}^{\sf T} {\mathbb X})^{-1}}$

，ただし $\hat{\tau} = {\mathbb E}_{\tau | D}[ \tau] = \frac{a}{b}$ .

$\underline{\tau {\text の更新}\ }$

更新後の $\tau$ の分布を， $q {}^\ast(\tau | D)$ とすれば，

$\ln q{}^\ast (\tau | D)$

$= \displaystyle{ {\mathbb E}_{\alpha, \beta} \Bigl[ \ln p(D, \alpha , \beta, \tau) \Bigl] }$

$= \displaystyle{ {\mathbb E}_{\alpha, \beta} \Bigl[ \ln p(D | \alpha , \beta, \tau) + \ln p(\tau) \Bigl] }$

$= \displaystyle{ {\mathbb E}_{\alpha, \beta} \Bigl[ \frac{d}{2} \ln \tau + \frac{\tau}{2} \sum_{i=1}^d (z_i - \alpha - \beta x_i^{\sf T})^2 + (a_0 -1)\ln \tau - b_0 \tau\Bigl] }$

となる．ただし， $\tau$ と無関係な定数を無視して計算している．ここで $\alpha | D$ および $\beta | D$ の分布が次で与えられるとしよう．

$\alpha | D \sim {\rm Normal}(\hat{\alpha}, \theta ^{-1})$ ,

$\beta | D \sim {\rm Normal} (\hat{\beta}, \Phi ^{-1})$ .

すると $\beta | D$ の更新後の分布は次のように表すことができる

$\tau | D \sim {\rm Gamma}(a, b)$

$\displaystyle{a = a_0 + \frac{d}{2}}$

$\displaystyle{ b = b_0 + \sum_{i = 1}^d (z_i - \hat{\alpha} - \hat{\beta} x_i ^{\sf T})^2 + \frac{d}{\theta} + {\bf tr} ({\mathbb X} \Phi ^{-1} {\mathbb X}^{\sf T} ) }$