非负矩阵分解NMF的乘法更新收敛性与初始化策略

非负矩阵分解（NMF）是一种将一个大矩阵拆解为两个更小、非负矩阵的乘积的技术。它在数据压缩、特征提取和模式识别等领域应用广泛。本文将聚焦其最经典的优化算法——乘法更新规则，讲解它为何能稳定收敛，并提供一系列确保算法效果的初始化操作指南。

理解NMF与乘法更新算法

给定 一个非负矩阵 $V$，维度为 $n \times m$。我们的目标是找到两个非负矩阵 $W$ ($n \times k$) 和 $H$ ($k \times m$)，使得它们的乘积 $WH$ 能近似还原 $V$，即 $V \approx WH$。这里的 $k$ 是用户设定的低秩维度，通常 $k \ll n, m$。

衡量近似好坏的常用标准是最小化欧氏距离的平方（Frobenius范数）：
$$ \min_{W, H} \| V - WH \|_F^2 \quad \text{s.t.} \quad W \ge 0, H \ge 0 $$
其中，下标 $F$ 表示矩阵的Frobenius范数，$\ge 0$ 表示矩阵中的每个元素都大于等于零。

乘法更新算法通过交替更新 $W$ 和 $H$ 来逐步降低上述目标函数值。其更新规则如下：

更新 $H$ 的规则：
$$ H_{aj} \leftarrow H_{aj} \frac{(W^T V)_{aj}}{(W^T W H)_{aj}} $$

更新 $W$ 的规则：
$$ W_{ia} \leftarrow W_{ia} \frac{(V H^T)_{ia}}{(W H H^T)_{ia}} $$

在编程实现时，$W^T V$、$W^T W H$、$V H^T$ 和 $W H H^T$ 这些矩阵运算可以一次性计算出整个矩阵，从而避免低效的逐元素循环。例如，更新 $H$ 的整体操作是：

1. 计算 分子矩阵：Num = W^T @ V
2. 计算 分母矩阵：Den = W^T @ W @ H
3. 执行 元素级乘法更新：H = H * (Num / Den)

乘法更新为何能收敛？

乘法更新规则并非凭空设计，它有坚实的优化理论基础。其收敛性主要基于以下两点：

1. 作为非负约束下的梯度下降变体：该规则可以看作是在非负约束 $W \ge 0, H \ge 0$ 下，对目标函数 $\|V - WH\|_F^2$ 进行梯度下降的一种自适应步长实现。传统的梯度下降更新公式为 $H \leftarrow H - \eta \nabla_H$，其中 $\nabla_H$ 是目标函数对 $H$ 的梯度。乘法更新通过构造一个与 $H$ 当前值和梯度相关的步长 $\eta$，巧妙地保证了更新后 $H$ 的所有元素依然为正（因为分子分母均由非负矩阵运算得到）。这种自适应步长通常比固定步长的梯度下降更高效。

2. 单调下降性证明：数学上可以严格证明，在每次更新 $H$ 或 $W$ 后，目标函数 $\|V - WH\|_F^2$ 的值是单调不增的。由于目标函数有下界（非负矩阵，且范数平方非负），单调下降的序列必然收敛。因此，乘法更新算法保证了目标函数值沿着迭代路径稳定下降，最终趋于一个局部最小值或驻点。

注意：收敛性保证的是目标函数值下降并趋于稳定，但它并不能保证找到的是全局最优解，因为目标函数通常是非凸的。算法的最终结果高度依赖于初始矩阵 $W$ 和 $H$ 的选择。

初始化策略：决定结果质量的关键第一步

由于NMF问题的非凸性，不同的初始化会导致算法收敛到不同的局部最优解。一个糟糕的初始化可能导致结果非常差。以下是几种主流的初始化策略及其具体操作步骤。

1. 随机初始化（最基础）

这是最简单的方法，但结果方差大，可能效果不佳。

1. 生成 两个随机矩阵：$W_0$ 和 $H_0$。
2. 设定 随机数的分布。通常使用均匀分布或正态分布，例如在 [0, 1] 区间内均匀采样。
3. 确保 矩阵中不出现零值。一个实用的技巧是添加一个很小的正数，例如 1e-6，以避免后续更新中出现除以零的风险。
4. 示例 (伪代码)：

   import numpy as np
   # 假设 V 是 (n, m)， k 是目标秩
   W0 = np.random.rand(n, k) + 1e-6
   H0 = np.random.rand(k, m) + 1e-6

2. 基于SVD的初始化（推荐，更稳定）

利用数据矩阵 $V$ 本身的低秩结构信息进行初始化，通常能获得更快的收敛速度和更好的最终结果。

1. 对矩阵 $V$ 执行 截断奇异值分解（SVD）。保留前 $k$ 个最大的奇异值及其对应的奇异向量。
2. 提取 左奇异向量矩阵 $U_k$ ($n \times k$) 和右奇异向量矩阵 $V_k^T$ ($k \times m$)，以及奇异值对角阵 $S_k$ ($k \times k$)。
3. 构造 初始的 $W$ 和 $H$。因为 $V \approx U_k S_k V_k^T$，我们希望 $WH \approx U_k S_k V_k^T$。
4. 设定 $W_0 = U_k$， $H_0 = S_k V_k^T$。
5. 处理负值：SVD得到的 $U_k$ 和 $V_k^T$ 可能包含负数，违反非负约束。一个常用策略是取其绝对值。
6. 进行缩放：直接使用绝对值的矩阵可能导致 $WH$ 与原始 $V$ 的数值尺度不匹配。一个改进是：计算 $W_0 = |U_k|$， 然后 求解 一个简单的最小二乘问题 $H_0 = \arg\min_{H\ge0} \|V - W_0 H\|_F^2$。由于 $W_0$ 已固定，这个问题有闭式解（非负最小二乘），或者直接使用 H0 = np.linalg.lstsq(W0, V) 后再将负值置零。

3. 平均值初始化（简单有效）

一个简单但往往比纯随机更稳健的方法是使用数据本身的统计特性。

1. 计算 整个数据矩阵 $V$ 的元素平均值，记为 $\mu$。
2. 初始化 $W_0$ 和 $H_0$ 的所有元素为 $\mu$。
3. 添加 一个小的随机扰动，例如 W0 = mu * np.ones((n,k)) + 0.1 * np.random.rand(n, k)，以打破对称性。

4. 随机凸组合初始化

这种方法通过模拟数据点的“混合”过程来初始化，更具解释性。

1. 随机选择 $k$ 个数据点（$V$ 的列）作为初始的“聚类中心”。这些中心构成 $W_0$ 的列，$W_0 = [v_{i1}, v_{i2}, ..., v_{ik}]$。
2. 计算 每个数据点 $v_j$ 与这 $k$ 个中心的“相似度”（例如，使用余弦相似度或高斯核函数）。将相似度值归一化，使得对于每个点 $j$，其对应 $k$ 个相似度值之和为 1。
3. 将 这 $m$ 组归一化的相似度值，按列排列 形成初始矩阵 $H_0$。这样，$WH$ 就近似于用 $k$ 个中心的凸组合来重构每个数据点。

实施初始化策略的操作流程

无论你选择哪种策略，都应遵循以下操作流程以确保成功：

1. 分析数据：观察输入矩阵 $V$ 的尺寸、稀疏性和数值范围。这对选择初始化方法有指导意义。
2. 选择方法：
   • 对数据结构不了解或需要快速实验时，尝试 平均值初始化。
   • 追求稳定且较好的初始结果时，首选 基于SVD的初始化。
   • 对计算成本极度敏感且数据简单时，可使用 随机初始化并进行多次尝试。
3. 生成初始矩阵：严格按所选方法的步骤生成 $W_0$ 和 $H_0$。
4. 执行预检查：
   • 检查 $W_0$ 和 $H_0$ 的维度是否正确：$W_0$ 为 $n \times k$，$H_0$ 为 $k \times m$。
   • 检查 所有元素是否为非负数。如果发现负数（如来自SVD），执行 绝对值操作或非负最小二乘解。
   • 可选：计算 初始近似误差 $\|V - W_0 H_0\|_F$，用于后续评估算法进展。
5. 进入主循环：将 $W_0$ 和 $H_0$ 作为乘法更新迭代的起点，设定 最大迭代次数和收敛阈值（例如，当目标函数值变化小于 $1e-4$ 时停止）。
6. 监控与调优：记录 每次迭代后的目标函数值。如果收敛速度过慢或陷入平台期，考虑 重新初始化或尝试更高级的初始化策略。