非负矩阵分解(NMF)详解
引言
非负矩阵分解(Non-negative Matrix Factorization,简称NMF)是一种在机器学习和数据挖掘领域广泛应用的矩阵分解技术。与传统的矩阵分解方法(如奇异值分解SVD)不同,NMF的最大特点是对分解后的矩阵元素施加了非负约束,这使得分解结果具有更好的可解释性,特别适合处理自然界中的非负数据,如图像、文本和音频等。
基本原理
假设我们有一个非负矩阵 V ∈ R + m × n V \in \mathbb{R}^{m \times n}_{+} V∈R+m×n(其中 R + m × n \mathbb{R}^{m \times n}_{+} R+m×n 表示所有元素非负的 m × n m \times n m×n 实矩阵空间),NMF的目标是找到两个非负矩阵 W ∈ R + m × r W \in \mathbb{R}^{m \times r}_{+} W∈R+m×r 和 H ∈ R + r × n H \in \mathbb{R}^{r \times n}_{+} H∈R+r×n,使得:
V ≈ W × H V \approx W \times H V≈W×H
其中 r r r 通常选择为远小于 m m m 和 n n n 的值(即 r ≪ min ( m , n ) r \ll \min(m, n) r≪min(m,n)),这样可以实现数据的降维和压缩。
从直观上理解,如果将 V V V 视为包含 n n n 个数据样本(每个样本有 m m m 个特征)的数据集,那么 W W W 可以被视为 r r r 个基向量(或主题),而 H H H 则表示这些基向量的组合系数。由于所有元素都是非负的,我们可以将这种分解理解为一种"加性模型",即每个数据样本都是由基本组件的非负组合构成的。
目标函数与优化
为了找到最佳的 W W W 和 H H H,我们需要定义一个目标函数来衡量 V V V 与 W × H W \times H W×H 之间的差异。常用的目标函数包括:
1. Frobenius范数
最常用的是基于Frobenius范数的平方误差:
min W , H ≥ 0 ∥ V − W H ∥ F 2 = ∑ i = 1 m ∑ j = 1 n ( V i j − ( W H ) i j ) 2 \min_{W, H \geq 0} \|V - WH\|_F^2 = \sum_{i=1}^{m}\sum_{j=1}^{n}(V_{ij} - (WH)_{ij})^2 W,H≥0min∥V−WH∥F2=i=1∑mj=1∑n(Vij−(WH)ij)2
展开这个表达式,我们可以得到:
∥ V − W H ∥ F 2 = Tr ( ( V − W H ) T ( V − W H ) ) = Tr ( V T V − V T W H − H T W T V + H T W T W H ) \|V - WH\|_F^2 = \text{Tr}((V-WH)^T(V-WH)) = \text{Tr}(V^TV - V^TWH - H^TW^TV + H^TW^TWH) ∥V−WH∥F2=Tr((V−WH)T(V−WH))=Tr(VTV−VTWH−HTWTV+HTWTWH)
其中 Tr ( ⋅ ) \text{Tr}(\cdot) Tr(⋅) 表示矩阵的迹。
2. 散度(Divergence)
另一种常用的目标函数是基于Kullback-Leibler散度的:
min W , H ≥ 0 D K L ( V ∥ W H ) = ∑ i = 1 m ∑ j = 1 n ( V i j log V i j ( W H ) i j − V i j + ( W H ) i j ) \min_{W, H \geq 0} D_{KL}(V\|WH) = \sum_{i=1}^{m}\sum_{j=1}^{n}\left(V_{ij}\log\frac{V_{ij}}{(WH)_{ij}} - V_{ij} + (WH)_{ij}\right) W,H≥0minDKL(V∥WH)=i=1∑mj=1∑n(Vijlog(WH)ijVij−Vij+(WH)ij)
这个散度函数在信息论中有深刻的意义,它衡量了两个概率分布之间的差异。当我们将矩阵 V V V 和 W H WH WH 归一化为概率分布时,KL散度提供了一种自然的度量方式。无论选择哪种目标函数,NMF的优化问题都是非凸的,这意味着可能存在多个局部最优解。因此,通常采用迭代算法来求解。
求解算法
乘法更新规则(Multiplicative Update Rules)
Lee和Seung在2001年提出的乘法更新规则是最经典的NMF求解算法。对于Frobenius范数目标函数,更新规则为:
H k j ← H k j ( W T V ) k j ( W T W H ) k j H_{kj} \leftarrow H_{kj} \frac{(W^T V)_{kj}}{(W^T W H)_{kj}} Hkj←Hkj(WTWH)kj(WTV)kj
W i k ← W i k ( V H T ) i k ( W H H T ) i k W_{ik} \leftarrow W_{ik} \frac{(V H^T)_{ik}}{(W H H^T)_{ik}} Wik←Wik(WHHT)ik(VHT)ik
这些更新规则可以通过梯度下降法推导出来。例如,对于 H H H,目标函数关于 H k j H_{kj} Hkj 的梯度为:
∂ ∂ H k j ∥ V − W H ∥ F 2 = − 2 ( W T V ) k j + 2 ( W T W H ) k j \frac{\partial}{\partial H_{kj}} \|V - WH\|_F^2 = -2(W^TV)_{kj} + 2(W^TWH)_{kj} ∂Hkj∂∥V−WH∥F2=−2(WTV)kj+2(WTWH)kj
设定学习率为 η k j = H k j 2 ( W T W H ) k j \eta_{kj} = \frac{H_{kj}}{2(W^TWH)_{kj}} ηkj=2(WTWH)kjHkj,则梯度下降更新为:
H k j ← H k j − η k j ⋅ ∂ ∂ H k j ∥ V − W H ∥ F 2 = H k j ( W T V ) k j ( W T W H ) k j H_{kj} \leftarrow H_{kj} - \eta_{kj} \cdot \frac{\partial}{\partial H_{kj}} \|V - WH\|_F^2 = H_{kj} \frac{(W^T V)_{kj}}{(W^T W H)_{kj}} Hkj←Hkj−ηkj⋅∂Hkj∂∥V−WH∥F2=Hkj(WTWH)kj(WTV)kj
对于KL散度目标函数,更新规则为:
H k j ← H k j ∑ i = 1 m W i k V i j ( W H ) i j ∑ i = 1 m W i k H_{kj} \leftarrow H_{kj} \frac{\sum_{i=1}^{m} W_{ik} \frac{V_{ij}}{(WH)_{ij}}}{\sum_{i=1}^{m} W_{ik}} Hkj←Hkj∑i=1mWik∑i=1mWik(WH)ijVij
W i k ← W i k ∑ j = 1 n H k j V i j ( W H ) i j ∑ j = 1 n H k j W_{ik} \leftarrow W_{ik} \frac{\sum_{j=1}^{n} H_{kj} \frac{V_{ij}}{(WH)_{ij}}}{\sum_{j=1}^{n} H_{kj}} Wik←Wik∑j=1nHkj∑j=1nHkj(WH)ijVij
这些更新规则保证了在每次迭代中目标函数单调递减,且自动满足非负约束(因为是乘法更新,只要初始值非负,更新后的值也一定非负)。
交替最小二乘法(Alternating Least Squares, ALS)
另一种常用的方法是交替最小二乘法,其基本思路是:
-
固定 W W W,求解 H H H:
min H ≥ 0 ∥ V − W H ∥ F 2 \min_{H \geq 0} \|V - WH\|_F^2 H≥0min∥V−WH∥F2 -
固定 H H H,求解 W W W:
min W ≥ 0 ∥ V − W H ∥ F 2 \min_{W \geq 0} \|V - WH\|_F^2 W≥0min∥V−WH∥F2
在每一步中,我们都将问题转化为带非负约束的最小二乘问题。例如,当固定 W W W 时,对于 H H H 的每一列 h j h_j hj,我们需要解决:
min h j ≥ 0 ∥ v j − W h j ∥ 2 2 \min_{h_j \geq 0} \|v_j - Wh_j\|_2^2 hj≥0min∥vj−Whj∥22
其中 v j v_j vj 是 V V V 的第 j j j 列。这是一个标准的非负最小二乘问题,可以使用投影梯度下降等方法求解。
NMF的几何解释
从几何角度看,NMF可以理解为在非负象限中寻找一个由 r r r 个向量生成的锥体,使得原始数据点尽可能地落在这个锥体内或附近。这与主成分分析(PCA)寻找最佳线性子空间的思路不同,NMF寻找的是非负锥体,这使得分解结果具有更好的可解释性。具体来说,如果将 V V V 的列向量视为数据点,那么 W W W 的列向量可以视为生成锥体的基向量,而 H H H 的每一列则表示将对应的数据点表示为这些基向量的非负线性组合的系数。
从数学上讲,NMF寻找的是一个凸锥(convex cone):
C = { ∑ k = 1 r α k w k : α k ≥ 0 , k = 1 , 2 , … , r } \mathcal{C} = \left\{ \sum_{k=1}^{r} \alpha_k w_k : \alpha_k \geq 0, k = 1,2,\ldots,r \right\} C={k=1∑rαkwk:αk≥0,k=1,2,…,r}
其中 w k w_k wk 是 W W W 的第 k k k 列。原始数据点被近似为这个凸锥中的点。
NMF的变种与扩展
稀疏NMF
为了获得更加稀疏的解,可以在目标函数中添加L1正则化项:
min W , H ≥ 0 ∥ V − W H ∥ F 2 + α ∥ W ∥ 1 + β ∥ H ∥ 1 \min_{W, H \geq 0} \|V - WH\|_F^2 + \alpha\|W\|_1 + \beta\|H\|_1 W,H≥0min∥V−WH∥F2+α∥W∥1+β∥H∥1
其中 ∥ W ∥ 1 = ∑ i , k ∣ W i k ∣ \|W\|_1 = \sum_{i,k}|W_{ik}| ∥W∥1=∑i,k∣Wik∣ 和 ∥ H ∥ 1 = ∑ k , j ∣ H k j ∣ \|H\|_1 = \sum_{k,j}|H_{kj}| ∥H∥1=∑k,j∣Hkj∣ 分别是 W W W 和 H H H 的L1范数, α \alpha α 和 β \beta β 是正则化参数。
对于这种正则化的NMF,更新规则变为:
H k j ← H k j ( W T V ) k j ( W T W H ) k j + β H_{kj} \leftarrow H_{kj} \frac{(W^T V)_{kj}}{(W^T W H)_{kj} + \beta} Hkj←Hkj(WTWH)kj+β(WTV)kj
W i k ← W i k ( V H T ) i k ( W H H T ) i k + α W_{ik} \leftarrow W_{ik} \frac{(V H^T)_{ik}}{(W H H^T)_{ik} + \alpha} Wik←Wik(WHHT)ik+α(VHT)ik
正交NMF
为了使基向量之间更加正交,可以添加正交约束:
min W , H ≥ 0 ∥ V − W H ∥ F 2 s.t. W T W = I \min_{W, H \geq 0} \|V - WH\|_F^2 \quad \text{s.t.} \quad W^TW = I W,H≥0min∥V−WH∥F2s.t.WTW=I
这里 I I I 是单位矩阵。正交约束使得基向量之间相互独立,有助于提高解的唯一性和可解释性。实际上,完全正交的约束在非负条件下可能过于严格,因此通常采用软约束形式:
min W , H ≥ 0 ∥ V − W H ∥ F 2 + λ ∥ W T W − I ∥ F 2 \min_{W, H \geq 0} \|V - WH\|_F^2 + \lambda\|W^TW - I\|_F^2 W,H≥0min∥V−WH∥F2+λ∥WTW−I∥F2
其中 λ \lambda λ 是权衡参数。
卷积NMF
对于时序数据,可以使用卷积NMF,其中 V V V 被建模为 W W W 和 H H H 的卷积:
V ≈ ∑ k = 1 r W k ∗ H k V \approx \sum_{k=1}^{r} W^k * H^k V≈k=1∑rWk∗Hk
其中 ∗ * ∗ 表示卷积操作, W k W^k Wk 是第 k k k 个基本模式, H k H^k Hk 是其激活。
在频域中,卷积可以表示为:
F ( V ) ≈ ∑ k = 1 r F ( W k ) ⊙ F ( H k ) \mathcal{F}(V) \approx \sum_{k=1}^{r} \mathcal{F}(W^k) \odot \mathcal{F}(H^k) F(V)≈k=1∑rF(Wk)⊙F(Hk)
其中 F \mathcal{F} F 表示傅里叶变换, ⊙ \odot ⊙ 表示元素级乘法。
深入理解NMF的数学性质
收敛性分析
Lee和Seung证明了乘法更新规则保证目标函数单调递减。对于Frobenius范数目标函数,我们可以构造一个辅助函数 G ( H , H ′ ) G(H, H') G(H,H′),满足:
F
(
H
)
≤
G
(
H
,
H
′
)
F(H) \leq G(H, H')
F(H)≤G(H,H′)
F
(
H
)
=
G
(
H
,
H
)
F(H) = G(H, H)
F(H)=G(H,H)
其中 F ( H ) = ∥ V − W H ∥ F 2 F(H) = \|V - WH\|_F^2 F(H)=∥V−WH∥F2, H ′ H' H′ 是当前估计。通过最小化 G ( H , H ′ ) G(H, H') G(H,H′) 关于 H H H 的值,我们可以得到更新规则:
H k j ← H k j ′ ( W T V ) k j ( W T W H ′ ) k j H_{kj} \leftarrow H'_{kj} \frac{(W^T V)_{kj}}{(W^T W H')_{kj}} Hkj←Hkj′(WTWH′)kj(WTV)kj
可以证明,这个更新规则保证 F ( H ) ≤ F ( H ′ ) F(H) \leq F(H') F(H)≤F(H′),即目标函数单调递减。
唯一性问题
与SVD不同,NMF的解通常不是唯一的。具体来说,对于任意可逆矩阵 S S S(满足 S S S 和 S − 1 S^{-1} S−1 都是非负的),如果 V ≈ W H V \approx WH V≈WH,那么 V ≈ ( W S ) ( S − 1 H ) V \approx (WS)(S^{-1}H) V≈(WS)(S−1H) 也是一个有效的分解。然而,在实践中,由于问题的非凸性和算法的初始化条件,NMF通常会收敛到特定的局部最优解。此外,添加额外的约束(如稀疏性或正交性)可以减少解的不确定性。
复杂度和扩展性
对于大规模数据,NMF的计算复杂度可能成为一个挑战。标准乘法更新规则的每次迭代复杂度为 O ( m n r ) O(mnr) O(mnr),其中 m m m 和 n n n 是矩阵 V V V 的维度, r r r 是潜在维度。为了处理大规模数据,研究人员提出了各种加速技术,如随机梯度下降、分布式计算和在线学习等。例如,在线NMF算法可以逐块处理数据,更新公式为:
W t = arg min W ≥ 0 ∥ V t − W H t ∥ F 2 + λ ∥ W − W t − 1 ∥ F 2 W_t = \arg\min_{W \geq 0} \|V_t - WH_t\|_F^2 + \lambda\|W - W_{t-1}\|_F^2 Wt=argW≥0min∥Vt−WHt∥F2+λ∥W−Wt−1∥F2
H t = arg min H ≥ 0 ∥ V t − W t H ∥ F 2 H_t = \arg\min_{H \geq 0} \|V_t - W_tH\|_F^2 Ht=argH≥0min∥Vt−WtH∥F2
其中 V t V_t Vt 是第 t t t 个数据块, W t W_t Wt 和 H t H_t Ht 是相应的更新。
NMF的理论基础与连接
与概率模型的联系
NMF可以从概率角度理解。例如,当使用KL散度作为目标函数时,NMF等价于假设数据服从泊松分布的概率模型:
V i j ∼ Poisson ( ( W H ) i j ) V_{ij} \sim \text{Poisson}((WH)_{ij}) Vij∼Poisson((WH)ij)
最大化似然函数等价于最小化KL散度:
max W , H ≥ 0 log P ( V ∣ W , H ) ≡ min W , H ≥ 0 D K L ( V ∥ W H ) \max_{W, H \geq 0} \log P(V|W,H) \equiv \min_{W, H \geq 0} D_{KL}(V\|WH) W,H≥0maxlogP(V∣W,H)≡W,H≥0minDKL(V∥WH)
这种概率解释使得NMF可以自然地扩展到贝叶斯框架中,引入先验分布和不确定性估计。
与其他矩阵分解方法的比较
NMF与其他矩阵分解方法(如SVD、PCA)有着密切的联系。例如,当我们放松非负约束时,Frobenius范数下的NMF问题退化为截断SVD问题。
然而,非负约束导致了本质上不同的解。SVD寻找的是正交基向量,这些基向量可能包含正负值,并且通常是全局的(影响整个数据空间);而NMF寻找的是非负基向量,这些基向量通常是局部的、部分的表示。从数学上看,SVD解决的是:
min U , Σ , V ∥ X − U Σ V T ∥ F 2 s.t. U T U = I , V T V = I \min_{U, \Sigma, V} \|X - U\Sigma V^T\|_F^2 \quad \text{s.t.} \quad U^TU = I, V^TV = I U,Σ,Vmin∥X−UΣVT∥F2s.t.UTU=I,VTV=I
而NMF解决的是:
min W , H ∥ X − W H ∥ F 2 s.t. W ≥ 0 , H ≥ 0 \min_{W, H} \|X - WH\|_F^2 \quad \text{s.t.} \quad W \geq 0, H \geq 0 W,Hmin∥X−WH∥F2s.t.W≥0,H≥0
NMF变体
半监督NMF
在有标签信息的情况下,可以将监督信息整合到NMF中:
min W , H ≥ 0 ∥ V − W H ∥ F 2 + γ ∥ Y − Q H ∥ F 2 \min_{W, H \geq 0} \|V - WH\|_F^2 + \gamma\|Y - QH\|_F^2 W,H≥0min∥V−WH∥F2+γ∥Y−QH∥F2
其中 Y Y Y 是标签矩阵, Q Q Q 是标签投影矩阵, γ \gamma γ 是权衡参数。这样,分解不仅捕捉数据的内在结构,还利用了标签信息。
鲁棒NMF
为了处理异常值和噪声,可以使用L1范数代替Frobenius范数:
min W , H ≥ 0 ∥ V − W H ∥ 1 = ∑ i , j ∣ V i j − ( W H ) i j ∣ \min_{W, H \geq 0} \|V - WH\|_1 = \sum_{i,j} |V_{ij} - (WH)_{ij}| W,H≥0min∥V−WH∥1=i,j∑∣Vij−(WH)ij∣
L1范数对异常值不如L2范数敏感,因此能提供更鲁棒的分解。