python - 非负矩阵分解无法收敛

Question

我正在尝试使用 Kullback-Liebler 散度作为相似性度量来实现非负矩阵分解。该算法描述于: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 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 的最后一行是这样的，

By reversing the roles of H and W, the update rule for W can similarly be shown to be nonincreasing.

因此我们可以推测更新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/WH
H *= (np.dot(V_over_WH.T, W) / W.sum(axis=0)).T

---

解释:

如果您查看 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 和 mu 按顺序索引的矩阵。这与 H 的索引相同。所以我们想将 H 乘以这个比率元素。完美的。 NumPy 默认按元素相乘。

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

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

---

所以，把它们放在一起:

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)