この記事では，研究のサーベイをまとめていきたいと思います。ただし，全ての論文が網羅されている訳ではありません。また，分かりやすいように多少意訳した部分もあります。ですので，参考程度におさめていただければ幸いです。間違えている箇所がございましたらご指摘ください。随時更新予定です。
読みたい場所へジャンプ！
NMFの問題設定
詳しくはこちらの記事をご覧ください。
NMFの問題設定は、ある観測行列Yを２つの非負値行列の積に分解するということです。
\begin{eqnarray}
\boldsymbol{Y} &\simeq& \boldsymbol{H}\boldsymbol{U} (\equiv \boldsymbol{F})
\end{eqnarray}
Y：観測値
H：基底行列
U：係数行列
F：モデル化行列
それぞれの行列のサイズは以下の通りとし、要素は実数を取るものとします。
\begin{eqnarray}
\boldsymbol{Y} &\in& \mathbb{R} (K \times N)\\
\boldsymbol{F} &\in& \mathbb{R} (K \times N)\\
\boldsymbol{H} &\in& \mathbb{R} (K \times M)\\
\boldsymbol{U} &\in& \mathbb{R} (M \times N)
\end{eqnarray}
これをベクトル表記で書くと、以下のようになります。
\boldsymbol{ y_n } &\simeq& \displaystyle \sum_{ m = 1 }^{ M } \boldsymbol{ h_m }u_{m,n} (n = 1, 2, \ldots, N)\\
\end{eqnarray}
\begin{eqnarray}
\boldsymbol{Y} &=& [\boldsymbol{ y_1 }, \ldots, \boldsymbol{ y_N }] = (y_{k,n})_{K \times N}\\
\boldsymbol{H} &=& [\boldsymbol{ h_1 }, \ldots, \boldsymbol{ h_M }] = (h_{k,m})_{K \times M} \\
\boldsymbol{U} &=& [\boldsymbol{ u_1 }, \ldots, \boldsymbol{ u_N }] = (h_{m,n})_{M \times N}\\
\end{eqnarray}
我々が目指すのは、いかに$\boldsymbol{HU}$が$\boldsymbol{Y}$と「近く」なるかです。そこで、数学的な距離として「乖離度」を導入します。
「数学的な距離」
●ユークリッド距離（二乗距離）
●Iダイバージェンス
●板倉斎藤擬距離
●βダイバージェンス
上から順番に、定義は以下の通りです。
\mathcal{D}_{EU}(y|x) &=& (y – x)^2\\
\mathcal{D}_{KL}(y|x) &=& y\log\frac{y}{x} – y + x\\
\mathcal{D}_{IS}(y|x) &=& \frac{y}{x} – \log\frac{y}{x} – 1\\
\mathcal{D}_{β}(y|x) &=& \frac{y^{\beta}}{\beta(\beta – 1)} + \frac{x^\beta}{\beta} – \frac{yx^{\beta – 1}}{\beta – 1}
\end{eqnarray}
上記乖離度を目的関数として、補助関数法を利用して乗法更新式を導出したものが、NMFのアルゴリズムになります。この記事では、それぞれの乖離度が最尤推定問題として捉えられることを示していきたいと思います。
生成プロセスを考慮した統計的アプローチ
以下では
\begin{eqnarray}
\sum_{m} h_{k,m}u_{m,n} = x_{k,n}
\end{eqnarray}
なる$x_{k,n}$を使って説明していきます。結論からお伝えすると、以下のようになります。
①ユークリッド距離基準
②Iダイバージェンス基準
③ISダイバージェンス基準
④βダイバージェンス基準
のNMFはそれぞれ、観測行列$\boldsymbol{Y}$の各行列要素$y_{k,n}$が$x_{k,n}$を平均とした
①正規分布
②ポアソン分布
③指数分布
④tweedie分布
に従って独立に生成されたと仮定した場合の$\boldsymbol{H}$と$\boldsymbol{U}$の最尤推定問題と等価になります。以下、その証明をお伝えしていきます。
NMF⇔最尤推定の証明
まずは、NMF側から見ていきましょう。それぞれの乖離度を改めて確認します。
\mathcal{D}_{EU}(y|x) &=& (y – x)^2\\
\mathcal{D}_{KL}(y|x) &=& y\log\frac{y}{x} – y + x\\
\mathcal{D}_{IS}(y|x) &=& \frac{y}{x} – \log\frac{y}{x} – 1\\
\mathcal{D}_{β}(y|x) &=& \frac{y^{\beta}}{\beta(\beta – 1)} + \frac{x^\beta}{\beta} – \frac{yx^{\beta – 1}}{\beta – 1}
\end{eqnarray}
以下では、各分布から導出される対数尤度差が上記乖離度と等価になることを示していきます。そのために、各分布の確率密度関数を確認しておきましょう。
\mathrm{Normal}(z| \mu, \sigma^2) &=& \frac{1}{\sqrt{2\pi} \sigma} e^{- \frac{(z – \mu)^2}{2 \sigma^2}}\\
\mathrm{Poisson}(z| \mu ) &=& \frac{\mu^z e^{- \mu}}{z!} (z \geq 0)\\
\mathrm{Exponential}(z| \mu) &=& \frac{1}{\mu} e^{- \frac{z}{\mu}} (z \geq 0)\\
\mathrm{Tweedie}(z| \phi) &=& a(z, \phi) e^{\frac{1}{\phi} (z \rho(\mu) – \kappa(\mu))}
\end{eqnarray}
ただし、$\rho(\mu)$と$\kappa(\mu)$は以下を満たす関数とします。
\begin{eqnarray}
\rho(\mu) &=&
\begin{cases}
\frac{\mu^{\beta – 1} – 1}{\beta – 1} & ( \beta \neq 1 )\\
\log \mu & ( \beta = 1 )
\end{cases}\\
\kappa(\mu) &=&
\begin{cases}
\frac{\mu^{\beta} – 1}{\beta} & ( \beta \neq 0 )\\
\log \mu & ( \beta = 0 )
\end{cases}
\end{eqnarray}
βダイバージェンスはβ≠0,1で定義される関数でしたが、β→0,1でそれぞれISダイバージェンスとIダイバージェンスに収束する性質を利用するためにも、式（23）（24）のようにβ=0,1の値を定めてしまいます。
Tweedie分布は、正規分布・ポアソン分布・指数分布の一般化です。$\phi$が分散パラメータ、$\beta$がTweedie分布を決定づけるパラメータになります。また、$a(z, \phi)$は$\beta$によって異なりるので閉形式で記述することができません。
さて、それぞれについて対数尤度$\mathrm{L}$を最大にする$x_{k,n}$を計算していきたいと思います。また、簡単のためインデックスの$k$と$n$は省略します。偏微分して最大値の候補が実際に最大値を取るかは、$x$と$y$に非負値の条件があることに注意して考えていきましょう。
正規分布
対数尤度は
\begin{eqnarray}
\mathrm{L}_{N} &=& \log p(y|x)\\
&=& \log \frac{1}{\sqrt{2\pi} \sigma} e^{- \frac{(y – x)^2}{2 \sigma^2}}\\
&=& – \frac{(y – x)^2}{2 \sigma^2} – \frac{1}{2} \log 2\pi – \log \sigma
\end{eqnarray}
となります。この式を$x$で偏微分して、対数尤度を最大にする$x$を求めます。
\begin{eqnarray}
\frac{ \partial \mathrm{L}_{N} }{ \partial x }
&=& \frac{1}{\sigma^2} (y – x)
\end{eqnarray}
したがって、対数尤度$\mathrm{L}_{N}$は$x = y$のときに最大になります。
ポアソン分布
対数尤度は
\begin{eqnarray}
\mathrm{L}_{P} &=& \log p(y|x)\\
&=& \log \frac{x^y e^{- x}}{y!}\\
&=& y \log x – x – \log y!
\end{eqnarray}
となります。この式を$x$で偏微分して、対数尤度を最大にする$x$を求めます。
\begin{eqnarray}
\frac{ \partial \mathrm{L}_{P} }{ \partial x }
&=& \frac{y}{x} – 1 \\
&=& \frac{y – x}{x}
\end{eqnarray}
したがって、対数尤度$\mathrm{L}_{P}$は$x = y$のときに最大になります。
指数分布
対数尤度は
\begin{eqnarray}
\mathrm{L}_{E} &=& \log p(y|x)\\
&=& \log \frac{1}{x} e^{- \frac{y}{x}} \\
&=& – \frac{y}{x} – \log x \\
\end{eqnarray}
となります。この式を$x$で偏微分して、対数尤度を最大にする$x$を求めます。
\begin{eqnarray}
\frac{ \partial \mathrm{L}_{E} }{ \partial x }
&=& \frac{y}{x^2} – \frac{1}{x}\\
&=& \frac{y – x}{x^2}
\end{eqnarray}
したがって、対数尤度$\mathrm{L}_{E}$は$x = y$のときに最大になります。
Tweedie分布
対数尤度は
\mathrm{L}_{T} &=& \log p(y|x)\\
&=& \log a(y, \phi) e^{\frac{1}{\phi} (y \rho(x) – \kappa(x))}\\
&=& \log a(y, \phi) + \frac{1}{\phi}(y \rho(x) – \kappa(x))
\end{eqnarray}
となります。この式を$x$で偏微分して、対数尤度を最大にする$x$を求めます。
\begin{eqnarray}
\frac{ \partial \mathrm{L}_{T} }{ \partial x }
&=& \frac{1}{\phi}(y \rho'(x) – \kappa'(x))\\
\end{eqnarray}
βの値によって確率密度関数が変わるので、場合分けしていきましょう。
★β=0のとき
\begin{eqnarray}
\frac{ \partial \mathrm{L}_{T} }{ \partial x }
&=& \frac{1}{\phi}(yx^{- 2} – \frac{1}{x})\\
&=& \frac{1}{\phi x}\frac{y – x}{x}
\end{eqnarray}
したがって、対数尤度$\mathrm{L}_{T}$は$x = y$のときに最大になります。
ここでもβ=0のとき、指数分布に基づく尤度に等しくなることが分かります。
★β=1のとき
\begin{eqnarray}
\frac{ \partial \mathrm{L}_{T} }{ \partial x }
&=& \frac{1}{\phi}(yx^{- 1} – 1)\\
&=& \frac{1}{\phi}(\frac{y – x}{x})
\end{eqnarray}
したがって、対数尤度$\mathrm{L}_{T}$は$x = y$のときに最大になります。
ここでもβ=1のとき、ポアソン分布に基づく尤度に等しくなることが分かります。
★β≠0,1のとき
\begin{eqnarray}
\frac{ \partial \mathrm{L}_{T} }{ \partial x }
&=& \frac{1}{\phi}(yx^{\beta – 2} – x^{\beta – 1})\\
&=& \frac{1}{\phi}x^{\beta – 2}(y – x)
\end{eqnarray}
したがって、対数尤度$\mathrm{L}_{T}$は$x = y$のときに最大になります。
以上より、βの値に関わらず対数尤度$\mathrm{L}_{T}$は$x = y$のときに最大になることが分かりました。
対数尤度と乖離度の関係
さて、上までの流れをまとめてみます。$y_{k,n}$が$x_{k,n}$を平均とした各分布より独立に生成されたと仮定すると、各分布の確率密度関数の対数を取った対数尤度は分布の種類にかかわらず$y_{k,n} = x_{k,n}$のときに最大になりました。これを定式化してみます。
\begin{eqnarray}
&&\mathrm{L}(y) \geq \mathrm{L}(x)\\
&\Leftrightarrow& \mathrm{L}(y) – \mathrm{L}(x) \geq 0
\end{eqnarray}
上のように表すことができます。ここで、$\mathrm{L}(y) – \mathrm{L}(x)$に注目してみましょう。$\mathrm{L}(y) – \mathrm{L}(x)$は「$x = y$のときに最小となり」「常に非負値の」$x$と$y$の近さを表す尺度であると解釈できます。
そしてなんと、実際に各分布に対して$\mathrm{L}(y) – \mathrm{L}(x)$を計算すると、乖離度の尺度と一致するのです。この事実が、NMFが観測値の生成プロセスを統計的に仮定した最尤推定値と等価であることを裏付けているのです。以下、計算をしていきます。
正規分布
\mathrm{L}(y) – \mathrm{L}(x)
&=& (- \frac{(y – y)^2}{2\sigma^2} – \frac{1}{2} \log 2\pi – \log \sigma) – (- \frac{(y – x)^2}{2\sigma^2} – \frac{1}{2} \log 2\pi – \log \sigma)\\
&=& \frac{(y – x)^2}{2\sigma^2}\\
&\propto& (y – x)^2 \\
&=& \mathcal{D}_{EU}(y|x) \\
\end{eqnarray}
見事、ユークリッド距離基準の乖離度と一致しました。
ポアソン分布
\mathrm{L}(y) – \mathrm{L}(x)
&=& (y \log y – y – \log y!) – (y \log x – x – \log y!)\\
&=& y \log \frac{y}{x} – y + x\\
&=& \mathcal{D}_{I}(y|x)
\end{eqnarray}
見事、Iダイバージェンス基準の乖離度と一致しました。
指数分布
\mathrm{L}(y) – \mathrm{L}(x)
&=& (- 1 – \log y) – (- \frac{y}{x} – \log x)\\
&=& \frac{y}{x} – \log \frac{y}{x} – 1\\
&=& \mathcal{D}_{IS}(y|x)
\end{eqnarray}
見事、ISダイバージェンス基準の乖離度と一致しました。
Tweedie分布
\mathrm{L}(y) – \mathrm{L}(x)
&=& (\log a(y, \phi) + \frac{1}{\phi} (y\rho(y) – \kappa(y))) – (\log a(y, \phi) + \frac{1}{\phi} (y\rho(x) – \kappa(x)))\\
&=& \frac{1}{\phi} (y(\rho(y) – \rho(x)) – (\kappa(y) – \kappa(x)))\\
\end{eqnarray}
★β=0のとき
\frac{1}{\phi} (y(\rho(y) – \rho(x)) – (\kappa(y) – \kappa(x)))
&=& \frac{1}{\phi} (y(\frac{y^{-1} – 1}{-1}) – (\log y – \log x))\\
&\propto& \frac{y}{x} – \log{y}{x} – 1 \\
&=& \mathcal{D}_{IS}(y|x)
\end{eqnarray}
★β=1のとき
\frac{1}{\phi} (y(\rho(y) – \rho(x)) – (\kappa(y) – \kappa(x)))
&=& \frac{1}{\phi} (y(\log y – \log x) – ((y – 1) – (x – 1)))\\
&\propto& y \log \frac{y}{x} – y + x\\
&=& \mathcal{D}_{I}(y|x)
\end{eqnarray}
★β≠0,1のとき
\frac{1}{\phi} (y(\rho(y) – \rho(x)) – (\kappa(y) – \kappa(x)))
&=& \frac{1}{\phi} (y \frac{(y^{\beta – 1}- 1) – (x_{\beta – 1} – 1)}{\beta – 1}) – (\frac{(y^{\beta}- 1) – (x_{\beta} – 1)}{\beta})\\
&\propto& \frac{y^{\beta}}{\beta(\beta – 1)} + \frac{x^{\beta}}{\beta} – \frac{yx^{\beta – 1}}{\beta – 1}\\
&=& \mathcal{D}_{\beta}(y|x)
\end{eqnarray}
見事、βダイバージェンス基準の乖離度と一致しました。
注意事項
今回の生成モデルの仮定は、「非負の観測値」に対して行っていました。「複素データの観測値」に対して分散変動ガウス分布を仮定すると、ISダイバージェンス基準のNMFに帰着することが示されています。
●「Nonnegative Matrix Factorization with the Itakura-Saito
Divergence: With Application to Music Analysis」（ F´evotte et al.,2009,Neural Computation）
さらなる拡張
βダイバージェンスの拡張として、Bregmanダイバージェンスが使われています。Bregmanダイバージェンスは指数分布族と対応することが示されています。詳しくは以下の論文をご覧ください。
●「非負値行列因子分解とその音響信号処理への応用」（亀岡,2015,日本統計学会誌）
●「非負値行列因子分解」（亀岡,2012,NTTコミュニケーション科学基礎研究所 解説：特集 計測・センシングのアルゴリズム）
●「零過剰 Tweedie 分布に基づく非負値行列因子分解について」（阿部/宿久,2015,日本計算機統計学会大会論文集）
●「非負値行列因子分解とその音響信号処理への応用」（亀岡,2015,日本統計学会誌）
●「Nonnegative Matrix Factorization with the Itakura-Saito
Divergence: With Application to Music Analysis」（ F´evotte et al.,2009,Neural Computation）
