非负矩阵分解不收敛

Question

我正在尝试使用Kullback-Liebler散度作为相似性度量来实现非负矩阵分解。该算法在：http://hebb.mit.edu/people/seung/papers/nmfconverge.pdf中描述。下面是我的python / numpy实现，以及一个运行它的示例矩阵。

简而言之，该算法应该学习矩阵W(n×r)和H(r×m)，使得V(n×m)近似为WH。你从W和H的随机值开始，通过遵循Seung和Lee论文中描述的更新规则，你应该越来越接近W和H的良好近似值。

事实证明，该算法可以单调地减少发散度，但在我的实现中不会发生这种情况。相反，它会在两个散度值之间进行交替。如果你看一下W和H，你会发现由此产生的因子分解并不是特别好。

我想知道在计算W的更新时是使用更新的H还是旧的H。我尝试了这两种方法，但它不会改变实现的行为。

我已经根据这篇论文检查了我的实现很多次，我看不出我做错了什么。有没有人能解释一下这个问题？

import numpy as np

def 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, H


def 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, H


V = np.arange(0.01,1.01,0.01).reshape(10,10)

W, H = factor(V, 6)

Answer 1

如何消除交替效应：

定理证明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)，交互效果就会消失。

简化代码的：

在处理Python数组时，使用它们的mean和sum方法，并避免使用NumPy 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/WH
H *= (np.dot(V_over_WH.T, W) / W.sum(axis=0)).T

Explanation:

如果您查看H的公式5更新规则，首先请注意V和(W H)的索引是相同的。因此，您可以将V / (W H)替换为

V_over_WH = V/WH

接下来，请注意，在分子中，我们对索引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和µ按该顺序索引。这与H的索引相同。所以我们想要用H乘以这个比率。太完美了。默认情况下，NumPy按元素对数组进行乘法。

因此，我们可以将H的整个更新规则表示为

H *= (np.dot(V_over_WH.T, W) / W.sum(axis=0)).T

So，把所有这些放在一起：

import numpy as np
np.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_WH


def 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, H

V = 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)