协同过滤与隐语义模型

在机器学习问题中，我们见到的数据集通常是如下的格式：

input target

$X_{1}$ $y_{1}$

$X_{2}$ $y_{2}$

... ...

$X_{N}$ $y_{N}$

，一个输入向量的集合 $\{X_{i}\}$ 以及对应的数据集合 $\{y_{i}\}$ , $\{y_{i}\}$ 就是我们想要去预测的值。

对于这样的数据集，非常适合使用在像回归模型这样标准的机器学习算法中。

但是在这里，我们的数据集用如下的形式来表示会更加合适：

user item rating

$\vec{u_{1}}$ $\vec_{v_{1}}$ $r(\vec{u_{1}},\vec{v_{1}})$

$\vec{u_{2}}$ $\vec_{v_{2}}$ $r(\vec{u_{2}},\vec{v_{2}})$

... ... ...

$\vec{u_{N}}$ $\vec_{v_{N}}$ $r(\vec{u_{N}},\vec{v_{N}})$

,作为一个用户、主题和得分的三元组的集合，从中我们可以获得确定的用户给确定的主题的打分。对于所有的用户其做出评分的主题都不可能完全相同，因此我们总是会缺失某个用户对于主题的评分。

我们可以用矩阵来形象化地表示：

也就是说，在矩阵中的许许多多的评分是缺失的（不是0分，是完全没有打过分）。

有了对待数据的抽象方式以后，我们就可以继续下面的环节了。

首先是加载我们的数据

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import time
import warnings
warnings.filterwarnings('ignore')
np.random.seed(1)
%matplotlib inline
plt.style.use('ggplot')data = pd.read_csv('ml-latest-small/ratings.csv')
movies = pd.read_csv('ml-latest-small/movies.csv')
movies = movies.set_index('movieId')[['title', 'genres']]

这里加载了两组数据，一个是评分数据一个是电影的数据，数据内容如下

评分数据

电影的数据

继续对数据集进行进一步的了解

# How many users?
print(data.userId.nunique(), 'users')# How many movies?
print(data.movieId.nunique(), 'movies')# How possible ratings?
print(data.userId.nunique() * data.movieId.nunique(), 'possible ratings')# How many do we have?
print(len(data), 'ratings')
print(100 * (float(len(data)) / 
(data.userId.nunique() * data.movieId.nunique())), '% of possible ratings')

运行结果

668 users
10325 movies
6897100 possible ratings
105339 ratings
1.5272940801206305 % of possible ratings

我们拥有：

700左右的用户

10000左右的电影

可能的打分组合数为7000000种左右

但实际上我们拥有的打分组合数为100000左右，只占整个打分可能组合数的1.5%左右；由此可见我们的矩阵是相当稀疏的。

# Number of ratings per users
fig = plt.figure(figsize=(10, 10))
ax = plt.hist(data.groupby('userId').apply(lambda x: len(x)).values, bins=50)
plt.title("Number of ratings per user")

可以看到，用户的评分积极性并不是很高。

# Number of ratings per movie
fig = plt.figure(figsize=(10, 10))
ax = plt.hist(data.groupby('movieId').apply(lambda x: len(x)).values, bins=50)
plt.title('Number of ratings per movie')

由于用户的打分积极性不高，所以每部电影被打分的次数也是屈指可数。

# Ratings distribution
fig = plt.figure(figsize=(10, 10))
ax = plt.hist(data.rating.values, bins=5)
plt.title("Distribution of ratings")

用户打分的范围是0~5分，这里查看一下打分的分布情况；总体来说是符合事实分布的，因为大多数的事物都会得到一个中等的评价，好的和坏的都是少数。

# Average rating per user
fig = plt.figure(figsize=(10, 10))
ax = plt.hist(data.groupby('userId').rating.mean().values, bins=10)
plt.title("Average rating per user")

这里展现的是用户的平均打分的分布；个人理解这个分布是比较正常的，大部分的人的平均打分都会比较中庸。

# Average rating per movie
fig = plt.figure(figsize=(10, 10))
ax = plt.hist(data.groupby('movieId').rating.mean().values, bins=10)
plt.title('Average rating per movie')

电影的平均得分的分布，个人感觉也很正常，还是那句话，事物的好的和坏的都是少数的，大部分都是中庸的。

# Top Movies
average_movie_rating = data.groupby('movieId').mean()
top_movies = average_movie_rating.sort_values('rating', ascending=False).head(10)
pd.concat([movies.loc[top_movies.index.values],average_movie_rating.loc[top_movies.index.values].rating], axis=1)

这里是根据电影的平均分来排出电影的top10；但是大家（爱看外国电影的朋友）会发现这些电影几乎都没有听说过！其原因很简单，有很多冷门电影都只有很少的观众，很少的评分，所以平均分就会比较高。一个极限情况下的例子就是一个电影只有一次评分（只有一个人看过），而这个人评了5分。那么这个电影就是排名很高了。所以这种评定top10的方法并不合理。

# Robust Top Movies - Lets weight the average
#rating by the square root of number of ratings
top_movies = data.groupby('movieId').apply(lambda x:len(x)**0.5 * x.mean())
.sort_values('rating', ascending=False).head(10)
pd.concat([movies.loc[top_movies.index.values], average_movie_rating.loc[top_movies.index.values].rating], axis=1)

这次的top10是用sqrt(评分数量)来加持每部电影的平均得分，也就是根据平均得分和评分数量来综合排序。这次的排序结果就好多了，起码《肖申克的救赎》我是看过很多次的，《辛德勒的名单》、星球大战系列也是都看过的。

controversial_movies = data.groupby('movieId').apply(lambda x:len(x)**0.25 * x.std())
.sort_values('rating', ascending=False).head(10)
pd.concat([movies.loc[controversial_movies.index.values], average_movie_rating.loc[controversial_movies.index.values].rating], axis=1)

这是一个最具争议的电影（也就是评分的方差最大的）。

下面我们开始预测所有的缺失的 $\tilde{\gamma }$ 。

协同过滤试图通过分享行(users)和列(items)之间的信息来补全uset-item矩阵。

一些想法：

（1）用每一个item的平均值来作为预测值

这是一个高效而简单的想法，结果是所有在该item上缺失评分的用户，都会得到一个相同的补全评分。同时对于那些评分本来就较少的item的预测结果也不会很好。

（2）利用item的属性来为每一个user建立一个模型，这就是基于内容的过滤。

可以很好的工作，但是有时候属性有时候并不是那么好（我们如何才能知道好动作片与坏片之间的区别？）。同时我们也需要训练成千上万的分类器（每个用户训练一个），这也就造成另外一个问题————用户之间没有信息的共享。

更大的一个问题是，有时候我们都没有属性，就是只有一个itemId和userId决定的评分的情况。

（3）利用相近的users或者items来进行预测，也就是把相近的users看好的items推荐给该user或者把相近的items上评分较高的users归属到该item。

也就是k-Nearest-Neighbours算法，一个标准的机器学习算法。这可以很好地工作，但需要一些工程化的工作才能很好地运行并且能够很好地扩展（因为我们需要查询或预先计算很多信息）。

（4）发掘items的隐藏的属性以及每个用户对这些属性的偏好（基于users和items的learning "embeddings"）

假设每一个item都拥有一组隐藏的属性，同时每一个user都item的每一个属性都有自己的偏好。那么评分结果就可以看做是两者的内积，如下图所示：

如果从整个矩阵的角度来看，那么是针对这个矩阵的分解：

有了这个思想以后，最重要的问题就是如何来找到这些隐藏的属性。我们将建立一个机器学习模型并使用数学优化方法来找到item和user的隐藏的属性。

警告！数学！警告！数学！警告！数学！

我们定义如下的函数：

$\tilde{r}=f(u,i)$

我们期望我们的模型输出的 $\tilde{r}$ 与真实的r之间的误差能够尽量地减少；最基础的一个误差衡量标准就是mean squared error：

$\wp =\frac{1}{N}\sum_{k=1...N}^{ }(\tilde{r}_{k}-r_{k})^2$

定义user u的隐藏属性为向量 $\vec_{p_{u}}$ ,定义item i的隐藏属性为向量 $\vec_{q_{i}}$ ，带入向量后的模型如下：

$f(u,i)={\vec_{q_{i}}}^{T}\vec_{p_{u}}$

在建立模型的过程中我们还需要考虑一些奇淫巧技的东西；因为我们要面对极端的user和item，比如一个user就是喜欢打高分，而一个item总是得到很低的分数。所以我们在建立模型的时候为每一个user和每一个item引入自己特有的偏移，同时再引入一个总的偏移。这样能够使的两个属性向量能够更加专心地来表示属性，而不是还要表示偏移。

我们定义主偏移： $b$

user u的偏移： $b_{u}$

item i的偏移： $b_{i}$

最终我们的模型变成了如下的样子：

$f(u,i)=b+b_{i}+b_{u}+{\vec_{q_{i}}}^T\vec_{p_{u}}$

数学、数学、还是数学！

有了模型以后就需要进行优化来找到偏移和属性向量的值；存在着许许多多的优化方法，这里使用随机梯度下降法。

通过对数据集中的数据进行迭代，每次将参数朝着适当的方向移动来降低均方误差；直到达到等待时长或者误差不在降低为止。也就是沿着参数（通常用 $\theta$ 来表示）的负梯度方向来进行移动，表示为：

$\frac{\partial \wp }{\partial \theta}$ $\theta = \{b, b_{i},b_{u},{\vec_{q}}_{i},{\vec_{p}}_{u}\}$ 针对于每一个user和item

梯度下降的迭代公式是：

$\theta \leftarrow \theta - \alpha\frac{\partial\wp}{\partial\theta}$

其中 $\alpha$ 表示每次迭代的步长，称之为学习率。

这里使用的是随机梯度下降，也就是每次从数据集中随机抽取一个数据进行梯度下降迭代；而标准的梯度下降是在整个数据集上完成计算后，进行一次梯度下降迭代。从总体上来说两者的效果相差不多，但是随机梯度下降法要快的多。

概括：

》在训练集上进行迭代

》更新参数来最小化误差

》一直迭代到没有提升的空间或者达到迭代上限

一图展示梯度下降（--其实展示的不怎么好）：

现在可以开始训练我们的模型了：

把数据分为训练集和测试集

ratings = data[['userId', 'movieId', 'rating']].values# Shuffle training examples so that movies and users are evenly distributed
np.random.shuffle(ratings)n_users, n_items, _ = ratings.max(axis=0) + 1
n = len(ratings)split_ratios = [0, 0.7, 0.85, 1]
train_ratings, valid_ratings, test_ratings = 
[ratings[int(n*lo):int(n*up)] for (lo, up) in zip(split_ratios[:-1], split_ratios[1:])]

定义model对象，实现的代码内容有：

》初始化参数

》保存于加载参数

》对给定的user-item对进行预测

》更新参数来减小误差

代码如下：

class Model(object):gradients = ['dL_db', 'dL_dbu', 'dL_dbv', 'dL_dU', 'dL_dV']def __init__(self, latent_factors_size, users, items):self.model_parameters = []self.n_users = usersself.n_items = itemsfor (name, value) in self.initialize_parameters(latent_factors_size):setattr(self, name, value)self.model_parameters.append(name)# Used to save parameters during the optimizationdef save_parameters(self):return [(name, np.copy(getattr(self, name))) for name in self.model_parameters]# Used to reload the best parameters once the optimization is finisheddef load_parameters(self,  parameters):for (name, value) in parameters:setattr(self, name, value)# Random embedding generation from normal distribution, given a size and variancedef initialize_parameters(self, latent_factors_size=100, std=0.05):U = np.random.normal(0., std, size=(int(self.n_users) + 1, latent_factors_size))V = np.random.normal(0., std, size=(int(self.n_items) + 1, latent_factors_size))u = np.zeros(int(self.n_users) + 1)v = np.zeros(int(self.n_items) + 1)return zip(('b', 'u', 'v', 'U', 'V'), (0, u, v, U, V))# Compute the gradient of the biases and embedings, given the user-itemdef compute_gradient(self, user_ids, item_ids, loc_ratings):predicted_ratings = self.predict(user_ids, item_ids)redidual = loc_ratings - predicted_ratings# biasesdL_db = -2 * redidualdL_dbu = -2 * redidualdL_dbv = -2 * redidual# embeddingseu = self.U[int(user_ids)]ev = self.V[int(item_ids)]dL_dU = -2 * redidual * evdL_dV = -2 * redidual * eu# Regularizationl2 = 0.1dl2eu_dU = l2 * 2 * eudl2ev_dV = l2 * 2 * evdl2bu_dbu = l2 * 2 * self.u[int(user_ids)]dl2bv_dbv = l2 * 2 * self.v[int(item_ids)]dL_dbu = dL_dbu + dl2bu_dbudL_dbv = dL_dbv + dl2bv_dbvdL_dU = dL_dU + dl2eu_dUdL_dV = dL_dV + dl2ev_dVresult = []for x in Model.gradients:result.append((x, eval(x)))return dict(result)# Sum of the biases and dot product of the embeddingsdef predict(self, user_ids, item_ids):user_ids = user_ids.astype('int')item_ids = item_ids.astype('int')return sum([self.b, self.u[user_ids], self.v[item_ids],(self.U[user_ids] * self.V[item_ids]).sum(axis=-1)])# Preform a gradient descent stepdef update_parameters(self, luser, litem, lrating, learning_rate=0.005):gradients = self.compute_gradient(luser, litem, lrating)self.b = self.b - learning_rate * gradients['dL_db']int_luser = int(luser)int_litem = int(litem)self.u[int_luser] = self.u[int_luser] - learning_rate * gradients['dL_dbu']self.v[int_litem] = self.v[int_litem] - learning_rate * gradients['dL_dbv']self.U[int_luser] = self.U[int_luser] - learning_rate * gradients['dL_dU']self.V[int_litem] = self.V[int_litem] - learning_rate * gradients['dL_dV']

几个有用的函数：

def sample_random_training_index():"""Generate a random number"""return np.random.randint(0, len(train_ratings))def compute_rmse(x, y):"""Compute root mean squared error between x and y"""return ((x - y) ** 2).mean() ** 0.5def get_rmse(loc_ratings, loc_model):return compute_rmse(loc_model.predict(*loc_ratings.T[:2]), loc_ratings.T[2])def get_trainset_rmse(loc_model):return get_rmse(train_ratings, loc_model)def get_validset_rmse(loc_model):return get_rmse(valid_ratings, loc_model)def get_testset_rmse(loc_model):return get_rmse(test_ratings, loc_model)

初始化模型和优化参数

model = Model(latent_factors_size=100, users=n_users, items=n_items)
model.b = train_ratings[:, 2].mean()sgd_iteration_count = 0
best_validation_rmse = 9999
patience = 0
update_frequency = 10000train_errors = []
valid_errors = []
test_errors = []best_parameters = None

进行梯度下降优化

一些和优化相关的概念的补充说明：

》我们会用验证集来衡量模型的表现。如果经过N轮的迭代，模型在验证集上的表现都不在提升，那么训练就会终止。这个N就是（patience）；

》我们会每隔10000次迭代，就保存一次训练集、验证集、测试集的误差；

》每当遇到能使模型在验证集上的误差是当前最小的模型参数，就保存该组参数。

start_time = time.time()while True:try:if sgd_iteration_count % update_frequency == 0:train_set_rmse = get_trainset_rmse(model)valid_set_rmse = get_validset_rmse(model)test_set_rmse = get_testset_rmse(model)train_errors.append(train_set_rmse)valid_errors.append(valid_set_rmse)test_errors.append(test_set_rmse)print("Iteration:   ", sgd_iteration_count)print("Validation RMSE:", valid_set_rmse)if valid_set_rmse < best_validation_rmse:print('Test RMSE    :', test_set_rmse)print('Best validation error up to now !')patience = 0best_validation_rmse = valid_set_rmsebest_parameters = model.save_parameters()else:patience += 1if patience >= 20:print("Exceed patience for optimization, stopping !")breakprint()training_idx = sample_random_training_index()user, item, rating = train_ratings[training_idx]model.update_parameters(user, item, rating)sgd_iteration_count += 1except KeyboardInterrupt:print('Stopped Optimization')print('Current valid set performance=%s'% compute_rmse(model.predict(*valid_ratings.T[:2]),valid_ratings[:, 2]))print('Current test set performance=%s' % compute_rmse(model.predict(*test_ratings.T[:2]),test_ratings[:, 2]))breakmodel.load_parameters(best_parameters)stop_time = time.time()print('Optimization time: ', (stop_time - start_time) / 60., 'minutes')

输出结果，略！

我们应该等待多少次迭代？

模型建立的怎么样？

优化算法的表现如何？

通过学习曲线，我们尝试回答以上问题中的部分问题。

x = update_frequency * np.arange(len(train_errors))
fig = plt.figure(num=None, figsize=(10, 10), dpi=500)
plt.plot(x, train_errors, 'r--', x, valid_errors, 'bs', x, test_errors, 'g^')
plt.legend(['training error', 'validation error', 'testing error'])plt.show()

可以看到大部分的学习在优化的早期就已经完成了，但是多等一会以后就能够得到更好的结果。

Training error总是比validation error 和 testing error 要低，当差别过大的时候，就会产生过拟合（overfitting）。

》较大的embedding size就会引起较严重的过拟合

》较小的embedding size则过拟合现象比较轻，但是在测试集上的表现上就比较差了。最好的embedding size在中间的某处！

还有很多其他的方法来提高模型的表现和泛化能力：

》正则化：使得embedding向量的值较小，从而减轻过拟合。给 user/item bias 和 embedding不同的正则化参数。

》学习率迭代下降：随着优化的迭代次数增加，降低学习率。

》Implicit feedback : use who rated what information, without the rating value（这个没看懂是做什么）。

》调整超参数 (learning rate, embedding size, regularization, ...)

》等等

评价

现在我们有了一个模型，那我们如何来评价这个模型，更重要的是如何来使用这个模型？

# Evaluation !!!
test_predictions = model.predict(*test_ratings.T[:2])
print(test_predictions)
test_df = pd.DataFrame({'userId': test_ratings[:, 0],'movieId': test_ratings[:, 1],'rating': test_ratings[:, 2],'prediction': test_predictions})
print(test_df.head())

测试集的全局性能的标准错误度量标准

# Standard error metrics for global performance of test set
print('Root Mean Squared Error\t\t', 
((test_df.rating - test_df.prediction) ** 2).mean() ** 0.5)
print('Mean Absolute Error\t\t', (test_df.rating - test_df.prediction).abs().mean())
print('Mean Absolute Percentile Error\t\t',100 * ((test_df.rating - test_df.prediction).abs() / test_df.rating).mean(),'%')

输出：

Root Mean Squared Error		0.856369783573
Mean Absolute Error		0.661424603843
Mean Absolute Percentile Error	28.5558836737 %

查看一下每个真实的评分对应的模型评分的平均值

plt.plot(*test_df.groupby('rating').prediction.mean().reset_index().values.T)
plt.plot(np.arange(0, 6), np.arange(0, 6), '--')
plt.xlim([0.5, 5])
plt.ylim([0.5, 5])
plt.xlabel('True Rating')
plt.ylabel('Mean Predicted Rating')
plt.show()

作者说模型的评分比较保守，意思应该是对于一个电影，不会给出过高或者过低的分数；这是正常现象！

每个用户实际给出的最高分与预测的最高分的比较：

best_predicted_rating_per_user = test_df.groupby('userId')
.apply(lambda x: x.sort_values('prediction', ascending=False).head(1).rating)
best_rating_per_user = test_df.groupby('userId')
.apply(lambda x: x.sort_values('rating', ascending=False).head(1).rating)
#求平均值
print('Best rating per user\t\t', best_rating_per_user.mean())
print('Best predicted rating per user\t', best_predicted_rating_per_user.mean())

输出：

Best rating per user 4.726386806596702
Best predicted rating per user 4.149175412293853

可以看到真实的平均用户最高分和预测的平均用户最高分的差别还是有的。

你还可以自己想出很多很多的性能指标来衡量！

但是最终，最重要的还是在实际应用中的表现！看看在实际应用中我们的模型表现如何！

一种衡量隐藏属性的方法就是放在一起比较！这里拿电影为例，计算两个电影的隐藏属性之间的相似度，然后人来判断这两个电影是否很相似，从而判断模型的表现！

提取电影的影藏属性为<id,embedding>的dict，然后定义两种相似度的度量

movies_embeddings = dict([(i, model.V[i]) for i in movies.index.values])def compute_cosine_similarity(movieId):movie_embedding = movies_embeddings[movieId]movie_embedding_norm = (movie_embedding ** 2).sum() ** 0.5similarity = dict([(movie,((movie_embedding * emb).sum())/ (((emb ** 2).sum() ** 0.5) * movie_embedding_norm))for (movie, emb) in movies_embeddings.items()])return similaritydef compute_euclidian_similarity(movieId):movie_embedding = movies_embeddings[movieId]similarity = dict([(movie, -((movie_embedding - emb) ** 2).sum() ** 0.5)for (movie, emb) in movies_embeddings.items()])return similarity

给定某个电影，查看距离最近的和距离最远的某几个电影

movie_id = 79132  # inception
sorted_movies = sorted(compute_cosine_similarity(movie_id).items(), key=lambda xx: xx[1])print('Closet')
for i in range(1, 16):l_id, sim = sorted_movies[-i]print(i, movies.loc[l_id].title, '\t', movies.loc[l_id].genres, sim)print()print('Farthest')
for i in range(15, -1, -1):l_id, sim = sorted_movies[i]print(i, movies.loc[l_id].title)print('\t', movies.loc[l_id].genres, sim)print()

输出结果：略。总体看来，还是靠谱的!