推荐算法——非负矩阵分解(NMF)

一、矩阵分解回顾

在博文推荐算法——基于矩阵分解的推荐算法中，提到了将用户-商品矩阵进行分解，从而实现对未打分项进行打分。矩阵分解是指将一个矩阵分解成两个或者多个矩阵的乘积。对于上述的用户-商品矩阵(评分矩阵)，记为Vm×n_V_{m\times n}，可以将其分解成两个或者多个矩阵的乘积，假设分解成两个矩阵_Wm×k_W_{m\times k}和_Hk×n_H_{k\times n}，我们要使得矩阵_Wm×k_W_{m\times k}和_Hk×n_H_{k\times n}的乘积能够还原原始的矩阵_Vm×_n_V_{m\times n}：

Vm×n≈Wm×k×Hk×n=V^m×n

V_{m\times n}\approx W_{m\times k}\times H_{k\times n}=\hat{V}_{m\times n}

其中，矩阵Wm×k_W_{m\times k}表示的是_m_m个用户与_k_k个主题之间的关系，而矩阵_Hk×_n_H_{k\times n}表示的是_k_k个主题与_n_n个商品之间的关系。

通常在用户对商品进行打分的过程中，打分是非负的，这就要求：

Wm×k⩾0

W_{m\times k}\geqslant 0

Hk×n⩾0

H_{k\times n}\geqslant 0

这便是非负矩阵分解(Non-negtive Matrix Factorization, NMF)的来源。

二、非负矩阵分解

2.1、非负矩阵分解的形式化定义

上面简单介绍了非负矩阵分解的基本含义，简单来讲，非负矩阵分解是在矩阵分解的基础上对分解完成的矩阵加上非负的限制条件，即对于用户-商品矩阵Vm×n_V_{m\times n}，找到两个矩阵_Wm×k_W_{m\times k}和_Hk×_n_H_{k\times n}，使得：

Vm×n≈Wm×k×Hk×n=V^m×n

V_{m\times n}\approx W_{m\times k}\times H_{k\times n}=\hat{V}_{m\times n}

同时要求：

Wm×k⩾0

W_{m\times k}\geqslant 0

Hk×n⩾0

H_{k\times n}\geqslant 0

2.2、损失函数

为了能够定量的比较矩阵Vm×n_V_{m\times n}和矩阵_V^m×n\hat{V}_{m\times n}的近似程度，在参考文献1中作者提出了两种损失函数的定义方式：

• 平方距离

∥A−B∥2=∑i,j(Ai,j−Bi,j)2

\left | A-B \right |^2=\sum_{i,j}\left ( A_{i,j}-B_{i,j} \right )^2

• KL散度

D(A∥B)=∑i,j(Ai,jlogAi,jBi,j−Ai,j+Bi,j)

D\left ( A\parallel B \right )=\sum_{i,j}\left ( A_{i,j}log\frac{A_{i,j}}{B_{i,j}}-A_{i,j}+B_{i,j} \right )

在KL散度的定义中，D(A∥B)⩾0D\left ( A\parallel B \right )\geqslant 0，当且仅当A=_B_A=B时取得等号。

当定义好损失函数后，需要求解的问题就变成了如下的形式，对应于不同的损失函数：

求解如下的最小化问题：

• minimize∥V−WH∥2s.t.W⩾0,H⩾0 \begin{matrix} minimize\; \left | V-WH \right |^2\ s.t.\; W\geqslant 0,H\geqslant 0 \end{matrix}minimizeD(V∥WH)s.t.W⩾0,H⩾0 \begin{matrix} minimize\; D\left ( V\parallel WH \right )\ s.t.\; W\geqslant 0,H\geqslant 0 \end{matrix}

2.3、优化问题的求解

在参考文献1中，作者提出了乘法更新规则(multiplicative update rules)，具体的操作如下：

对于平方距离的损失函数：

Wi,k=Wi,k(VHT)i,k(WHHT)i,k

W_{i,k}=W_{i,k}\frac{\left ( VH^T \right )_{i,k}}{\left ( WHH^T \right )_{i,k}}

Hk,j=Hk,j(WTV)k,j(WTWH)k,j

H_{k,j}=H_{k,j}\frac{\left ( W^TV \right )_{k,j}}{\left ( W^TWH \right )_{k,j}}

对于KL散度的损失函数：

Wi,k=Wi,k∑uHk,uVi,u/(WH)i,u∑vHk,v

W_{i,k}=W_{i,k}\frac{\sum_{u}H_{k,u}V_{i,u}/\left ( WH \right )_{i,u}}{\sum_{v}H_{k,v}}

Hk,j=Hk,j∑uWu,kVu,j/(WH)u,j)∑vWv,k

H_{k,j}=H_{k,j}\frac{\sum_{u}W_{u,k}V_{u,j}/\left ( WH \right )_{u,j})}{\sum_{v}W_{v,k}}

上述的乘法规则主要是为了在计算的过程中保证非负，而基于梯度下降的方法中，加减运算无法保证非负，其实上述的乘法更新规则与基于梯度下降的算法是等价的，下面以平方距离为损失函数说明上述过程的等价性：

平方损失函数可以写成：

l=∑i=1m∑j=1nVi,j−(∑k=1rWi,k⋅Hk,j)2

l=\sum_{i=1}^{m}\sum_{j=1}^{n}\left V_{i,j}-\left ( \sum_{k=1}^{r}W_{i,k}\cdot H_{k,j} \right ) \right ^2

使用损失函数对Hk,_j_H_{k,j}求偏导数：

∂l∂Hk,j=∑i=1m∑j=1n2(Vi,j−(∑k=1rWi,k⋅Hk,j))⋅(−Wi,k)=−2(WTV)k,j−(WTWH)k,j

\begin{align*} \frac{\partial l}{\partial H_{k,j}}&= \sum_{i=1}^{m}\sum_{j=1}^{n}\left 2\left ( V_{i,j}-\left ( \sum_{k=1}^{r}W_{i,k}\cdot H_{k,j} \right ) \right )\cdot \left ( -W_{i,k} \right ) \right \ &= -2\left \left ( W^TV \right )_{k,j}-\left ( W^TWH \right )_{k,j} \right \end{align*}

则按照梯度下降法的思路：

Hk,j=Hk,j−ηk,j∂l∂Hk,j

H_{k,j}=H_{k,j}-\eta _{k,j}\frac{\partial l}{\partial H_{k,j}}

即为：

Hk,j=Hk,j+ηk,j(WTV)k,j−(WTWH)k,j

H_{k,j}=H_{k,j}+\eta _{k,j}\left \left ( W^TV \right )_{k,j}-\left ( W^TWH \right )_{k,j} \right

令ηk,j=Hk,j(WTWH)k,j\eta _{k,j}=\frac{H_{k,j}}{\left ( W^TWH \right )_{k,j}}，即可以得到上述的乘法更新规则的形式。

2.4、非负矩阵分解的实现

对于如下的矩阵：

通过非负矩阵分解，得到如下的两个矩阵：

对原始矩阵的还原为：

实现的代码

#!/bin/python

from numpy import * 

def load_data(file_path):
    f = open(file_path)
    V = []
    for line in f.readlines():
        lines = line.strip().split("\t")
        data = []
        for x in lines:
            data.append(float(x))
        V.append(data)
    return mat(V)

def train(V, r, k, e):
    m, n = shape(V)
    W = mat(random.random((m, r)))
    H = mat(random.random((r, n)))

    for x in xrange(k):
        #error 
        V_pre = W * H
        E = V - V_pre
        #print E
        err = 0.0
        for i in xrange(m):
            for j in xrange(n):
                err += E[i,j] * E[i,j]
        print err

        if err < e:
            break

        a = W.T * V
        b = W.T * W * H
        #c = V * H.T
        #d = W * H * H.T
        for i_1 in xrange(r):
            for j_1 in xrange(n):
                if b[i_1,j_1] != 0:
                    H[i_1,j_1] = H[i_1,j_1] * a[i_1,j_1] / b[i_1,j_1]

        c = V * H.T
        d = W * H * H.T
        for i_2 in xrange(m):
            for j_2 in xrange(r):
                if d[i_2, j_2] != 0:
                    W[i_2,j_2] = W[i_2,j_2] * c[i_2,j_2] / d[i_2, j_2]

    return W,H 


if __name__ == "__main__":
    #file_path = "./data_nmf"
    file_path = "./data1"

    V = load_data(file_path)
    W, H = train(V, 2, 100, 1e-5 )

    print V
    print W
    print H
    print W * H

收敛曲线如下图所示：

'''
Date:20160411
@author: zhaozhiyong
'''

from pylab import *
from numpy import *

data = []

f = open("result_nmf")
for line in f.readlines():
    lines = line.strip()
    data.append(lines)

n = len(data)
x = range(n)
plot(x, data, color='r',linewidth=3)
plt.title('Convergence curve')
plt.xlabel('generation')
plt.ylabel('loss')
show()

参考文献

• Algorithm for Non-negative Matrix Factorization
• 白话NMF（Non-negative Matrix Factorization）——Matlab 实现