Algorithms for Non-negative Matrix Factorization 非負矩陣分解

阿新 • • 發佈：2018-11-21

NIPS 2000 經典論文翻譯。

摘要

非負矩陣分解（NMF）是一種可以有效處理多變數資料的方法。本文介紹、分析了兩種不同的 NMF 演算法，這兩種演算法僅在更新規則（update rule）中使用的乘性因子（multiplicative factor）有所區別。其中一種可以對傳統的最小二乘誤差進行最小化（minimize），而另一種可以對廣義 Kullback-Leibler 散度（KL 散度）進行最小化。可以使用與證明最大化期望演算法收斂性類似的輔助函式來證明這兩種演算法的單調收斂性。這兩種演算法均可理解為用斜向最陡下降法（diagonally rescaled gradient descent）對因子進行最優化，以保證演算法收斂。

簡介

PCA、向量量化（Vector Quantization）等無監督學習演算法可以理解為在不同約束條件下對資料矩陣進行分解。根據其約束的不同，分解所得的因子的會表現出大相徑庭的性質。比如，PCA 僅使用了弱正交約束，從而得到非常分散的表示，對這些表示使用消去法來產生多樣性；向量量化使用一種嚴格的全域性最優型約束，最終會得到互斥的資料聚類原型。

我們之前已經證明過，在矩陣分解用於學習資料的部分表示中，非負性（non-negative）是一種非常有用的約束。學習得到的非負基向量是分散的，但仍可通過稀疏的組合，在重建時得到效果良好的表達向量。在本文中，我們詳細分析了這兩種用於在資料中學習最優的非負因子的數值演算法。

非負矩陣分解

下面我們正式開始分析如何用演算法解決以下問題：

在非負矩陣分解（NMF）中，給定非負矩陣V，找到非負矩陣因子W和H，使得：

NMF 可以應用下面的方法用於對多變數資料進行統計分析。給定一組多變數的 n 維資料向量，其向量位於一個 $n\times x$ 矩陣 V 的列中（m 表示資料集中的示例數）。然後將此矩陣近似分解為 $n\times r$ 的 W 矩陣與 $r\times m$ 的 H 矩陣。通常 r 要小於 n 或 m，以使 W 和 H 小於原始矩陣 V。最終得到的是原始資料矩陣的壓縮形態。

公式(1)中約等於的意義在於它可以將公式逐列用 $v\approx Wh$ 來表示，其中 v 和 h 是矩陣 V 和矩陣 H 的對應的列。也就是說，每個資料向量 v 近似地由矩陣 W 的各列線性組合而成，同時用 h 的分量進行加權。因此可以被認為 W 包含了對 V 中的資料的線性近似優化的基向量。由於要使用少量的基向量來表示大量的資料向量，因此只有在基向量發現數據中的潛在結構時才能實現較好的近似。

本文不會涉及關於 NMF 的應用，而會側重於在技術方面探討非負矩陣分解的技術。當然，已經有許多其它的矩陣分解方式在數值線性代數中得到了廣泛的研究，但是以前的大多數工作都不適用於非負性約束情況。

在此，我們討論了基於迭代更新 W 和 H 的兩種 NMF 演算法。由於這兩種演算法易於實現，同時能保證其收斂性，因此它們在現實情況中非常實用。其他演算法可能在總計算時間方面更有效率，但是更難實現，並且很難推廣到不同的代價函式（cost function）。因子與我們類似的演算法，已經被用於對發射斷層掃描和天文影象進行反捲積（deconvolution）。

在我們演算法的每次迭代中，會用當前值乘某些取決於公式(1)中的“近似程度”的因數，來找到 W 或 H 的新值。我們可以證明“近似程度”會隨著不斷應用這些乘法更新規則而單調減小。這正意味著更新規則的重複迭代可以保證矩陣分解演算法收斂到區域性最優。

代價函式

為了找到 $V\approx WH$ 的近似解，我們首先需要定義一個代價函式，用以量化近似的程度。可以使用兩個非負矩陣 A 和 B 的距離來構造此代價函式。一種使用的距離度量方法為：計算 A 和 B 之間的歐幾里得距離（Euclidean distance）的平方值。

(2): ||A-B||^2 = \sum_{ij}(A_{ij} - B_{ij})^2

此公式下界為 0，僅當 A=B 時距離消失。

另一種實用的度量方式為：

(3): D(A||B) = \sum_{ij}(A_{ij} \log{\frac{A_{ij}}{B_{ij}}} - A_{ij}+B_{ij})

與歐幾里得距離相同，它的下界也為 0，且在 A=B 時距離消失。但它不能被稱為“距離”，因為這個式子在 A 與 B 中並不對稱，因此我們將其稱為 A 對於 B 的“散度”（divergence）。它可以歸納為 KL 散度或者相對熵，當 $\sum_{ij}A_{ij}=\sum_{ij}B_{ij}=1$ 時，A 與 B 可以看做是標準化的概率分佈。

現在，我們可以按照以下兩種公式來將 NMF 化為最優化問題：

最優化問題1：在約束條件 $W, H \geq 0$ 下，以 W 和 H 作為引數，最小化 $||V - WH||^2$ 。

最優化問題2：在約束條件 $W, H \geq 0$ 下，以 W 和 H 作為引數，最小化 $D(V||WH)$ 。

雖然方程 $||V - WH||^2$ 和 $D(V||WH)$ 在只考慮 W 或 H 之一時為凸，但在同時考慮 WH 兩個變數時不為凸。因此，尋找一種可以找到全域性最小值的演算法去解決以上兩個最優化問題是不切實際的。但是，還有許多數值優化方法可以用於尋找區域性最小值。

雖然梯度下降法（Gradient descent）的收斂速度很慢，但它的實現最為簡單。其它方法（比如共軛梯度法）可以更快地收斂（至少在區域性最小值附近會更快），但是它們比梯度下降更復雜。此外，梯度下降方法的收斂對步長的選擇非常敏感，這對於大規模應用十分不利。

乘法更新規則

我們發現在解決上述兩個最優化問題時，在速度與實現難度中權衡，“乘法更新規則”是一種綜合性能很好方法。

定理1：歐幾里得距離 $||V-WH||$ 在下面的更新規則中呈非增：

(4): H_{a\mu} \leftarrow H_{a\mu}\frac{(W^T V)_{a\mu}}{(W^T W H)_{a\mu}}

(4): W_{ia} \leftarrow W_{ia}\frac{(V H^T)_{ia}}{(W H H^T)_{ia}}

在上述更新規則中，W 與 H 在距離公式的駐點上時，歐幾里得距離將固定不動。

定理2：散度 $D(V|WH)$ 在下面的更新規則中呈非增：

(5): H_{a\mu} \leftarrow H_{a\mu}\frac{\frac{\sum_{i}W_{ia}V_{i\mu}}{WH_{i\mu}}}{\sum_k W_{ka}}

(5): W_{ia} \leftarrow W_{ia}\frac{\frac{\sum_{\mu}H_{a\mu}V_{i\mu}}{WH_{i\mu}}}{\sum_v H_{av}}

在上述更新規則中，W 和 H 在散度公式的駐點上時，散度將不再更新。

上述定理的證明將在後面給出。我們可以發現，每次更新都是乘以一個因子。特別地，當V = WH時，可以直觀地看出這個乘數因子是一樣的，當更新規則固定時，才會得到完美的分解。

乘法與加法更新規則

可以將乘法更新與梯度下降更新進行對比。特別的，對 H 進行更新以減小平方距離可以記為：

(6):H_{a\mu} \leftarrow H_{a\mu} + \eta_{a\mu}[(W^TV)_{a\mu} - (W^T WH)_{a\mu}]

如果將 $\eta_{a\mu}$ 設定為較小的正數，上式就和常規的梯度下降等價。只要數字充分小，就能在更新時減小 $||V-WH||$ 。

如果我們按照斜向最陡調整變數，並設定：

(7): \eta_{a\mu}=\frac{H_{a\mu}}{(W^TWH)_{a\mu}}

就能得到定理 1 中給出的 H 更新規則。注意，該調整可得出乘性因子（分母中的梯度的正分量和因子的分子中的負分量的絕對值）。

對於散度公式，我們按照下述公式調整斜向最陡梯度下降：

(8): H_{a\mu} \leftarrow H_{a\mu} + \eta_{a\mu}[\sum_{i}W_{ia}\frac{V_{i\mu}}{(WH)_{i\mu}}-\sum_{i}W_{ia}]

同樣的，如果將 $\eta_{a\mu}$ 設定為較小的正數，上式就和常規的梯度下降等價。只要數字充分小，就能在更新時減小 $D(V||WH)$ 。如果設定：

(9): \eta_{a\mu} = \frac{H_{a\mu}}{\sum_i W_{ia}}

那麼就能得到定理 2 中給出的 H 更新規則。該調整可得出乘性因子（分母中的梯度的正分量和因子的分子中的負分量的絕對值）。

由於我們對 $\eta_{a\mu}$ 的取值並不夠小，看起來不能保證這種調整過後的梯度下降的代價函式減小。不過讓人驚訝的是，如下節所示，上述假設是事實。

收斂證明

我們將使用一個類似於 EM 演算法的輔助函式來證明定理 1 與定理 2。

定義 1： $G(h,h')$ 是 $F(h)$ 的輔助函式，滿足以下條件成立：

根據下面的引理，此輔助函式是一個有用的概念。（在圖1中的插圖也顯示了這一點）

引理 1：如果 G 為輔助函式，則 F 在下述更新時為非增：

證明： $F(h^{t+1}) \leq G(h^{t+1}, h^t) \leq G(h^t,h^t) = F(h^t)$

請注意，只有在 $h^t$ 為 $G(h,h^t)$ 的全域性最小值時滿足 $F(h^{t+1})=F(h^t)$ 。如果 F 的導數存在，且在 $h^{t}$ 的鄰域連續，也就是說 $\nabla F(h^t) = 0$ 。因此通過公式11反覆更新，我們就能得到目標函式收斂的區域性最小值 $h_{min} = \arg\min_h F(h)$