Python实现基于KL散度的NMF算法及两种KL散度对比
非负矩阵分解 (NMF, Non-negative Matrix Factorization) 是一种常用的数据降维和特征提取技术，在图像处理、文本挖掘、推荐系统等领域有着广泛的应用。NMF 的核心思想是将一个非负矩阵分解为两个非负矩阵的乘积，从而得到数据的低维表示。而KL散度 (Kullback-Leibler Divergence) 是一种衡量两个概率分布差异的指标，常用于NMF的目标函数。
今天，咱们就来聊聊如何用Python实现基于KL散度的NMF算法，并且对比两种不同的KL散度在不同数据集上的表现。别担心，我会尽量讲得通俗易懂，保证有一定Python基础的朋友都能跟上。
1. NMF和KL散度：理论基础
在深入代码之前，咱们先来简单回顾一下NMF和KL散度的理论知识。这部分对于理解后面的代码实现至关重要。
1.1 非负矩阵分解 (NMF)
给定一个非负矩阵 V (n × m)，NMF的目标是找到两个非负矩阵 W (n × k) 和 H (k × m)，使得：
V ≈ WH
其中，k 通常远小于 n 和 m，从而实现降维的目的。W 可以看作是基矩阵，H 可以看作是系数矩阵。NMF 的关键在于“非负”约束，这使得分解后的矩阵具有更好的可解释性。
举个例子，假设 V 是一个包含多个文档的词频矩阵，每一行代表一个文档，每一列代表一个词语。那么，NMF 可以将 V 分解为 W 和 H。W 可以看作是主题矩阵，每一行代表一个主题，每一列代表该主题下各个词语的权重；H 可以看作是文档-主题矩阵，每一行代表一个文档，每一列代表该文档属于各个主题的概率。
1.2 KL散度 (Kullback-Leibler Divergence)
KL散度，也叫相对熵，是用来衡量两个概率分布 P 和 Q 之间差异的非对称性度量。其定义如下：
KL(P || Q) = ∑ P(i) * log(P(i) / Q(i))
其中，P(i) 和 Q(i) 分别表示两个概率分布在第 i 个事件上的概率。KL散度越小，表示两个分布越接近。
在NMF中，我们可以用KL散度来衡量原始矩阵 V 和重构矩阵 WH 之间的差异。常用的KL散度有两种：
原始KL散度 (Standard KL Divergence)：
D_KL(V || WH) = ∑ (V_ij * log(V_ij / (WH)_ij) - V_ij + (WH)_ij)改进的KL散度 (Alternative KL Divergence): 有时候也称为 I-divergence
D_altKL(V || WH) = ∑( V_ij *( (WH)_ij/V_ij - log((WH)_ij/V_ij) -1))
这两种KL散度在数学上并不等价，它们在NMF中的表现也会有所不同。我们将在后面的实验中对比它们的性能。
2. Python代码实现
有了理论基础，咱们就可以开始动手写代码了。我们将使用Python的NumPy库来进行矩阵运算，使用Scikit-learn库来加载数据集和评估结果。
2.1 导入必要的库
import numpy as np
from sklearn.decomposition import NMF
from sklearn.datasets import fetch_20newsgroups
from sklearn.feature_extraction.text import TfidfVectorizer
from time import time
2.2 定义NMF函数 (使用原始KL散度)
def nmf_kl(V, k, max_iter=200, tol=1e-4):
"""基于原始KL散度的NMF算法
参数：
V: 输入的非负矩阵 (n × m)
k: 降维后的维度
max_iter: 最大迭代次数
tol: 收敛阈值
返回值：
W: 基矩阵 (n × k)
H: 系数矩阵 (k × m)
"""
n, m = V.shape
W = np.random.rand(n, k)
H = np.random.rand(k, m)
for i in range(max_iter):
# 更新 H
H = H * (W.T @ (V / (W @ H + 1e-9))) / (W.T @ np.ones((n, m)) + 1e-9)
# 更新 W
W = W * (((V / (W @ H + 1e-9)) @ H.T)) / (np.ones((n, m)) @ H.T + 1e-9)
# 计算损失函数
WH = W @ H
loss = np.sum(V * np.log(V / (WH + 1e-9) + 1e-9) - V + WH)
if i > 0 and abs(prev_loss - loss) < tol:
break
prev_loss = loss
if (i + 1) % 10 == 0:
print(f'Iteration: {i + 1}, Loss: {loss:.4f}')
return W, H
2.3 定义NMF函数 (使用改进的KL散度)
def nmf_alt_kl(V, k, max_iter=200, tol=1e-4):
"""基于改进的KL散度的NMF算法
参数和返回值同上
"""
n, m = V.shape
W = np.random.rand(n, k)
H = np.random.rand(k, m)
ones_matrix = np.ones_like(V)
for i in range(max_iter):
# 更新H
numerator_H = W.T @ ( (W @ H) / (V+1e-9) )
denominator_H = W.T @ (ones_matrix)
H = H * numerator_H / (denominator_H + 1e-9)
#更新W
numerator_W = (((W@H)/ (V+1e-9)) @ H.T)
denominator_W = ones_matrix @ H.T
W = W* numerator_W / (denominator_W + 1e-9)
# 计算损失函数
WH = W @ H
loss = np.sum(V * (WH/(V+1e-9) - np.log(WH/(V+1e-9)) -1))
if i > 0 and abs(prev_loss - loss) < tol:
break
prev_loss = loss
if (i + 1) % 10 == 0:
print(f'Iteration: {i + 1}, Loss: {loss:.4f}')
return W, H
2.4 代码解释
这两个函数的核心都是基于乘法更新规则 (Multiplicative Update Rule) 来迭代更新 W 和 H。乘法更新规则能够保证 W 和 H 的非负性，并且在理论上可以证明其收敛性。
1e-9是为了防止除零错误。max_iter和tol控制迭代次数和收敛条件。- 代码中计算了每次迭代后的损失函数值 (loss)，并打印出来，方便我们观察算法的收敛情况。
3. 实验对比
现在，我们来对比一下两种KL散度在不同数据集上的表现。我们将使用20 Newsgroups数据集，这是一个常用的文本分类数据集，包含20个不同主题的新闻组帖子。
3.1 加载数据
# 加载20 Newsgroups数据集
datasets = fetch_20newsgroups(shuffle=True, random_state=1, remove=('headers', 'footers', 'quotes'))
# 使用TF-IDF向量化文本
vectorizer = TfidfVectorizer(max_df=0.95, min_df=2, stop_words='english')
V = vectorizer.fit_transform(datasets.data).T # 注意这里转置了，使得每一列代表一个文档
print(f'Shape of V: {V.shape}')
3.2 运行NMF算法并计时
k = 10 # 设置降维后的维度
print("\n--- Standard KL Divergence ---")
t0 = time()
W_kl, H_kl = nmf_kl(V.toarray(), k)
print(f"Done in {time() - t0:.3f}s.")
print("\n--- Alternative KL Divergence ---")
t0 = time()
W_altkl, H_altkl = nmf_alt_kl(V.toarray(), k)
print(f"Done in {time() - t0:.3f}s.")
3.3 结果分析
运行上述代码后，你可以观察到两种KL散度下NMF算法的迭代次数、损失函数值和运行时间。一般来说，改进的KL散度 (nmf_alt_kl)可能会比原始KL散度收敛得更快，损失函数值更低，但实际情况会因数据集和参数设置而异。你可以尝试不同的 k 值、max_iter 和 tol 值，观察它们对结果的影响。
为了更直观地比较两种KL散度的性能，我们可以计算重构误差 (Reconstruction Error)，即原始矩阵 V 和重构矩阵 WH 之间的差异。我们可以使用 Frobenius 范数 (Frobenius Norm) 来衡量这种差异：
def reconstruction_error(V, W, H):
"""计算重构误差"""
return np.linalg.norm(V - W @ H)
print(f"\nReconstruction Error (Standard KL): {reconstruction_error(V.toarray(), W_kl, H_kl):.4f}")
print(f"Reconstruction Error (Alternative KL): {reconstruction_error(V.toarray(), W_altkl, H_altkl):.4f}")
4. 总结与拓展
今天，我们学习了如何用Python实现基于KL散度的NMF算法，并对比了两种不同的KL散度在20 Newsgroups数据集上的表现。通过实验，我们可以发现，不同的KL散度在NMF中的性能会有所差异。改进的KL散度在某些情况下可能会有更好的表现。
当然，NMF算法还有很多可以改进和拓展的地方，比如：
- 不同的初始化方法： 我们可以尝试不同的 W 和 H 初始化方法，比如随机初始化、SVD初始化等，看看它们对结果的影响。
- 不同的更新规则： 除了乘法更新规则，还有其他一些更新规则，比如梯度下降法、投影梯度法等。
- 正则化： 我们可以在目标函数中加入正则化项，比如 L1 正则化或 L2 正则化，来控制 W 和 H 的稀疏性。
- 稀疏NMF: 如果需要更稀疏的表示，可以尝试稀疏NMF算法。
- 在线NMF: 如果数据量很大，可以尝试在线NMF算法，它可以分批处理数据，减少内存消耗。
希望这篇文章能帮助你更好地理解NMF算法和KL散度，并在实际应用中发挥作用。如果你有任何问题或者想法，欢迎在评论区留言交流!
