第13章 NMF 与非负约束分解：可解释表示与分解优化

关键词：非负矩阵分解 (NMF)、乘法更新 (MU)、HALS、非负最小二乘 (NNLS)、单纯形几何、高光谱解混

学习目标：

  1. 直觉：理解为什么“非负”约束能带来“基于部分（Parts-based）”的表示，以及其与 PCA 的本质几何区别。
  2. 建模：掌握 Frobenius 范数、KL 散度及 $\beta$-Divergence 对应的噪声物理意义。
  3. 算法：超越经典的乘法更新，掌握更高效的 HALS 和投影梯度法。
  4. 应用：深入理解 NMF 在高光谱解混（Unmixing）中的核心地位。

13.1 动机：从“全脸叠加”到“五官拼图”

在第12章中，我们讨论了低秩模型。奇异值分解（SVD）和主成分分析（PCA）是寻找低秩表示的黄金标准，数学性质完美（正交性、解析解）。但在图像处理的某些领域，PCA 的结果令人费解。

13.1.1 符号的物理困境

PCA 的基向量（Eigenfaces）看起来像什么？它们通常是“幽灵脸”——包含正像素和负像素。要重建一张人脸，PCA 可能会说：“取一张平均脸，加上 0.5 倍的特征脸A，再减去 0.3 倍的特征脸B”。

然而，在物理世界中：

  • 光强不能是负的。
  • 化学物质的浓度不能是负的。
  • 出现的频率不能是负的。

13.1.2 心理学与加性模型

Lee 和 Seung 在 1999 年的 Nature 论文中指出，人类认知倾向于“基于部分（Parts-based）”的表示。非负矩阵分解（NMF） 强制所有元素非负，这迫使算法不能通过“减法”来消除误差，只能通过“加法”来组合特征。

  • 结果：NMF 学习到的基向量通常对应有意义的局部结构（如鼻翼、眼角、嘴唇）。
  • 本质：NMF 是一种加性组合（Additive Combination）模型。

13.2 几何视角：单纯形锥（Simplicial Cone）

13.2.1 数学定义

给定非负数据矩阵 $X \in \mathbb{R}^{m \times n}_{\ge 0}$（$n$ 个样本，特征维度 $m$），寻找两个非负矩阵 $W \in \mathbb{R}^{m \times r}_{\ge 0}$ 和 $H \in \mathbb{R}^{r \times n}_{\ge 0}$，使得： $$ X \approx WH $$ 通常秩 $r \ll \min(m, n)$。

13.2.2 PCA vs NMF 的几何对比

想象数据点分布在第一象限的一个扇形区域内。

      Feature 2
         ^
         |      PCA 寻找方差最大的轴
         |      (可以是负方向，正交)
         |            /
         |          /   数据点云
         |        / . . .
         |      / . . . . .
      W1 \    / . . . . . . .
 (NMF基)  \ /   . . . . . .     / W2 (NMF基)
           O -----------------------
            \     . . . .     /
             \      . .     /
              \           /
               \        /
                \     /
                 \  /
                  v
  • PCA：寻找数据方差最大的正交轴（椭球的主轴）。数据点是基向量的线性组合，系数可正可负。基向量本身通常都不在第一象限。
  • NMF：寻找一个锥（Cone），使得所有数据点都包含在这个锥内部（或尽可能靠近）。
    • 锥的边缘就是基向量 $W$ 的列。
    • 数据点 $x_j$ 是 $W$ 中列的非负线性组合（即 $x_j$ 落在 $W$ 张成的凸锥内）。

直觉：如果数据确实是由几个纯粹的端元混合而成（如几个颜色的油漆混合），NMF 试图找到这几个“极点”。

13.3 目标函数与噪声模型

选择什么样的 $D(X, WH)$ 取决于你假设数据中蕴含什么样的噪声。这是最大似然估计（MLE）的直接推论。

13.3.1 Frobenius 范数（高斯噪声）

$$ \min_{W,H \ge 0} \frac{1}{2} |X - WH|_F^2 = \frac{1}{2} \sum_{i,j} (X_{ij} - (WH)_{ij})^2 $$

  • 统计假设：$X_{ij} = (WH)_{ij} + \epsilon_{ij}, \quad \epsilon \sim \mathcal{N}(0, \sigma^2)$。
  • 特点：对大误差敏感，适合一般自然图像处理、去噪。
  • 梯度： $$ \nabla_W = (WH - X)H^T, \quad \nabla_H = W^T(WH - X) $$

13.3.2 广义 KL 散度（泊松噪声）

$$ \min_{W,H \ge 0} D_{KL}(X || WH) = \sum_{i,j} \left( X_{ij} \log \frac{X_{ij}}{(WH)_{ij}} - X_{ij} + (WH)_{ij} \right) $$

  • 统计假设：$X_{ij} \sim \text{Poisson}((WH)_{ij})$。
  • 适用场景：光子计数成像（如天文图像、低光照显微镜）、文本挖掘（词频计数）。
  • 注意：当 $(WH)_{ij} \to 0$ 且 $X_{ij}>0$ 时，损失趋于无穷大。这迫使模型不能轻易预测 0。

13.3.3 $\beta$-Divergence（统一族）

通过标量 $\beta$ 统一上述度量：

  • $\beta=2$: Euclid (Frobenius)
  • $\beta=1$: KL Divergence
  • $\beta=0$: Itakura-Saito (常用于音频谱分解)

13.4 优化算法详解

NMF 的目标函数是非凸的（Non-convex），可能有多个局部极小值。但它是一个双凸（Bi-convex）问题：固定 $W$，对 $H$ 是凸的；固定 $H$，对 $W$ 是凸的。

13.4.1 乘法更新 (Multiplicative Update, MU)

这是最经典（也是代码库中最常见）的算法。它本质上是一种梯度下降法，但步长 $\eta$ 被巧妙构造成乘法形式，从而自然保证非负性。

以 Frobenius 范数为例，梯度下降形式：$H \leftarrow H - \eta \nabla_H$。 令步长 $\eta_{kj} = \frac{H_{kj}}{(W^TWH)_{kj}}$，代入整理得：

算法流程：

  1. 初始化 $W, H$ 为非负随机矩阵。
  2. 重复直到收敛： $$ H_{kj} \leftarrow H_{kj} \frac{(W^T X)_{kj}}{(W^T W H)_{kj} + \epsilon} $$ $$ W_{ik} \leftarrow W_{ik} \frac{(X H^T)_{ik}}{(W H H^T)_{ik} + \epsilon} $$

    优点：代码极其简单，无需调参，自动保持非负。 缺点：收敛慢（一阶方法），容易陷入“零锁死”（一旦某元素变为0，乘法无法将其救回）。

13.4.2 分层交替最小二乘 (HALS / Rank-1 Residue)

这是目前处理标准 NMF 最高效的算法之一，比 MU 快得多。 思路：我们不一次更新整个 $H$，而是一次只更新 $H$ 的一行（对应一个特征）。这利用了残差的秩-1 结构。

对于第 $k$ 个特征（$W$ 的第 $k$ 列 $w_k$，$H$ 的第 $k$ 行 $h^{(k)}$）：

  1. 计算除去第 $k$ 个特征后的残差：$R^{(k)} = X - \sum_{j \neq k} w_j h^{(j)}$。
  2. 我们希望 $w_k h^{(k)} \approx R^{(k)}$。
  3. 这是一个简单的秩-1 近似问题，有闭式解，然后投影到非负空间： $$ h^{(k)} \leftarrow \max\left(0, \frac{w_k^T R^{(k)}}{|w_k|^2}\right) $$ $$ w_k \leftarrow \max\left(0, \frac{R^{(k)} (h^{(k)})^T}{|h^{(k)}|^2}\right) $$

13.4.3 投影梯度法 (Projected Gradient)

利用 Lin & More 的框架： $$ x_{k+1} = P_{\ge 0} [ x_k - \alpha_k \nabla f(x_k) ] $$ 其中关键在于步长 $\alpha_k$ 的搜索（通常使用 Armijo 规则）。这比 MU 有更好的收敛理论保证。

13.5 解的非唯一性与正则化

NMF 最大的理论痛点是不唯一性（Non-uniqueness）。 $$ WH = (WQ)(Q^{-1}H) $$ 只要 $Q$ 及其逆矩阵是非负的（例如排列矩阵和对角缩放矩阵），解就是等价的。更糟糕的是，在这个锥内部，基向量可以稍微“晃动”而不影响对数据的包围。

为了获得物理上有意义的解，我们需要先验（Prior）。

13.5.1 稀疏正则 (Sparse NMF)

如果你认为图像是由“少量”部件组成的，或者每个像素只包含“少量”物质： $$ \min |X - WH|_F^2 + \lambda |H|_1 + \gamma |W|_1 $$

  • $L_1$ on $H$：每个样本只激活很少的基（稀疏编码）。
  • $L_1$ on $W$：基向量本身是稀疏的（例如基向量大部分是黑色的，只有局部有值）。

13.5.2 空间/流形正则

在图像中，相邻像素通常具有相似的物质构成。 $$ \dots + \lambda \sum_{i,j} S_{ij} |h_i - h_j|^2 $$ 其中 $S$ 是像素邻接矩阵。这迫使系数矩阵 $H$ 在空间上平滑（Smooth NMF）。

13.6 核心应用：高光谱解混 (Hyperspectral Unmixing)

这是 NMF 在遥感和计算成像中最直观的应用。

场景：卫星拍摄地面，每个像素包含数百个波段的光谱信息。一个像素实际上是覆盖范围内草地、水体、路面的混合。

线性混合模型 (LMM)：

  • $X \in \mathbb{R}^{B \times P}$：观测数据（$B$ 个波段，$P$ 个像素）。
  • $W \in \mathbb{R}^{B \times r}$ (Endmembers/端元)：纯物质的光谱签名（如纯净水的反射光谱）。
  • $H \in \mathbb{R}^{r \times P}$ (Abundances/丰度)：物质在每个像素中的占比。

物理约束：

  1. 非负性：光强和占比非负。
  2. 和为一 (ASC)：一个像素内所有物质占比之和为 100%。 $$ \sum_{k=1}^r H_{kj} = 1, \quad \forall j $$

求解策略： 带有 ASC 约束的 NMF。通常通过在 $X$ 和 $W$ 下方拼接一行全为 $\delta$ 的值，或者使用投影到单纯形（Simplex Projection）的步骤来求解。

13.7 常见变体速览

| 变体名称 | 目标/约束 | 典型应用 |

变体名称 目标/约束 典型应用
Weighted NMF $|W \odot (X-WH)|_F^2$ 处理缺失数据（Inpainting）、置信度加权
Orthogonal NMF $H H^T = I$ 或 $W^T W = I$ 聚类（等价于 K-Means 的软化版本）
Convex NMF 约束 $W$ 必须是数据 $X$ 的凸组合 质心必须在数据簇内部，可解释性更强
Deep NMF $X \approx W_1 W_2 \dots W_k H$ 层级特征学习，深度学习的可解释替代

13.8 本章小结

  1. 物理意义：NMF 通过非负约束实现了“部分构成整体”的加性模型，比 PCA 更具物理可解释性。
  2. 几何本质：NMF 是在寻找一个包裹数据的单纯形锥（Simplicial Cone）。
  3. 优化求解：
    • 乘法更新（MU）最简单但慢。
    • HALS（分层交替最小二乘）通常是最高效的基准算法。
    • 问题是非凸的，初始化很关键。
  4. 正则化：为了克服旋转不确定性，稀疏（$L_1$）和空间平滑先验是必须的。
  5. 杀手级应用：高光谱解混完全依赖于 NMF 及其变体。

13.9 练习题

基础题

  1. 乘法更新的单调性： 虽然我们没有推导，但请验证：如果在 $W, H$ 的更新公式中，如果重构误差为 0（即 $X = WH$），那么更新因子是多少？$W$ 和 $H$ 还会变化吗？

    Hint观察公式中的分子分母。如果 $X=WH$，则 $W^TX = W^TWH$，分数项为1。

  2. K-Means 也是 NMF? 如果限制矩阵 $H$ 的每一列必须是 One-hot 向量（即只有一个元素为1，其余为0），NMF 等价于什么算法？此时 $W$ 的列代表什么？

    Hint这退化为 K-Means 聚类。$H$ 是聚类标签，$W$ 的列是聚类中心。

  3. 维度计算： 有一批 1000 张 $64 \times 64$ 的人脸图像。我们将每张图拉成向量。

    • $X$ 的维度是多少？
    • 如果设定秩 $r=49$，则 $W$ 和 $H$ 的参数量总和是多少？相比原始数据压缩了多少倍？
      Hint$X: 4096 \times 1000$。参数量：$4096 \times 49 + 49 \times 1000$。计算比值。

挑战题

  1. KL 散度的梯度推导： 对于 $D_{KL} = \sum (X \log \frac{X}{Y} - X + Y)$，其中 $Y=WH$。 请证明 $\frac{\partial D_{KL}}{\partial W_{ik}} = \sum_j (1 - \frac{X_{ij}}{(WH)_{ij}}) H_{kj}$。 并由此推导出 KL 散度下的乘法更新公式。

    Hint利用链式法则 $\frac{\partial D}{\partial Y} \cdot \frac{\partial Y}{\partial W}$。注意对数求导项。

  2. 初始化与 SVD： 既然 SVD 给了最好的低秩近似，我们能否用 SVD 结果初始化 NMF？但是 SVD 的 $U, V$ 包含负数。 请设计一种策略，利用 SVD 的结果构造非负的 $W_{init}, H_{init}$，使其比随机初始化更接近最优解。

    Hint参考 NNDSVD (Non-negative Double SVD) 算法。一种简单的思路是分离正负部：$x = x_+ - x_-$，但更严谨的方法是解秩-1 的非负近似。

  3. 缺失值处理： 如果在成像时，某些像素被死点掩膜（Mask）覆盖了，即 $X$ 中有未观测项。如何修改 NMF 的目标函数和乘法更新公式来处理这种情况？

    Hint引入二值掩膜矩阵 $M$，$ \min |M \odot (X - WH)|_F^2 $。更新公式中所有涉及 $X$ 和 $WH$ 的地方都需要乘上 $M$。

13.10 常见陷阱与错误 (Gotchas)

  1. 数据预处理错误（最常见！）： 很多习惯了 PCA 的人会先对数据做 中心化（Subtract Mean）。

    • 错误：$X \leftarrow X - \text{mean}(X)$。这会引入负数，直接导致 NMF 算法失效或物理意义丧失。
    • 正确：NMF 不需要中心化。如果必须处理背景，应将其作为 $W$ 的一列单独建模。

  2. 尺度模糊（Scaling Ambiguity）： 算法跑完后，你可能会发现 $W$ 的某些列数值极大（e.g., $10^5$），而对应的 $H$ 行数值极小（e.g., $10^{-5}$）。这虽然数学上正确，但在可视化和分析时非常头疼。

    • Fix：每次迭代或迭代结束后，手动归一化：$W_{\cdot k} \leftarrow W_{\cdot k} / |W_{\cdot k}|$, $H_{k \cdot} \leftarrow H_{k \cdot} \times |W_{\cdot k}|$。

  3. 零锁死（Zero-locking）： 在乘法更新中，如果初始化 $W_{ik}=0$ 或迭代中某一步变成了 0，它将永远保持为 0。

    • Fix：初始化不要包含 0。在分母和分子中加极小的 $\epsilon$。或者使用梯度类方法（如 ALS/PG）允许从 0 变回非 0。

  4. 秩 $r$ 的选择： 这没有解析解。通常通过观察残差曲线的“拐点”，或者使用交叉验证（Mask 掉一部分数据看预测能力）来确定。不要盲目相信 $r$ 越大越好，过大的 $r$ 会导致基向量重复或包含噪声（Overfitting）。