基于矩阵分解的CF算法实现

所用数据集：链接：https://pan.baidu.com/s/1OLQE7mpefXGRpADyVEkpVQ
提取码：7x5c

一、矩阵分解发展史

1.1 Traditional SVD

通常SVD矩阵分解指的是SVD（奇异值）分解技术，在这我们姑且将其命名为Traditional SVD（传统并经典着）其公式如下：

Traditional SVD分解的形式为3个矩阵相乘，中间矩阵为奇异值矩阵。如果想运用SVD分解的话，有一个前提是要求矩阵是稠密的，即矩阵里的元素要非空，否则就不能运用SVD分解。

很显然我们的数据其实绝大多数情况下都是稀疏的，因此如果要使用Traditional SVD，一般的做法是先用均值或者其他统计学方法来填充矩阵，然后再运用Traditional SVD分解降维，但这样做明显对数据的原始性造成一定影响。

1.2 FunkSVD（LFM）

刚才提到的Traditional SVD首先需要填充矩阵，然后再进行分解降维，同时存在计算复杂度高的问题，因为要分解成3个矩阵，所以后来提出了Funk SVD的方法，它不在将矩阵分解为3个矩阵，而是分解为2个用户-隐含特征，项目-隐含特征的矩阵，Funk SVD也被称为最原始的LFM模型

借鉴线性回归的思想，通过最小化观察数据的平方来寻求最优的用户和项目的隐含向量表示。同时为了避免过度拟合（Overfitting）观测数据，又提出了带有L2正则项的FunkSVD，上公式：

以上两种最优化函数都可以通过梯度下降或者随机梯度下降法来寻求最优解。

1.3 BiasSVD

在FunkSVD提出来之后，出现了很多变形版本，其中一个相对成功的方法是BiasSVD，顾名思义，即带有偏置项的SVD分解：

它基于的假设和Baseline基准预测是一样的，但这里将Baseline的偏置引入到了矩阵分解中

SVD++:

人们后来又提出了改进的BiasSVD，被称为SVD++，该算法是在BiasSVD的基础上添加了用户的隐式反馈信息：

显示反馈指的用户的评分这样的行为，隐式反馈指用户的浏览记录、购买记录、收听记录等。

SVD++是基于这样的假设：在BiasSVD基础上，认为用户对于项目的历史浏览记录、购买记录、收听记录等可以从侧面反映用户的偏好。

二、LFM

LFM也就是前面提到的Funk SVD矩阵分解

用隐语义模型进行协同过滤的目标

揭示隐藏的特征，这些特征能够解释为什么给出对应的预测评分。
这些特征可能是无法用语言描述的，事实上我们并不知道，”玄学“

我们可以认为，用户之所以给电影打出这样的分数，是有内在原因的，我们可以挖掘出影响用户打分的隐藏因素，进而根据未评分电影与这些隐藏因素的关联度，决定此未评分电影的预测评分

应该有一些隐藏的因素，影响用户的打分，比如电影：演员、题材、年代…甚至不一定是人直接可以理解的隐藏因子

找到隐藏因子，可以对 user 和 item 进行关联（找到是由于什么使得 user 喜欢/不喜欢此 item, 什么会决定 user 喜欢/不喜欢此 item） , 就可以推测用户是否会喜欢某一部未看过的电影

2.1 LFM原理解析

LFM(latent factor model)隐语义模型核心思想是通过隐含特征联系用户和物品，如下图：

P矩阵是User-LF矩阵，即用户和隐含特征矩阵。可以理解成用户对各物品类别的一个偏好信息。LF有三个，表示共总有三个隐含特征。
Q矩阵是LF-Item矩阵，即隐含特征和物品的矩阵，可以理解成各个物品归属到各个类别的信息。
R矩阵是User-Item矩阵，表示用户对物品的偏好信息，由P*Q得来
能处理稀疏评分矩阵

利用矩阵分解技术，将原始User-Item的评分矩阵（稠密/稀疏）分解为P和Q矩阵，然后利用\(P*Q\)还原出User-Item评分矩阵\(R\)。整个过程相当于降维处理，其中：

矩阵值\(P_{11}\)表示用户1对隐含特征1的权重值
矩阵值\(Q_{11}\)表示隐含特征1在物品1上的权重值
矩阵值\(R_{11}\)就表示预测的用户1对物品1的评分，且\(R_{11}=\vec{P_{1,k}}\cdot \vec{Q_{k,1}}\)

利用LFM预测用户对物品的评分，\(k\)表示隐含特征数量：

\[\begin{split} \hat {r}_{ui} &=\vec {p_{uk}}\cdot \vec {q_{ik}} \\&={\sum_{k=1}}^k p_{uk}q_{ik} \end{split} \]

因此最终，我们的目标也就是要求出P矩阵和Q矩阵及其当中的每一个值，然后再对用户-物品的评分进行预测。

2.2 损失函数

矩阵分解得到的预测评分矩阵，与原评分矩阵R在已知评分项上可能会有误差，我们的目标是找到一个最好的分解方式，让分解之后的预测评分矩阵误差最小。

同样对于评分预测我们利用平方差来构建损失函数：

\[\begin{split} Cost &= \sum_{u,i\in R} (r_{ui}-\hat{r}_{ui})^2 \\&=\sum_{u,i\in R} (r_{ui}-{\sum_{k=1}}^k p_{uk}q_{ik})^2 \end{split} \]

加入L2正则化，以防过拟合，λ则需要在实际场景中反复进行实验来得到合适的值。

\[Cost = \sum_{u,i\in R} (r_{ui}-{\sum_{k=1}}^k p_{uk}q_{ik})^2 + \lambda(\sum_U{p_{uk}}^2+\sum_I{q_{ik}}^2) \]

对损失函数求偏导：

\[\begin{split} \cfrac {\partial}{\partial p_{uk}}Cost &= \cfrac {\partial}{\partial p_{uk}}[\sum_{u,i\in R} (r_{ui}-{\sum_{k=1}}^k p_{uk}q_{ik})^2 + \lambda(\sum_U{p_{uk}}^2+\sum_I{q_{ik}}^2)] \\&=2\sum_{u,i\in R} (r_{ui}-{\sum_{k=1}}^k p_{uk}q_{ik})(-q_{ik}) + 2\lambda p_{uk} \\\\ \cfrac {\partial}{\partial q_{ik}}Cost &= \cfrac {\partial}{\partial q_{ik}}[\sum_{u,i\in R} (r_{ui}-{\sum_{k=1}}^k p_{uk}q_{ik})^2 + \lambda(\sum_U{p_{uk}}^2+\sum_I{q_{ik}}^2)] \\&=2\sum_{u,i\in R} (r_{ui}-{\sum_{k=1}}^k p_{uk}q_{ik})(-p_{uk}) + 2\lambda q_{ik} \end{split} \]

2.3 随机梯度下降法优化

梯度下降更新参数\(p_{uk}\)：

\[\begin{split} p_{uk}&:=p_{uk} - \alpha\cfrac {\partial}{\partial p_{uk}}Cost \\&:=p_{uk}-\alpha [2\sum_{u,i\in R} (r_{ui}-{\sum_{k=1}}^k p_{uk}q_{ik})(-q_{ik}) + 2\lambda p_{uk}] \\&:=p_{uk}+\alpha [\sum_{u,i\in R} (r_{ui}-{\sum_{k=1}}^k p_{uk}q_{ik})q_{ik} - \lambda p_{uk}] \end{split} \]

同理：

\[\begin{split} q_{ik}&:=q_{ik} + \alpha[\sum_{u,i\in R} (r_{ui}-{\sum_{k=1}}^k p_{uk}q_{ik})p_{uk} - \lambda q_{ik}] \end{split} \]

随机梯度下降： 向量乘法每一个分量相乘求和

\[\begin{split} &p_{uk}:=p_{uk}+\alpha [(r_{ui}-{\sum_{k=1}}^k p_{uk}q_{ik})q_{ik} - \lambda_1 p_{uk}] \\&q_{ik}:=q_{ik} + \alpha[(r_{ui}-{\sum_{k=1}}^k p_{uk}q_{ik})p_{uk} - \lambda_2 q_{ik}] \end{split} \]

由于P矩阵和Q矩阵是两个不同的矩阵，通常分别采取不同的正则参数，如\(\lambda_1\)和\(\lambda_2\)

算法实现

'''
LFM Model
'''
import pandas as pd
import numpy as np

# 评分预测    1-5
class LFM(object):

    def __init__(self, alpha, reg_p, reg_q, number_LatentFactors=10, number_epochs=10, columns=["uid", "iid", "rating"]):
        self.alpha = alpha # 学习率
        self.reg_p = reg_p    # P矩阵正则
        self.reg_q = reg_q    # Q矩阵正则
        self.number_LatentFactors = number_LatentFactors  # 隐式类别数量
        self.number_epochs = number_epochs    # 最大迭代次数
        self.columns = columns

    def fit(self, dataset):
        '''
        fit dataset
        :param dataset: uid, iid, rating
        :return:
        '''

        self.dataset = pd.DataFrame(dataset)

        self.users_ratings = dataset.groupby(self.columns[0]).agg([list])[[self.columns[1], self.columns[2]]]
        self.items_ratings = dataset.groupby(self.columns[1]).agg([list])[[self.columns[0], self.columns[2]]]

        self.globalMean = self.dataset[self.columns[2]].mean()

        self.P, self.Q = self.sgd()

    def _init_matrix(self):
        '''
        初始化P和Q矩阵，同时为设置0，1之间的随机值作为初始值
        :return:
        '''
        # User-LF
        P = dict(zip(
            self.users_ratings.index,
            np.random.rand(len(self.users_ratings), self.number_LatentFactors).astype(np.float32)
        ))
        # Item-LF
        Q = dict(zip(
            self.items_ratings.index,
            np.random.rand(len(self.items_ratings), self.number_LatentFactors).astype(np.float32)
        ))
        return P, Q

    def sgd(self):
        '''
        使用随机梯度下降，优化结果
        :return:
        '''
        P, Q = self._init_matrix()

        for i in range(self.number_epochs):
            print("iter%d"%i)
            error_list = []
            for uid, iid, r_ui in self.dataset.itertuples(index=False):
                # User-LF P
                ## Item-LF Q
                v_pu = P[uid] #用户向量
                v_qi = Q[iid] #物品向量
                err = np.float32(r_ui - np.dot(v_pu, v_qi))

                v_pu += self.alpha * (err * v_qi - self.reg_p * v_pu)
                v_qi += self.alpha * (err * v_pu - self.reg_q * v_qi)
                
                P[uid] = v_pu 
                Q[iid] = v_qi

                # for k in range(self.number_of_LatentFactors):
                #     v_pu[k] += self.alpha*(err*v_qi[k] - self.reg_p*v_pu[k])
                #     v_qi[k] += self.alpha*(err*v_pu[k] - self.reg_q*v_qi[k])

                error_list.append(err ** 2)
            print(np.sqrt(np.mean(error_list)))
        return P, Q

    def predict(self, uid, iid):
        # 如果uid或iid不在，我们使用全剧平均分作为预测结果返回
        if uid not in self.users_ratings.index or iid not in self.items_ratings.index:
            return self.globalMean

        p_u = self.P[uid]
        q_i = self.Q[iid]

        return np.dot(p_u, q_i)

    def test(self,testset):
        '''预测测试集数据'''
        for uid, iid, real_rating in testset.itertuples(index=False):
            try:
                pred_rating = self.predict(uid, iid)
            except Exception as e:
                print(e)
            else:
                yield uid, iid, real_rating, pred_rating

if __name__ == '__main__':
    dtype = [("userId", np.int32), ("movieId", np.int32), ("rating", np.float32)]
    dataset = pd.read_csv("datasets/ml-latest-small/ratings.csv", usecols=range(3), dtype=dict(dtype))

    lfm = LFM(0.02, 0.01, 0.01, 10, 100, ["userId", "movieId", "rating"])
    lfm.fit(dataset)

    while True:
        uid = input("uid: ")
        iid = input("iid: ")
        print(lfm.predict(int(uid), int(iid)))

三、BiasSvd

BiasSvd其实就是前面提到的Funk SVD矩阵分解基础上加上了偏置项（Bias）也就是衡量商品和用户自身的因素。原因：经过观测，评分数据大部分都是和用户或物品无关的因素产生的效果，即有很大一部分因素是和用户对物品的喜好无关而只取决于用户或物品本身特性。这些独立于用户或独立于物品的因素称为偏置(Bias)部分。

3.1 BiasSvd 思想

Funk-SVD方法通过学习用户和物品的特征向量进行预测，即用户和物品的交互信息。用户的特征向量代表了用户的兴趣，物品的特征向量代表了物品的特点，且每一个维度相互对应，两个向量的内积表示用户对该物品的喜好程度。但是我们观测到的评分数据大部分都是都是和用户或物品无关的因素产生的效果，即有很大一部分因素是和用户对物品的喜好无关而只取决于用户或物品本身特性的。例如，对于乐观的用户来说，它的评分行为普遍偏高，而对批判性用户来说，他的评分记录普遍偏低，即使他们对同一物品的评分相同，但是他们对该物品的喜好程度却并不一样。同理，对物品来说，以电影为例，受大众欢迎的电影得到的评分普遍偏高，而一些烂片的评分普遍偏低，这些因素都是独立于用户或产品的因素，而和用户对产品的的喜好无关。

我们把这些独立于用户或独立于物品的因素称为偏置(Bias)部分，将用户和物品的交互即用户对物品的喜好部分称为个性化部分。

事实上，在矩阵分解模型中偏好部分对提高评分预测准确率起的作用要大大高于个性化部分所起的作用，以Netflix Prize推荐比赛数据集为例为例，Yehuda Koren仅使用偏置部分可以将评分误差降低32%，而加入个性化部分能降低42%，也就是说只有10%是个性化部分起到的作用，这也充分说明了偏置部分所起的重要性，剩下的58%的误差Yehuda Koren将称之为模型不可解释部分，包括数据噪音等因素。

3.2 目标函数

偏置部分主要由三个子部分组成，分别是

训练集中所有评分记录的全局平均数μ，表示了训练数据的总体评分情况，对于固定的数据集，它是一个常数。
用户偏置b\(_{u}\)，独立于物品特征的因素，表示某一特定用户的打分习惯。例如，对于批判性用户对于自己的评分比较苛刻，倾向于打低分；而乐观型用户则打分比较保守，总体打分要偏高。
物品偏置b\(_{i}\)，特立于用户兴趣的因素，表示某一特定物品得到的打分情况。以电影为例，好片获得的总体评分偏高，而烂片获得的评分普遍偏低，物品偏置捕获的就是这样的特征。

则偏置部分表示为 u + b\(_{u}\) + b\(_{i}\)

利用BiasSvd预测用户对物品的评分，\(k\)表示隐含特征数量：

\[\begin{split} \hat {r}_{ui} &=\mu + b_u + b_i + \vec {p_{uk}}\cdot \vec {q_{ki}} \\&=\mu + b_u + b_i + {\sum_{k=1}}^k p_{uk}q_{ik} \end{split} \]

3.2 损失函数

同样对于评分预测我们利用平方差来构建损失函数：

\[\begin{split} Cost &= \sum_{u,i\in R} (r_{ui}-\hat{r}_{ui})^2 \\&=\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i -{\sum_{k=1}}^k p_{uk}q_{ik})^2 \end{split} \]

加入L2正则化：

\[Cost = \sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})^2 + \lambda(\sum_U{b_u}^2+\sum_I{b_i}^2+\sum_U{p_{uk}}^2+\sum_I{q_{ik}}^2) \]

对损失函数求偏导：

\[\begin{split} \cfrac {\partial}{\partial p_{uk}}Cost &= \cfrac {\partial}{\partial p_{uk}}[\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})^2 + \lambda(\sum_U{b_u}^2+\sum_I{b_i}^2+\sum_U{p_{uk}}^2+\sum_I{q_{ik}}^2)] \\&=2\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})(-q_{ik}) + 2\lambda p_{uk} \\\\ \cfrac {\partial}{\partial q_{ik}}Cost &= \cfrac {\partial}{\partial q_{ik}}[\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})^2 + \lambda(\sum_U{b_u}^2+\sum_I{b_i}^2+\sum_U{p_{uk}}^2+\sum_I{q_{ik}}^2)] \\&=2\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})(-p_{uk}) + 2\lambda q_{ik} \end{split} \]

\[\begin{split} \cfrac {\partial}{\partial b_u}Cost &= \cfrac {\partial}{\partial b_u}[\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})^2 + \lambda(\sum_U{b_u}^2+\sum_I{b_i}^2+\sum_U{p_{uk}}^2+\sum_I{q_{ik}}^2)] \\&=2\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})(-1) + 2\lambda b_u \\\\ \cfrac {\partial}{\partial b_i}Cost &= \cfrac {\partial}{\partial b_i}[\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})^2 + \lambda(\sum_U{b_u}^2+\sum_I{b_i}^2+\sum_U{p_{uk}}^2+\sum_I{q_{ik}}^2)] \\&=2\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})(-1) + 2\lambda b_i \end{split} \]

3.3 随机梯度下降法优化

梯度下降更新参数\(p_{uk}\)：

\[\begin{split} p_{uk}&:=p_{uk} - \alpha\cfrac {\partial}{\partial p_{uk}}Cost \\&:=p_{uk}-\alpha [2\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})(-q_{ik}) + 2\lambda p_{uk}] \\&:=p_{uk}+\alpha [\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})q_{ik} - \lambda p_{uk}] \end{split} \]

同理：

\[\begin{split} q_{ik}&:=q_{ik} + \alpha[\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})p_{uk} - \lambda q_{ik}] \end{split} \]

\[b_u:=b_u + \alpha[\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik}) - \lambda b_u] \]

\[b_i:=b_i + \alpha[\sum_{u,i\in R} (r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik}) - \lambda b_i] \]

随机梯度下降：

\[\begin{split} &p_{uk}:=p_{uk}+\alpha [(r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})q_{ik} - \lambda_1 p_{uk}] \\&q_{ik}:=q_{ik} + \alpha[(r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik})p_{uk} - \lambda_2 q_{ik}] \end{split} \]

\[b_u:=b_u + \alpha[(r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik}) - \lambda_3 b_u] \]

\[b_i:=b_i + \alpha[(r_{ui}-\mu - b_u - b_i-{\sum_{k=1}}^k p_{uk}q_{ik}) - \lambda_4 b_i] \]