非负矩阵分解无法收敛

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.我尝试了两种方式,并且它不会改变实现的行为.

我已经多次检查了我对文件的实现情况,而且我没有看到我做错了什么.任何人都可以对这个问题有所了解吗？

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

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

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

这两个更新都减少了分歧.

试一试:WH = W.dot(H)在计算之前坚持下去W_coeff,交替效果消失.

简化代码:

处理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方法那样通用.(它不允许你指定要求和的轴.我们设法sum通过一次调用NumPy sum或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默认情况下以元素方式乘以数组.

因此,我们可以表达整个更新规则Has

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)