MCMCと変分近似

今回は代表的な２つの確率分布の近似推定手法であるMCMCと変分近似を比較します。変分近似に関しては複数回にわけて記事にしているのでそちらを参照されるとよいです。

変分近似（Variational Approximation）の基本（１）

変分近似（Variational Approximation）の基本（２）

変分近似（Variational Approximation）の基本（３）

さて、MCMC（Markov Chain Monte Carlo）は、サンプリング手法の一種です。サンプリングでは、解析的に計算できない事後分布の統計量などを、データをサンプリングすることによって近似的に求めます。今回はMCMCの中でも一番シンプルで便利なギブスサンプリング（Gibbs sampling）と呼ばれる手法を紹介します。

前回取り上げた2次元ガウス分布の近似問題をまた例に取り上げます。今回はこの分布を近似推定するギブスサンプリングのアルゴリズムを導き、過去に導出した変分近似による推定アルゴリズムと簡単な比較してみたいと思います。

[必要な知識]

下記をさらっとだけ確認しておくといいです。

• 前回の変分近似の内容
• 多次元ガウス分布
• KL divergence

ギブスサンプリングでは、推定したい多変数の事後分布に対して、下記のような手続きで個目のサンプル値とを拾ってきます。

他の変数のサンプル値で条件付けされた分布から新たにサンプル取る、というような手続きですね。これ以上の多変数の場合も同様です。事後分布から一度にすべての変数のサンプルを取るのが難しいのであれば1つずつ変数をサンプルしてしまえばいい、という発想です。この方法が真の事後分布からサンプルしていることを証明することができるのですが、今回は割愛します。*1

さて、変分近似の更新則と比較してみると面白いです。 

<変分近似>

- より単純な事後分布の関数形を仮定し、積分を解析的に行えるようにする。

- 注目している確率分布以外の確率分布に関して期待値を取ることにより、を得る。

<ギブスサンプリング>

- より単純な条件付き分布からサンプルし、積分をサンプルで近似する。

- 注目している確率変数以外の確率変数をサンプルされた値で固定し、をサンプルする。

こうして比較してみると、ギブスサンプリングと変分近似が兄弟のような関係になっていることがわかります。あとで見るように、実際にアルゴリズムを導いてみると同じような式が出てきます。

さて、前回と同様、今回ギブスサンプリングで推定したいのは下記のような2次元ガウス分布です。

前にも触れましたが、推定したい分布が2次元ガウス分布だとわかっているので、実はあんまり意味のない問題を解いていることになります。しかしこの単純な例を使うと、真の分布との間の近似誤差を解析的に計算できるので、２つの手法の近似精度を定量的に比較できるというメリットがあります。というわけで、今回は「1次元のガウス分布は簡単にサンプルできるけど、2次元になるとサンプルできない」という架空の想定で進めたいと思います。

では、実際に条件付き分布を求めてみましょう。を固定されたサンプル値であるとして、の確率分布を求めてみます。

これを正規化することにより1次元のガウス分布が得られます。

ただし、

です。

1次元のガウス分布によるサンプリングは「簡単」なんでしたね。変分近似と似ていますが、分布を細かくした結果が十分に簡単にサンプリングできるような形になっている必要があります*2。また実装でこの式からをサンプルするには、お使いのプログラミング言語にたぶんあるrandn()のような1次元ガウス分布からサンプルを取るライブラリーを使えばOKです。 

のサンプルも全く同様なので省略します。

さて、目的は2次元のガウス分布を近似したかったんでした。上の手続きで個のサンプルを得た後、求めたい分布の統計量（平均、分散）は次のように推定することができます。

では実際に動作させてみましょう。 

 

 青い楕円が推定したい真の2次元ガウス分布で、赤い楕円が推定中の分布です。ピンク色の点が実際のサンプルです。最終的に形状（共分散）も含めてきちんと正解の分布に収束していっていることが分かります。これは分解してしまった変数間の相関を補足できない変分近似とは対照的です。（前回の記事参照）

次にギブスサンプリングと前回導いた変分近似の近似能力を比較してみたいと思います。評価は、真の分布（2次元ガウス分布）とのKL divergenceを直接計算することにより定量的に求めたいと思います。繰り返しになりますが、今回は既知の2次元ガウス分布を使っているためこのような比較ができるのですが、実際は真の確率分布に対するKL divergenceは解析的には求まりません。

 本当は2つの手法を比較するためには横軸は実計算時間を取らないとだめなのですが、どっちもこの例では高速すぎるので繰り返し回数にしました。一般的には変分近似の方がギブスサンプリングよりも早く収束するのですが、この例だと分布が十分に簡単なので、ギブスサンプリングの方が少ない繰り返し数であっという間に変分近似（最適値.46）を抜き去っています。

以上、今回の記事からわかる範囲と個人的な見解にもとづいて変分近似とギブスサンプリングの利点欠点をまとめます。

<変分近似>

- 良い点

計算が早い。収束判定がしやすく、バグも発見しやすい（目的関数があるので）。

- 悪い点

導出がちょっと大変。近似の仕方によっては重要な相関が取れなくなる。

<ギブスサンプリング>

- 良い点

理論的には真の事後分布からサンプルしていることになる。導出が比較的簡単。並列化しやすい。

- 悪い点

大規模なモデルでは収束が遅い。またどれくらいサンプルすればいいのか判断がしにくい。

一般的には小規模問題ならギブスサンプリング、大規模問題なら変分近似、という風に解釈もできますが、計算環境や問題によっても変わってくるので出来れば両方試すのが理想です。

発展的な話題として、2つの技術を組み合わせたような手法を作ることも可能です。他にも、崩壊ギブスサンプリング（collapsed Gibbs sampling）という、日本語だと少し中2病こじらせたような名前の手法があり、これに関してもまた別の機会で書きたいと思います。

今回の記事がよくわからん！という方には，次のような入門書もあります．

books.rakuten.co.jp

*1:サンプリングの手続きが不変（invariant）であることと、エルゴディック（ergodic）であることを示す必要があります。ギブスサンプリングは証明済みであるので安心して使えます。

*2:今回の例のように必ずしも正規化できる必要はありません。よくわからない確率分布が出てきても、例えば棄却サンプリングのような方法を使えばうまく点がサンプルできる可能性があります。

変分近似（Variational Approximation）の基本（３）

「作って遊ぶ」を題目として掲げておきながらまだ作っても遊んでもいなかったので、今回はそろそろ何か動くものを載せたいと思います。

さて、前回得られた変分近似のアルゴリズムを導出するための手引きを使って、今回は世界で一番簡単だと思われる２次元ガウス分布に対して近似推定をやってみたいと思います。*1

[必要な知識]

下記をさらっとだけ確認しておくといいです。

• 前回の記事の内容
• 多次元ガウス分布

今回は２次元のガウス分布の近似推定を例として行いますが、実を言うと、多次元ガウス分布は積分も解析的にできますしサンプリングも簡単にできるような単純な分布なので、近似分布をわざわざ求める意味は皆無です。しかし、この例は計算がとてもシンプルで変分近似の導出手順を説明しやすいのと、変分近似が「近似してしまっているもの」が何なのか明確化することができるので、基本を説明するには十分な例だと思います。

では、次のような２次元のガウス分布を変分近似を使って近似推定してみましょう。

\[ p(x_1,x_2|\mu_1,\mu_2, \Lambda) = \mathcal{N} \Bigl(\left[\begin{array}{r} x_1 \\ x_2 \end{array}\right] \bigg| \left[ \begin{array}{r} \mu_1 \\ \mu_2 \end{array} \right], \left[ \begin{array}{cc} \Lambda_{1,1} & \Lambda_{1,2} \\ \Lambda_{2,1} & \Lambda_{2,2} \end{array} \right]^{-1} \Bigr) \\ \propto exp\Bigl\{ \Bigl(\left[\begin{array}{r} x_1 \\ x_2 \end{array}\right] - \left[ \begin{array}{r} \mu_1 \\ \mu_2 \end{array} \right] \Bigr)' \left[ \begin{array}{cc} \Lambda_{1,1} & \Lambda_{1,2} \\ \Lambda_{2,1} & \Lambda_{2,2} \end{array} \right] \Bigl(\left[\begin{array}{r} x_1 \\ x_2 \end{array}\right] - \left[ \begin{array}{r} \mu_1 \\ \mu_2 \end{array} \right]\Bigr) \Bigr\} \]

$x_1$,$x_2$はガウス分布からサンプルされる確率変数です。$\mu_1$,$\mu_2$は平均値、$\Lambda$は$x$の精度行列*2で、２つとも今回は適当な値で固定されたパラメータです。

さて、これをある近似分布$q(x_1,x_2)$で推定しようと思います。２変数の確率分布を分解して推定しようとしているので、次のように2つの分布に分解するしか今回は選択がないです。

\[q(x_1, x_2) = q(x_1)q(x_2) \]

さて、前回紹介した公式

\[ \ln q(z_1) = \langle \ln p(z_1, z_2| x) \rangle_{q(z_2)} + c \]

を適用してみましょう。*3

\[ \ln q(x_1) = \langle \ln p(x_1, x_2|\mu_1,\mu_2, \Lambda) \rangle_{q(x_2)} + c \\ = -\frac{1}{2} \langle x_1^{2}\Lambda_{1,1} - 2 x_1 (\Lambda_{1,1} \mu_1 - \Lambda_{1,2}(x_2 - \mu_2)) \rangle + c \\ = -\frac{1}{2} \{x_1^{2}\Lambda_{1,1} -2 x_1 (\Lambda_{1,1} \mu_1 - \Lambda_{1,2}(\langle x_2 \rangle - \mu_2)) \} + c \]

求めたいのは$x_1$に関する確率分布です。なので$x_1$にだけ注目し、無関係な項をすべて定数$c$に吸収させてしまうのが計算上のポイントです。さらにブラケット$\langle \cdot \rangle$を使って表現した期待値計算ですが、ここでは$x_2$のみに適用してあげればOKで、$x_2$に無関係な項たちはブラケットをするりと抜け出すことができます。

さて、この式をよく見てみると、$x_1$の「上に凸の2次関数」になっていることがわかります。対数計算の結果が上に凸の2次関数になっているということは、この確率分布は1次元のガウス分布であることを表しています。したがってこの式から、平均と分散を求めてあげれば近似分布が求まります。*4

\[ q(x_1) = \mathcal{N}(x_1| m_1, \Lambda_{1,1}^{-1}) \]

ただし、

\[ m_1 = \langle x_1 \rangle = \mu_1 - \Lambda_{1,1}^{-1}\Lambda_{1,2}(\langle x_2 \rangle - \mu_2) \]

です。

$q(x_2)$の計算も同様に計算でき、下記のようになります。*5。

\[ q(x_2) = \mathcal{N}(x_2| m_2, \Lambda_{2,2}^{-1}) \]

\[ m_2 = \langle x_2 \rangle = \mu_2 - \Lambda_{2,2}^{-1}\Lambda_{2,1}(\langle x_1 \rangle - \mu_1) \]

以上から、得られる疑似コードは次のようになります。

1. $m_2$をランダムに初期化する。
2. $m_1$を更新する。
3. $m_2$を更新する。
4. 以上、2と3を十分な回数まで繰り返す。

結果的に近似分布の精度は$\Lambda_{1,1}$と$\Lambda_{2,2}$のまま更新されず、平均値だけ更新されていくようなアルゴリズムになりましたね。

さて、これを実装して動かした結果が次のものです。

上図では、青い楕円が推定したい真のガウス分布で、赤い楕円が推定中の近似分布です（σ=1にあたるところで楕円を描いています）。繰り返し回数は50回にしています。ランダムな平均値から出発して、だいたい15回目くらいの更新で収束しているように見えますね。下図では対応する真の分布と近似分布との間のKL divergenceを繰り返しごとにプロットしてみました。

要点を少し挙げてみます。

１、近似分布が共分散を表現できていない

上図の赤い楕円（近似分布）は常に軸に平行になっており、絶対に斜め向きにはなってくれません。これは近似分布の独立性（分解）を仮定しているため、x1とx2の相関が表現できなくなっているためです。このように変分近似では、最初に独立性を仮定してしまった変数間の相関は取ることができません。

２、KL divergenceが単調に減少している

下図を見ればわかるように、真の事後分布と近似分布との間のKL divergenceが単調に減少していることが分かります。これは、毎回の更新でKL divergenceを最小化する方向に近似分布を修正しているからです。もしこれが増加する場合はバグなので、更新式やソースコードを見直してみる必要があります。

今回の実験に関するソースコードは時間があるときにGitHubに上げようと思います。

すっかり忘れていましたが、Juliaで実装したものをGitHubに公開しました。

MLBlog/demo_simpleVI.jl at master · sammy-suyama/MLBlog · GitHub

基本的にはJulia上で

julia> include("demo_simpleVI.jl")

と叩けば今回のような図が色々出てくるかと思います。

次回以降は、もうちょっと複雑な、でも現実的なモデルに対して変分近似を適用してみたいと思います。

[続き・関連]

MCMCと変分近似 - 作って遊ぶ機械学習。

*1:まったく同じ例がPRMLにもあります。世界で一番簡単な例なので許してください。

*2:精度行列は共分散行列の逆行列なのですが、今回の例では共分散行列を使うよりもこの精度行列の方が導出がスッキリします。

*3:前回導いた公式では右辺は観測データによって条件づけされていますが、今回はそれに該当するのはガウス分布のパラメータ（$\mu$、$\Lambda$）です。ベイズ学習における推論では、既知の固定パラメータと観測データには数学的な区別はありません。

*4:高校数学で教わった平方完成を使います。2次式を$-\frac{1}{2}(x - m)' \Lambda (x - m)$と置いて展開してあげれば、逆からどのように平均と精度を求めたらいいかがわかるかと思います。

*5:実は2つの更新式を連立方程式として解くと2つの平均値$\langle x_1 \rangle$と$\langle x_2 \rangle$が解析に求まってしまいます。今回は実際の多くの応用のように、あえて繰り返しアルゴリズムを使って解を求めます。

Googleの人工知能が囲碁のトッププレイヤーを撃破

オセロ、チェス、将棋に続き、ついに囲碁までもコンピュータが人間を上回るようになったみたいですね。

www.wired.com

概略は以下の通りです。 

・GoogleのDeepMindの研究員が開発したAlphaGoという囲碁の人工知能システムが、欧州のチャンピオンを倒した。

・深層学習（ディープラーニング）を使って過去の対局データ（棋譜）を大量に学習し、人間のように局面を「直観的に」判断するような技術を実現した。

・他にも強化学習や従来の探索アルゴリズムを組み合わせたり、AI同士を競わせるなどして戦略を向上させた。

んーなるほど。どれだけ過去に同じアプローチで研究が行われていたかはちょっと調べていないのですが、私はもう数年はかかるんじゃないかと思っていました。

・囲碁の難しさ

囲碁はオセロやチェスと比べると格段に難しい問題です。IBMが開発したディープブルーというシステムがチェスのチャンピオンを倒したのは1997年ですが、当時の技術はしらみつぶし探索に近い手法を利用して最善の手を選択していました。

しかし囲碁ではボードのサイズが19x19と非常に大きく、単純計算をするとたった7手先を読むだけで80京通り（800000000000000000通り）の組み合わせを探索しなければならず、現状ではどんなにパワフルなコンピュータであってもしらみつぶし探索は使えません。

・深層学習（ディープラーニング）を利用

今回囲碁チャンピオンを倒したアルゴリズムに関して特筆すべき点は、パターン認識のアルゴリズムである深層学習を用いているところでしょう。これは極めて自然なアプローチだと思います。囲碁のプロ棋士たちはよく「模様」という言葉を使って盤上の戦況を判断します。これは本当に文字通り「模様」で、特に対局の序盤から中盤では「白石はこのあたりが厚いから（集まっているから）有利そうだ」とか「このあたりに石を打っておけば黒の陣地拡大を抑えられる」といった、大局的な「見た目」から次の一手を決めます。終盤になってくると、可能な手数が限られてくるので、しらみつぶし探索のようなロジカルな読みに次第に移行していきます。大局的な（直観的な）局面の判断と、局所的な（精緻な）読みの組み合わせを要するところが囲碁の難しいところであり、面白いところでもあります。今回のアルゴリズムでは、このような人間の棋士が普通に行っているような「直観的な判断」を、画像認識などで使われているパターン認識の技術を用いてコンピュータに学習させたわけです。

・大量の棋譜データで学習

さて、学習の仕組みがあるだけではアルゴリズムは強くなりません。特筆すべきもう1つのポイントは、充実した過去のプロ棋士の対局データ（棋譜）を利活用したことです。これをアルゴリズムに与えることにより、どのような局面でどのように手を打ったら良いのかを学習することができます。このアプローチは将棋のAIでも同じで、電脳戦でプロ棋士に勝ったコンピュータは大量の棋譜データから指し手の方針を学習しています。

・現時点で「最良」のアプローチではあるが「最高」ではない

ここからは私の個人的な考えです。私はこの深層学習を使ったアプローチが現時点では「最良」ではありますが、囲碁を解くための「最高」のアプローチではないと思っています。強い囲碁AIを作るという課題は、深層学習が得意とする画像や音声などの認識技術とは決定的に違う点があります。それは目標の設定の仕方です。

例えば、画像認識は、写真に写っている物体が何なのかを当てるといった問題ですが、これは人間が正しいラベルを与えているので、人間の認識能力を実現することが最終的なゴールになります。したがって深層学習のような人間の認識プロセスを模倣するようなアプローチがおそらく一番素直かつ最適なアプローチであると言えるでしょう。

しかし囲碁はゴールが異なります。囲碁は厳格なルールなもとに作られた、極めて人工的な問題です。人間の指し手を模倣すること自体がゴールではありません。囲碁では、究極的にはしらみつぶし探索が最強であり、人間は過去の長い歴史の経験をもとに「そこそこ効率的な」探索手法を発見し、それを棋譜や定石にまとめたに過ぎません。極端な話になりますが、例えば量子コンピュータなどが実現して極めて困難な探索問題をしらみつぶしに近い手法で高速に解けるようになったとしたら、今のプロ棋士や今回のテクノロジーとは全く異なる指し筋になるでしょうし、比べ物にならないほど強いものになると思われます。

人間のプレイヤーを倒してしまった以上、今後この研究が進むかどうかはわかりません。しかし、人間もAI技術も囲碁という問題に対してある程度の近似解を得たに過ぎず、我々はまだ囲碁を制していないと言えると思います。

最尤推定、MAP推定、ベイズ推論

今回は、最尤推定、MAP推定（事後確率最大化推定、正則化）、ベイズ推論*1の関係性を見ていきたいと思います。結論から言うと、最尤推定とMAP推定はベイズ推論の特殊な近似方法であると見ることができます。

[必要な知識]

下記をさらっとだけ確認しておくといいです。

• KL divergence
• 確率の加法定理、乗法定理

\[ \newcommand{\argmax}{\mathop{\rm arg~max}\limits} \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \]

$x$を観測データ、$\theta$をパラメータとした確率モデル$p(x, \theta)$に対して、それぞれの推定方法は一般的には下記のように認識されているようです。

・最尤推定

\[ \theta_{ML} = \argmax_{\theta} \{ p(x|\theta) \} \tag{1} \]

・MAP推定（正則化）

\[ \theta_{MAP} = \argmax_{\theta} \{ p(x|\theta)p(\theta) \} \tag{2} \]

・ベイズ推論

\[ p(\theta|x)=\frac{p(x|\theta)p(\theta)}{p(x)} \tag{3} \]

ベイズ推論では基本的にベイズの定理を用いて事後分布を求めているだけであり、argmaxのような関数の「最適化」は含まれていません。*2

まずはベイズ推論に比較的「近い」MAP推定から見ていきましょう。MAP推定をベイズ推論の近似であるとすると、事後分布$p(\theta|x)$に次のようなパラメトリックな分布を仮定していることになります。

\[ p(\theta | x) \approx q(\theta|\hat{\theta}) = \delta(\theta - \hat{\theta}) \]

ここで$\delta(\theta)$はデルタ分布であり、次のような性質を持ちます。

\[ \delta(\theta) = \begin{cases} +\infty & (\theta=0) \\ 0 & (otherwise) \end{cases} \]

\[ \int_{-\infty}^{+\infty} \delta(\theta) d\theta = 1 \]

\[ \int_{-\infty}^{+\infty} \delta(\theta)f(\theta) d\theta = \langle f(\theta) \rangle_{\delta(\theta)} = f(0) \]

簡単に言ってしまうと、$\theta=0$の点で無限大に尖った形状を持つような確率分布です。したがって$\delta(\theta - \hat{\theta})$は$\theta = \hat{\theta}$の時に無限大の値を持つような確率分布です。

さて、この仮定のもとで真の事後分布とのKL divergenceを最小化してみましょう。

\[ KL(q(\theta|\hat{\theta})||p(\theta|x)) \\ = \langle \ln q(\theta|\hat{\theta}) \rangle_{q(\theta|\hat{\theta})} - \langle \ln p(\theta|x) \rangle_{q(\theta|\hat{\theta})} \\ = \ln q(\hat{\theta}|\hat{\theta}) - \ln p(x|\hat{\theta}) - \ln p(\hat{\theta}) \\ = - \{ \ln p(x|\hat{\theta}) + \ln p(\hat{\theta}) \} + c \]

$\ln q(\hat{\theta}|\hat{\theta})$は$\hat{\theta}$の値にかかわらず無限大の値を取るため、 定数項$c$に吸収させました。さて、この式を$\hat{\theta}$に関して最小化するので、

\[ \theta_{MAP} = \argmin_{\hat{\theta}} \{ KL(q(\theta|\hat{\theta})||p(\theta|x)) \} \\ = \argmax_{\theta} \{ \ln p(x|\theta) + \ln p(\theta) \} = \argmax_{\theta} \{ p(x|\theta)p(\theta) \} \]

となり、(2)が得られます。

次に最尤推定を見てみましょう。最尤推定をベイズ推論の近似であるとすると、事後分布$p(\theta | x)$と事前分布$p(\theta)$に次のような仮定を置いていることになります。

\[ p(\theta | x) \approx q(\theta|\hat{\theta}) = \delta(\theta - \hat{\theta}) \]

\[ p(\theta) = c \]

事後分布が無限大に尖っている分布に対して、事前分布は無限に平坦な無情報事前分布です*3。さて、MAP推定の場合と同様、真の事後分布とのKL divergenceを最小化してみます。

\[ KL(q(\theta|\hat{\theta})||p(\theta|x)) \\ = \langle \ln q(\theta|\hat{\theta}) \rangle_{q(\theta|\hat{\theta})} - \langle \ln p(\theta|x) \rangle_{q(\theta|\hat{\theta})} \\ = \ln q(\hat{\theta}|\hat{\theta}) - \ln p(x|\hat{\theta}) - \ln p(\hat{\theta}) \\ = - \{ \ln p(x|\hat{\theta}) \} + c \]

ここでも$\hat{\theta}$に無関係な項は$c$にまとめました。この式を$\hat{\theta}$に関して最小化したいので、

\[ \theta_{ML} = \argmin_{\hat{\theta}} \{KL(q(\theta|\hat{\theta})||p(\theta|x)) \} \\ = \argmax_{\theta} \{ \ln p(x|\theta) \} = \argmax_{\theta} \{ p(x|\theta) \} \]

となり、(1)が得られます。

以上のように、ベイズ推論の観点からすると、MAP推定や最尤推定は、事前分布や事後分布に対して非常に極端な分布を仮定して推定していることになります。

・MAP推定：

　事後分布を無限に尖った確率分布であると仮定

・最尤推定：

　事後分布は無限に尖った確率分布であり、事前分布は無限に平坦な確率分布であると仮定

事実、KL divergenceに基づく近似の観点からすると、どちらの推定方法も真の事後分布から「無限大に遠い」ような近似事後分布を計算していることになってしまっています。
ただし2つの手法にも利点はあります。それはデルタ関数を使うことによって難しい積分計算（期待値計算）を回避することができるという点です。積分さえ回避してしまえば、あとはKL divergenceを勾配法などを用いて最小化すればOKです。勾配法を使うには微分を要しますが、一般的に積分をするよりははるかに簡単です。

*1:ここでは推定（estimation）と推論（inference）を分けています。推定はある特定の値を何らかの方法で求めるのに対して、ここでは推論は確率計算によりある確率分布を求めることを意味しています。ちなみにMAP推定のことをベイズ推定と呼ぶ場合もあり、非常にややこしいことになっています。

*2:事後分布が解析的に求められない場合に限り、変分近似やラプラス近似のような最適化手法を使います。

*3:厳密に言うと、連続値を取る確率分布の場合、定数cを$-\infty$から$+\infty$まで積分すると無限大に発散してしまうため、$p(\theta) = c$は正しい確率分布ではないです。ここでは便宜上、$\theta$に事前に偏りがないという意味で定数$c$を置いています。

変分近似（Variational Approximation）の基本（２）

さて、前回は変分近似の目的（複雑過ぎて解析解が得られないような確率分布の近似）と、近似のための指標（KL divergence）に関して解説しました。

今回は、変分近似の「公式」を導いてみたいと思います。近似分布$q(z)$に関して「分解の仮定」を置くことにより、数値最適化における勾配法のような繰り返しの更新手続きが得られ、近似分布が数値的に計算できることを示します。

[必要な知識]

下記をさらっとだけ確認しておくといいです。

• 前回の記事の内容
• 勾配法

前回は、近似分布と真の事後分布の間のKL divergenceをなるべく小さくすることで、うまく近似分布を得ようという方針にしました。

\[ KL(q(z)||p(z|x)) = - \int q(z) \ln \frac{p(z|x)}{q(z)} dz \]

で、ここが一番のミソなのですが、変分近似では、近似したい分布を次のように複数の分布に分解されるという仮定を置きます。

\[ p(z|x) \approx q(z) = \prod_{i} q(z_i) \]

今回は簡単のために、２つの分布に分解するということにしておきましょう。*1

\[ q(z) = q(z_1)q(z_2) \]

変分近似におけるKL divergenceの最小化方法は、数値最適化における偏微分を使った勾配法と似たような繰り返し更新の戦略を取ります。つまり、$q(z_1)$と$q(z_2)$を一度にえいっと更新するのではなく、一方を固定して一方を更新、というのを繰り返します。

1. $q(z_2)$をランダムに初期化する。
2. $q(z_2)$を固定した上で、KL divergenceを最小化する$q(z_1)$を求める。
3. $q(z_1)$を固定した上で、KL divergenceを最小化する$q(z_2)$を求める。
4. 以上の2，3を十分な回数まで繰り返す。*2

 では、実際に上記のステップ2のための更新式を導いてみましょう。（ステップ3はステップ2と同じです。）

ここからちょっと数学になります。$q(z_1)$を更新するためには、$q(z_2)$を固定し、KL divergenceを$q(z_1)$のみの関数であると考えます*3。 式を式に代入し、KL divergeneの式から$z_1$に無関係な項を$c$にまとめてしまうと、

\[ KL(q(z)||p(z|x)) = - \int \int q(z_1)q(z_2) \ln \frac{p(z_1, z_2|x)}{q(z_1)q(z_2)} dz_1dz_2 \\ = - \int \int q(z_1)q(z_2) \{ \ln p(z_1, z_2|x) - \ln q(z_1) - \ln q(z_2) \} dz_1dz_2 \\ = - \int q(z_1) \bigl\{ \int q(z_2) \ln p(z_1, z_2|x) dz_2 - \int q(z_2) \ln q(z_1)dz_2 \bigr\}dz_1 \\ + \int q(z_1)q(z_2)\ln q(z_2) dz_1dz_2 \\ = - \int q(z_1) \{\int q(z_2) \ln p(z_1, z_2|x) dz_2 - \ln q(z_1) \}dz_1 + c \;\;\;\;(1) \]

積分計算が入れ子になっているので注意してください。さらに計算を進めると、この式は結局次のような開始地点とはまた別のKL divergenceに落ち着くことが分かります。

\[ (1) = - \int q(z_1) \ln \frac{exp \{\int q(z_2) \ln p(z_1, z_2|x) dz_2 \}}{q(z_1)}dz_1 + c \\ =KL(q(z_1) || exp \{ \int q(z_2) \ln p(z_1, z_2|x) dz_2 \} ) +c \]

KL divergenceは2つの分布が一致するとき、最小値0を取ります。したがって式の最小値は次のようにして得ることができます。

\[ \ln q(z_1) = \int q(z_2) \ln p(z_1, z_2|x) dz_2 + c \]

ここでの定数項$c$は分布を正規化（積分が1になる）するための定数です。上記の期待値計算（積分）はちょっと長ったらしいので、ブラケット$\langle \cdot \rangle$を使って書き直すのが便利です。

\[ \ln q(z_1) = \langle \ln p(z_1, z_2|x) \rangle_{q(z_2)} + c \]

一応この式が変分近似の公式になります。「自分以外の分布で期待値を取る」と覚えておけばOKです。

実際の応用では、具体的な確率モデルを設定し、上の期待値計算を頑張って手計算します。一つ大事なポイントは、この式を評価した結果が簡単に正規化できる（積分できる）確率分布になっていることです。これに関しては実際、モデルや分解の仮定を色々変えて手計算してみるという、地道な試行錯誤が必要な場合もあります。一方で、線形ガウスモデルや混合モデル、パラメータに対する共役事前分布の設定など、うまく計算できることが良く知られているモデルの構築方法があるので、よっぽどのことがない限り、そういった便利な部品を組み合わせるのが良いと思います。*4

さて、ひとまず確率モデルに基づく変分近似の適用方法をまとめておきましょう。

1. 確率モデルを定義する。（混合ガウス分布、HMMなど）
2. 事後分布に関する分解の仮定をする。
3. 公式を使って更新式を導出する。
4. 実装して動かす。*5

一応変分近似を使った近似アルゴリズムを動かすのに必要な知識はこれだけなのですが、なんだか抽象的で、なんだかしっくりこないなぁという感じがするんじゃないかと思います。

次回は、この方法を使って簡単な変分近似の使用例と、簡単ではない（けど現実的な）変分近似の使用例を紹介したいと思います。

[続き・関連]

変分近似（Variational Approximation）の基本（３） - 作って遊ぶ機械学習。

*1:実際の応用では、具体的なモデルに合わせて適切な分解を思いつく必要があります。例えば標準的な混合モデルでは、パラメータと隠れ変数の２つに分解すれば効率の良い近似アルゴリズムが得られます（EMアルゴリズム）。もっと複雑なモデルだと、２つ以上の細かい分解が必要になってくる場合や、すべての変数レベルでバラバラにする場合（平均場近似、mean field approximation）もあります。基本的には計算ができる最低限の分解数に収めていくのが最適化の観点では良いのですが、より多くのメモリを使用するという欠点もあります。

*2:繰り返し計算を行う回数は、"変分下限"と呼ばれる実数値を毎ステップごとに評価して決めるのが一般的です。変分下限は変分近似における目的関数で、具体的な計算方法はまた別の機会に解説します。

*3:この場合KL divergenceは関数$q(z_1)$の関数になっているので、汎関数と呼ばれます。汎関数を最適化する問題なので"変分法"という名前がついています。今回の記事ではKL divergenceの一般的な性質を使って更新則を導きますが、数学的な変分を使って同じ更新則を導くことも可能です。

*4:公式の右辺が正規化できるようにモデルを組むのが理想ですが、どうしても正規化できないような複雑なケースにはさらなる近似を考えるしかありません。例えば、さらに細かい分解の仮定を考えてみる、$q(z)$に適当なパラメトリックな分布（ガウス分布など）をおいて勾配法などの数値最適化法を適用してパラメータを求める、サンプリングを用いて必要な$q(z)$の統計量を求める、ラプラス近似を用いる、などです。

*5:ここでも変分下限の計算も実装しておくと、デバッグに便利です。また別の機会に紹介します。

変分近似（Variational Approximation）の基本（１）

初回の記事で変分近似はけっこう重たいのですが、今後ここで頻繁に使っていこうと考えているのでとりあえずご紹介です。

変分近似（variational approximation）とは、確率分布を近似的に求める方法のひとつです*1。一般的には確率分布を求めるには正規化（積分して１になるようにする）しなければならないのですが、複雑な分布（例えば潜在変数モデルの事後分布）になってくると、どうしても解析的に積分ができなくなってしまいます。変分近似ではこのような複雑すぎて正規化できないような確率分布を、もっとシンプルな確率分布たちの積に分解する（＝独立性を仮定する）ことにより近似します。分解を仮定することによって変数の依存関係を簡略化し、数値最適化でいうところの偏微分を使った勾配法と似たようなことが確率分布の推論に対しても行えるようになります。

これが使えるようになると、様々なデータサイエンスの課題に合わせて確率モデルを作り*2、自分で自由に分布推定ができるようになります。実際に、画像や音声、金融データ、生命情報、自然言語、各種センサーデータなど、現在まででほぼすべての機械学習の問題に適用されてきています。

[必要な知識]

下記をさらっとだけ確認しておくといいです。

• 確率の加法定理（sum rule）と乗法定理（product rule）、ベイズの定理（Bayes' theorem）
• KL divergence 

今、次のような確率モデルを考えたいと思います。

\[ p(x,z) \]

$x$は観測データで、$z$は推定したい未知の変数（欠損値やパラメータ、未来の予測値など）で、ともに多次元ベクトルってことにしておいてください。今回は連続値を取る変数を仮定しますが、離散値でもまったく同じ議論になります。

機械学習の目的はzの事後分布$p(z|x)$を下記のようにベイズの定理を用いて推定することです。

\[ p(z|x) = \frac{p(x|z)p(z)}{p(x)} = \frac{p(x|z)p(z)}{\int p(x,z) dz} \tag{1} \]

例えば普通に$x$と$z$がともにガウス分布に従うようなモデルでは、式(1)の分母の積分が公式を使えば簡単に行えるので、事後分布は普通に手計算で一発で解けます。これを解析的に解けるとか、closed formで解けるとかって言います。

ただし、今回はこれがどうしてもできないと仮定します。つまり式(1)の積分計算がめちゃくちゃ複雑で、解析解が得られない状態にあるとします。

こういうときに登場するのが変分近似のような近似推論法です。事後分布を次のような別の関数形で近似します。

\[ p(z|x) \approx q(z) \]

$q(z)$の具体的な関数（ガウス分布だとか）は仮定していないことに注意してください。

今やりたいことは、$q(z)$が$p(z|x)$と、なるべく「似る」ようにしたいということです。

２つの確率分布がどれだけ「似ていないか」を表す指標の１つとして、KL divergenceがあります。例えば、混乱を避けるために確率変数$w$を一時的に使うと、確率分布$p(w)$と$q(w)$の間の（q(w)から見た*3）KL divergenceは

\[ KL(q(w)||p(w)) = - \int q(w) \ln \frac{p(w)}{q(w)} dw \]

のように定義されます。$q(w)=p(w)$が成り立つときこの式は0になります。

今回は二つの確率分布$p(z|x)$と$q(z)$をなるべく「似せ」たいので、この2つの確率分布の間のKL divergenceを最小化することにより目的を達成したいと思います。つまり、

\[ KL(q(z)||p(z|x)) = - \int q(z) \ln \frac{p(z|x)}{q(z)} dz \]

を最小にするような$q(z)$を求めることが目標になります。

しかしここで疑問が残ります。

$p(z|x)$は、積分こそできないものの、確かに何らかの形状が存在するような確率分布です。しかし最初の仮定の通り、この分布は直接手計算をして求めることはできない。直接計算できない分布と、近似分布$q(z)$の間の距離を、いったいどうやって縮めるのか？

ちょっと長くなったのでここでいったん切ります。

[続き・関連]

変分近似（Variational Approximation）の基本（２） - 作って遊ぶ機械学習。

今回の記事がよくわからん！という方には，次のような入門書もあります．

books.rakuten.co.jp

*1:他にも、変分推論（variational inference）とかただ単に変分法（variational method）とかって呼んだりもします。ベイズモデルであることを強調する場合には、変分ベイズ（variational Bayes）と呼ぶこともあります。

*2:確率モデルの作り方に関してはグラフィカルモデルの記事をご参考ください。

http://machine-learning.hatenablog.com/entry/2016/02/10/184755

*3:一般に$KL(q||p)$と$KL(p||q)$は一致しません。

開設しました。

機械学習に関するブログを開設しました。

確率モデルに基づく機械学習の基本的なテクニックの紹介から、データサイエンスに関する一般的な話題まで取り上げたいと思っています。

・なぜ機械学習のブログ？

このブログを始めるに至った理由は、確率モデルを使ったシンプルで便利な機械学習の技術を、日本のあらゆる分野の技術者や学生にも使ってもらいたいと思ったからです。

ビッグデータやIoTといった言葉に代表されるように、21世紀は蓄積された大規模データをもとに予測や推定などの解析を行う、いわゆる「データサイエンス」の時代になると言われています。ただ残念ながら現時点では、データの量やバリエーションの急速な増加に対して、効率の良いデータ解析環境（知識、人材、ツール）が整っていないという状況にあります。データという広大な未知の新大陸がすぐ近くにあるにも関わらず、それを開拓するためのノウハウを持っている技術者が非常に少ないのが現状です。

このブログでは、こうした未知の課題を前にして呆然と立ちすくんでしまうのではなく、確率モデルという道具を使いこなすことにより極めて「フォーマルに」取り組めることを示したいと思っています。

・作って遊ぶ！

このブログでは「作って遊ぶ」をモットーに、解きたい課題に合わせた確率モデルの構築方法とその推論手段（近似アルゴリズムなど）に焦点を当てたいと思っています。確率モデルは非常に便利なツールです。レゴブロックのように部品と部品を自由に組み合わせることによってありとあらゆる課題に対してアプローチができる柔軟性を持っています。作曲家はコードやリズムの知識を利用して、自由に表現したい音楽を紡ぎだすことができます。優れたプログラマーは仕様さえはっきりすれば、絵を描くように軽快にソースコードを書き上げてしまいます。それと同じようなことがデータサイエンスでも可能であり、データの特性を精査し課題（何を予測・推定したいか）を決めれば、あとは確率モデルという道具を使って自由自在に解析アルゴリズムを作ることができます。確率モデルはデータサイエンスの世界で「自由」を手にするための手段です。

また、「作って遊ぶ」を体現するためのもう１つの手段として、ブログではなるべくソースコードと実行結果を載せていこうと思っています。数式に関する導出がいまいちピンと来ない場合でも、ソースコードを見ればアルゴリズムが一体何をやっているのか直感だけでもつかむことができるはずです。

・他の話題も

確率モデルの他にも、ディープラーニングなど最近よく話題になる機械学習技術に関しても取り上げるかもしれません。ただ、ディープラーニングや人工知能と呼ばれる技術は応用先が非常に限られているという欠点があります。これからのIoTの時代では、人間が経験したことのないデータ（例えばある産業用機器に取り付けられた大量のセンサーデータなど）に対する解析技術が必要になってきます。こういった課題に対してわざわざ「人間の脳」をシミュレートする人工知能のような技術を使うのは少しナンセンスでしょう。とはいうものの、こういった技術は伝統的な画像や音声の識別問題に対しては非常に自然なアプローチであり、実際に良い性能も出しているので、機会があれば触れたいと思います。

また小難しい理論だけでなく、経営者やリーダーの方が組織のデータサイエンスの取り組みに関してちょっと頭の片隅に入れておいた方がいいような注意点なんかもお話しできればなぁと思っています。データを使って嘘をつくのはものすごく簡単です。これからの５年１０年、大量の「なんちゃって」データサイエンティストが登場し、あの手この手で判断者を騙そうとしてくると思います（本人が騙していると気づいてない場合も多いので非常にやっかいです！）。そういった状況の中、正しい判断を行うためには多少のデータに対する「感覚」が必要になってくるでしょう。

・キーワード

参考として、ブログの話題の中心になりそうなキーワードをまとめておきます。

機械学習、パターン認識、人工知能、データサイエンス、ビッグデータ、確率分布、ベイズ学習、潜在変数モデル、ベイジアンノンパラメトリクス、確率測度、ガウス過程、ディープラーニング、最適化、変分近似、MCMC、時系列解析、ネットワークモデル、強化学習、etc

・参考文献

このブログでは紹介しきれなかった基本事項や，発展的な応用モデルなどは下記の書籍に載っていますので是非ご参考ください．books.rakuten.co.jp

より発展的な内容や，ディープラーニングとの関連性に関しては下記の書籍が詳しいです．

books.rakuten.co.jp

よろしくお願いいたします。

Sammy