矩阵分解之梯度下降及其python代码实现

分解矩阵

矩阵分解是矩阵分析的重要内容，这种技术就是将一个矩阵分解成若干个矩阵的乘积，通常是两个矩阵或者是三个矩阵的乘积。

例如：有如下R(5,4)的打分矩阵：（“-”表示用户没有打分）

其中打分矩阵R(n,m)是n行和m列，n表示user个数，m行表示item个数

那么，如何根据目前的矩阵R（5,4）如何对未打分的商品进行评分的预测（如何得到分值为0的用户的打分值）？

矩阵分解的思想可以解决这个问题，其实这种思想可以看作是有监督的机器学习问题（回归问题）。

矩阵R可以近似表示为P与Q的乘积：R（n,m）≈ P(n,K)*Q(K,m)

矩阵分解的过程中，将原始的评分矩阵R（n,m)分解成P(n,K)和Q(K,m)两个矩阵和的乘积：

矩阵P(m,k)表示n个user和K个特征之间的关系矩阵，这K个特征是一个中间变量，矩阵Q(n,k)的转置是矩阵Q(m,K)，矩阵Q(m,K)表示m个item和K个特征之间的关系矩阵，这里的K值是自己控制的，可以使用交叉验证的方法获得最佳的K值。为了得到近似的R(n,m)，必须求出矩阵P(m,k)和Q(m,K)，如何求它们呢？

构造损失函数(loss function)

损失函数(loss function)或代价函数(cost function)来度量预测错误程度

也就是说，对于一个机器学习算法，如何评价一个算法是否是比较好的算法，需要提前定义一个损失函数，来判断这个算法是否是最优的，而下面不断的优化求梯度下降，使得损失函数最小，也是为了让一个算法达到意义上的最优。

下面是我们的假设目标函数的公式，我们需要通过costfunction求得参数θ0 和 θ1，其中m为样本个数

梯度下降

梯度下降法,又称最速下降法

1847年由著名的数学家柯西Cauchy给出。
基本思想:假设我们爬山，如果想最快的上到山顶，那么我们应该从山势最陡的地方上山。也就是山势变化最快的地方上山。

同样，如果从任意一点出发，需要最快搜索到函数最大值，那么我们也应该从函数变化最快的方向搜索。函数变化最快的方向是什么呢？
------函数的梯度。

如果函数为一元函数，梯度就是该函数的导数。
如果需要找的是函数极小点，那么应该从负梯度的方向寻找，该方法称之为梯度下降法。
要搜索极小值C点，在A点必须向x增加方向搜索，此时与A点梯度方向相反；在B点必须向x减小方向搜索，此时与B点梯度方向相反。总之，搜索极小值，必须向负梯度方向搜索。

针对文章开始提到的问题我们可以用梯度下降的方法解决

1. 首先令
2. 损失函数：使用原始的评分矩阵与重新构建的评分矩阵之间的误差的平方作为损失函数

即：如果R(i,j)已知，则R(i,j)的误差平方和为：

最终，需要求解所有的非“-”项的损失之和的最小值：
　　

3. 使用梯度下降法获得修正的p和q分量：

• 求解损失函数的负梯度：
  
• 根据负梯度的方向更新变量：
  

4. 不停迭代直到算法最终收敛（直到sum(e^2) <=阈值）

加入正则项的损失函数求解

1. 首先令
2. 通常在求解的过程中，为了能够有较好的泛化能力，会在损失函数中加入正则项，以对参数进行约束，加入
   正则的损失函数为：
   
   即：
   
3. 使用梯度下降法获得修正的p和q分量：

• 求解损失函数的负梯度：
  

• 根据负梯度的方向更新变量：
  

4. 不停迭代直到算法最终收敛（直到sum(e^2) <=阈值）

Python代码实现

# !/usr/bin/env python
# encoding: utf-8
__author__ = 'Scarlett'
#矩阵分解在打分预估系统中得到了成熟的发展和应用
# from pylab import *
import matplotlib.pyplot as plt
from math import pow
import numpy


def matrix_factorization(R,P,Q,K,steps=5000,alpha=0.0002,beta=0.02):
    Q=Q.T  # .T操作表示矩阵的转置
    result=[]
    for step in range(steps):
        for i in range(len(R)):
            for j in range(len(R[i])):
                if R[i][j]>0:
                    eij=R[i][j]-numpy.dot(P[i,:],Q[:,j]) # .dot(P,Q) 表示矩阵内积
                    for k in range(K):
                        P[i][k]=P[i][k]+alpha*(2*eij*Q[k][j]-beta*P[i][k])
                        Q[k][j]=Q[k][j]+alpha*(2*eij*P[i][k]-beta*Q[k][j])
        eR=numpy.dot(P,Q)
        e=0
        for i in range(len(R)):
            for j in range(len(R[i])):
                if R[i][j]>0:
                    e=e+pow(R[i][j]-numpy.dot(P[i,:],Q[:,j]),2)
                    for k in range(K):
                        e=e+(beta/2)*(pow(P[i][k],2)+pow(Q[k][j],2))
        result.append(e)
        if e<0.001:
            break
    return P,Q.T,result

if __name__ == '__main__':
    R=[
        [5,3,0,1],
        [4,0,0,1],
        [1,1,0,5],
        [1,0,0,4],
        [0,1,5,4]
    ]

    R=numpy.array(R)

    N=len(R)
    M=len(R[0])
    K=2

    P=numpy.random.rand(N,K) #随机生成一个 N行 K列的矩阵
    Q=numpy.random.rand(M,K) #随机生成一个 M行 K列的矩阵

    nP,nQ,result=matrix_factorization(R,P,Q,K)
    print("原始的评分矩阵R为：\n",R)
    R_MF=numpy.dot(nP,nQ.T)
    print("经过MF算法填充0处评分值后的评分矩阵R_MF为：\n",R_MF)

#-------------损失函数的收敛曲线图---------------

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

运行结果如图：

原始的评分矩阵R为：
[5 3 0 1]
[4 0 0 1]
[1 1 0 5]
[1 0 0 4]
[0 1 5 4]

经过MF算法填充0处评分值后的评分矩阵R_MF为：
[5.05216756 2.75396636 6.14752912 0.99926854]
[3.92049773 2.13729219 4.86314053 1.00146253]
[1.13246882 0.62163262 3.32199015 4.96708256]
[0.9313801 0.51115749 2.68983381 3.98193247]
[2.8126172 1.53632276 4.8401713 4.01543862]