14.6 非负矩阵分解

非负矩阵的分解（Lee and Seung, 1999）是最近出现的类似主成分分析的另一个可选方法，它适用于数据和成分都是非负数。所以它可用于对例如图片这样的非负数据建模。

$N\times p$ 的数据矩阵 $\mathbf{X}$ 可近似为：

其中 $\mathbf{W}$ 的维度为 $N\times r$，$\mathbf{H}$ 的维度为 $r\times p$，$r\leq\max(N,p)$。假设 $x_{ij},w_{ik},h_{kj}\geq0$。

对下式最大化可得出矩阵 $\mathbf{W}$ 和 $\mathbf{H}$：

这其实就是假设 $x_{ij}$ 服从均值为 $(\mathbf{W}\mathbf{H})_{ij}$ 的泊松分布的模型的对数几率函数。对于正数数据来说，这个假设一般是合理的。

下式中的交替算法（Lee and Seung, 2001）可收敛至 $L(\mathbf{W},\mathbf{H})$ 的一个局部最大化解：

这个算法可从对 $L(\mathbf{W},\mathbf{H})$ 的最大化的最小化过程推导得出（MM 算法，练习 14.23），并且也与对数线型模型的 迭代比例尺度（iterative proportional scaling） 算法有关联（练习 14.24）。

图 14.33 展示了来自 Lee and Seung (1999) 的一个示例[^8]，其中对比了非负矩阵分解（NMF）、向量量化（VQ，等价于 K 均值聚类）、和主成分分析（PCA）。将这三个学习方法应用于一个 $N=2429$ 人脸图片的数据库中，每个图片包含了 $19\times19$ 个像素点，得到一个 $2429\times381$ 的数据矩阵 $\mathbf{X}$。如图中的 $7\times7$ 排列的拼接图所示（每个图有 $19\times19$ 个像素），每个方法都得到了一组 $r=49$ 个基图片。其中正数取值用黑色像素点表示，负数取值用红色像素点表示。右上展示了人脸图片的某个特定示例，将它近似为这些基图片的一个线性叠加。线性叠加的参数被展示在每个拼接图的旁边，为一个 $7\times7$ 的矩阵[^2]，等号右侧展示了叠加的结果。文章作者指出 NMF 与 VQ 和 PCA 不同，它学习出的表达人脸的是一组与人脸不同部分相像的基图片。

Donoho and Stodden (2004) 指出了非负矩阵分解的一个可能比较严重的问题。即使是在 $\mathbf{X}=\mathbf{W}\mathbf{H}$ 严格整理的情况下，这个分解也可能不是唯一的。图 14.34 演示了这个问题。数据点是在 $p=2$ 维度的空间上，在数据点和坐标之间有一片“空白区域”。因此可将基向量 $h_1$ 和 $h_2$ 选为这个区域上的任意两个向量，然后每个数据点都可被这两个向量的非负线性组合准确地表达出来。这个不唯一性意味着上述算法得到的解依赖于初始点，这可能会妨碍了这个分解的可解释性。不过尽管有这个解释性上的缺点，非负矩阵分解和它的应用仍然吸引到很多的关注。

14.6.1 原型分析

这个方法来自 Cutler and Breiman (1994)，使用本身就是数据点线性组合的原型点来近似数据点。从这点来看它与 K 均值聚类有相似之处。不过 原型分析（archetypal analysis） 并不是用一个最近的原型点来近似每个数据点，而是用一组原型点的一个凸性组合来近似每个数据点。凸性组合的使用导致了原型点会分布在数据云的一个凸包（convex hull）上。这样来看，这些原型点是“纯粹的（pure）”，或“典型的（archetypal）”。

与式 14.72 一样，将 $N\times p$ 的数据矩阵 $\mathbf{X}$ 近似为：

其中 $\mathbf{W}$ 的维度为 $N\times r$，$\mathbf{H}$ 的维度为 $r\times p$。假设 $w_{ik}\geq0$，并且 $\sum_{k=1}^rw_{ik}=1,\forall i$。所以 $p$ 维空间上的 $N$ 个数据点（$\mathbf{X}$ 的行）可表达为 $r$ 个原型点（$\mathbf{H}$ 的行）的凸性组合。同时假设：

其中 $\mathbf{B}$ 的维度为 $r\times N$，$b_{ki}\geq0$，并且 $\sum_{i=1}^Nb_{ki}=1,\forall k$。所以原型点本身就是数据点的凸性组合。基于式 14.75 和 14.76，我们求解下式最小化的权重矩阵 $\mathbf{W}$ 和 $\mathbf{B}$：

这个函数可以通过交替最小化的方式求解，每个最小化分别都是一个凸优化问题。不过这个问题整体上不是凸性的，所以这个算法收敛到的是准则函数的一个局部最小值解。

图 14.35 展示了二维空间上模拟数据的例子。上面三图展示的是原型分析的结果，而下面三图展示的是 K 均值聚类的结果。为了最优地通过原型点的凸性组合来重构数据点，就应该将原型点分布在数据点的凸包上。这在图 14.35 的上三图中能看得出，而且 Cutler and Breiman (1994) 也证明了在普遍情况下也的确如此。下三图中展示的 K 均值聚类则是将原型点置于数据点云的中间。

也可以将 K 均值聚类理解为原型分析的一个特例，$\mathbf{W}$ 的每一行中只有一个位置为一，而其他位置都是零。

同时也注意到原型分析模型（式 14.75）与非负矩阵分解模型（式 14.72）的表达式完全一样。尽管如此，两个模型所适用的场景和目的都有些不同。非负矩阵分解的目标是对数据矩阵 $\mathbf{X}$ 的列的近似，而其主要关注的输出结果是代表着数据中主要的非负成分的 $\mathbf{W}$ 列向量。原型分析关注的则是通过代表着原型数据点的 $\mathbf{H}$ 的行来对 $\mathbf{X}$ 的行的近似。同时，非负矩阵分解假设 $r\leq p$。如果 $r=p$，我们可令 $\mathbf{W}$ 等于将列向量缩放到和为一的数据矩阵 $\mathbf{X}$，就可得到一个完全的重构。与之相反，原型分析要求 $r\leq N$，但允许 $r>p$。例如在图 14.35 中，$p=2$、$N=50$，而 $r$ 等于 2、4、或 8。附加约束条件式 14.76 意味着，即使 $r>p$，原型分析的近似不会是完全准确的。

图 14.36 展示了在图 14.22 中数字“3”的数据库上应用了原型分析的结果。图 14.36 中的三行是分别指定了两个、三个、和四个原型点，三次使用原型分析所得出的原型点结果。与预期的相符，这个算法得出的原型点是大小和形状都比较极端的“3”。

本节练习

练习 14.23

Algorithm for non-negative matrix factorization (Wu and Lange, 2007). A function $g(x,y)$ to said to minorize a function $f(x)$ if

for all x, y in the domain. This is useful for maximizing $f(x)$ since it is easy to show that $f(x)$ is nondecreasing under the update

There are analogous definitions for majorization, for minimizing a function $f(x)$. The resulting algorithms are known as MM algorithms, for “minorize- maximize” or “majorize-minimize” (Lange, 2004). It also can be shown that the EM algorithm (8.5) is an example of an MM algorithm: see Sec- tion 8.5.3 and Exercise 8.2 for details.

1. Consider maximization of the function L(W, H) in (14.73), written here without the matrix notation $$L(\mathbf{W}, \mathbf{H}) = \sum_{i}^N \sum_{j=1}^p \left[ x_{ij} \log \left( \sum_{k=1}^r w_{ik} h_{kj} \right) -\sum_{k=1}^r w_{ik} h_{kj} \right]$$ Using the concavity of P log(x), show that for any set of r values $y_k\geq0$ and $0\leq c_k\leq1$ with $\sum_{k=1}^rc_k=1$, $$\log \left( \sum_{k=1}^r y_k \right) \geq \sum_{k=1}^r c_k \log(y_k / c_k)$$ Hence $$\log \left( \sum_{k=1}^r w_{ik} h_{kj} \right) \geq \sum_{k=1}^r \frac{a_{ikj}^s}{b_{ij}^s} \log \left( \frac{b_{ij}^s}{a_{ikj}^s} w_{ik} h_{kj}\right)$$ where $$a_{ikj}^s = w_{ik}^s h_{kj}^s \text{ and } b_{ij}^s = \sum_{k=1}^r w_{ik}^s h_{kj}^s$$ and s indicates the current iteration.

2. Hence show that, ignoring constants, the function $$\begin{align} g(\mathbf{W}, \mathbf{H} | \mathbf{W}^s, \mathbf{H}^s) =& \sum_{i=1}^N \sum_{j=1}^p \sum_{k=1}^r x_{ij} \frac{a_{ikj}^s}{b_{ij}^s} \left( \log w_{ik} + \log h_{kj} \right) \\ &- \sum_{i=1}^N \sum_{j=1}^p \sum_{k=1}^r w_{ik} h_{kj} \end{align}$$ minorizes $L(\mathbf{W},\mathbf{H})$.

3. Set the partial derivatives of $g(\mathbf{W},\mathbf{H}|\mathbf{W}^s,\mathbf{H}^s)$ to zero and hence derive the updating steps (14.74).

练习 14.24

Consider the non-negative matrix factorization (14.72) in the rank one case (r = 1).

1. Show that the updates (14.74) reduce to $$\begin{align} w_i &\leftarrow w_i \frac{\sum_{j=1}^p x_{ij}}{\sum_{j=1}^p w_i h_j} \\ h_j &\leftarrow h_j \frac{\sum_{i=1}^N x_{ij}}{\sum_{i=1}^N w_i h_j} \end{align}\tag{14.121}$$ where $w_i=w_{i1}$, $h_j=h_{1j}$ . This is an example of the iterative proportional scaling procedure, applied to the independence model for a two-way contingency table (Fienberg, 1977, for example).
2. Show that the final iterates have the explicit form $$w_i = c \cdot \frac{\sum_{j=1}^p x_{ij}}{\sum_{i=1}^N\sum_{j=1}^p x_{ij}}, h_k = \frac{1}{c} \cdot \frac{\sum_{i=1}^N x_{ij}}{\sum_{i=1}^N\sum_{j=1}^p x_{ij}} \tag{14.122}$$ for any constant c > 0. These are equivalent to the usual row and column estimates for a two-way independence model.