推荐系统之矩阵分解，梯度下降

为什么要在推荐系统中使用矩阵分解？

如果从通信原理的视角来看机器学习的话，可以将机器学习看做一门研究如何从含有冗余的数据中提取出对我们有用的数据最本质的特征的学科。传统的推荐方法通过对User-Item评分矩阵采用基于近邻的协同过滤来发现与用户具有相似偏好的其他用户，进而产生推荐。但这种方法存在着数据稀疏和信息冗余大的缺陷，即矩阵中评分元素稀疏，以及信息量并非随着向量维度增加而线性增加。

为了解决传统方法存在的以上缺陷，有学者利用数学中的奇异值分解相关概念提出了应用于推荐系统的矩阵分解方法。矩阵分解就是把原来的大评分矩阵，近似的分解成小矩阵的乘积，在实际推荐计算时不再使用大矩阵，而是使用分解得到的两个小矩阵。

具体来说就是，假设用户物品的评分矩阵A是m乘n维，即一共有m个用户，n个物品.通过一套算法转化为两个矩阵U和V,矩阵U的维度是m乘k，矩阵V的维度是n乘k。

该方法基于这样一个假设，用户对一个产品的喜爱程度主要由K个因素（比如价格、外观等）决定，

表示第n个用户对第i个因素的偏好程度，而
表示第x个产品满足第i个因素上的程度，
即为用户n对产品x最终的喜好程度。

例子

有如下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在这里插入图片描述
)
矩阵分解的过程中，将原始的评分矩阵分解成两个矩阵和的乘积：
矩阵P(n,K)表示n个user和K个特征之间的关系矩阵，这K个特征是一个中间变量，矩阵Q(K,m)的转置是矩阵Q(m,K)，矩阵Q(m,K)表示m个item和K个特征之间的关系矩阵，这里的K值是自己控制的，可以使用交叉验证的方法获得最佳的K值。为了得到近似的R(n,m)，必须求出矩阵P和Q，如何求它们呢？
【方法】

1. 首先

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

   如果R(i,j)已知，则R(i,j)的误差平方和为：　
   最终，需要求解所有的非“-”项的损失之和的最小值
   3.使用梯度下降法获得修正的p和q分量：
   　　求解损失函数的负梯度：
   
   根据负梯度的方向更新变量

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

【预测】利用上述的过程，我们可以得到矩阵和，这样便可以为用户 i 对商品 j 进行打分：

手动推导(忽略丑字^_)

注：带^符号的值是预测值

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()

运行结果：

参考：矩阵分解

梯度下降

梯度下降就是给一个函数找最小值，机器学习恰恰就是给损失函数找最小值，值最小就代表着误差最小，也就是情况最好，所以梯度下降让机器看上去训练好了自己，说的文艺一点，就是让机器学会了学习。梯度下降是个啥？在一个漆黑的夜晚，黑的啥也看不见，只能用脚探，一是感受下地面，看看哪陡朝哪走，向下最陡的方向，就是梯度的负方向，梯度下面加个负号，二是走的快，就是步子迈的大，也叫学习率，（步长=学习率），值越小代表你越苟，小步走，怕掉进沟，大步走小心伸到腰，所以步子不能太大也不能太小，选择一个合适的步长刚刚好，合在一起，从现在位置朝梯度的负方向走一个合适的步长就能到一个更接近山底的新位置。（新位置=现在位置-步长*梯度），重复这个过程到下山，就能一步一步最快朝下走。这就是机器学习中的梯度下降算法。

如果需要找的是函数极小点，那么应该从负梯度的方向寻找，该方法称之为梯度下降法。
要搜索极小值C点，在A点必须向x增加方向搜索，此时与A点梯度方向相反；在B点必须向x减小方向搜索，此时与B点梯度方向相反。总之，搜索极小值，必须向负梯度方向搜索。
我们可以用一阶泰勒公式来推倒

梯度下降法-步骤

假设函数y=f(x1,x2,…,xn) 只有一个极小点。
初始给定参数为 X0=(x10,x20,…,xn0)。从这个点如何搜索才能找到原函数的极小值点？
方法：
1.首先设定一个较小的正数η,ε.
2.求当前位置处的各个偏导数：
假设函数y=f(x1,x2,…,xn) 只有一个极小点。
初始给定参数为 X0=(x10,x20,…,xn0)。从这个点如何搜索才能找到原函数的极小值点？
方法：
1.首先设定一个较小的正数η,ε.
2.求当前位置处的各个偏导数：

3.修改当前函数的参数值，公式如下：

4. 如果参数变化量小于ε，退出；否则返回2。
例：任给一个初始出发点，设为x0=-4，利用梯度下降法求函数y=x2/2-2x的极小值。


则x’=1.994时，y取得最小值（极小值），若要找山脊那种较多山尖的最小值，建议多找几个初始点。
η是步长，ε是变化量
梯度下降终止条件：1.迭代次数2.变化量小于规定的ε

Python代码实现

一维问题
f(x)=x*x

```python
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
def f(x):
    return np.power(x, 2)
def d_f_1(x):
    return 2.0 * x
def d_f_2(f, x, delta=1e-4):
    return (f(x+delta) - f(x-delta)) / (2 * delta)
# plot the function
xs = np.arange(-10, 11)
plt.plot(xs, f(xs))
plt.show()
learning_rate = 0.1
max_loop = 30
x_init = 10.0
x = x_init
lr = 0.1
for i in range(max_loop):
    # d_f_x = d_f_1(x)
    d_f_x = d_f_2(f, x)
    x = x - learning_rate * d_f_x
    print(x)
print('initial x =', x_init)
print('arg min f(x) of x =', x)
print('f(x) =', f(x))

运行结果：