一.矩阵分解是什么？

度娘定义：矩阵分解(Matrix Factorization, MF)是将矩阵拆解为数个矩阵的乘积。

狐仙定义：矩阵分解(Matrix Factorization, MF)就是将一个大矩阵分解为两个或多个小矩阵。

二.矩阵分解干什么？

如上图所示，矩阵算法就是将用户和产品矩阵中的数据，分解成两个矩阵（用User矩阵和Item矩阵），两个矩阵相乘得到的结果就是预测评分。当我们要计算第i 个用户对第j 个item的预测评分时），我们就可以用User矩阵的第i行和Item矩阵的第j 列做内积，这个内积的值就是预测评分了。对于某个用户进行推荐时，我们把他的用户向量和所有item向量做内积，然后按内积从大到小排序，取出前K 个item，过滤掉历史item后推荐给用户。
那MF是如何从评分矩阵中分解出User矩阵和Item矩阵的呢？简单来说，MF把User矩阵和Item矩阵作为未知量，用它们表示出每个用户对每个item的预测评分，然后通过最小化预测评分跟已知实际评分的差异，学习出User矩阵和Item矩阵。

也就是说：

我们的目的就是利用已知的评分来预测出未知的评分。

三.矩阵算法推导

1.首先我们要定义一个类似于上图的评分矩阵，用 $R$ 表示，其维度为 $N ? M$ ，也就是 $R$ 为 $N$ 行 $M$ 列矩阵。
然后我们将其分解 $P$ 矩阵与 $Q$ 矩阵，其中 $P$ 矩阵维度为 $N ? K$ ， $Q$ 矩阵维度为 $K ? M$
于是我们可以得到
$R\approx$ $\hat{R}=P*Q$
对于 $P$ ， $Q$ 矩阵的解释，直观上， $P$ 矩阵是 $N$ 个用户对 $K$ 个主题的关系， $Q$ 矩阵是 $K$ 个主题跟 $M$ 个物品的关系，至于 $K$ 个主题具体是什么，在算法里面 $K$ 是一个参数，需要调节的，通常 $10\sim100$ 之间。

2.对于 $\hat{R}$ 矩阵：
$\hat{r}_{ij}=p_{i} ^{T}q_{j}=\sum_{k=1}^{K}p_{ik}q_{kj}$
$\hat{R}$ 与 $R$ 的维度相同，其中 $\hat{r}_{ij}$ 是 $\hat{R}$ 第 $i$ 行第 $j$ 列的元素值。

3.求损失函数并更新变量：
使用原始的评分矩阵 $R$ 与重新构建的评分矩阵 $\hat{R}$ 之间的误差的平方作为损失函数，即：
$e_{ij}^{2}=(r_{ij}-\hat{r}_{ij})^{2}=(r_{ij}-\sum_{k=1}^{K}p_{ik}q_{kj})^{2}$
通过梯度下降法，更新变量：
求导： $\frac{?}{?_{p{ik}}}e_{ij}^{2}=-2(r_{ij}-\sum_{k=1}^{K}p_{ik}q_{kj})q_{kj}=-2e_{ij}q_{kj}$ $\frac{?}{?_{q{kj}}}e_{ij}^{2}=-2(r_{ij}-\sum_{k=1}^{K}p_{ik}q_{kj})p_{ik}=-2e_{ij}p_{ik}$
根据负梯度的方向更新变量： $p_{ik}'=p_{ik}-α\frac{?}{?{p_{ik}}}e_{ij}^{2}=p_{ik}+2αe_{ij}q_{kj}$ $q_{kj}'=q_{kj}-α\frac{?}{?{q_{kj}}}e_{ij}^{2}=q_{kj}+2αe_{ij}p_{ik}$

4.在损失函数中加入正则化惩罚项：通常在求解的过程中，为了能够有较好的泛化能力，会在损失函数中加入正则项，以对参数进行约束。加入正则项后的计算过程如下： $E_{ij}^{2}=(r_{ij}-\sum_{k=1}^{K}p_{ik}q_{kj})^{2}+\frac{β}{2}\sum_{k=1}^{K}(p_{ik}^{2}+q_{kj}^{2})$ 通过梯度下降法，更新变量：求导： $\frac{?}{?{p_{ik}}}E_{ij}^{2}=-2(r_{ij}-\sum_{k=1}^{K}p_{ik}q_{kj})q_{kj}+βp_{ik}=-2e_{ij}q_{kj}+βp_{ik}$ $\frac{?}{?{q_{kj}}}E_{ij}^{2}=-2(r_{ij}-\sum_{k=1}^{K}p_{ik}q_{kj})p_{ik}+βq_{kj}=-2e_{ij}p_{ik}+βq_{kj}$ 根据负梯度的方向更新变量： $p_{ik}'=p_{ik}-α(\frac{?}{?{p_{ik}}}e_{ij}^{2}+βp_{ik})=p_{ik}+α(2e_{ij}q_{kj}-βp_{ik})$ $q_{kj}'=q_{kj}-α(\frac{?}{?{q_{kj}}}e_{ij}^{2}+βq_{kj})=q_{kj}+α(2e_{ij}p_{ik}-βq_{kj})$

算法终止：每次更新完 $\hat{R}$ 后，计算一次 $l o s s$ 值，若 $l o s s$ 值非常小或者到达最大迭代次数，结束算法。于是就得到了我们最终的预测矩阵 $\hat{R}$ 。

四.Python代码实现

import numpy as np
import math
import matplotlib.pyplot as plt

#定义矩阵分解函数
def Matrix_decomposition(R,P,Q,N,M,K,alpha=0.0002,beta=0.02):
    Q = Q.T #Q 矩阵转置
    loss_list = [] #存储每次迭代计算的 loss 值
    for step in range(5000):
        #更新 R^
        for i in range(N):
            for j in range(M):
                if R[i][j] != 0:
                    #计算损失函数
                    error = R[i][j]
                    for k in range(K):
                        error -= P[i][k]*Q[k][j]
                    #优化 P,Q 矩阵的元素
                    for k in range(K):
                        P[i][k] = P[i][k] + alpha*(2*error*Q[k][j]-beta*P[i][k])
                        Q[k][j] = Q[k][j] + alpha*(2*error*P[i][k]-beta*Q[k][j])
                    
        loss = 0.0
        #计算每一次迭代后的 loss 大小，就是原来 R 矩阵里面每个非缺失值跟预测值的平方损失
        for i in range(N):
            for j in range(M):
                if R[i][j] != 0:
                    #计算 loss 公式加号的左边
                    data = 0
                    for k in range(K):
                        data = data + P[i][k]*Q[k][j]
                    loss = loss + math.pow(R[i][j]-data,2)
                    #得到完整 loss 值
                    for k in range(K):
                        loss = loss + beta/2*(P[i][k]*P[i][k]+Q[k][j]*Q[k][j])
                    loss_list.append(loss)
        plt.scatter(step,loss)
        #输出 loss 值
        if (step+1) % 1000 == 0:
            print("loss={:}".format(loss))
        #判断
        if loss < 0.001:
            print(loss)
            break
    plt.show()    
    return P,Q
        
if __name__ == "__main__":
    N = 5
    M = 4
    K = 5
    R = np.array([[5,3,0,1],
                [4,0,0,1],
                [1,1,0,5],
                [1,0,0,4],
                [0,1,5,4]]) #N=5,M=4
    print("初始评分矩阵：")
    print(R)
    #定义 P 和 Q 矩阵
    P = np.random.rand(N,K) #N=5,K=2
    Q = np.random.rand(M,K) #M=4,K=2
    
    print("开始矩阵分解：")
    P,Q = Matrix_decomposition(R,P,Q,N,M,K)
    print("矩阵分解结束。")
    print("得到的预测矩阵：")
    print(np.dot(P,Q))

运行结果分析：
初始评分矩阵：
[[5 3 0 1]
 [4 0 0 1]
 [1 1 0 5]
 [1 0 0 4]
 [0 1 5 4]]
开始矩阵分解：
loss=8.727735279273917
loss=1.5983412898789944
loss=1.2664582185054736
loss=1.2241666331627166
loss=1.2166150584503614

矩阵分解结束。
得到的预测矩阵：
[[4.97903684 2.97131217 3.41300234 1.00412214]
 [3.97892547 0.6736204  3.33839296 1.00341136]
 [1.00662457 0.97898601 3.36790504 4.97254172]
 [0.99649186 0.5035927  2.9660996  3.98258083]
 [3.95311908 1.02446502 4.9821697  3.98566116]]