矩阵分解

MF(矩阵分解)算法概述

在之前我们介绍了协调过滤算法，但由于协调过滤存在以下问题：

• 协调过滤处理稀疏矩阵的能力较弱
• 协调过滤中，相似度矩阵维护难度大(物品和用户数量都比较大)
• 协调过滤算法的头部效应强

为了解决这些问题，我们提出了矩阵分解算法。其解决问题的大致思路是：希望为每个用户和物品生成一个隐向量，将用户和物品定位到隐向量的表示空间上，距离相近的用户和物品表面兴趣特点接近，应该把距离相近的视频推荐给目标用户。下面举例说明：

如上图所示，左边是一个用户物品的共现矩阵，我们利用矩阵分解把矩阵分成了$3*5(令k=5)$，和$5(k)*3$的两个矩阵相乘，我们可以从图中看出，两个矩阵的行和列的标签是一样的，我们用隐向量给每个用户和每首音乐都打上标签，如果这个人喜欢小清新的歌，那么就应该给他推荐符合小清新风格的音乐。
然而在实际应用中要注意：

• 这里的隐向量是不可解释的，并不是如上图所示可以给用户和物品打上具体的标签，这是需要模型自己去学习的。
• 隐向量个数k决定了隐向量的表达能力的强弱，k越大，表达的信息就越强，用户的兴趣和物品的分类划分的越具体。
  
  所以基于用户矩阵U和物品矩阵V，用户u对物品i的预测评分就可以表示为：

$${\hat{r}}_{ui}=q_{i}^{T}p_{u}$$

其中$p_{u}$是用户u在用户矩阵U中对应的行向量，$q_{i}$是物品i在物品矩阵V中对应的列向量。

矩阵分解的求解

特征值分解

首先我们引入特征值，特征向量的概念

$$Av=\lambda v$$

其中A是方阵，v是A的特征向量，$\lambda$是A的特征值。
则A可以分解为

$$A=Q\Sigma Q^{-1}$$

其中Q是该矩阵的特征向量组成的矩阵，$\Sigma$是一个对角阵，每个对角线上元素就是该矩阵的特征值。
但由于其只适用于方阵，所以这种方法便被淘汰了。

SVD 奇异值分解

其计算步骤为(了解即可)：

从上图可知，SVD分解的计算过程比较复杂，计算复杂度达到$O(mn^{2})$级别，这对于商品数量动辄上百万的场景显然是不可接受的。

BaseSVD

基于我们上文的评分函数构造出损失函数:

对该函数进行梯度下降，最小化损失函数，即可求解出用户矩阵和物品矩阵。
详细求解过程参考下文。

LFM

由于不同用户的打分体系不同，有的人认为打3分已经是很低的分数了，二而有的人认为打1分才是比较差的评价，再加上不同物品的衡量标准也有所区别，比如电子产品和日用品的评分标准就回会很大，为了消除用户和物品打分的偏差，我们在君子分解时加入用户和物品的偏差向量,则我们的预测评分变成了:

$$r_{ui}=u+b_{i}+b_{u}+q_{i}^{T}p_{u}$$

其中u是全局偏差常数，可使用评分的全局平均分，$b_{i}$是物品偏差系数，可使用物品i收到的所以评分均值，$b_{u}$是用户偏差系数，可使用用户u给出的所有评分均值。
因此我们的目标函数也变成了:

SVD++

由于物品之间可能存在某种关联性，这种关联性可能会影响物品的评分，(比如之前看过一部电影，就可能对他的续作也产生不错的印象，从而给打较高的分)，所以SVD++就引入了用户评过分的历史物品，将这些物品间的相似度也考虑了进去。
期预测函数如下图所示：

1.问题引入

 有如下R(5,4)矩阵:("-"表示用户没有打分)
其中打分矩阵R(n,m)是n行和m列，n表示user个数,m表示item个数。

 想要根据目前的矩阵R对未打分的商品进行评分的预测。为了解决这个问题，就用到了矩阵分解的方法。

2.问题分析

2.1构造损失函数

 为了求出未打分的值，可以将矩阵R(n,m)分解为P(n,k)*Q(K,m),所以可以得到一个预测值$\hat{R}=P*Q$,来预测矩阵R，那么此时我们的问题就转换成了如何令$\hat{R}$与R最为接近，由此我们引入了损失函数的概念，即让$\hat{R}-R$变小，当小到一定程度后我们即可认为$\hat{R}\approx{R}$,为了简化计算(令该式子始终大于0)，我们给其加入了一个平方变为

$$Loss(P,Q)=\frac{1}{2}(\hat{R-R})^2=\frac{1}{2}(\sum_{i=1}^{m}\sum_{j=1}^{n}({R_{ij}-\sum_{s=1}^{k}P_{is}Q_{js}}))^{2}$$

2.1梯度下降

 梯度下降法的核心思想即是沿函数的导数方向的不断改变自身参数，从到最终能到达最低点。梯度下降法的一些注意事项在多元线性回归中做了一些介绍，有兴趣可以查看。
 故运用梯度下降法,先求其偏导数：

$\frac{\partial}{\partial{P_{ik}}}L(P,Q)$
$=\frac{1}{2}(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})*2*(-1)*(\frac{\partial{\sum_{k=1}^{k}P_{ik}Q_{kj}})}{\partial{P_{ik}}})$
$=-(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})Q_{kj}$
$\frac{\partial}{\partial{Q_{kj}}}L(P,Q)$
$=\frac{1}{2}(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})*2*(-1)*(\frac{\partial{\sum_{k=1}^{k}P_{ik}Q_{kj}})}{\partial{P_{kj}}})$
$=-(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})P_{ik}$

应用梯度下降法不断修改当前参数值：
$P_{ik}=P_{ik}-\alpha*(-(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})Q_{kj})$
$=P_{ik}+\alpha(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})Q_{kj})$

$Q_{kj}=P_{ik}-\alpha*(-(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})P_{ik})$
$=P_{ik}+\alpha(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})P_{ik})$
($\alpha$为learning rate ，即学习率)

2.3正则化

 当训练的数据不够时或特征数量大于样本数量时，就会出现过拟合的情况，我们为了解决过拟合，就会采用正则化的手段或者减少一些非必要的样本特征。常用的有L1范数和L2范数。
①LP范数不是一个范数，而是一组范数
$||x||_{p}=(\sum_{i}^{n}|x_{i}|^{p})^{\frac{1}{p}}$
②p=1时，几位L1范数
$||x_{1}||=\sum_{i=1}^{n}|x_{i}|$
②L2范数,平方再开方，表示向量的距离。
$||x_{2}||=(\sum_{i=1}^{n}|x_{i}|^{2})^{\frac{1}{2}}$
L1正则化：

$$J({\theta})=\frac{1}{2m}[\sum_{i=1}^{m}{(y^{i}-h_{\theta}(x^{i}))^{2}+\lambda\sum_{j=1}{n}|\theta_{j}|}]$$

L2正则化：

$$J({\theta})=\frac{1}{2m}[\sum_{i=1}^{m}{(y^{i}-h_{\theta}(x^{i}))^{2}+\lambda\sum_{j=1}{n}\theta_{j}^{2}}]$$

式子中m为样本数量，n为特征个数。

3.代码实现

import numpy as np
import matplotlib.pyplot as plt

R=np.array([[5,3,0,1],
            [4,0,0,1],
            [1,1,0,5],
            [1,0,0,4],
            [0,1,5,4]])

M=R.shape[0]
N=R.shape[1]
K=2
P=np.random.rand(M,K)
Q=np.random.rand(K,N)
print(P)
print(Q)

[[0.65360868 0.08726216]
 [0.69762591 0.03076026]
 [0.88861613 0.75595254]
 [0.31559948 0.49104854]
 [0.33077934 0.40914307]]
[[0.68124148 0.49563369 0.14058849 0.73546283]
 [0.31494592 0.86798718 0.48578754 0.95628464]]

$$Loss=\frac{1}{2}(R-PQ)^{2}=\frac{1}{2}(\sum_{i=1}^{m}\sum_{j=1}^{n}({R_{ij}-\sum_{s=1}^{k}P_{is}Q_{js}}))^{2}$$

#损失函数
def cost(R,P,Q):
    e=0
    for i in range(R.shape[0]):
        for j in range(R.shape[1]):
            if(R[i,j]>0):
                e+=(R[i,j]-np.dot(P[i,:],Q[:,j]))**2
                
    return e/2

L2正则化：

$$Loss=\frac{1}{2}(R-PQ)^{2}=\frac{1}{2}(\sum_{i=1}^{m}\sum_{j=1}^{n}){R_{ij}-\sum_{s=1}^{k}P_{is}Q_{js}})^{2}+\frac{\lambda}{2}\sum_{s=1}{k}({P_{ik}^{2}}+Q_{kj}^{2})$$

#损失函数正则
def cost_re(R,P,Q,lambda1):
    e=0
    for i in range(R.shape[0]):
        for j in range(R.shape[1]):
            if(R[i][j]>0):
                e+=(R[i,j]-np.dot(P[i,:],Q[:,j]))**2
                for k in range(P.shape[1]):
                    e+=lambda1*(P[i,k]**2+Q[k,j]**2)/2    
    return e/2

对损失函数求偏导：
$\frac{\partial}{\partial{P_{ik}}}L(P,Q)$
$=\frac{1}{2}(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})*2*(-1)*(\frac{\partial{\sum_{k=1}^{k}P_{ik}Q_{kj}})}{\partial{P_{ik}}})$
$=-(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})Q_{kj}$
$\frac{\partial}{\partial{Q_{kj}}}L(P,Q)$
$=\frac{1}{2}(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})*2*(-1)*(\frac{\partial{\sum_{k=1}^{k}P_{ik}Q_{kj}})}{\partial{P_{kj}}})$
$=-(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})P_{ik}$
梯度下降法：
$P_{ik}=P_{ik}-\alpha*(-(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})Q_{kj})$
$=P_{ik}+\alpha(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})Q_{kj})$

$Q_{kj}=P_{ik}-\alpha*(-(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})P_{ik})$
$=P_{ik}+\alpha(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})P_{ik})$

#梯度下降
def grad(R,P,Q,lr,epochs):
    costList=[]
    for s in range(epochs+1):
        for i in range(R.shape[0]):
            for j in range(R.shape[1]):
                if(R[i][j]>0):
                    e=R[i,j]-np.dot(P[i,:],Q[:,j])
                    for k in range(P.shape[1]):
                        grad_p=e*Q[k][j]
                        grad_q=e*P[i][k]

                        P[i][k]=P[i][k]+lr*grad_p
                        Q[k][j]=Q[k][j]+lr*grad_q
            
        if s%50==0:
            e=cost(R,P,Q)
            costList.append(e) 
            #print(e)
    return P,Q,costList
lr=0.0001
epochs=10000
p,q,costList=grad(R,P,Q,lr,epochs)
print(np.dot(p,q))

[[5.06116901 2.87426131 2.43375669 1.00254423]
 [3.96425112 2.26334177 2.08974137 0.96504113]
 [1.0366021  0.90360134 5.30277467 4.91333664]
 [0.99043153 0.81279639 4.2952649  3.93863408]
 [1.64541544 1.19616657 4.78417292 4.23883696]]

x=np.linspace(0,10000,201)
plt.plot(x,costList,'r')
plt.xlabel('epochs')
plt.ylabel('cost')
plt.show()

​

加入正则项后：
$P_{ik}=P_{ik}-\alpha*(-(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})P_{ik})+\lambda{P_{ik}}$
$=P_{ik}+\alpha(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})P_{ik})-\lambda{P_{ik}}$
$Q_{kj}=Q_{kj}-\alpha*(-(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})Q_{kj})+\lambda{Q_{kj}}$
$=P_{ik}+\alpha(R_{ij}-\sum_{k=1}^{k}P_{ik}Q_{kj})Q_{kj})-\lambda{Q_{kj}}$

#梯度下降（正则）
def grad_re(R,P,Q,lr,epochs,lambda1):
    costList=[]
    for s in range(epochs+1):
        for i in range(R.shape[0]):
            for j in range(R.shape[1]):
                if(R[i][j]>0):
                    e=R[i,j]-np.dot(P[i,:],Q[:,j])
                    for k in range(P.shape[1]):
                        grad_p=e*Q[k][j]#求梯度
                        grad_q=e*P[i][k]

                        #对PQ同时进行梯度下降，改变其值
                        P[i][k]=P[i][k]+lr*grad_p-lambda1*P[i][k]
                        Q[k][j]=Q[k][j]+lr*grad_q-lambda1*Q[k][j]
            
        if s%50==0:#求损失
            e=cost_re(R,P,Q,lambda1)
            costList.append(e) 
            #print(e)
    return P,Q,costList
lambda1=0.0001
p,q,costList=grad_re(R,P,Q,0.003,epochs,lambda1)
print(np.dot(p,q))

[[4.92262848 2.9460198  2.02962105 1.00548811]
 [3.9420148  2.37520018 1.85600337 0.99804788]
 [1.00407789 0.99148638 6.03037773 4.90035828]
 [0.99660486 0.90775439 4.8875762  3.94603309]
 [1.15092516 1.00009806 4.95103649 3.97741499]]

x=np.linspace(0,10000,201)
plt.plot(x,costList,'r')
plt.xlabel('epochs')
plt.ylabel('cost')
plt.show()

​

 通过观察代价函数随迭代步数的变化图可以看出：进行正则化后的代价函数的图像比没有进行正则化的图像较为圆滑了些许，可以得出正则化确实是一种解决过拟合问题的一种手段。
 此外，除了正则化以外，还可以通过增加样本数据集的数量或者减少样本特征，抛弃一些不必要的样本特征来避免过拟合。