３１．GaP-NMFでベイズ推定を理解しよう - keep learning blog（キープラーニングブログ）

――ふるさとは遠きにありて思ふもの。そして悲しくうたふもの。――

室生犀星「抒情小曲集」

GaP-NMFでベイズ推定を理解しよう

GaP-NMFでベイズ推定を理解しよう

1年間ご愛読ありがとうございました的なことを言ってから、舌の根も乾かぬうちに更新します。実は私、もう日本に帰国して自宅で研究活動を続けています。

日本に帰ってきて最初の感想は、湿気が凄いことと、物価が安いことでした。PCR検査は陰性だったので、食事を調達するための外出は許されているのですが、食品がとにかく安い。え？これが100円？という謎の感動がありました。

ちなみに、政府のCOVID-19拡大防止策により、海外から帰国した人は原則二週間の自主隔離が要請されています。そのため、7月1日から元の職場に復帰するためには、二週間早く帰国せねばならなかったわけです。

さて、研究活動といっても、論文は私の手を離れましたし、新しいことを始めるには短すぎる期間だし、会社への報告書もほぼ作成済みなので、やれることと言えば論文を読んで勉強することくらいしかありません。

今回の記事は、前々から興味があった確率論的アプローチの機械学習を実践的に理解するために、GaP-NMFというアルゴリズムの背景を整理して、Pythonで実装してみよう、という内容です。

f:id:yuki0718:20200621143207p:plain — マンション広告風の北欧留学チラシ

NMFの確率論的解釈

以前の記事で、音声処理やコンピュータビジョンの分野で発展した、非負値行列因子分解（NMF: Non-negative matrix factorization）という手法があることをご紹介しました。

keep-learning.hatenablog.jp

NMFは、対象となるデータ行列S（音声ならSpectrogram）が、行数K×列数Mの基底行列Hと行数M×列数Nの活性化行列Uとの行列積 $S \approx H \times U$ で近似できることを仮定するのでした。

ここで、任意の基底数Mは、Sの行列サイズよりも小さい（M＜min(K, N)）ものとされ、HとUの更新アルゴリズムを回す前に決めておくべき固定値となります。

このとき、HとUの更新アルゴリズムは、以下のようなデータ行列Sと近似値 $H \times U$ との乖離度D（Deviation）を最小化するように決定されるのでした。

【二乗和誤差（Euclidean distance）】

$\displaystyle D(S \| H \times U)=\sum_{k=1}^{K} \sum_{n=1}^{N} \{ S_{k,n}-(H \times U)_{k,n} \} ^{2}$

【一般化カルバック・ライブラー情報量（Kullback–Leibler divergence）】

$\displaystyle D_{KL}(S \| H \times U)=\sum_{k=1}^{K} \sum_{n=1}^{N} \left\{ S_{k,n} \log{\frac{S_{k,n}}{(H \times U)_{k,n}}}-S_{k,n}+(H \times U)_{k,n} \right\}$

【板倉・斎藤擬距離（Itakura–Saito divergence）】

$\displaystyle D_{IS}(S \| H \times U)=\sum_{k=1}^{K} \sum_{n=1}^{N} \left\{ \frac{S_{k,n}}{(H \times U)_{k,n}}-\log{\frac{S_{k,n}}{(H \times U)_{k,n}}}-1 \right\}$

ここまでは過去記事の復習です。今回は、これを確率論の立場から解釈するため、以下のような条件付き確率を考えてみます。

$\displaystyle P(S \| H, U)=\prod_{k=1}^{K}\prod_{n=1}^{N} Gauss \bigl(S_{k,n} \| (H \times U)_{k,n}, \sigma^{2} \bigr) \propto \prod_{k=1}^{K}\prod_{n=1}^{N} e^{-\frac{1}{2\sigma^{2}} \bigl\{ S_{k,n}-(H \times U)_{k,n} \bigr\} ^{2} }$

$\displaystyle P(S \| H, U)=\prod_{k=1}^{K}\prod_{n=1}^{N} Poisson \bigl(S_{k,n} \| (H \times U)_{k,n} \bigr) \propto \prod_{k=1}^{K}\prod_{n=1}^{N} (H \times U)_{k,n}^{S_{k,n}} e^{-(H \times U)_{k,n}}$

$\displaystyle P(S \| H, U)=\prod_{k=1}^{K}\prod_{n=1}^{N} Exponential \bigl(S_{k,n} \| (H \times U)_{k,n} \bigr) \propto \prod_{k=1}^{K}\prod_{n=1}^{N} \frac{1}{(H \times U)_{k,n}} e^{-\frac{S_{k,n}}{(H \times U)_{k,n}}}$

これらは、データSの各行列要素 $S_{k, n}$ が、確率分布関数（ガウス分布、ポアソン分布、指数分布）から独立に生成されると仮定した場合の生成モデル（Generative model）を表しています*1。

以前の記事でご紹介したとおり、生成モデルの確率分布関数 $P(S \| \theta)$ が定義されたとき、そのモデルを特徴付けるパラメータ $\theta$ は最尤推定法という手法で計算可能なのでした。

keep-learning.hatenablog.jp

最尤推定法とは、データSの尤もらしさを表す尤度（Likelihood） $P(S \| \theta)$ の対数をとり、これを最大化するように $\theta$ を決定する手法でした。今の場合、パラメータ $\theta$ には $H$ と $U$ の各成分が対応します。

実は、この対数尤度を計算してみると、上式で天下り的に定義しておいた生成モデル3種類の対数尤度は、（HとUに依らない定数項を除いて）3種類の乖離度Dに負号を付けたものと一致することが分かります。

したがって、NMFで典型的に用いられる3種類の乖離度Dを最小化するようにHとUを決定することは、それぞれに対応する生成モデルの対数尤度を最大化するようにHとUを決定すること（最尤推定法）と数学的に同値だと分かります。

単に確率論の立場から同じことを解釈し直しているだけなので、得られるHとUの更新アルゴリズムは全く同じです。しかし、この確率論的なアプローチは、ベイズ推定の手法をNMFに適用するための入り口になります。

事後分布と変分ベイズ法

NMFを確率論の立場から解釈できたので、続いてはベイズ推定（Bayesian inference）をする際に便利な手法の一例*2として、変分ベイズ法（Variational Bayesian methods）というアプローチをご紹介します。

まず最初に、ベイズ推定の話をする上では、条件付き確率の性質（周辺化、同時確率、ベイズの定理等）を理解しておくことが基本となります。条件付き確率の性質については以下の記事で説明しているので、ここでは説明を割愛します。

keep-learning.hatenablog.jp

まずはゴールからお話します。ベイズ推定の最終目標は、あるデータ $S$ が与えられた（観測された）状態で、モデルのパラメータ $\theta$ の確率分布関数 $P(\theta \| S)$ を計算することです。

一般に、この確率分布関数 $P(\theta \| S)$ のことを、事後分布（Posterior probability distribution）と呼びます。前節でご説明した生成モデルの確率分布関数は、 $P(S \| \theta)$ （条件変数の順序が逆）であることに注意してください。

さて、今回は（そしてほとんどのケースでは）生成モデルを採用しているので、以下のようにベイズの定理を用いて事後確率を間接的に計算することになります。

$\displaystyle Bayes' theorem: P(\theta \| S)=\frac{P(S, \theta)}{P(S)}=\frac{P(S \| \theta)P(\theta)}{P(S)}$

右辺の分子に出てくる $P(\theta)$ を事後確率と対比して事前分布（Prior probability distribution）と呼びます。事前分布は、多くの場合にパラメータ $\theta$ の性質を考慮し、かつ生成モデルとの相性が良いものが選択されます*3。

問題は右辺の分母 $P(S)$ です。これは確率の和が1になるように定められる正規化因子なので、分子の和さえ計算できれば解析的に答えが出ます。計算できるならそれで終わりです。今回は、 $P(S)$ の計算ができない場合の近似アプローチです。

その近似アプローチが変分ベイズ法です。ちょっと導入が唐突になりますが、計算できないと申し上げた $P(S)$ の対数から出発して、任意の確率分布関数 $Q(\theta)$ による下限を設定するところから始めます。すなわち、

$\displaystyle \log{P(S)}=\log{\sum_{\theta} P(S, \theta)}=\log{\sum_{\theta} \frac{Q(\theta)P(S, \theta)}{Q(\theta)}} \geqq L(Q(\theta))$

ここに出てきた下限 $L(Q(\theta))$ はいわゆる汎関数（関数の関数）であるため、変分下限（Variational lower bound）と呼ばれます。変分下限は、Jensenの不等式*4という有名な数学公式によって以下のように定義されます。

$\displaystyle L(Q(\theta)) \stackrel{\mathrm{def}}{\equiv} \sum_{\theta} Q(\theta)\log{\frac{P(S, \theta)}{Q(\theta)}}$

この変分下限が何を意味しているかを正しく理解するには、不等式の左辺と右辺の差を取ってみると分かります。つまり、解析的に計算できない $\log{P(S)}$ と変分下限 $L(Q(\theta))$ との差を計算します。

$\displaystyle \log{P(S)}-L(Q(\theta))$

まず、Jensenの不等式で導入された任意の関数 $Q(\theta)$ は、確率分布関数なので $1=\sum_{\theta} Q(\theta)$ を常に満たします。よって、1を第一項に挿入しても値は変わらず、

$\displaystyle =\sum_{\theta} Q(\theta)\log{P(S)}-\sum_{\theta} Q(\theta)\log{\frac{P(S, \theta)}{Q(\theta)}}$

対数の性質から、

$\displaystyle =\sum_{\theta} Q(\theta)\log{\frac{Q(\theta)}{P(S, \theta)/P(S)}}$

同時確率を周辺確率で除算したものは条件付き確率に等しいので、

$\displaystyle =\sum_{\theta} Q(\theta)\log{\frac{Q(\theta)}{P(\theta \| S)}}$

この式は、事後分布 $P(\theta \| S)$ とJensenの不等式で導入された任意関数 $Q(\theta)$ との間の差異を測る相対エントロピー*5の形をしており、両者が完全に一致したとき最小値0をとります。

したがって、相対エントロピーを最小化するように任意関数 $Q(\theta)$ を選ぶことができれば、 $Q(\theta)$ はベイズ推定の最終目標である事後分布 $P(\theta \| S)$ の近似値になることが期待されます。

また、計算できないとはいえ $P(S)$ は任意関数 $Q(\theta)$ に依存しない定数なので、相対エントロピーの最小化問題は、そのまま変分下限の最大化問題と等価です。

以上をまとめると、「変分下限 $L(Q(\theta))$ を最大化するような $Q(\theta)$ を求めれば、その $Q(\theta)$ が求めたい事後分布 $P(\theta \| S)$ の近似値になっている」。これが変分ベイズ法の要点です。

平均場近似による変分下限の交互最適化

変分ベイズ法では、変分下限をパラメータ変数ごとの交互最適化によって局所最大値へ近づける戦略を取ります。すなわち、物理の磁気スピンモデル等でもよく用いられる平均場近似（Mean field approximation）を行います。

平均場近似は、パラメータ変数 $\theta=\{ \theta_{1},\theta_{2},\theta_{3} \}$ （ここでは3つとします）同士で確率分布が独立であるという仮定に基づいています。つまり、

$Q(\theta)=Q(\theta_{1})Q(\theta_{2})Q(\theta_{3})$

のような直積に解を限定して変数分離するわけです。この形を仮定すると、前節でご説明した変分下限は以下のように整理されます。

$\displaystyle L(Q(\theta))=\sum_{\theta} Q(\theta_{1})Q(\theta_{2})Q(\theta_{3})\log{\frac{P(S, \theta)}{Q(\theta_{1})Q(\theta_{2})Q(\theta_{3})}}$

このように変数分離ができると、変分下限はパラメータ変数ごとに偏微分が計算可能になります。それぞれの偏微分を0とおけば変分下限を最大化するための方程式が以下のとおり得られます。

$\displaystyle Q(\theta_{1}) \propto exp\biggl\{ \sum_{\theta_{2}} \sum_{\theta_{3}} Q(\theta_{2}) Q(\theta_{3}) \log{P(S, \theta)} \biggr\}$

$\displaystyle Q(\theta_{2}) \propto exp\biggl\{ \sum_{\theta_{1}} \sum_{\theta_{3}} Q(\theta_{1}) Q(\theta_{3}) \log{P(S, \theta)} \biggr\}$

$\displaystyle Q(\theta_{3}) \propto exp\biggl\{ \sum_{\theta_{1}} \sum_{\theta_{2}} Q(\theta_{1}) Q(\theta_{2}) \log{P(S, \theta)} \biggr\}$

この式は、求めたいパラメータ変数以外の変数に関して、同時確率の対数の期待値を計算するような形になっています。

GaP-NMFの定式化

変分ベイズ法で近似関数 $Q(\theta)$ を得るための方程式が導けたので、NMFを無限基底に拡張したガンマ過程NMF（Gamma process-NMF）、通称GaP-NMFを定式化できます。

GaP-NMFは別名ノンパラメトリックベイズモデルとも呼ばれ、2010年頃にM. D. Hoffmanらが発表した論文「Bayesian nonparametric matrix factorization for recorded music」の中で初めて提案されました。

まずは、最初の節でご紹介した、乖離度Dとして板倉・斎藤擬距離を採用した場合に対応する確率分布関数、

を参考にして、以下のような生成モデルを定義します。

$\displaystyle P(S \| w, H, U) \propto \prod_{k=1}^{K}\prod_{n=1}^{N} \frac{1}{\sum_{m=1}^{M} w_{m}H_{k,m}U_{m,n}} e^{-\frac{S_{k,n}}{\sum_{m=1}^{M} w_{m}H_{k,m}U_{m,n}}}$

元のモデルと異なる点は、各基底番号 $m$ に対応した重み係数 $w_{m}$ が導入されている点です。

基底数の上限 $M$ は依然として事前に決めておくべき固定値ですが、これを十分に大きくとることによって実質的に無限な基底数を想定できることになります。

ここで、重み $w_{m}$ としては、限られた数の $m$ でのみ値を持ち、それ以外の $m$ ではほぼ0になるような、スカスカ（Sparce）な分布が望ましいです。

なぜなら、データを近似するために理想的な基底の候補がいくつもあるのは不自然ですし、実効的な基底数が減るほど計算が軽くなるからです（重みがほぼ0の基底番号 $m$ は計算上無視できる）。

このような性質を満たす分布の候補の一つとして、ガンマ分布（Gamma distribution）があります。ここでは $w_{m}$ の確率分布関数 $P(w_{m})$ として以下のガンマ分布を採用します。

$\displaystyle P(w_{m})=Gamma \biggl( w_{m} \| \frac{\alpha}{M}, \alpha c \biggr)=\frac{(\alpha c)^{\frac{\alpha}{M}}}{\Gamma(\frac{\alpha}{M})} w_{m}^{(\frac{\alpha}{M}-1)} e^{-\alpha c w_{m}}$

ガンマ分布の第一引数 $\frac{\alpha}{M}$ は形状母数（Shape parameter）と呼ばれ、小さくすると急峻な分布になります。他方、第二引数 $\alpha c$ は逆尺度母数（Inverse scale parameter）と呼ばれ、小さくすると分布の幅が広がります。

基底数の上限 $M$ を $\alpha$ よりも十分大きい値にセットしておけば、形状母数が小さくなって0に向かって急峻な分布が得られるため、重み係数 $w_{m}$ は多くの $m$ で0になることが期待されます。

同様に、NMFで分解される基底行列Hと活性化行列Uについてもスカスカな分布が理想的なので、HとUの事前分布にも同じく適当なガンマ分布を採用します（単純化のため形状母数と逆尺度母数は同じ値a,bを採用）。

$\displaystyle P(H_{k,m})=Gamma(H_{k,m} \| a, a)=\frac{a^{a}}{\Gamma(a)} H_{k,m}^{(a-1)} e^{-aH_{k,m}}$

$\displaystyle P(U_{m,n})=Gamma(U_{m,n} \| b, b)=\frac{b^{b}}{\Gamma(b)} U_{m,n}^{(b-1)} e^{-bU_{m,n}}$

以上のモデル定義により、同時確率の対数 $\log{P(S,\theta)}=\log{P(S,w,H,U)}=\log{P(S \| w,H,U)}+\log{P(w)}+\log{P(H)}+\log{P(U)}$ を解析的に書き下せるので、前節で導出した変分ベイズ法の変分下限最大化方程式に当てはめることができます。

ただし、生成モデル $\log{P(S \| w,H,U)}$ の部分はそのままだと計算の都合が悪いので、原著論文では $\displaystyle -\frac{1}{x}$ に対するJensenの不等式や $-\log{x}$ に対する1次のTeylor展開を利用して更なる下限を定義しています。

計算の詳細は原著論文に譲るとして、結果的に事後分布はパラメータ３つで特徴付けられる、一般化逆ガウス分布（GIG: Generalized inverse Gaussian distribution）という形になるそうです。

$\displaystyle Q(H_{k,m})=GIG\biggl(a, \: a+E_{Q}[w_{m}]\sum_{n} \frac{E_{Q}[U_{m,n}]}{\omega_{k,n}}, \: E_{Q}[w_{m}^{-1}]\sum_{n} S_{k,n}\phi_{k,m,n}^{2}E_{Q}[U_{m,n}^{-1}] \biggr)$

$\displaystyle Q(U_{m,n})=GIG\biggl(b, \: b+E_{Q}[w_{m}]\sum_{k} \frac{E_{Q}[H_{k,m}]}{\omega_{k,n}}, \: E_{Q}[w_{m}^{-1}]\sum_{k} S_{k,n}\phi_{k,m,n}^{2}E_{Q}[H_{k,m}^{-1}] \biggr)$

$\displaystyle Q(w_{m})=GIG\biggl(\frac{\alpha}{L}, \: \alpha c+\sum_{k,n} \frac{E_{Q}[H_{k,m}U_{m,n}]}{\omega_{k,n}}, \: \sum_{k,n} S_{k,n}\phi_{k,m,n}^{2}E_{Q}[H_{k,m}^{-1}U_{m,n}^{-1}] \biggr)$

ここに出てくる $E_{Q}[x]$ という表記は期待値計算を表します。つまり、 $x$ に確率分布関数 $Q$ を掛けて変数について和を取る操作です。GIGの具体的な形や期待値の公式は原著論文に記載されているのでご参照ください。

また、 $\omega$ と $\phi$ はTeylor展開とJensenの不等式で下限を設定するために必要な補助変数です。これらについても、以下の等号成立条件を満たすように更新すると効率よく収束します。

$\displaystyle \omega_{k,n}=\sum_{m} E_{Q}[w_{m}H_{k,m}U_{m,n}]$

$\displaystyle \left. \phi_{k,m,n}=\biggl( E_{Q}[w_{m}^{-1}H_{k,m}^{-1}U_{m,n}^{-1}] \biggr)^{-1} \middle/ \sum_{m} \biggl( E_{Q}[w_{m}^{-1}H_{k,m}^{-1}U_{m,n}^{-1}] \biggr)^{-1} \right.$

Pythonで実装

早速Pythonで実装しました。最近、コードはGithubに上げた方が管理しやすいと（今さら）気が付いたので、ブログに直接書き込むことはしません。ご興味がおありの方は以下のURLをご参照ください。

github.com

これまでのCOVID-19自主隔離期間のうち、1/3くらいをこの実装に捧げました。ベイズ推定系のアルゴリズムを組むのは初めてで勝手が分からなかったのと、変分下限が突如発散するバグが何度も出たので、途中泣きそうになりました。

調べてみたところ、どうやらこのアルゴリズムはGIGの三番目の引数の値が小さくなると計算が不安定になって発散しやすいらしく、そこでハマっていたようです。まさに理想（理論）と現実（実装）の差です。例外処理を一行入れただけで大人しく収束しました。

今回はおもしろいデモ音声とかはないのですが、以前使ったピアノ音のサンプルをGaP-NMFで分解してみたところ、基底数は9くらいに落ち着きました。収束の様子は以下のとおりです。

f:id:yuki0718:20200623000444p:plain

上述したとおり、ベイズ変分法は変分下限を最大化するアプローチなので、収束の目安は変分下限です。スペクトログラムのサイズが大きいと変分下限はかなり大きい値になるので、収束条件には変化率を採用しました。

f:id:yuki0718:20200623000457p:plain

また、GaP-NMFのユニークな点は、「get_valid_M」という関数で重み $w_{m}$ の値が有効な（相対的に大きい）基底番号 $m$ を特定し、それ以外は計算から落としてしまうところにあります。

そのため、プログラム開始直後（初期基底数100）は一連のパラメータ更新を終えるまで時間がかかるのですが、収束してくると基底数が最適な数に向かって減ってくる（行列サイズが小さくなる）ので、更新速度は急激に速くなります。

NMFを実装したことがある人なら、基底数 $M$ というパラメータがパフォーマンスに与える影響の大きさを理解できると思います。決め打ちで試すしかなかった基底数自体を最適化できるこのアルゴリズムは画期的です。

また、スカスカな分布としてはガンマ分布以外にもベータ分布やディリクレ分布などを使ったタイプも提案されているみたいなので、そのうち暇があったら実装してみたいと思います。

今日のところは以上です。やってみて感じたのは、ベイズ推定系のアルゴリズムは求められる数学の知識が比較的高度なので*6、そもそも論文を読み解くところからしてキツいということです。

ただ、機械学習を1年間やった経験からすると、ディープラーニングに対抗できそうなモデルベースの手法は、ベイズ推定のような確率論的なアプローチくらいなんじゃないでしょうか。今後は（たぶんもうすでに）両者の融合領域が活発になるはずです。

そういう意味では、今回の実装は正直けっこう辛かったですけどやって良かったですし、今後もこの手のアルゴリズムを勉強していきたいなあと思いました。

稚文をお読みいただきありがとうございました。

*1:なお、今のところはこの関数の形を採用することに特段根拠があるわけではなく、「とりあえず解析的に知られている分布を使っていくつか生成モデルを定義してみました」という程度の話です。

*2:他の例としては、ギブスサンプリングに代表されるMCMC（マルコフ連鎖モンテカルロ法）などが用いられます。MCMCはベイズ変分法と違ってグローバルな最適解を得られることが保証されています。ただ、MCMCは論文を読むとけっこう「うげっ」となる数式が並んでいるので、まだまともに勉強できていません。いつかは取り組まないとダメだと思ってますけど。

*3:共役事前分布（Conjugate prior distribution）と呼ばれる、生成モデルに掛け算しても関数のタイプが変わらないものがよく用いられます。

*4:簡単に言うと、「確率の積分後に対数を取った値は、確率の対数を取ってから積分した値よりも必ず大きい」という内容の不等式です。

*5:カルバック・ライブラー情報量（Kullback–Leibler divergence）と呼ばれることもあります。これは乖離度Dに出てくる一般化KLダイバージェンスとは異なるのでご注意ください。

*6:高度と言っても、基礎を丁寧に押さえればやってることは（たぶん）シンプルです。要は何事も慣れなんでしょうね。