4

我正在尝试使用 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)
4

1 回答 1

8

如何消除交替效应：

定理 2 证明的最后一行是：

通过颠倒 H 和 W 的角色，W 的更新规则可以类似地被证明是非递增的。

因此我们可以推测更新H可以独立于更新来完成 W。这意味着更新后H：

H = H * H_coeff

我们还应该WH在更新之前更新中间值W：

WH = W.dot(H)
W = W * W_coeff

两个更新都减少了分歧。

试试看：WH = W.dot(H)在计算 之前坚持W_coeff，交替效应消失。

简化代码：

在处理 NumPy 数组时，使用它们的meanandsum方法，避免使用 Pythonsum函数：

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函数，因为

  • 它比 NumPysum方法慢，并且
  • 它不像 NumPysum方法那样通用。（它不允许您指定求和的轴。我们设法sum通过一次调用 NumPysum或mean方法来消除对 Python 的两次调用。）

消除三重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

解释：

如果您查看 的等式 5 更新规则H，首先注意到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)
于 2013-05-15T20:38:12.527 回答