欠損値のある行列の非負値行列因子分解

非負値行列因子分解（NMF）はQiitaでもいくつか紹介されていますが、
欠損値のある場合についての情報が見当たらなかったので、紹介します。

参考文献

Lee, D. D., & Seung, H. S. (n.d.). Algorithms for Non-negative Matrix Factorization.

問題設定

NMFでは、二次元の計測データYについて以下のコスト関数を最小化する非負の行列 W および H を求めます。

• 二乗距離
  $$
  \left|| \mathrm{Y} - \mathrm{WH} \right||^2
  = \frac{1}{2}\sum_{i j} \left\{
  y_{i j} - \sum_k w_{i k} h_{k j}
  \right\}^2
  \tag{1}
  $$

• 一般化KL距離
  $$
  \mathrm{KL}\left[\mathrm{Y} || \mathrm{WH} \right]
  = \sum_{i, j} \left\{
  y_{i j} \log\frac{y_{i j}}{\sum_k w_{i k} h_{k j}}
  - y_{i j} + \sum_k w_{i k}h_{k j}
  \right\}
  \tag{2}
  $$

欠損値があるデータYの場合を考えます。
上式の総和の部分を欠損でない要素の組
$\mathcal{I}=\{(i_0, j_0), (i_1, j_1)... \}$
についての総和に置き換えます。

• 二乗距離
  $$
  \begin{aligned}
  \left|| \mathrm{Y} - \mathrm{WH} \right||^2
  &= \frac{1}{2}\sum_{i, j \in \mathcal{I}} \left\{
  y_{i j} - \sum_k w_{i k} h_{k j}
  \right\}^2 \\
  &= \frac{1}{2}\sum_{i, j}\left\{
  \left(y_{i j} - \sum_k w_{i k} h_{k j}\right)m_{i,j}
  \right\}^2
  \end{aligned}
  \tag{3}
  $$

• 一般化KL距離
  $$
  \begin{aligned}
  \mathrm{KL}\left[\mathrm{Y} || \mathrm{WH} \right]
  &= \sum_{i, j\in \mathcal{I}} \left\{
  y_{i j} \log\frac{y_{i j}}{\sum_k w_{i k} h_{k j}}
  - y_{i j} + \sum_k w_{i k}h_{k j}
  \right\}\\
  &= \sum_{i, j} \left\{
  \left(y_{i j} \log\frac{y_{i j}}{\sum_k w_{i k} h_{k j}}
  - y_{i j} + \sum_k w_{i k}h_{k j}\right)m_{i j}
  \right\}
  \end{aligned}
  \tag{4}
  $$

ここでそれぞれの二行目の式は、有効な要素に関する総和をマスク行列
$\mathrm{M}=\{m_{i j}\}$ で置き換えたものです。
マスク行列M は0と1からなる行列で、Yの欠損値に対応する要素には0、それ以外には1 が入っています。

更新式

二乗距離、一般化KL距離を採用した時のそれぞれの更新式は以下のようになります。
なお、これらの導出については後ほどで述べます。

二乗距離（欠損値なし）

$$
w_{i, k} \gets w_{i k}
\frac{
\sum_j y_{i j} h_{k j}
}
{
\sum_{j k'} w_{i k'} h_{k' j} h_{k j}
}
$$

$$
h_{k j} \gets h_{k j}
\frac{
\sum_i w_{i k} y_{i j}
}
{
\sum_{i k'} w_{i k} w_{i k'} h_{k' j}
}
$$

Python では、行列積と要素演算を使うと非常に簡単に記述できます。

def update_squared_loss(Y, W, H):
  """  Update W and H with squared loss  """
  W = W * np.dot(Y, H.T) / np.dot(np.dot(W, H), H.T)
  H = H * np.dot(W.T, Y) / np.dot(W.T, np.dot(W, H))
  return W, H

二乗距離（欠損値あり）

$$
w_{i k} \gets w_{i k}
\frac{
\sum_j y_{i j} m_{i j} h_{k j}
}
{
\sum_{j k'} w_{i k'} h_{k' j} h_{k j} m_{i j}
}
$$

$$
h_{k j} \gets h_{k j}
\frac{
\sum_i w_{i k} y_{i j} m_{i j}
}
{
\sum_{i k'} w_{i k} w_{i k'} h_{k' j} m_{i j}
}
$$

Python 擬似コード

def update_squared_loss_w_mask(Y, W, H, M):
  """  Update W and H with squared loss with missing values """
  W = W * np.dot(Y * M, H.T) / np.dot(np.dot(W, H) * M, H.T)
  H = H * np.dot(W.T, Y * M) / np.dot(W.T, np.dot(W, H) * M)
  return W, H

一般化KL距離（欠損値なし）

$$
w_{i k} \gets w_{i k}
\frac{
\sum_j \frac{y_{i j}}{\sum_{k'}w_{i k'} h_{k' j}} h_{k j}
}
{
\sum_{j} h_{k j}
}
$$

$$
h_{k j} \gets h_{k j}
\frac{
\sum_i \frac{y_{i j}}{\sum_{k'}w_{i k'} h_{k' j}} w_{i k}
}
{
\sum_{i} w_{i k}
}
$$

Python 擬似コードは以下のとおりです。

def update_kl_loss(Y, W, H):
  """  Update W and H with squared loss  """
  W = W * np.dot(Y / np.dot(W, H), H.T) / np.sum(H, axis=1, keepdims=True)
  H = H * np.dot(W.T, Y / np.dot(W, H)) / np.sum(W, axis=0, keepdims=True)
  return W, H

一般化KL距離（欠損値あり）

$$
w_{i k} \gets w_{i k}
\frac{
\sum_j \frac{y_{i j} m_{i j}}{\sum_{k'}w_{i k'} h_{k' j}} h_{k j}
}
{
\sum_{j} h_{k j} m_{i j}
}
$$

$$
h_{k j} \gets h_{k j}
\frac{
\sum_i \frac{y_{i j} m_{i j}}{\sum_{k'}w_{i k'} h_{k' j}} w_{i k}
}
{
\sum_{i} w_{i k} m_{i j}
}
$$

Python 擬似コードは以下のとおりです。

def update_kl_loss_w_mask(Y, W, H, M):
  """  Update W and H with squared loss  """
  W = W * np.dot(Y * M / np.dot(W, H), H.T) / np.dot(M, H.T)
  H = H * np.dot(W.T, Y * M / np.dot(W, H)) / np.sum(W.T, M)
  return W, H

更新式の導出

二乗距離（欠損値なし）

まず、コスト関数(1)を W に関して最小化することにします。
Wに関するコスト関数を F(W) とおくと、

$$
\begin{aligned}
\mathrm{F(W)} &= \frac{1}{2}\sum_{i j} \left\{
y_{i j} - \sum_k w_{i k} h_{k j}
\right\}^2 \\
&= \frac{1}{2}\sum_{i j} \left\{
y_{i j}^2 - 2 y_{i j} \sum_k w_{i k} h_{k j}
+ \left( \sum_k w_{i k} h_{k j} \right)^2
\right\}
\end{aligned}
\tag{5}
$$
のように書けます。

方針としては、F(W) を上から抑える単純な関数Gを求め、その関数について W を最適化します。
具体的には、以下の性質を持つ関数 G(W; W') を探します。

$$
\mathrm{
F(W) = G(W; W), ;;; F(W) \leq G(W; W')
}
\tag{6}
$$

上から抑える関数のイメージについては、
チュートリアル：非負値行列因子分解
などを参考にしてください。

ここで、Jensenの不等式を思い出します。
Jensenの不等式では 下に凸の関数 $f(x)$ に対して、 $\sum_k r_k = 1$ となる $r_k$ を用いて

$$
f\left(\sum_k r_k x_k\right)
\leq
\sum_k r_k f(x_k)
$$

の関係が成り立ちます。
これを式5の第3項に適用すると $x^2$ が下に凸であることから、

$$
\begin{aligned}
\left( \sum_k w_{i k} h_{k j} \right)^2
&= \left( \sum_k r_{i j k} \frac{w_{i k} h_{k j}}{r_{i j k}} \right)^2 \\
&\leq \sum_k r_{i j k} \left(\frac{w_{i k} h_{k j}}{r_{i j k}} \right)^2
\end{aligned}
$$

を示せます。
なお $r_{i j k}$ は $\sum_k r_{i j k} = 1$ を満たす定数です。
天下り的ですが、
$$
r_{i j k} = \frac{w_{i k'} h_{k j}}{\sum_{k'} w_{i k'} h_{k' j}}
$$
としたときに、以下が式6を満たすことを示します。
$$
\begin{aligned}
\mathrm{G(W; W')}=
\frac{1}{2}\sum_{i, j} \left\{
y_{i j}^2 - 2 y_{i j} \sum_k w_{i k} h_{k j}

• \sum_k r_{i j k} \left(\frac{w_{i k} h_{k j}}{r_{i j k}}\right)^2
  \right\}
  \end{aligned}
  \tag{7}
  $$

まず、 $\mathrm{F(W) \leq G(W; W')}$ は上に述べた通り、Jensenの不等式から示すことができます。
また $\mathrm{F(W) = G(W; W)}$ も第3項を約分すると示せます。
そのため式7を最小化するWを用いると、$\mathrm{F(W)}$ は必ず $\mathrm{F(W')}$ より小さくなることがわかります。

式7を最小化するWは、勾配がゼロの点から求めます。

$$
\begin{aligned}
\frac{\partial \mathrm{G(W; W')}}{w_{i k}} = 0
;;\to;;
w_{i k} &= \frac{\sum_j h_{k j}}{\sum_{k, j} \frac{h_{k j}^2}{r_{i j k}}}
= w'_{i k} \frac{\sum_j y_{i j} h_{k j}}
{\sum_{j, k'} w'_{i k} h_{k' j} h_{k j}}
\end{aligned}
$$
となって、上で述べた更新式が示せます。

また、式の対称性を考えるとHに関する更新式も同様に示せます。

二乗距離（欠損値あり）

欠損値ありの場合も、マスク行列を用いることで同様に更新式を導くことができます。

$$
\begin{aligned}
\mathrm{G(W; W')} =
\frac{1}{2}\sum_{i, j} \left\{
y_{i j}^2 m_{i j}^2 - 2 y_{i j} \sum_k w_{i k} h_{k j} m_{i j}

• \sum_k {r_{i j k}}\left(
  \frac{w_{i k} h_{k j}}{r_{i j k}}
  \right)^2 m_{i j}
  \right\}
  \end{aligned}
  $$

上式が式6の関係を満たすため、勾配がゼロになる点を探すと

$$
\begin{aligned}
\frac{\partial \mathrm{G(W; W')}}{w_{i k}} = 0
;;\to;;
w_{i k} &= \frac{\sum_j h_{k j} m_{i j}}
{\sum_{k j} \frac{h_{k j}^2}{r_{i j k}} m_{i j}}
= w'_{i k} \frac{\sum_j y_{i j} m_{i j} h_{k j}}
{\sum_{j k'} w'_{i, k} h_{k' j} h_{k j} m_{i j}}
\end{aligned}
$$

となり、上の更新式を示せます。

一般化KL距離（欠損値なし）

同様に、コスト関数(2)を W に関して最小化することにします。
Wに関するコスト関数を F(W) とおくと、

$$
\begin{aligned}
\mathrm{F(W)} &=
\left\{
y_{i j} \log y_{i j} - y_{i j} \log \sum_k w_{i k} h_{k j}

• \sum_k w_{i k} h_{k j}
  \right\}
  \end{aligned}
  $$

右辺括弧内第二項を上から抑える関数を探すことにします。
$-\log x$ が下に凸であることから、
$$
-\log \sum_k w_{i k} h_{k j} \leq
-\sum_k r_{i j k} \log \frac{w_{i k} h_{k j}}{r_{i j k}}
$$
を満たします。
なお $r_{i j k}$ は $\sum_k r_{i j k} = 1$ を満たす定数であり、
$$
r_{i j k} = \frac{w_{i k} h_{k j}}{\sum_{k'} w_{i k'} h_{k' j}}
$$
としたときに等号が成立します。

そのため、以下の $\mathrm{G(W; W')}$ は式6の関係をを満たします。
$$
\begin{aligned}
\mathrm{G(W; W')} &=
\left\{
y_{i j} \log y_{i j}

• y_{i j} \sum_k r_{i j k} \log \frac{w_{i k} h_{k j}}{r_{i j k}}

• \sum_k w_{i k} h_{k j}
  \right\}
  \end{aligned}
  $$

上式の勾配をゼロにする W とすることが、コスト関数を減少させる更新式を求めることになります。
$$
\begin{aligned}
\frac{\partial \mathrm{G(W; W')}}{w_{i k}} = 0
;;\to;;
w_{i k} &= \frac{\sum_j y_{i j} r_{i j k}}
{\sum_j h_{k j}}
= w'_{i k} \frac{\sum_j \frac{y_{i j}}{\sum_{k'} w'_{i k'} h_{k' j}} h_{k j}}
{\sum_j h_{k j}}
\end{aligned}
$$

一般化KL距離（欠損値あり）

これも上と同様に、マスク行列を用いることで更新式を得ることができます。

$$
\begin{aligned}
\mathrm{G(W; W')} &=
\left\{
y_{i j} m_{i j} \log y_{i j}

• y_{i j} m_{i j} \sum_k r_{i j k} \log \frac{w_{i k} h_{k j}}{r_{i j k}}

• m_{i j} \sum_k w_{i k} h_{k j}
  \right\}
  \end{aligned}
  $$

上式が式6の関係を満たすため、勾配がゼロになる点を探すと

$$
\begin{aligned}
\frac{\partial \mathrm{G(W; W')}}{w_{i k}} = 0
;;\to;;
w_{i k} &= \frac{\sum_j y_{i j} m_{i j} r_{i j k}}
{\sum_j h_{k j} m_{i j}}
= w'_{i k} \frac{\sum_j \frac{y_{i j} m_{i j}}
{\sum_{k'} w'_{i k'} h_{k' j}} h_{k j}}
{\sum_j h_{k j} m_{i j}}
\end{aligned}
$$

となり、上の更新式を示せます。

まとめ

NMFを協調フィルタリングに用いるためには、欠損値を扱える必要が有りますが、
簡単にアクセスできる情報がなかったので、ここにまとめました。

欠損値を扱うことができると、交差検証法も使えることになります。
因子の数を決める手段の1つになるのではないでしょうか。