我正在尝试使用库尔巴克-莱布勒散度作为相似性度量来实现非负矩阵分解。该算法在以下文献中有描述：http://hebb.mit.edu/people/seung/papers/nmfconverge.pdf。下面是我在Python/Numpy中的实现，以及一个用于运行的示例矩阵。
简而言之，该算法应该学习矩阵W(n x r)和H(r x m)，使得V(n x m)近似等于WH。你从W和H中的随机值开始，通过遵循Seung和Lee论文中描述的更新规则，你应该越来越接近W和H的良好近似值。
该算法被证明可以单调地减少散度度量，但在我的实现中却不是这样。相反，它会在两个散度值之间交替。如果你查看W和H，你会发现得到的分解效果并不好。
我在计算W的更新时，考虑过是使用更新后的H还是旧的H。我尝试了两种方式，但这并没有改变实现的行为。
我已经多次将我的实现与论文进行了对比，我没有发现自己做错了什么。谁能帮我解答这个问题？
import numpy as npdef update(V, W, H, r, n, m): n,m = V.shape WH = W.dot(H) # equation (5) H_coeff = np.zeros(H.shape) for a in range(r): for mu in range(m): for i in range(n): H_coeff[a, mu] += W[i, a] * V[i, mu] / WH[i, mu] H_coeff[a, mu] /= sum(W)[a] H = H * H_coeff W_coeff = np.zeros(W.shape) for i in range(n): for a in range(r): for mu in range(m): W_coeff[i, a] += H[a, mu] * V[i, mu] / WH[i, mu] W_coeff[i, a] /= sum(H.T)[a] W = W * W_coeff return W, Hdef factor(V, r, iterations=100): n, m = V.shape avg_V = sum(sum(V))/n/m W = np.random.random(n*r).reshape(n,r)*avg_V H = np.random.random(r*m).reshape(r,m)*avg_V for i in range(iterations): WH = W.dot(H) divergence = sum(sum(V * np.log(V/WH) - V + WH)) # equation (3) print "At iteration " + str(i) + ", the Kullback-Liebler divergence is", divergence W,H = update(V, W, H, r, n, m) return W, HV = np.arange(0.01,1.01,0.01).reshape(10,10)W, H = factor(V, 6)回答：
如何消除交替效应：
定理2证明的最后一行写道，
通过交换H和W的角色，可以同样证明W的更新规则是非增的。
因此我们可以推测，更新H可以独立于更新W进行。这意味着在更新H之后：
H = H * H_coeff我们还应该在更新W之前更新中间值WH：
WH = W.dot(H)W = W * W_coeff这两个更新都会减少散度。
尝试一下：只需在计算W_coeff之前插入WH = W.dot(H)，交替效应就会消失。
简化代码：
在处理NumPy数组时，使用它们的mean和sum方法，避免使用Python的sum函数：
avg_V = sum(sum(V))/n/m可以写成
avg_V = V.mean()以及
divergence = sum(sum(V * np.log(V/WH) - V + WH)) # equation (3)可以写成
divergence = ((V * np.log(V_over_WH)) - V + WH).sum() 避免使用Python内置的sum函数，因为
- 它比NumPy的
sum方法慢，并且 - 它不如NumPy的
sum方法那样多功能。（它不允许你指定要求和的轴。我们设法用一次对NumPy的sum或mean方法的调用，消除了对Python的sum的两次调用。）
消除三重for循环：
但在速度和可读性方面更大的改进可以通过替换以下代码来实现
H_coeff = np.zeros(H.shape)for a in range(r): for mu in range(m): for i in range(n): H_coeff[a, mu] += W[i, a] * V[i, mu] / WH[i, mu] H_coeff[a, mu] /= sum(W)[a]H = H * H_coeff替换为
V_over_WH = V/WHH *= (np.dot(V_over_WH.T, W) / W.sum(axis=0)).T解释：
如果你查看方程5中H的更新规则，首先注意到V和(W H)的索引是相同的。所以你可以用
V_over_WH = V/WH来替换V / (W H)。接下来，注意在分子中我们是对索引i求和，i是W和V_over_WH的第一个索引。我们可以将其表示为矩阵乘法：
np.dot(V_over_WH.T, W).T而分母只是：
W.sum(axis=0).T如果我们将分子和分母相除
(np.dot(V_over_WH.T, W) / W.sum(axis=0)).T我们得到一个由剩余的两个索引alpha和mu按此顺序索引的矩阵。这与H的索引相同。所以我们希望按元素将H乘以这个比率。完美。NumPy默认按元素乘法数组。
因此，我们可以将H的整个更新规则表示为
H *= (np.dot(V_over_WH.T, W) / W.sum(axis=0)).T所以，总结起来：
import numpy as npnp.random.seed(1)def update(V, W, H, WH, V_over_WH): # equation (5) H *= (np.dot(V_over_WH.T, W) / W.sum(axis=0)).T WH = W.dot(H) V_over_WH = V / WH W *= np.dot(V_over_WH, H.T) / H.sum(axis=1) WH = W.dot(H) V_over_WH = V / WH return W, H, WH, V_over_WHdef factor(V, r, iterations=100): n, m = V.shape avg_V = V.mean() W = np.random.random(n * r).reshape(n, r) * avg_V H = np.random.random(r * m).reshape(r, m) * avg_V WH = W.dot(H) V_over_WH = V / WH for i in range(iterations): W, H, WH, V_over_WH = update(V, W, H, WH, V_over_WH) # equation (3) divergence = ((V * np.log(V_over_WH)) - V + WH).sum() print("At iteration {i}, the Kullback-Liebler divergence is {d}".format( i=i, d=divergence)) return W, HV = np.arange(0.01, 1.01, 0.01).reshape(10, 10)# V = np.arange(1,101).reshape(10,10).astype('float')W, H = factor(V, 6)