４.2 最適化に基づく推論手法 斎藤惇

Table of contents 1. マルコフ連鎖モンテカルロ法の問題 2. 4.2.1 変分推論法 3. 4.2.2.1 線形次元削減への適用 4. 4.2.2.2 混合ガウス分布への適用 5. 4.2.3 　ラプラス近似 6. 4.2.4 　モーメントマッチングによる近似 7. 4.2.4.3 ガウス分布の例 8. 4.2.4.4 ガンマ分布の例 9. 4.2.5 例：モーメントマッチングによるプロビット回帰モデルの学 習 10. 4.2.6 期待値伝搬法 1

マルコフ連鎖モンテカルロ法の問 題

マルコフ連鎖モンテカルロ法の問題 • 長所： • 理論的には、無限回サンプルすれば真の分布を取得 • 短所： • 実用上何回サンプルすれば、よいのかが不明 • 計算コストがかかる 2

4.2.1 変分推論法

変分推論 • 変分推論　 variational inference method • p(X) = ∫ p(X, Z)dZ で潜在変数を周辺化して消去する必要 • 事後分布計算に必要な解析不可能な積分を近似的に数値計算 • 近似的に計算するために、勾配を使用した最適化の問題に置き 換える • 　エビデンス下界 ELBO evidence lower bound ln p(X) ≥ L(ξ) (1) • 上記の式の右辺を最大化することで、所望の左辺に近づける • F = −L(ξ) を変分エネルギーと呼ぶ 3

変分推論 • ELBO の方法は複数あるが、事後分布 q(Z; ξ) を近似分布 q ( Z; ξopt. ) によって近似する • 近似分布と事後分布の KLdivergence を最小化 q ( Z; ξopt. ) = argmin ξ D KL[q(Z; ξ)∥p(Z|X)] (2) • 対数周辺尤度は、ELBO を用いて、 • 4を3へ代入して KLD の定義式から計算すると Z を含む後が確率 分布の積分が 1 になるので消える。 ln p(X) = L(ξ) + D KL[q(Z; ξ)∥p(Z|X)] (3) L(ξ) = ∫ q(Z; ξ) ln p(X,Z) q(Z;ξ) dZ (4) • 3は一定の値なので、3の第２項を最小化すると4が最大化される。 4

変分推論 (近似分布) • 平均場近似 • 潜在変数の集合 {Z1, ..., ZM} の各々の要素について、事後分布の 独立性を仮定 q(Z) = M ∏ i=1 q (Zi) (5) • 長所：メモリーが少なく計算可能 • 短所：依存関係は無視しているので、精度に限界 5

例：平均場近似による潜在変数モデルの学習 • 次元削減手法 • 大きな次元 N の高次元観測データ集合 X = {x1, ..., xN} を低次元 の潜在的な集合 X = {x1, ..., xm} で表現 • 次元削減手法のその他の例 • 主成分分析、独立成分分析、行列分解、K 平均方法 6

4.2.2.1 線形次元削減への適用

4.2.2.1 線形次元削減への適用 • 事後分布を平均場近似に基づいて近似 • 以降で解説する手法は変分 EM アルゴリズム • のちに関連する内容 • 6.1.1 の変分事故符号化器 7.5.1 のガウス過程潜在変数モデル 7

4.2.2.1 線形次元削減への適用 • 線形回帰モデルと同じように観測データ X は Z の線型結合と固 定ノイズ σ2 x で記述される。 p(X|Z, W) = N ∏ n=1 p (xn|zn, W) = N ∏ n=1 N ( xn|Wzn, σ2 x I ) (6) • Z は潜在変数なので、未観測の変数の集合であり、観測された値 ではない。 • 潜在変数の生成は独立なガウス分布に従うと仮定 p(Z) = N ∏ n=1 N (zn|0, I) (7) • 上記の仮定の生成分布のパラメータもガウス分布にしたがう と仮定 p(Z, W|X) ≈ q(Z)q(W) (8) 8

4.2.2.1 線形次元削減への適用 • 真の事後分布を近似分布の積によって分解近似 p(Z, W|X) ≈ q(Z)q(W) (9) • 周辺尤度の対数の下界を計算 p(X) = ∫ p(X, Z, W)dZdW (10) • L = ∫ q(Z)q(W) ln p(X, Z, W) q(Z)q(W) dZdW = Eq(Z)q(W) [ln p(X|Z, W)] − D KL[q(Z)∥p(Z)] − D KL[q(W)∥p(W)] (11) • qi(W)i , qi(Z) がそれぞれ、ループカウンタ i ごとに更新されてい く • 片方を固定して、もう片方を更新することを繰り返す。 • q(Z)i (変分 M ステップ), q(W)i （変分 E ステップ）, を固定した場 合を計算する。 9

変分 M ステップ • q(Z)i を固定した場合を計算する。 Lqi(Z) = Eqi(Z)qi+1 (W) [ln p(X|Z, W)] − D KL [qi+1 (W)∥p(W)] + c = Eqi+1 (W) [ ln exp ( Eqi(Z) [ln p(X|Z, W)] ) p(W) qi+1 (W) ] + c = −D KL [qi+1 (W)∥ri(W)] + c (12) • ただし、 ri(W) ∝ exp ( Eqi(Z) [ln p(X|Z, W)] ) p(W) (13) ∫ ri(W)dW = 1 (14) • KL div の最小値を求めると、最適解は、以下 qi+1 (W) = ri(W) (15) 10

変分 E ステップ • L は,Z, W について対称式なので、最適解は同様の計算により求 まる。 qi+1 (Z) = ri+1 (Z) ri+1 (Z) ∝ exp ( Eqi+1 (W) [ln p(X|Z, W)] ) p(Z) (16) ∫ ri+1 (Z)dZ = 1 (17) • 線形次元削減では、ri(W), ri(Z) が解析的に計算可能 • 共役分布の計算をするはず。 11

4.2.2.2 混合ガウス分布への適用

4.2.2.2 混合ガウス分布への適用 • 離散的な潜在変数 S を用いた場合は、クラスタイリングのアル ゴリズムを導出可能 • X = {x1 , ..., xN} を K 個のグループに分ける。 • 離散の潜在変数 S = {s1 , ..., sK} • 各々の潜在変数は,sn ∈ {0, 1}K, and ∑K k=1 sn,k = 1 • 混合ガウス分布＝各々のデータ点が、K 個の異なるガウス分布に したがう分布 p(X|S, W) = N ∏ n=1 p (xn|sn, W) = N ∏ n=1 N ( xn|Wsn, σ2 x I ) (18) • 潜在変数はカテゴリ分布にしたがい、そのパラメータはガウス 事前分布にしたがう。事後分布は、平均場近似により分解する。 p(S) = N ∏ n=1 Cat (sn|π) , where p(S, W|X) ≈ q(S)q(W) (19) 12

4.2.3 　ラプラス近似

4.2.3 　ラプラス近似 • 事後分布を対数事後分布の形状を点推定の最大値のまわりのテ イラー展開で２次の後まで近似する p(Z|X) ≈ N ( Z|ZMAP, {Λ (ZMAP)}−1 ) (20) Λ(Z) = −∇2 Z ln p(Z|X) (21) Λ(Z) = −∇2 Z ln p(Z|X) (22) • ラプラス近似の手順 • 勾配降下法やニュートン・ラフソン法で対数事後分布の最大値 ZMAP を計算 • 最大値の周りで、ヘッセ行列の値を計算して、精度行列を求める。 13

4.2.3 　ラプラス近似 • テイラー展開の式 ln p(Z|X) ≈ ln p (ZMAP|X) + (Z − ZMAP)⊤ ∇Z ln p(Z|X)|z = ZMAP + (Z − ZMAP)⊤ ∇2 Z ln p(Z|X)|z = ZMAP (Z − ZMAP) = ln p (ZMAP|X) + (Z − ZMAP)⊤ ∇2 Z ln p(Z|X)|z = ZMAP (Z − ZMAP) (23) • 指数をとることにより、事後分布が近似的に導出できる p(Z|X) ∝ exp ( − (Z − ZMAP)⊤ (ZMAP) (Z − ZMAP) ) (24) 　　 • 長所：既存の計算手法（MAP 推定ベース）から、拡張が可能　　 • 短所：メモリーが必要（ヘッセ行列） 14

4.2.4 　モーメントマッチングによ る近似

4.2.4 　モーメントマッチングによる近似 モーメントマッチング：仮定密度フィルタリングや期待値伝搬 法の基礎となる概念 15

4.2.4.1 モーメントマッチング • 近似分布が簡単な指数型分布族を用いて表現可能とする。 q(z; η) = h(z) exp ( η⊤t(z) − a(η) ) (25) • 最小化したい KL D KL[p(z)∥q(z; η)] (26) • 微分して、最小値の導出 • 指数型分布族の式を KL の式に代入して計算 D KL[p(z)∥q(z; η)] = −Ep[ln q(z; η)] + Ep[ln p(z)] = −η⊤Ep[t(z)] + a(η) + c (27) • 勾配を計算して、ゼロとすると、対数分配関数の勾配が十分統計 量の期待値であるので、 Eq[t(z)] = Ep[t(z)] (28) 16

4.2.4.2 仮定密度フィルタリング • データ集合 D1 を観測後の事後分布 p (θ|D1 ) ∝ p (D1 |θ) p(θ) (29) • 尤度関数 p (θ|D1 ) と、事前分布 p(θ) の間に共役性が成立すれば p (θ|D1 , D2 ) , p (θ|D1 , D2 , D3 ) もすべて、解析的（条件つける前 と、後でのパラーメタの比較によって）に事後分布に取り込む ことが可能 • しかし、共役性が成り立たない場合もある。 q 1 (θ) ≈ r 1 (θ) = 1 Z 1 p (D1 |θ) p(θ) (30) Z 1 = ∫ p (D1 |θ) p(θ)dθ (31) 17

4.2.4.2 仮定密度フィルタリング • p (Di+1 |θ) に対する近似分布 qi (θ) を設定 (正規化項 Z 1 = ∫ p (D1 |θ) p (θ)) して、qi (θ) のパラメータを決定 • パラメータの決定は右辺のモーメントを計算して、それと一致 するようなモーメントをもつように q 1 (θ) をのパラメタ決定 q 1 (θ) ≈ r 1 (θ) = 1 Z 1 p (D1 |θ) p(θ) (32) • 追加の観測データ Di+1 に対しては、同じように近似計算 qi+1 (θ) ≈ ri+1 (θ) = 1 Zi+1 p(Di+1 |θ)qi(θ) (33) • 以降では、パラメータ θ ∈ R に関する分布 q(θ) がガウスあるい はガンマ分布の場合を計算する。 qi+1 (θ) ≈ ri+1 (θ) = 1 Zi+1 fi+1 (θ)qi(θ) (34) • 上の式の qi+1 のモーメントを ri+1 のモーメントに合わせること により、近似分布を qi から qi+1 へ更新 18

4.2.4.3 ガウス分布の例

4.2.4.3 ガウス分布の例 ・近似分布 qi(θ) = N (θ|µi, vi) (35) ・正規化項 Zi+1 = ∫ fi+1 (θ)qi(θ)dθ = ∫ fi+1 (θ) 1 √ 2πvi exp ( − (θ − µi)2 2vi ) dθ (36) 19

4.2.4.3 ガウス分布の例 ・正規化項の対数を平均パラメタで偏微分して、 ∂ ∂µi ln Zi+1 = 1 Zi+1 ∫ fi+1 (θ)N (θ|µi, vi) θ − µi vi dθ = Eri+1 [θ] − µi vi (37) ・分布 ri+1 (θ) の１次モーメントがもとまる。 Eri+1 [θ] = µi + vi ∂ ∂µi ln Zi+1 (38) 20

4.2.4.3 ガウス分布の例 ・正規化項目の対数を標準偏差パラメタで偏微分 ∂ ∂vi ln Zi+1 = 1 Zi+1 ∫ fi+1 (θ)N (θ|µi, vi) { − 1 2vi + (θ − µi)2 2v2 i } dθ = − 1 2vi + 1 2v2 i { Eri+1 [ θ2 ] − 2µiEri+1 [θ] + µ2 i } (39) ・分布 ri+1 (θ) の２次のモーメントが計算される Eri+1 [ θ2 ] = 2v2 i ∂ ∂vi ln Zi+1 + v + 2µiEri+1 [θ] − µ2 i (40) 21

4.2.4.3 ガウス分布の例 ・分布 ri+1 (θ) の１次と２次のモーメントからパラメタの更新式を導 出する ・平均のパラメタ更新 (１次のモーメントを代入) µi+1 = Eri+1 [θ] = µi + vi ∂ ∂µi ln Zi+1 (41) ・分散パラメタ更新 (１次と２次のモーメントを代入) vi+1 =Eri+1 [ θ2 ] − Eri+1 [θ]2 = vi − v2 i {( ∂ ∂µi ln Zi+1 )2 − 2 ∂ ∂vi ln Zi+1 } (42) 22

4.2.4.4 ガンマ分布の例

4.2.4.4 ガンマ分布の例 ・正規化定数 Zi+1 = Z (ai, bi) = ∫ fi+1 (θ) Gam (θ|ai, bi) dθ (43) ・１次と２次のモーメント 1 Z(ai,bi) ∫ θfi+1 (θ) Gam (θ|ai, bi) dθ = Z(ai+1,bi)ai Z(ai,bi)bi 1 Z(ai,bi) ∫ θ2fi+1 (θ) Gam (θ|ai, bi) dθ = Z(ai+2,bi)ai(ai+1) Z(ai,bi)b2 i (44) 23

4.2.4.4 ガンマ分布の例 ・１次と２次のモーメントから平均と分散を計算すると、 ai+1 bi+1 = Z(ai+1,bi)ai Z(ai,bi)bi ai+1 b2 i+1 = Z(ai+2,bi)ai(ai+1) Z(ai,bi)b2 i − { Z(ai+1,bi)ai Z(ai,bi)bi }2 (45) ・ガンマ分布のパラメータの更新式 ai+1 = { Z (ai, bi) Z (ai + 2, bi) Z (ai + 1, bi)−2 ai+1 ai − 1 }−1 bi+1 = { Z (ai + 2, bi) Z (ai + 1, bi)−1 ai+1 bi − Z (ai + 1, bi) Z (ai, bi)−1 ai bi }−1 (46) 24

4.2.5 例：モーメントマッチングに よるプロビット回帰モデルの学習

プロビット回帰モデルの学習 ・プロビット回帰：２値 yn ∈ {−1. + 1} を予測 ・プロビット回帰の尤度関数は、標準正規分布の累積分布関数 Φ を 用いて p(Y|X, w) = N ∏ n=1 p (yn|xn, w) = N ∏ n=1 Φ (yn wxn) (47) ・パラメータの事前分布は、固定の文 s なんをもつガウス分布を仮定 する。 p(w) = N (w|0, v 0 ) (48) ・以下の周辺尤度は解析的に計算できないので、モーメントマッチ ングによる近似計算を行う。 Z = ∫ p(Y|X, w)p(w)dw (49) 25

プロビット回帰モデルの学習 ・更新式 q 1 (w) ≈ p (w|y 1 , x 1 ) = 1 Z 1 p (y 1 |x 1 , w) p(w), qi+1 (θ) ≈ ri+1 (θ) = 1 Zi+1 fi+1 (θ)qi(θ (50) ・正規化定数は以下のように解析的に計算できる。 Zi+1 = ∫ p (yi+1 |xi+1 , w) N (w|µi, vi) dw =Φ (ai+1 ) (51) ai+1 = yi+1 xi+1 µi √ 1 + vi y2 i+1 x2 i+1 (52) 26

4.2.6 期待値伝搬法

4.2.6 期待値伝搬法 ・仮定密度フィルタリングでは、逐次的入力の入力順番に近似結果 が依存する ・期待値伝搬方法では、バッチ学習で仮定密度フィルタリングがで きるように一般化した。 ・確率モデルを因子の積で表記する。 p(X, θ) = p(θ) N ∏ n=1 p (xn|θ) = N ∏ n=0 fn(θ) (53) fn(θ) = { p(θ), if n = 0 p (xn|θ) , if n > 0 (54) 27

4.2.6 期待値伝搬法 ・この因子で、事後分布を書き直す。 p(θ|X) = p(θ) ∏N n=1 p (xn|θ) p(X) = ∏N n=0 fn(θ) p(X) (55) ・この事後分布に対する近似分布を、近似因子の積で表現する q(θ) = 1 Z N ∏ n=0 ˜ fn(θ) (56) ・近似因子としてガウス分布の確率密度関数をとる。 ・近似因子の積を正規化した分布もガウス r 分布となり、以下の式 q(θ) = N  θ| ( N ∑ n=0 Σ−1 n )−1 N ∑ n=0 Σ−1 n µn, ( N ∑ n=0 Σ−1 n )−1   (57) 28

4.2.6 期待値伝搬法 ・現在の近似分布 qold (θ) q λi(θ) = ∏ j̸=i ˜ fj(θ) = qold (θ) ˜ fi(θ) (58) ・この近似分布から、現在の近似因子を取り除く。そして、モデルの 因子をかけて、正規化する。 r(θ) = 1 Zi fi(θ)q \i(θ) (59) ・この r(θ) のモーメントを計算して、新しい近似分布 q new(θ) のモ ーメントとする。 ・新しく変更した近似分布を使って近似因子の更新を行う。 ・˜ fi(θ) を すべての i = 0, 1, 2, ... で実行すると、近似分布が更新される。 ˜ fi(θ) ← Zi q new(θ) q \i(θ) (60) 29