ML学习笔记 #13 协同过滤推荐算法

本文将介绍机器学习在推荐系统（Recommender Systems）上的应用。推荐系统在互联网环境中的重要性不亚于搜索引擎，具有极高的落地价值。推荐系统的发展十分迅速，基于机器学习的方法在现在看来早已过时，但依然有学习的必要。本文将以电影推荐为例展开介绍。

推荐系统概述

传统推荐算法可以根据数据源的不同分为三类：

• 基于人口统计学的推荐（Demographic-based Recommendation）：依据用户间的相似性；
• 基于内容的推荐（Content-based Recommendation）：依据内容间的相似性；
• 基于协同过滤的推荐（Collaborative Filtering-based Recommendation）：结合前两者。

假设我们是一个电影供应商，现在需要根据用户的观影历史以及评分，来预测每个人可能会给他们没看过的电影打多少分，并以此作为推荐的依据。首先定义如下符号：

• \(n_u\) 代表用户的数量；
• \(n_m\) 代表电影的数量；
• \(r(i, j)\) 代表用户 \(j\) 是否给电影 \(i\) 评过分，如果评过则 \(r(i,j)=1\)；
• \(y^{(i, j)}\) 代表用户 \(j\) 给电影 \(i\) 的评分，前提是 \(r(i,j)=1\)；
• \(m^{(j)}\) 代表用户 \(j\) 评过分的电影的总数。

基于人口统计学的推荐算法

假设我们经过人口统计学的方法，获得了每个用户的属性作为其特征向量 \(u\)，那么我们可以找出 \(u\) 相近的若干用户，则其中某个用户喜欢的电影，其他用户也可能喜欢。此时我们不需要知道关于电影的任何信息。

这种方法的本质是根据「用户间的相似度」，将相似用户喜爱的内容推荐给当前用户。因此，其最重要的数据源就是用户画像（User Profile）：通过收集与分析用户的社会属性、生活习惯等主要信息的数据之后，将用户的信息抽象成标签。

原理上有点类似聚类算法，这里不展开介绍。但需要注意的是，这种方法存在物品冷启动问题（Item Cold Start），即当新物品出现时，由于未与任何用户有过数据交互历史，无法被主动推荐给与相关性较高的目标用户。

然而，获得用户的基本属性并非一件容易的事情（涉及隐私），包括地理位置、设备信息、社交关系、浏览器等等，因此，更好的办法是学习一个用户的偏好 \(\theta\)。一旦获得了 \(\theta\)，我们就可以根据不同电影的属性来预测 \(\theta\rightarrow y\)，其中 \(y\) 就是用户对某部电影的可能评分。

基于内容的推荐算法

现在假设我们对每部电影构建了一个特征向量 \(x\)，例如 \(x_1\) 代表浪漫程度、$ x_2$ 代表动作程度等等。那么对于某个用户而言，我们将其对电影的评分设为 \(y\)，目标就是做线性回归预测 \(x\rightarrow y\)。

具体地，我们记第 \(i\) 部电影的特征向量为 \(x^{(i)}\in\mathbb R^{n+1}\)（包含偏置项），那么对于第 \(j\) 个用户，线性回归的目标就是学习一个参数 \(\theta^{(j)}\)（代表用户对电影特征的偏好），使得： \[ \min_{\theta ^{\left( j \right)}} \frac{1}{2}\sum_{i:r\left( i,j \right) =1}{\left( \left( \theta ^{\left( j \right)} \right) ^Tx^{\left( i \right)}-y^{\left( i,j \right)} \right) ^2}+\frac{\lambda}{2}\sum_{k=1}^n{\left( \theta _{k}^{\left( j \right)} \right) ^2} \] 注意此处去掉了分母中的 \(m^{(j)}\)，因为这个常数不会影响结果。上面的代价函数只是针对一个用户的，由于每个用户的线性回归都是独立的，我们可以将所有用户的代价函数求和： \[ \min_{\theta} J\left( \theta ^{\left( 1 \right)},\cdots ,\theta ^{\left( n_u \right)} \right) :=\frac{1}{2}\sum_{j=1}^{n_u}{\sum_{i:r(i,j)=1}{\left( \left( \theta ^{\left( j \right)} \right) ^Tx^{\left( i \right)}-y^{\left( i,j \right)} \right) ^2}}+\frac{\lambda}{2}\sum_{j=1}^{n_u}{\sum_{k=1}^n{\left( \theta _{k}^{\left( j \right)} \right) ^2}} \] 使用梯度下降法来求解最优解，计算导函数： \[ \frac{\partial J}{\partial \theta _{k}^{\left( j \right)}}=\sum_{i:r(i,j)=1}{\left( \left( \theta ^{\left( j \right)} \right) ^Tx^{\left( i \right)}-y^{\left( i,j \right)} \right) \theta _{k}^{\left( j \right)}}+\lambda \theta _{k}^{\left( j \right)}\;\;[k\ne 0] \] 那么对于一个特征向量为 \(x\) 的电影，用户 \(j\) 对它的评分的预测值就是：\(\left( \theta ^{\left( j \right)} \right) ^Tx\)。

协同过滤 | Collaborative Filtering

基于内容的推荐算法的缺点在于，我们需要知道描述每部电影的特征向量，然而这一点通常很难做到。所以我们需要不是基于内容的推荐算法——协同过滤。协同过滤的本质是一种「特征学习」算法，即自行学习模型所需要的特征，这在深度学习领域是非常经典的思想（尤其是对于复杂特征的建立）。

初始版本

现在我们不知道每部电影的特征向量，但是我们可以询问用户以得到用户的参数 \(\theta^{(j)}\)（譬如让用户选择对不同类型电影的偏好），然后反过来，用 \(\theta^{(j)}\) 去训练出 \(x^{(i)}\)，得到每部电影的特征。具体地，对于第 \(i\) 部电影，我们可以学习它的特征 \(x^{(i)}\)，使得： \[ \min_{x^{\left( i \right)}} \frac{1}{2}\sum_{j:r\left( i,j \right) =1}{\left( \left( \theta ^{\left( j \right)} \right) ^Tx^{\left( i \right)}-y^{\left( i,j \right)} \right) ^2}+\frac{\lambda}{2}\sum_{k=1}^n{\left( x_{k}^{\left( i \right)} \right) ^2} \] 由于每部电影的线性回归是独立的，所以我们可以放在一起训练： \[ \min_x J(x^{\left( 1 \right)},\cdots ,x^{\left( n_m \right)}):=\frac{1}{2}\sum_{i=1}^{n_m}{\sum_{j:r\left( i,j \right) =1}{\left( \left( \theta ^{\left( j \right)} \right) ^Tx^{\left( i \right)}-y^{\left( i,j \right)} \right) ^2}}+\frac{\lambda}{2}\sum_{i=1}^{n_m}{\sum_{k=1}^n{\left( x_{k}^{\left( i \right)} \right) ^2}} \] 训练过程可能用到导函数： \[ \frac{\partial J}{\partial x_{k}^{\left( i \right)}}=\sum_{j:r\left( i,j \right) =1}{\left( \left( \theta ^{\left( j \right)} \right) ^Tx^{\left( i \right)}-y^{\left( i,j \right)} \right)}x_{k}^{\left( i \right)}+\lambda x_{k}^{\left( i \right)}\,\,[k\ne 0] \] 看到这里我们已经发现，这个就是基于人口统计学的推荐算法啊！也就是说「用户」和「内容」是相互的，已知 \(\theta^{(j)}\)，我们可以学习 \(x^{(i)}\)；已知 \(x^{(i)}\)，我们可以学习 \(\theta^{(j)}\)。

因此，我们可以采用迭代的想法（期望最大化算法的思想），随机初始化一个 \(\theta^{(j)}\)，学习出 \(x^{(i)}\)，再用学习出的 \(x^{(i)}\) 去学习 \(\theta^{(j)}\)，再用新的 \(\theta^{(j)}\) 去学习 \(x^{(i)}\)。如此反复迭代，最终得到稳定的电影特征和用户参数。这就是最初始版本的协同过滤算法。

改进版本

事实上，我们没有反复迭代的必要。观察用 \(\theta^{(j)}\) 训练 \(x^{(i)}\) 的优化目标和用 \(x^{(i)}\) 训练 \(\theta^{(j)}\) 的优化目标，我们可以发现，它们除了正则化项以外其实是相同的，都是预测值与实际值的平方误差损失，所以，我们将两个目标综合起来，优化以下函数即可： \[ \begin{aligned} &\min_{x,\theta} J(x^{\left( 1 \right)},\cdots ,x^{\left( n_m \right)},\theta ^{\left( 1 \right)},\cdots ,\theta ^{\left( n_u \right)})\\ &=\frac{1}{2}\sum_{\left( i,j \right) :r\left( i,j \right) =1}{\left( \left( \theta ^{\left( j \right)} \right) ^Tx^{\left( i \right)}-y^{\left( i,j \right)} \right) ^2}+\frac{\lambda}{2}\sum_{i=1}^{n_m}{\sum_{k=1}^n{\left( x_{k}^{\left( i \right)} \right) ^2}}+\frac{\lambda}{2}\sum_{j=1}^{n_u}{\sum_{k=1}^n{\left( \theta _{k}^{\left( j \right)} \right) ^2}}\\ \end{aligned} \] 值得注意的是，在综合起来之前，\(n\) 是我们人为选定的特征维度数，是一个定值；而现在，\(n\) 变成了一个超参数，因此我们也没有必要加上偏置项，所以这里 \(x^{(i)},\theta ^{(j)}\in \mathbb{R} ^n\).

上式的导函数为： \[ \begin{aligned} &\frac{\partial J}{\partial \theta _{k}^{\left( j \right)}}=\sum_{i:r\left( i,j \right) =1}{\left( \left( \theta ^{\left( j \right)} \right) ^Tx^{\left( i \right)}-y^{\left( i,j \right)} \right) \theta _{k}^{\left( j \right)}+\lambda \theta _{k}^{\left( j \right)}}\\ &\frac{\partial J}{\partial x_{k}^{\left( i \right)}}=\sum_{j:r\left( i,j \right) =1}{\left( \left( \theta ^{\left( j \right)} \right) ^Tx^{\left( i \right)}-y^{\left( i,j \right)} \right) x_{k}^{\left( i \right)}+\lambda x_{k}^{\left( i \right)}}\\ \end{aligned} \] 接下来就可以梯度下降算法完成优化了。最终收敛得到的向量 \(x^{(i)},\theta ^{(j)}\)，可以用于相似电影的推荐，也可以用于相似用户的发掘。

向量化版本

为了代码的运行效率，我们先将该算法向量化。

• 矩阵 \(Y=\left[ y^{\left( i,j \right)} \right] \in \mathbb{R} ^{n_m\times n_u}\)，即第 \(i\) 行第 \(j\) 列表示用户 \(j\) 对电影 \(i\) 的评分；
• 矩阵 \(X=\left[ \begin{array}{c}\left( x^{\left( 1 \right)} \right) ^T\\\vdots\\\left( x^{\left( n_m \right)} \right) ^T\\\end{array} \right] \in \mathbb{R} ^{n_m\times n}\)，即第 \(i\) 行表示电影 \(i\) 的特征向量；
• 矩阵 \(\Theta =\left[ \begin{array}{c}\left( \theta ^{\left( 1 \right)} \right) ^T\\\vdots\\\left( \theta ^{\left( n_u \right)} \right) ^T\\\end{array} \right] \in \mathbb{R} ^{n_u\times n}\)，即第 \(j\) 行表示用户 \(j\) 的偏好向量。

如此，线性回归的预测值可以构成矩阵： \[ Y=X\Theta^T\in\mathbb R^{n_m\times n_u} \] 上式可以视作对矩阵 \(Y\) 的一个低秩分解（Low Rank Factorization）。

均值归一化与冷启动

当一个新用户下载软件、新电影上市时，我们还是会遇到冷启动问题，即没有电影或用户与之产生过交互！此时我们的矩阵 \(Y\) 会新增一行或一列 N/A 值（没有评分），如果我们还是将其代入到目标函数中，则对其产生影响的只有正则化项： \[ \frac{\lambda}{2}\sum_{i=1}^{n_m}{\sum_{k=1}^n{\left( x_{k}^{\left( i \right)} \right) ^2}}+\frac{\lambda}{2}\sum_{j=1}^{n_u}{\sum_{k=1}^n{\left( \theta _{k}^{\left( j \right)} \right) ^2}} \] 而显然当 \(\theta _{k}^{\left( j \right)}\) 或 \(x_{k}^{\left( i \right)}\) 取全 \(0\) 时，正则化项取值最小。代入 \(y^{\left( i,j \right)}=\left( \theta ^{\left( j \right)} \right) ^Tx^{\left( i \right)}\) 进行预测，得到的所有值也都是 \(0\)。但实际上新用户不会给所有电影都打 \(0\) 分，新电影也不会收到所有人的 \(0\) 分差评。

为了能真正进行推荐，我们选择将 \(Y\) 矩阵按行进行均值归一化（Mean Normalization），即减去每部电影得分的平均值。此时我们再将 \(Y'\) 代入到目标函数中计算，得到的结果虽然 \(\theta _{k}^{\left( j \right)}\) 或 \(x_{k}^{\left( i \right)}\) 还是取全 \(0\)，且 \(Y'\) 中对应的行或列也是全 \(0\)，但当我们恢复均值时，原始的 \(Y\) 就有取值了。

不难发现，这里 \(Y\) 的取值就是该行或该列的其他值的平均值，也就是说我们朴素地假设新用户给电影评分就是电影现有平均分，新电影得到的评分也是用户的历史打分平均分。这也是我们平时启动软件看到的「热门电影」所做的事情。

代码实现

下面以 Coursera 上的数据集 ex8_movies.mat 为例，这是一个电影评分数据集 MovieLens 的子集。其中包含 \(Y\) 矩阵表示评分结果，\(R\) 矩阵表示用户是否给电影评过分。现在我们要在上面实现一个协同过滤算法，预测每个用户对所有未评分电影的打分。代码如下：

import numpy as np
from scipy.io import loadmat
import scipy.optimize as opt

# load data
data = loadmat('ex8_movies.mat')
Y = data['Y'] # (1682, 943)
R = data['R'] # (1682, 943)
n_movie, n_user = Y.shape # (1682, 943)
Ymean = Y.mean(axis=1, keepdims=True)

# parameter
n_features = 50 # numbers of feature
X = np.random.standard_normal((n_movie, n_features))
Theta = np.random.standard_normal((n_user, n_features))
lamb = 1 # regularization

def serialize(X, Theta):
    return np.concatenate((X.ravel(), Theta.ravel()))

def deserialize(param):
    return param[:n_movie * n_features].reshape(n_movie, n_features), \
           param[n_movie * n_features:].reshape(n_user, n_features)

# linear regression
def J(param, Y, R, lamb):
	# deserialize
	X, Theta = deserialize(param)
	# calculate
	inner = np.multiply(X @ Theta.T - Y, R)
	inner_term = np.power(inner, 2).sum() / 2
	reg_term = lamb / 2 * np.power(param, 2).sum()
	return inner_term + reg_term

def partJ(param, Y, R, lamb):
	# deserialize
	X, Theta = deserialize(param)
	# calculate
	inner = np.multiply(X @ Theta.T - Y, R)
	X_grad = inner @ Theta + lamb * X
	Theta_grad = inner.T @ X + lamb * Theta
	return serialize(X_grad, Theta_grad)

# train, about 10 mins
res = opt.minimize(fun=J, # cost function
                   x0=serialize(X, Theta), # serialize input param
                   args=(Y - Ymean, R, lamb), # input data
                   method='TNC',
                   jac=partJ) # serialize gradient
print(res)

# predict
X_trained, Theta_trained = deserialize(res.x)
Y_pred = X_trained @ Theta_trained.T + Ymean

# parse movie_id.txt
movie_list = []
with open('movie_ids.txt', encoding='latin-1') as file:
	for line in file:
		tokens = line.strip().split(' ')
		movie_list.append(' '.join(tokens[1:]))
movie_list = np.array(movie_list)

# sort rating and show result of first user
idx = np.argsort(Y_pred[:, 0])[::-1]
print('Top 10 movies for first user:')
for movie in movie_list[idx][:10]:
	print(movie)

使用 Scipy 第三方库优化结果为：

    fun: 11099.432068949513
    jac: array([-5.25906371e-06, -3.12407994e-07, -1.53592980e-06, ...,
      -5.59603399e-07, -3.29349047e-07, -3.10467495e-07])
message: 'Converged (|f_n-f_(n-1)| ~= 0)'
   nfev: 7740
    nit: 235
 status: 1
success: True
      x: array([-0.50327981,  0.29839561, -0.27672091, ...,  0.34445801,
       1.09017551, -0.34508789])

用该参数找到第一个用户预测评分最高的 \(10\) 部电影：

Top 10 movies for first user:
Titanic (1997)
Deer Hunter, The (1978)
Lawrence of Arabia (1962)
In the Name of the Father (1993)
Being There (1979)
Ran (1985)
Philadelphia (1993)
Tombstone (1993)
Great Escape, The (1963)
Manchurian Candidate, The (1962)