9. 吴恩达机器学习-推荐系统

1. 理论原理

1. 问题规划

  在机器学习中，针对一些问题，有算法可以为系统自动学习一套好的特征。因此，不要试图手动设计，但手写代码是我们目前为止常干的事。有一些设置，我们可以有一个算法，仅仅学习其使用的特征，推荐系统就是类型设置的一个例子。
  我们从一个例子开始定义推荐系统的问题。假设我们是一个电影供应商，我们有 5 部电影和 4 个用户，我们要求用户为电影打分。

  我们可以看出前三部是爱情片，后两部则是动作片，Alice 和 Bob 似乎更倾向于爱情片，而 Carol 和 Dave 则更倾向于动作片，并且没有一个用户给所有的电影都打过分。我们希望构建一个算法来预测他们每个人可能会给他们没看过的电影打多少分，并以此作为推荐的依据。
  通过推荐算法，预测结果可能如下图所示：

2. 基于内容的推荐算法

  在一个基于内容的推荐系统算法中，我们假设对于我们希望推荐的东西有一些数据，这些数据是有关这些东西的特征。在电影预测评分的例子中，我们假设每部电影都有两个特征，如 $x_1$ 代表电影的浪漫程度，$x_2$ 代表电影的动作程度，则每部电影都有一个特征向量，如 $x_1$ 是第一部电影的特征向量为 $[0.9,0]$。

  下面我们要基于这些特征来构建一个推荐系统算法。假设我们采用线性回归模型，我们可以针对每一个用户都训练一个线性回归模型，如 ${\theta }^{(1)}$ 是第一个用户的模型的参数。于是，我们有：${\theta }^{(j)}$ 表示用户 $j$ 的参数向量，$x^{(i)}$ 表示电影 $i$ 的特征向量，对于用户 $j$ 和电影 $i$，我们预测评分为：$({\theta }^{(j)}{)}^T{x}^{(i)}$。例如：$x^{(3)}=\left[\begin{array}{c}1\\ 0.99\\ 0\end{array}\right]$，${\theta }^{(1)}=\left[\begin{array}{c}0\\ 5\\ 0\end{array}\right]$，则第一个用户对第三部电影的评分预测为：$({\theta }^{(1)}{)}^T{x}^{(3)}=0\times 1+5\times 0.99+0\times 0=4.95$。针对用户 $j$，该线性回归模型的代价为预测误差的平方和，加上正则化项：$$\underset{{\theta }^{(j)}}{\min }\frac{1}{2}\sum _{i:r(i,j)=1}(({\theta }^{(j)}{)}^T{x}^{(i)}-y^{(i,j)}{)}^2+\frac{\lambda }{2}\sum _{k=1}^n({\theta }_k^{(j)}{)}^2$$其中，${\sum }_{i:r(i,j)=1}$表示我们只对那些用户 $j$ 评分过的电影进行求和。在一般的线性回归模型中，误差项和正则项应该都是乘以 $\frac{1}{2m}$，在这里我们将 $m$ 去掉，并且我们不对 ${\theta }_{(0)}$ 项进行正则化处理。
上面的代价函数只是针对一个用户，为了学习所有用户，我们将所有用户的代价函数求和：$$\underset{{\theta }^{(1)},{\theta }^{(2)},\dots ,{\theta }^{(n_u)}}{\min }\frac{1}{2}\sum _{j=1}^{n_u}\sum _{i:r(i,j)=1}(({\theta }^{(j)}{)}^T{x}^{(i)}-y^{(i,j)}{)}^2+\frac{\lambda }{2}\sum _{j=1}^{n_u}\sum _{k=1}^n({\theta }_k^{(j)}{)}^2$$如果我们要用梯度下降法来求解最优解，我们计算代价函数的偏导数后得到梯度下降的更新公式：$${\theta }_k^{(j)}:={\theta }_k^{(j)}-\alpha \sum _{k=1}^n(({\theta }^{(j)}{)}^T{x}^{(i)}-y^{(i,j)})x_k^{(i)}(fork=0)$$$${\theta }_k^{(j)}:={\theta }_k^{(j)}-\alpha \left(\sum _{k=1}^n\right(({\theta }^{(j)}{)}^T{x}^{(i)}-y^{(i,j)})x_k^{(i)}+\lambda {\theta }_k^{(j)})(fork\ne 0)$$

3. 协同过滤

  在之前的基于内容的推荐系统中，对于每一部电影，我们都掌握了可用的特征，使用这些特征训练出了每一个用户的参数。相反地，如果我们拥有用户的参数，我们可以学习得出电影的特征。但是如果当我们既没有用户的参数，又没有电影的特征时，这两种方法就不可行了。协同过滤算法可以同时学习这两者。
优化目标算法为：
  给定 ${\theta }^{(1)},\dots ,{\theta }^{(n_u)}$，去学习 $x^{(1)},\dots ,x^{(n_m)}$：$$\underset{x^{(1)},\dots ,x^{(n_m)}}{\min }\frac{1}{2}\sum _{i=1}^{n_m}\sum _{i:r(i,j)=1}(({\theta }^{(j)}{)}^T{x}^{(i)}-y^{(i,j)}{)}^2+\frac{\lambda }{2}\sum _{i=1}^{n_m}\sum _{k=1}^n(x_k^{(i)}{)}^2$$协同滤波算法为：
  ①给定 $x^{(1)},\dots ,x^{(n_m)}$，估计 ${\theta }^{(1)},\dots ,{\theta }^{(n_u)}$
  ②给定 ${\theta }^{(1)},\dots ,{\theta }^{(n_u)}$，估计 $x^{(1)},\dots ,x^{(n_m)}$
  ③通过迭代 $\theta \to x\to \theta \to x\to \theta \dots$ 从而得到更好的 $\theta ，x$，并且算法会收敛到一组好的用户参数及合理的特征。
协同过滤算法使用步骤：
  ①初始 $x^{(1)},\dots ,x^{(n_m)},{\theta }^{(1)},\dots ,{\theta }^{(n_u)}$ 为一些随机小值
  ②使用梯度下降算法最小化代价函数
  ③在训练完算法后，我们预测 $({\theta }^{(j)}{)}^T{x}^{(i)}$ 为用户 $j$ 给电影 $i$ 的评分
通过这个学习过程获得的特征矩阵包含了有关电影的重要数据，这些数据不总是人能读懂的，但是我们可以用这些数据作为给用户推荐电影的依据。例如：如果一位用户正在观看电影 $x^{(i)}$，我们可以寻找另一部电影 $x^{(j)}$，依据两部电影的特征向量之间的距离 $||x^{(i)}-x^{(j)}||$ 的大小。

4. 协同过滤算法

协同过滤的优化目标：
  ①给定 $x^{(1)},\dots ,x^{(n_m)}$，估计 ${\theta }^{(1)},\dots ,{\theta }^{(n_u)}$：$$\underset{{\theta }^{(1)},\dots ,{\theta }^{(n_u)}}{\min }\frac{1}{2}\sum _{j=1}^{n_u}\sum _{i:r(i,j)=1}(({\theta }^{(j)}{)}^T{x}^{(i)}-y^{(i,j)}{)}^2+\frac{\lambda }{2}\sum _{i=1}^{n_u}\sum _{k=1}^n({\theta }_k^{(i)}{)}^2$$  ②给定 ${\theta }^{(1)},\dots ,{\theta }^{(n_u)}$，估计 $x^{(1)},\dots ,x^{(n_m)}$：$$\underset{x^{(1)},\dots ,x^{(n_m)}}{\min }\frac{1}{2}\sum _{i=1}^{n_m}\sum _{i:r(i,j)=1}(({\theta }^{(j)}{)}^T{x}^{(i)}-y^{(i,j)}{)}^2+\frac{\lambda }{2}\sum _{i=1}^{n_m}\sum _{k=1}^n(x_k^{(i)}{)}^2$$  ③同时最小化 $x^{(1)},\dots ,x^{(n_m)},{\theta }^{(1)},\dots ,{\theta }^{(n_u)}$：$$J(x^{(1)},\dots ,x^{(n_m)},{\theta }^{(1)},\dots ,{\theta }^{(n_u)})=\frac{1}{2}\sum _{(i,j):r(i,j)=1}(({\theta }^{(j)}{)}^T{x}^{(i)}-y^{(i,j)}{)}^2+\frac{\lambda }{2}\sum _{i=1}^{n_u}\sum _{k=1}^n({\theta }_k^{(i)}{)}^2+\frac{\lambda }{2}\sum _{i=1}^{n_m}\sum _{k=1}^n(x_k^{(i)}{)}^2$$ $$\underset{x^{(1)},\dots ,x^{(n_m)},{\theta }^{(1)},\dots ,{\theta }^{(n_u)}}{\min }J(x^{(1)},\dots ,x^{(n_m)},{\theta }^{(1)},\dots ,{\theta }^{(n_u)})$$如果我们要用梯度下降法来求解最优解，我们计算代价函数的偏导数后得到梯度下降的更新公式为：$$x_k^{(i)}:=x_k^{(i)}-\alpha \left(\sum _{j:r(i,j)=1}\right(({\theta }^{(j)}{)}^T{x}_k^{(i)}-y^{(i,j)}){\theta }_k^{(j)}+\lambda x_k^{(i)})$$ $${\theta }_k^{(j)}:={\theta }_k^{(j)}-\alpha \left(\sum _{i:r(i,j)=1}\right(({\theta }^{(j)}{)}^T{x}_k^{(i)}-y^{(i,j)}){\theta }_k^{(j)}+\lambda {\theta }_k^{(j)})$$

5. 矢量化：低秩矩阵分解

  接下来介绍有关协同滤波算法的向量化实现，以及该算法可以做的其他事情。例如
    1. 当给出一件产品时，能否找到与之相关的其他产品
    2. 一位用户最近看上一件产品，有没有其他相关的产品可以推荐给他
  我们要做的是，实现另一种方法写出协同过滤算法的预测情况。我们有关于五部电影的数据集，将这些用户的电影评分，进行分组并存到一个矩阵中。

  我们有五部电影和四位用户，那么这个矩阵 $Y$ 就是一个 5 行 4 列的矩阵，它将这些电影的用户评分数据都存在矩阵里，如下所示：$$Y=\left[\begin{array}{cccc}5& 5& 0& 0\\ 5& ?& ?& 0\\ ?& 4& 0& ?\\ 0& 0& 5& 4\\ 0& 0& 5& 0\end{array}\right]$$  预测出评分为：$$\left[\begin{array}{cccc}({\theta }^{(1)}{)}^T{x}^{(1)}& ({\theta }^{(2)}{)}^T{x}^{(1)}& \dots & ({\theta }^{(n_u)}{)}^T{x}^{(1)}\\ ({\theta }^{(1)}{)}^T{x}^{(2)}& ({\theta }^{(2)}{)}^T{x}^{(2)}& \dots & ({\theta }^{(n_u)}{)}^T{x}^{(2)}\\ ⋮& ⋮& \ddots & ⋮\\ ({\theta }^{(1)}{)}^T{x}^{(n_m)}& ({\theta }^{(2)}{)}^T{x}^{(n_m)}& \dots & ({\theta }^{(n_u)}{)}^T{x}^{(n_m)}\end{array}\right]$$  找到相关推荐影片：当用户在看某部电影 $i$ 的时候，为了能给用户推荐 5 部新电影，我们需要做的是找出 5 部与电影 $i$ 非常相似的电影$j$。电影 $i$ 有一个特征向量 $x^{(i)}$，那么我们只要保证两部电影的特征向量 $x^{(i)}$ 和 $x^{(j)}$ 之间的距离 $$||x^{(i)}-x^{(j)}||$$ 很小，那就能很有力地表明电影 $i$ 和电影 $j$ 在某种程度上有相似，这样你就能给你的用户推荐几部不同的电影了。

6. 实施细节：均值规范化

  如果我们在上述预测电影评分例子的基础上新增一个用户 Eve，并且 Eve 没有为任何电影评分，那么我们以什么为依据为 Eve 推荐电影呢？

  我们首先需要对结果矩阵 $Y$ 进行均值归一化处理，将每一个用户对某一部电影的评分减去所有用户对该电影评分的平均值：$$Y=\left[\begin{array}{ccccc}5& 5& 0& 0& ?\\ 5& ?& ?& 0& ?\\ ?& 4& 0& ?& ?\\ 0& 0& 5& 4& ?\\ 0& 0& 5& 0& ?\end{array}\right],\mu =\left[\begin{array}{c}2.5\\ 2.5\\ 2\\ 2.25\\ 1.25\end{array}\right]\to Y=\left[\begin{array}{ccccc}2.5& 2.5& -2.5& -2.5& ?\\ 2.5& ?& ?& -2.5& ?\\ ?& 2& -2& ?& ?\\ -2.25& -2.25& 2.75& 1.75& ?\\ -2.25& -1.25& 3.75& -1.25& ?\end{array}\right]$$
  然后我们利用这个新的 $Y$ 矩阵来训练算法。 如果我们要用新训练出的算法来预测评分，则需要将平均值重新加回去，预测用户 $j$ 给电影 $i$ 的评分为 $(({\theta }^{(j)}){)}^T{x}^{(i)}+{\mu }_i$。对于 Eve，我们的新模型会认为她给每部电影的评分都是该电影的平均分。

2. 代码实现

import matplotlib.pyplot as plt
import seaborn as sns
sns.set(context="notebook", style="white", palette=sns.color_palette("RdBu"))
import numpy as np
import pandas as pd
import scipy.io as sio

加载数据并建立数据集
  $X$：电影特征数据，其维度为：(1682, 10)
  $Theta$：用户特征数据，其维度为：(943, 10)
  $Y$：用户对电影的评分
  $R$：若第 $i$ 部电影由第 $j$ 个用户评分，则：$R(i,j)=1$

movies_mat = sio.loadmat('../data/ex8_movies.mat')
Y, R = movies_mat.get('Y'), movies_mat.get('R')

m, u = Y.shape      # m：电影数    u:用户数
n = 10   # 电影的特征数

param_mat = sio.loadmat('../data/ex8_movieParams.mat')
theta, X = param_mat.get('Theta'), param_mat.get('X')

损失函数

# 将 X 和 theta 展开为一维数组并合并
def serialize(X, theta):
    return np.concatenate((X.ravel(), theta.ravel()))

# 将数据还原为原来的电影数据和用户数据
def deserialize(param, n_movie, n_user, n_features):
    return param[:n_movie * n_features].reshape(n_movie, n_features),param[n_movie * n_features:].reshape(n_user, n_features)
    
# 损失函数
def cost(param, Y, R, n_features):
    n_movie, n_user = Y.shape
    X, theta = deserialize(param, n_movie, n_user, n_features)
    
    inner = np.multiply(X @ theta.T - Y, R)
    
    return np.power(inner, 2).sum() / 2

# 梯度下降
def gradient(param, Y, R, n_features):
    n_movies, n_user = Y.shape
    X, theta = deserialize(param, n_movies, n_user, n_features)
    
    inner = np.multiply(X @ theta.T - Y, R)
    
    X_grad = inner @ theta
    
    theta_grad = inner.T @ X
    
    return serialize(X_grad, theta_grad)

# 正则化损失函数
def regularized_cost(param, Y, R, n_features, l=1):
    reg_term = np.power(param, 2).sum() * (l / 2)
    
    return cost(param, Y, R, n_features) + reg_term

# 正则化梯度下降
def regularized_gradient(param, Y, R, n_features, l=1):
    grad = gradient(param, Y, R, n_features)
    reg_term = l * param
    
    return grad + reg_term

# 使用子数据集计算损失函数
users = 4
movies = 5
features = 3

X_sub = X[:movies, :features]
theta_sub = theta[:users, :features]
Y_sub = Y[:movies, :users]
R_sub = R[:movies, :users]

param_sub = serialize(X_sub, theta_sub)
cost(param_sub, Y_sub, R_sub, features)         # 22.224603725685675

# 全部的真参数
param = serialize(X, theta)
# 真正的总损失
cost(serialize(X, theta), Y, R, 10)     # 27918.64012454421

梯度下降

n_movie, n_user = Y.shape

X_grad, theta_grad = deserialize(gradient(param, Y, R, 10), n_movie, n_user, 10)

# 判断一个表达式，在表达式条件为 false 的 时候触发异常
assert X_grad.shape == X.shape 
assert theta_grad.shape == theta.shape

正则化损失函数

# 在 ex8_config.m 数据里面，lambda=1.5，并且使用子数据集
regularized_cost(param_sub, Y_sub, R_sub, features, l=1.5)      # 31.34405624427422

# 总正则化损失
regularized_cost(param, Y, R, 10, l=1)       # 32520.682450229557

正则化梯度

n_movie, n_user = Y.shape

X_grad, theta_grad = deserialize(regularized_gradient(param, Y, R, 10), n_movie, n_user, 10)

assert X_grad.shape == X.shape
assert theta_grad.shape == theta.shape

解析 movie_id.txt

movie_list = []

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

重现评分

ratings = np.zeros(1682)

ratings[0] = 4
ratings[6] = 3
ratings[11] = 5
ratings[53] = 4
ratings[63] = 5
ratings[65] = 3
ratings[68] = 5
ratings[97] = 2
ratings[182] = 4
ratings[225] = 5
ratings[354] = 5

准备数据

Y, R = movies_mat.get('Y'), movies_mat.get('R')

# 现在我变成了user 0
Y = np.insert(Y, 0, ratings, axis=1)
Y.shape       # (1682, 944)

R = np.insert(R, 0, ratings != 0, axis=1)
R.shape   # (1682, 944)

n_features = 50
n_movie, n_user = Y.shape
l = 10

X = np.random.standard_normal((n_movie, n_features))
theta = np.random.standard_normal((n_user, n_features))

X.shape, theta.shape     # ((1682, 50), (944, 50))

param = serialize(X, theta)

# 对评分进行归一化处理
Y_norm = Y - Y.mean()
Y_norm.mean()        # 4.6862111343939375e-17

训练

import scipy.optimize as opt

res = opt.minimize(fun=regularized_cost,
                   x0=param,
                   args=(Y_norm, R, n_features, l),
                   method='TNC',
                   jac=regularized_gradient)

X_trained, theta_trained = deserialize(res.x, n_movie, n_user, n_features)
X_trained.shape, theta_trained.shape     # ((1682, 50), (944, 50))

X_trained[:10]

prediction = X_trained @ theta_trained.T

my_preds = prediction[:, 0] + Y.mean()
idx = np.argsort(my_preds)[::-1]
idx.shape      # (1682,)

my_preds[idx][:10]

for m in movie_list[idx][:10]:
    print(m)

# Titanic (1997)
# Star Wars (1977)
# Shawshank Redemption, The (1994)
# Forrest Gump (1994)
# Raiders of the Lost Ark (1981)
# Braveheart (1995)
# Return of the Jedi (1983)
# Usual Suspects, The (1995)
# Godfather, The (1972)
# Schindler's List (1993)