矩阵分解在协同过滤推荐算法中的应用以及简单的代码实现

转载出处
代码出处

---

1. 矩阵分解用于推荐算法要解决的问题

在推荐系统中，常常遇到用户-物品矩阵稀疏性问题，即：有很多用户和物品，也有少部分用户对少部分物品的评分，我们希望预测目标用户对其他未评分物品的评分，进而将评分高的物品推荐给目标用户。比如下面的用户物品评分表：

用户\物品	物品1	物品2	物品3	物品4	物品5	物品6	物品7
用户1	3	-	5	-	-	1	-
用户2	-	2	-	-	-	-	4
用户3	-	-	-	4	-	-	-
用户4	-	-	2	-	-	-	1
用户5	1	-	-	-	4	-	-

对于每个用户，我们希望较准确的预测出用户对未评分物品的评分。对于这个问题我们有很多解决方法，本文我们关注于用矩阵分解的方法来做。如果将m个用户和n个物品对应的评分看做一个矩阵M，我们希望通过矩、阵分解来解决这个问题。

2. 传统的奇异值分解SVD用于推荐

将这个用户物品对应的m×n矩阵M进行SVD分解，并通过选择部分较大的一些奇异值来同时进行降维，也就是说矩阵M此时分解为：$$M_{m\ast n}=P_{m\ast k}^T{Q}_{k\ast n}$$
我们知道SVD分解已经很成熟了，但是FunkSVD如何将矩阵M分解为P和Q呢？这里采用了线性回归的思想。我们的目标是让用户的评分和用矩阵乘积得到的评分残差尽可能的小，也就是说，可以用均方差作为损失函数，来寻找最终的P和Q。

对于某一个用户评分$m_{ij}$，如果用FunkSVD进行矩阵分解，则对应的表示为$q_j^T{p}_i$，采用均方差做为损失函数，则我们期望$(m_{ij}-q_j^T{p}_i{)}^2$尽可能的小，如果考虑所有的物品和样本的组合，则我们期望最小化下式：$$\sum _{ij}(m_{ij}-q_j^T{p}_i{)}^2$$只要能够最小化上式，并求出极值对应的$p_i,q_j$，则我们最终可以得到矩阵P和Q，那么对于任意矩阵M任意一个空白评分的位置，我们可以通过$q_j^T{p}_i$计算预测评分。很漂亮的方法！

当然，在实际应用中，我们为了防止过拟合，会加入一个L2的正则化项，因此正式的FunkSVD的优化目标函数$J(p,q)$是这样的：$$\underbrace{argmin}\sum (m_{ij}-q_j^T{p}_i{)}^2+\lambda ({\Vert p_i\Vert }_2^2+{\Vert q_j\Vert }_2^2)$$
其中$\lambda$为正则化系数，需要调参。对于这个优化问题，一般采用梯度下降进行优化得到结果。

将上式分别对$p_i,q_j$求导我们得到:$$\frac{\partial J}{\partial p_i}=-2(m_{ij}-q_j^T{p}_i)q_j+2\lambda p_i$$$$\frac{\partial J}{\partial q_j}=-2(m_{ij}-q_j^T{p}_i)q_j+2\lambda q_j$$
则梯度下降迭代时，$p_i,q_j$的迭代公式为：$$p_i=p_i+\alpha ((m_{ij}-q_j^T{p}_i)q_j-\lambda p_i)$$$$q_j=q_j+\alpha ((m_{ij}-q_j^T{p}_i)p_i-\lambda q_j)$$通过迭代我们最终可以得到P和Q矩阵，进而用于推荐。

4. BiasSVD算法用于推荐

在FunkSVD算法火爆之后，出现了很多FunkSVD的改进版算法。其中BiasSVD算是改进的比较成功的一种算法。BiasSVD假设评分系统包括三部分的偏置因素：一些和用户物品无关的评分因素，用户有一些和物品无关的评分因素，称为用户偏置项。而物品也有一些和用户无关的评分因素，称为物品偏置项。这其实很好理解。比如一个垃圾山寨货评分不可能高，自带这种烂属性的物品由于这个因素会直接导致用户评分低，与用户无关。

假设评分系统平均分为$\mu$，第i个用户的用户偏置为 $b_i$，而第j个物品的物品偏置为 $b_j$，则加入了偏置之后的优化目标为:$$\underbrace{argmin}\sum (m_{ij}-\mu -b_i-b_j-q_j^T{p}_i{)}^2+\lambda ({\Vert p_i\Vert }_2^2+{\Vert q_j\Vert }_2^2+{\Vert b_i\Vert }_2^2+{\Vert b_j\Vert }_2^2)$$

这个优化目标也可以采用梯度下降法求解。和FunkSVD不同的是，此时我们多了两个偏执项$b_i,b_j$，$p_i,q_j$的迭代公式和FunkSVD类似，只是每一步的梯度导数稍有不同而已，这里就不给出了。而$b_i,b_j$一般可以初始设置为0，然后参与迭代。这里给出$b_i,b_j$的迭代方法
$$b_i=b_i+\alpha (m_{ij}-\mu -b_i-b_j-q_j^T{p}_i-\lambda b_i)$$$$b_j=b_j+\alpha (m_{ij}-\mu -b_i-b_j-q_j^T{p}_i-\lambda b_j)$$
通过迭代我们最终可以得到P和Q，进而用于推荐。BiasSVD增加了一些额外因素的考虑，因此在某些场景会比FunkSVD表现好。

5. 矩阵分解推荐方法小结

FunkSVD将矩阵分解用于推荐方法推到了新的高度，在实际应用中使用也是非常广泛。当然矩阵分解方法也在不停的进步，目前张量分解和分解机方法是矩阵分解推荐方法今后的一个趋势。

对于矩阵分解用于推荐方法本身来说，它容易编程实现，实现复杂度低，预测效果也好，同时还能保持扩展性。这些都是它宝贵的优点。当然，矩阵分解方法有时候解释性还是没有基于概率的逻辑回归之类的推荐算法好，不过这也不影响它的流形程度。小的推荐系统用矩阵分解应该是一个不错的选择。大型的话，则矩阵分解比起现在的深度学习的一些方法不占优势。

6. Python实现

import numpy

def matrix_factorization(R, P, Q, K, steps=5000, alpha=2e-4, beta=2e-2):
    Q = Q.T
    n, m = R.shape
    ToF = (R != 0)
    # 这一部分不知道怎么简化，有简化版的麻烦告诉我 TAT
    for step in range(steps):
        for i in range(len(R)):
            for j in range(len(R[i])):
                if R[i][j] != 0:
                    mij = R[i][j] - numpy.dot(P[i,:], Q[:,j])
                    for k in range(K):
                        P[i][k] = P[i][k] + alpha * (2 * mij * Q[k][j] - beta * P[i][k])
                        Q[k][j] = Q[k][j] + alpha * (2 * mij * P[i][k] - beta * Q[k][j])
        mR = numpy.dot(P, Q)
        e = np.sum(np.multiply(ToF, np.power(R - mR, 2)))    #损函数第一项
        # 损失函数第二项
        a = np.multiply(np.sum(np.power(P, 2), axis=1).reshape(n,1), np.ones((n,m))) 
        b = np.multiply(np.sum(np.power(Q, 2), axis=0).reshape(m,1), np.ones((m,n))).T
        e += np.sum(0.02 * np.multiply(ToF, a+b))
        if e < 0.001:
            break
    return P, Q.T

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,M = R.shape
K = 2

P = numpy.random.rand(N,K)
Q = numpy.random.rand(M,K)

nP, nQ = matrix_factorization(R, P, Q, K)
nR = numpy.dot(nP, nQ.T)

# P
array([[ 2.22839121, -0.30868672],
       [ 1.78410159, -0.15645184],
       [ 0.95396851,  1.95232724],
       [ 0.8072353 ,  1.54787026],
       [ 1.20254888,  1.44023592]])
# Q
array([[ 2.17409209, -0.51189425],
       [ 1.2816283 , -0.21006825],
       [ 2.12220296,  1.63470871],
       [ 0.74945871,  2.17573691]])
# 计算后的R
array([[5.00274265, 2.92081451, 4.22448555, 0.99846612],
       [3.95888795, 2.31942065, 3.53047248, 0.99671242],
       [1.07463028, 0.81251106, 5.21600112, 4.96271043],
       [0.96265798, 0.7094172 , 4.24343413, 3.97274798],
       [1.87720352, 1.23867283, 4.90641901, 4.03483519]])