一、矩阵分解回想

在博文推荐算法——基于矩阵分解的推荐算法中，提到了将用户-商品矩阵进行分解。从而实现对未打分项进行打分。

矩阵分解是指将一个矩阵分解成两个或者多个矩阵的乘积。对于上述的用户-商品矩阵(评分矩阵)，记为$V_{m imes n}$。能够将其分解成两个或者多个矩阵的乘积，如果分解成两个矩阵$W_{m imes k}$和$H_{k imes n}$。我们要使得矩阵$W_{m imes k}$和$H_{k imes n}$的乘积能够还原原始的矩阵$V_{m imes n}$：

$$V_{m imes n}approx W_{m imes k} imes H_{k imes n}=hat{V}_{m imes n}$$

当中，矩阵$W_{m imes k}$表示的是$m$个用户与$k$个主题之间的关系，而矩阵$H_{k imes n}$表示的是$k$个主题与$n$个商品之间的关系。

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

$$W_{m imes k}geqslant 0$$

$$H_{k imes n}geqslant 0$$

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

二、非负矩阵分解

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

上面简介了非负矩阵分解的基本含义。简单来讲，非负矩阵分解是在矩阵分解的基础上对分解完毕的矩阵加上非负的限制条件。即对于用户-商品矩阵$V_{m imes n}$，找到两个矩阵$W_{m imes k}$和$H_{k imes n}$，使得：

$$V_{m imes n}approx W_{m imes k} imes H_{k imes n}=hat{V}_{m imes n}$$

同一时候要求：

$$W_{m imes k}geqslant 0$$

$$H_{k imes n}geqslant 0$$

2.2、损失函数

为了能够定量的比較矩阵$V_{m imes n}$和矩阵$hat{V}_{m imes n}$的近似程度。在參考文献1中作者提出了两种损失函数的定义方式：

• 平方距离

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

• KL散度

$$Dleft ( Aparallel B ight )=sum_{i,j}left ( A_{i,j}logfrac{A_{i,j}}{B_{i,j}}-A_{i,j}+B_{i,j} ight )$$

在KL散度的定义中，$Dleft ( Aparallel B ight )geqslant 0$。当且仅当$A=B$时取得等号。

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

求解例如以下的最小化问题：

• $$egin{matrix} minimize; left | V-WH ight |^2\ s.t.; Wgeqslant 0,Hgeqslant 0 end{matrix}$$

• $$egin{matrix} minimize; Dleft ( Vparallel WH ight )\ s.t.; Wgeqslant 0,Hgeqslant 0 end{matrix}$$

2.3、优化问题的求解

在參考文献1中，作者提出了乘法更新规则(multiplicative update rules)，详细的操作例如以下：

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

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

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

对于KL散度的损失函数：

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

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

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

平方损失函数能够写成：

$$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} ight ) ight ]^2$$

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

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

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

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

即为：

$$H_{k,j}=H_{k,j}+eta _{k,j}left [ left ( W^TV ight )_{k,j}-left ( W^TWH ight )_{k,j} ight ]$$

令$eta _{k,j}=frac{H_{k,j}}{left ( W^TWH ight )_{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("	")
        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 实现