基于矩阵分解的推荐算法

矩阵分解算法是学习协同过滤算法的基础，接下来就来讲一件什么是矩阵分解。
一，基于矩阵分解算法的相关理论介绍
很显然，我们要是想做推荐系统，最基本的一个数据就是用户对这个商品的特性的反应，就拿用户-物品的评分矩阵来说

矩阵中，描述了五个用户（U1,U2,U3,U4,U5）对四个物品（D1,D2,D3,D4）的评分，但是我们发现有些表格中没有填写评分
的数字，这个时候我们就需要把没有填写的评分根据已知的分数进行预测，然后根据评分的高低给用户进行推荐。
那么如何进行预测缺失的评分呢？那么我们可以把这些想象成点，然后可以根据梯度下降跟回归方程进行连续值的预测。假设是MN维矩阵（M行N列），然后我们可以分解成P矩阵和Q矩阵，其中P矩阵是MK,Q矩阵维度是K*N。（K是过渡数值，也就是过渡维度）

对于P,Q矩阵的解释，直观上是M个用户对K个主题的关系，Q矩阵是K个主题对N个物品的关系，具体K是多少，这就需要自己的调节了，K只是一个参数，通常取值在10到100之间

式子2

其中r代表 R原始矩阵，p，q代表P,Q矩阵。rij代表R矩阵第i行第j列的数字，后边的同理
我们对于矩阵的知识可以知道，这个矩阵i行j列的数值等于P矩阵i行的所有数值乘以q第j列的所有数值的和
然后我们可以知道回归方程的最小二乘法的公式如下
然后把式子2带入式子3中得到下图

式子3

对于损失函数的左边项，表示的是R^ 第i行，第j列的元素值，对于如何衡量，我们分解的好坏呢，式子3，给出了衡量标准，也就是损失函数，平方项损失，最后的目标，就是每一个元素(非缺失值)的e(i,j)的总和最小
然后我们需要用到L2正则化

因此最终的损失函数为：

然后可以使用梯度下降算法对其求解P矩阵和Q矩阵
最后结果为：

代码如下：

import numpy as np

#  original_matrix --> P*Q=R
def matrix_factorization(original_matrix, K, alpha, beta, epochs):
    '''
    :param original_matrix(mat): 原始矩阵
    :param K(int): 分解矩阵中间维度
    :param alpha(float): 学习率
    :param beta(float): 惩罚性系数
    :param epochs(int): 最大迭代次数
    :return: 分解后的两个矩阵P,Q
    '''
    original_matrix = np.mat(original_matrix)
    #np.mat 生成矩阵
    M, N = original_matrix.shape
    #np.shape 获取矩阵的行数和列数，然后赋值给M,N  即M代表行数 N代表列数
    P = np.mat(np.random.random((M, K)))
    Q = np.mat(np.random.random((K, N)))
    #np.random.random((M,K))  生成M行K列的浮点数 ，搭配上mat 也就是生成矩阵
    print(".................")
    print(P,Q)
    print(".................")
    loss = 1.0
    epoch = 0
    while loss >= 0.001 and epoch <= epochs:
        for m in range(M):
            for n in range(N):
                if original_matrix[m, n] > 0:  # 非缺失值
                    r = original_matrix[m, n]
                    r_ = 0  # R[m, n]
                    for k in range(K):  # 计算[m, n]位置的误差
                        r_ += P[m, k]*Q[k, n]
                    e = r - r_
                    for k in range(K):  # 更新P[m, :]与Q[:, n]的值
                        P[m, k] += alpha*(2*e*Q[k, n]-beta*P[m, k])
                        Q[k, n] += alpha*(2*e*P[m, k]-beta*Q[k, n])
        loss = 0.0
        for m in range(M):
            for n in range(N):
                if original_matrix[m, n] > 0:
                    r = original_matrix[m, n]
                    r_ = 0.0
                    regularization = 0.0
                    for k in range(K):
                        r_ += P[m, k]*Q[k, n]  # 偏差
                        regularization += P[m, k]**2+Q[k, n]**2  # L2正则化
                    e = r - r_
                    loss += e**2 + (beta/2)*regularization  # 总损失
        # if epoch % 200 == 0:
        #     print('epoch:{}, loss: {}'.format(epoch, loss))
        epoch += 1
    return P, Q


if __name__ == '__main__':
    low, high = 1, 10
    size = 10, 8
    original_matrix = np.random.randint(low=low, high=high, size=size)
    missing_rate = 0.3
    n_counts = int(size[0]*size[1]*missing_rate)
    while n_counts:
        pos = np.random.randint(size[0]*size[1])
        row, col = pos//size[1], pos%size[1]
        if original_matrix[row, col] != 0:
            original_matrix[row, col] = 0
            n_counts -= 1
    # print(original_matrix)
    P, Q = matrix_factorization(original_matrix, 6, 0.004, 0.02, 12000)
    print(np.dot(P, Q))

参考：
文章一
文章二