用于个性化推荐的矩阵分解概率模型#

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import xarray as xr

%config InlineBackend.figure_format = 'retina'
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")

动机#

所以你在Netflix上浏览想看的内容，但就是不喜欢这些推荐。你觉得自己可以做得更好。你只需要收集一些自己和朋友的评分数据，并构建一个推荐算法。这个笔记本将指导你完成这个过程！

我们将首先通过一些直觉来了解我们的模型将如何工作。然后我们将形式化我们的直觉。之后，我们将检查我们将使用的数据集。一旦我们对数据的样子有了一些概念，我们将定义一些基线方法来预测电影的偏好。接着，我们将研究概率矩阵分解（PMF），这是一种更复杂的贝叶斯方法，用于预测偏好。详细描述了PMF模型后，我们将使用PyMC进行MAP估计和MCMC推断。最后，我们将比较使用PMF获得的结果与基线方法获得的结果，并讨论结果。

直觉#

通常，如果我们想要某方面的推荐，我们会尝试找到与我们相似的人并询问他们的意见。如果Bob、Alice和Monty都与我相似，并且他们都喜欢犯罪剧，我可能会喜欢犯罪剧。当然，这并不总是正确的。这取决于我们如何定义“相似”。为了获得最佳的性价比，我们真正想要寻找的是那些口味最相似的人。口味是一个复杂的概念，我们可能希望将其分解为更易理解的部分。我们可能会尝试根据各种因素来描述每部电影。也许电影可以是忧郁的、轻松的、电影化的、对话密集的、大预算的等等。现在想象一下，我们浏览IMDB并为每部电影在每个类别中分配一个评分。它有多忧郁？它有多少对话？它的预算是什么？也许我们为每个类别使用0到1之间的数字。直观地说，我们可能会称这为电影的“轮廓”。

现在让我们假设我们回到那5部我们评分的电影。此时，通过查看我们喜欢和不喜欢的每部电影的影片简介，我们可以更全面地了解自己的偏好。也许我们对这5部电影的简介取平均值，并称之为我们理想的电影类型。换句话说，我们已经计算出了我们对各种类型电影的固有偏好。假设Bob、Alice和Monty都做了同样的事情。现在我们可以比较我们的偏好，并确定我们每个人之间的相似程度。我可能会发现Bob是最相似的，而其他两个人仍然比其他人更相似，但不如Bob那么相似。所以我希望从这三个人那里得到推荐，但当我做出最终决定时，我会更重视Bob的推荐，而不是我从Alice和Monty那里得到的推荐。

虽然上述过程听起来相当有效，但它也揭示了一个意想不到的额外信息来源。如果我们对某部电影评价很高，并且我们知道它的电影简介，我们可以与其他电影的简介进行比较。如果我们发现一部电影的数字非常接近，我们也很可能喜欢这部电影。这两种方法都被称为邻域方法。同时利用这两种方法的技术通常被称为协同过滤 [Koren 等, 2009]。我们讨论的第一种方法使用用户-用户相似性，而第二种方法使用项目-项目相似性。理想情况下，我们希望同时使用这两种信息来源。我们的想法是我们有很多项目可供选择，我们希望与其他人合作，将项目列表过滤到我们每个人最喜欢的项目。我的列表应该将我最喜欢的项目放在顶部，而将我最不喜欢的项目放在底部。其他人都希望同样的事情。如果我与其他人一起观看5部电影，并且我们有一个有效的计算过程来确定相似性，我们可以非常快速地按照我们的喜好对电影进行排序。

形式化#

让我们花点时间将我们一直在讨论的直观概念变得更加具体。我们有一组 \(M\) 部电影，或者称为项目（在上面的例子中 \(M = 100\)）。我们还有 \(N\) 个人，我们称他们为我们推荐系统的用户。对于每个项目，我们希望找到一个 \(D\) 维的因子组合（如上面的电影简介）来描述该项目。理想情况下，我们希望在不实际逐一手动标记所有电影的情况下完成这项工作。手动标记既慢又容易出错，因为不同的人可能会以不同的方式标记电影。因此，我们将每部电影建模为一个 \(D\) 维向量，即其潜在因子组合。此外，我们期望每个用户都有一些偏好，但在没有手动标记和平均过程的情况下，我们必须依赖潜在因子组合来学习每个用户的 \(D\) 维潜在偏好向量。我们唯一能观察到的是用户提供的 \(N \times M\) 评分矩阵 \(R\)。条目 \(R_{ij}\) 是用户 \(i\) 对项目 \(j\) 的评分。由于大多数用户不会对所有100部电影进行评分，因此这些条目中可能有许多是缺失的。我们的目标是用基于潜在变量 \(U\) 和 \(V\) 的预测评分来填补这些缺失值。我们用 \(R_{ij}^*\) 表示预测评分。我们还定义了一个指示矩阵 \(I\)，如果 \(R_{ij}\) 缺失，则条目 \(I_{ij} = 0\)，否则 \(I_{ij} = 1\)。

所以我们有一个\(N \times D\)的用户偏好矩阵，我们称之为\(U\)，还有一个\(M \times D\)的因子组合矩阵，我们称之为\(V\)。我们还有一个\(N \times M\)的评分矩阵，我们称之为\(R\)。我们可以将每一行\(U_i\)视为每个用户对每个\(D\)潜在因子的偏好程度。每一行\(V_j\)可以被视为每个项目由每个潜在因子描述的程度。为了做出推荐，我们需要一个合适的预测函数，该函数将用户偏好向量\(U_i\)和项目潜在因子向量\(V_j\)映射到一个预测的排名。这个预测函数的选择是一个重要的建模决策，并且已经使用了各种预测函数。也许最常见的是两个向量的点积，\(U_i \cdot V_j\) [Koren 等, 2009]。

为了更好地理解协同过滤技术，让我们探讨一个具体的例子。假设我们正在使用一个模型来推荐电影，该模型推断出五个潜在因素，\(V_j\)，对于\(j = 1,2,3,4,5\)。实际上，这些潜在因素通常无法以直接的方式解释，大多数模型也不会试图理解每个因素所捕获的信息。然而，为了解释的目的，我们假设这五个潜在因素可能最终捕获了我们上面讨论的电影特征。因此，我们的五个潜在因素是：情绪化、轻松愉快、电影感、对话和预算。然后，对于特定用户\(i\)，假设我们推断出一个偏好向量\(U_i = <0.5, 0.1, 1.5, 1.1, 0.3>\)。同样，对于特定项目\(j\)，我们推断出这些潜在因素的值：\(V_j = <0.5, 1.5, 1.25, 0.8, 0.9>\)。使用点积作为预测函数，我们将计算出3.425作为该项目的排名，这在我们的1到5评分尺度上大致是一个中性偏好。

\[ 0.5 \times 0.5 + 0.1 \times 1.5 + 1.5 \times 1.25 + 1.1 \times 0.8 + 0.3 \times 0.9 = 3.425 \]

数据#

MovieLens 100k 数据集 [Harper 和 Konstan, 2016] 由明尼苏达大学的 GroupLens 研究项目收集。该数据集包含 943 名用户对 1682 部电影的 100,000 条评分（1-5 分）。每个用户至少对 20 部电影进行了评分，并且我们拥有用户的基本信息（年龄、性别、职业、邮政编码）。每部电影包括基本信息，如标题、发行日期、视频发行日期和类型。我们将在该数据上实现一个适合协同过滤的模型，并根据均方根误差（RMSE）对其进行评估，以验证结果。

数据是通过MovieLens网站在1997年9月19日至1998年4月22日的七个月期间收集的。这些数据已经过清理——少于20条评分或没有完整人口统计信息的用户已从该数据集中删除。

让我们从探索数据开始。我们想要对数据的总体情况有一个大致的了解，并对其可能包含的模式有所感知。以下是用户评分数据：

data_kwargs = dict(sep="\t", names=["userid", "itemid", "rating", "timestamp"])
try:
    data = pd.read_csv("../data/ml_100k_u.data", **data_kwargs)
except FileNotFoundError:
    data = pd.read_csv(pm.get_data("ml_100k_u.data"), **data_kwargs)

data.head()

	userid	itemid	rating	timestamp
0	196	242	3	881250949
1	186	302	3	891717742
2	22	377	1	878887116
3	244	51	2	880606923
4	166	346	1	886397596

这里是电影详细数据：

# fmt: off
movie_columns  = ['movie id', 'movie title', 'release date', 'video release date', 'IMDb URL', 
                  'unknown','Action','Adventure', 'Animation',"Children's", 'Comedy', 'Crime',
                  'Documentary', 'Drama', 'Fantasy', 'Film-Noir', 'Horror', 'Musical', 'Mystery',
                  'Romance', 'Sci-Fi', 'Thriller', 'War', 'Western']
# fmt: on

item_kwargs = dict(sep="|", names=movie_columns, index_col="movie id", parse_dates=["release date"])
try:
    movies = pd.read_csv("../data/ml_100k_u.item", **item_kwargs)
except FileNotFoundError:
    movies = pd.read_csv(pm.get_data("ml_100k_u.item"), **item_kwargs)

movies.head()

	movie title	release date	video release date	IMDb URL	unknown	Action	Adventure	Animation	Children's	Comedy	...	Fantasy	Film-Noir	Horror	Musical	Mystery	Romance	Sci-Fi	Thriller	War	Western
movie id
1	Toy Story (1995)	1995-01-01	NaN	http://us.imdb.com/M/title-exact?Toy%20Story%2...	0	0	0	1	1	1	...	0	0	0	0	0	0	0	0	0	0
2	GoldenEye (1995)	1995-01-01	NaN	http://us.imdb.com/M/title-exact?GoldenEye%20(...	0	1	1	0	0	0	...	0	0	0	0	0	0	0	1	0	0
3	Four Rooms (1995)	1995-01-01	NaN	http://us.imdb.com/M/title-exact?Four%20Rooms%...	0	0	0	0	0	0	...	0	0	0	0	0	0	0	1	0	0
4	Get Shorty (1995)	1995-01-01	NaN	http://us.imdb.com/M/title-exact?Get%20Shorty%...	0	1	0	0	0	1	...	0	0	0	0	0	0	0	0	0	0
5	Copycat (1995)	1995-01-01	NaN	http://us.imdb.com/M/title-exact?Copycat%20(1995)	0	0	0	0	0	0	...	0	0	0	0	0	0	0	1	0	0

5行 × 23列

# Plot histogram of ratings
data.groupby("rating").size().plot(kind="bar");

../_images/8fcc5a578bcc8708b918dcd0daf99848a665caf1c75bd49ec47df4a55097dd2c.png

data.rating.describe()

count    100000.000000
mean          3.529860
std           1.125674
min           1.000000
25%           3.000000
50%           4.000000
75%           4.000000
max           5.000000
Name: rating, dtype: float64

这一定是一批不错的电影。从我们上面的探索中，我们知道大多数评分都在3到5之间，并且正面评分比负面评分更有可能。让我们看看每部电影的平均评分，看看我们是否有特别好的（或差的）电影。

movie_means = data.join(movies["movie title"], on="itemid").groupby("movie title").rating.mean()
movie_means[:50].plot(kind="bar", grid=False, figsize=(16, 6), title="Mean ratings for 50 movies");

/Users/severinhatt/miniconda3/envs/pymc-hack/lib/python3.9/site-packages/IPython/core/pylabtools.py:151: UserWarning: constrained_layout not applied because axes sizes collapsed to zero.  Try making figure larger or axes decorations smaller.
  fig.canvas.print_figure(bytes_io, **kw)

../_images/69929f6218fea99518c31e01f880f137b5cb02b7291a6597dbd6eb4f9f1cc714.png

虽然大多数电影通常都能获得用户的好评，但确实有几部电影表现得特别差。让我们来看看最差和最好的电影，只是为了娱乐一下：

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(16, 4), sharey=True)
movie_means.nlargest(30).plot(kind="bar", ax=ax1, title="Top 30 movies in data set")
movie_means.nsmallest(30).plot(kind="bar", ax=ax2, title="Bottom 30 movies in data set");

/Users/severinhatt/miniconda3/envs/pymc-hack/lib/python3.9/site-packages/IPython/core/pylabtools.py:151: UserWarning: constrained_layout not applied because axes sizes collapsed to zero.  Try making figure larger or axes decorations smaller.
  fig.canvas.print_figure(bytes_io, **kw)

../_images/f4c7b9e474727a83ee28a3f0920e60fe48a8d76950dc942dde8a4fb73432728d.png

我明白了。我们现在知道这些电影之间确实存在明显的受欢迎程度差异。其中一些电影比其他电影更好，而有些则非常糟糕。通过观察电影评分，我们能够发现这些普遍趋势。也许在用户之间也存在类似的趋势。可能有些用户比其他人更容易被娱乐。让我们来看一看。

user_means = data.groupby("userid").rating.mean().sort_values()
_, ax = plt.subplots(figsize=(16, 6))
ax.plot(np.arange(len(user_means)), user_means.values, "k-")

ax.fill_between(np.arange(len(user_means)), user_means.values, alpha=0.3)
ax.set_xticklabels("")
# 1000 labels is nonsensical
ax.set_ylabel("Rating")
ax.set_xlabel(f"{len(user_means)} average ratings per user")
ax.set_ylim(0, 5)
ax.set_xlim(0, len(user_means));

../_images/2a3f02d521f5f6a8d895f7e14d46885d5fe5ef26709d695dc230a05c9e5378fa.png

我们在这里看到了更加显著的趋势。一些用户对几乎所有内容都给予高分，而另一些（尽管数量不多）则几乎对所有内容都给予低分。在考虑用于预测用户对未观看电影偏好的模型时，这些观察将非常有用。

方法#

在探索了数据之后，我们现在准备深入研究并开始解决问题。我们希望预测每个用户对所有他或她尚未阅读的电影的喜爱程度。

基线#

每个好的分析都需要一些基准方法来进行比较。如果我们没有定义“好”的参考点，就很难声称我们取得了好的结果。我们将定义三种非常简单的基准方法，并使用这些方法计算均方根误差（RMSE）。我们的目标是通过我们构建的任何模型获得更低的RMSE分数。

均匀随机基线#

我们的第一个基线方法简单得不能再简单了。每当我们在\(R\)中看到缺失值时，我们就会用一个在范围[1, 5]内均匀随机抽取的数字来填充它。我们预计这种方法的表现会是最差的。

\[R_{ij}^* \sim Uniform\]

全球平均基线#

这种方法仅比上一种方法稍好一些。无论哪里有缺失值，我们都会用所有观察到的评分的平均值来填充它。

\[\text{global_mean} = \frac{1}{N \times M} \sum_{i=1}^N \sum_{j=1}^M I_{ij}(R_{ij})\]

\[R_{ij}^* = \text{global_mean}\]

均值基线#

现在我们要开始变得更聪明一些。我们设想一些用户可能很容易被娱乐，倾向于给所有电影更高的评分。其他用户可能则相反。此外，一些电影可能比其他电影更机智，因此所有用户可能会普遍给某些电影更高的评分。我们可以在上面的电影平均评分图中清楚地看到这一点。我们将尝试通过每个用户和每个电影的评分平均值来捕捉这些总体趋势。我们还将引入全局平均值来稍微平滑一下。因此，如果我们看到单元格 \(R_{ij}\) 中的缺失值，我们将全局平均值与 \(U_i\) 的平均值和 \(V_j\) 的平均值进行平均，并使用该值来填充它。

\[\text{user_means} = \frac{1}{M} \sum_{j=1}^M I_{ij}(R_{ij})\]

\[\text{movie_means} = \frac{1}{N} \sum_{i=1}^N I_{ij}(R_{ij})\]

\[R_{ij}^* = \frac{1}{3} \left(\text{user_means}_i + \text{ movie_means}_j + \text{ global_mean} \right)\]

# Create a base class with scaffolding for our 3 baselines.


def split_title(title):
    """Change "BaselineMethod" to "Baseline Method"."""
    words = []
    tmp = [title[0]]
    for c in title[1:]:
        if c.isupper():
            words.append("".join(tmp))
            tmp = [c]
        else:
            tmp.append(c)
    words.append("".join(tmp))
    return " ".join(words)


class Baseline:
    """Calculate baseline predictions."""

    def __init__(self, train_data):
        """Simple heuristic-based transductive learning to fill in missing
        values in data matrix."""
        self.predict(train_data.copy())

    def predict(self, train_data):
        raise NotImplementedError("baseline prediction not implemented for base class")

    def rmse(self, test_data):
        """Calculate root mean squared error for predictions on test data."""
        return rmse(test_data, self.predicted)

    def __str__(self):
        return split_title(self.__class__.__name__)


# Implement the 3 baselines.


class UniformRandomBaseline(Baseline):
    """Fill missing values with uniform random values."""

    def predict(self, train_data):
        nan_mask = np.isnan(train_data)
        masked_train = np.ma.masked_array(train_data, nan_mask)
        pmin, pmax = masked_train.min(), masked_train.max()
        N = nan_mask.sum()
        train_data[nan_mask] = rng.uniform(pmin, pmax, N)
        self.predicted = train_data


class GlobalMeanBaseline(Baseline):
    """Fill in missing values using the global mean."""

    def predict(self, train_data):
        nan_mask = np.isnan(train_data)
        train_data[nan_mask] = train_data[~nan_mask].mean()
        self.predicted = train_data


class MeanOfMeansBaseline(Baseline):
    """Fill in missing values using mean of user/item/global means."""

    def predict(self, train_data):
        nan_mask = np.isnan(train_data)
        masked_train = np.ma.masked_array(train_data, nan_mask)
        global_mean = masked_train.mean()
        user_means = masked_train.mean(axis=1)
        item_means = masked_train.mean(axis=0)
        self.predicted = train_data.copy()
        n, m = train_data.shape
        for i in range(n):
            for j in range(m):
                if np.ma.isMA(item_means[j]):
                    self.predicted[i, j] = np.mean((global_mean, user_means[i]))
                else:
                    self.predicted[i, j] = np.mean((global_mean, user_means[i], item_means[j]))


baseline_methods = {}
baseline_methods["ur"] = UniformRandomBaseline
baseline_methods["gm"] = GlobalMeanBaseline
baseline_methods["mom"] = MeanOfMeansBaseline

num_users = data.userid.unique().shape[0]
num_items = data.itemid.unique().shape[0]
sparsity = 1 - len(data) / (num_users * num_items)
print(f"Users: {num_users}\nMovies: {num_items}\nSparsity: {sparsity}")

dense_data = data.pivot(index="userid", columns="itemid", values="rating").values

Users: 943
Movies: 1682
Sparsity: 0.9369533063577546

概率矩阵分解#

概率矩阵分解 [Mnih 和 Salakhutdinov, 2008] 是一种从贝叶斯角度解决协同过滤问题的概率方法。评分 \(R\) 被建模为从高斯分布中抽取的样本。\(R_{ij}\) 的均值是 \(U_i V_j^T\)。精度 \(\alpha\) 是一个固定参数，反映了估计的不确定性；正态分布通常以精度（方差的倒数）重新参数化。通过在 \(U\) 和 \(V\) 上放置零均值球形高斯先验来控制复杂性。换句话说，\(U\) 的每一行都从均值 \(\mu = 0\) 且精度为单位矩阵 \(I\) 的某个倍数的多变量高斯分布中抽取。这些倍数分别是 \(U\) 的 \(\alpha_U\) 和 \(V\) 的 \(\alpha_V\)。因此，我们的模型定义为：

\(\newcommand\given[1][]{\:#1\vert\:}\)

\[ P(R \given U, V, \alpha^2) = \prod_{i=1}^N \prod_{j=1}^M \left[ \mathcal{N}(R_{ij} \given U_i V_j^T, \alpha^{-1}) \right]^{I_{ij}} \]

\[ P(U \given \alpha_U^2) = \prod_{i=1}^N \mathcal{N}(U_i \given 0, \alpha_U^{-1} \boldsymbol{I}) \]

\[ P(V \given \alpha_U^2) = \prod_{j=1}^M \mathcal{N}(V_j \given 0, \alpha_V^{-1} \boldsymbol{I}) \]

给定小的精度参数，\(U\) 和 \(V\) 的先验确保我们的潜在变量不会偏离 0 太远。这防止了过于强烈的用户偏好和项目因子组合被学习到。这通常被称为复杂度控制，其中模型的复杂度在这里通过潜在变量的大小来衡量。像这样控制复杂度有助于防止过拟合，从而使模型能够更好地泛化到未见过的数据。我们还必须为 \(R\) 的正态分布选择一个合适的 \(\alpha\) 值。因此，挑战变成了为 \(\alpha_U\)、\(\alpha_V\) 和 \(\alpha\) 选择合适的值。这个挑战可以通过 Nowlan 和 Hinton [1992] 讨论的软权重共享方法来解决。然而，出于本次分析的目的，我们将坚持使用从我们的数据中获得的点估计。

import logging
import time

import pytensor
import scipy as sp

# Enable on-the-fly graph computations, but ignore
# absence of intermediate test values.
pytensor.config.compute_test_value = "ignore"

# Set up logging.
logger = logging.getLogger()
logger.setLevel(logging.INFO)


class PMF:
    """Probabilistic Matrix Factorization model using pymc."""

    def __init__(self, train, dim, alpha=2, std=0.01, bounds=(1, 5)):
        """Build the Probabilistic Matrix Factorization model using pymc.

        :param np.ndarray train: The training data to use for learning the model.
        :param int dim: Dimensionality of the model; number of latent factors.
        :param int alpha: Fixed precision for the likelihood function.
        :param float std: Amount of noise to use for model initialization.
        :param (tuple of int) bounds: (lower, upper) bound of ratings.
            These bounds will simply be used to cap the estimates produced for R.

        """
        self.dim = dim
        self.alpha = alpha
        self.std = np.sqrt(1.0 / alpha)
        self.bounds = bounds
        self.data = train.copy()
        n, m = self.data.shape

        # Perform mean value imputation
        nan_mask = np.isnan(self.data)
        self.data[nan_mask] = self.data[~nan_mask].mean()

        # Low precision reflects uncertainty; prevents overfitting.
        # Set to the mean variance across users and items.
        self.alpha_u = 1 / self.data.var(axis=1).mean()
        self.alpha_v = 1 / self.data.var(axis=0).mean()

        # Specify the model.
        logging.info("building the PMF model")
        with pm.Model(
            coords={
                "users": np.arange(n),
                "movies": np.arange(m),
                "latent_factors": np.arange(dim),
                "obs_id": np.arange(self.data[~nan_mask].shape[0]),
            }
        ) as pmf:
            U = pm.MvNormal(
                "U",
                mu=0,
                tau=self.alpha_u * np.eye(dim),
                dims=("users", "latent_factors"),
                initval=rng.standard_normal(size=(n, dim)) * std,
            )
            V = pm.MvNormal(
                "V",
                mu=0,
                tau=self.alpha_v * np.eye(dim),
                dims=("movies", "latent_factors"),
                initval=rng.standard_normal(size=(m, dim)) * std,
            )
            R = pm.Normal(
                "R",
                mu=(U @ V.T)[~nan_mask],
                tau=self.alpha,
                dims="obs_id",
                observed=self.data[~nan_mask],
            )

        logging.info("done building the PMF model")
        self.model = pmf

    def __str__(self):
        return self.name

我们还需要函数来计算MAP并在我们的PMF模型上进行采样。当观测噪声方差\(\alpha\)和先验方差\(\alpha_U\)和\(\alpha_V\)都保持固定时，最大化对数后验等价于最小化带有二次正则化项的平方误差和目标函数。

\[ E = \frac{1}{2} \sum_{i=1}^N \sum_{j=1}^M I_{ij} (R_{ij} - U_i V_j^T)^2 + \frac{\lambda_U}{2} \sum_{i=1}^N \|U\|_{Fro}^2 + \frac{\lambda_V}{2} \sum_{j=1}^M \|V\|_{Fro}^2, \]

其中 \(\lambda_U = \alpha_U / \alpha\), \(\lambda_V = \alpha_V / \alpha\), 和 \(\|\cdot\|_{Fro}^2\) 表示 Frobenius 范数 [Mnih 和 Salakhutdinov, 2008]。最小化这个目标函数得到一个局部最小值，这本质上是一个最大后验概率（MAP）估计。虽然可以使用快速随机梯度下降程序来找到这个 MAP，但我们将使用 pymc 中内置的工具来找到它。特别是，我们将使用 Powell 优化（scipy.optimize.fmin_powell）的 find_MAP。找到这个 MAP 估计后，我们可以将其用作 MCMC 采样的起点。

由于这是一个相当复杂的模型，我们预计MAP估计会花费一些时间。因此，让我们在找到它之后保存它。请注意，我们在下面定义了一个用于查找MAP的函数，假设它将接收一个包含一些变量的命名空间。然后我们将该函数附加到PMF类中，在该类初始化后，它将具有这样一个命名空间。PMF类是这样分段定义的，因此我可以在每段之间说几句话，使其更清晰。

def _find_map(self):
    """Find mode of posterior using L-BFGS-B optimization."""
    tstart = time.time()
    with self.model:
        logging.info("finding PMF MAP using L-BFGS-B optimization...")
        self._map = pm.find_MAP(method="L-BFGS-B")

    elapsed = int(time.time() - tstart)
    logging.info("found PMF MAP in %d seconds" % elapsed)
    return self._map


def _map(self):
    try:
        return self._map
    except:
        return self.find_map()


# Update our class with the new MAP infrastructure.
PMF.find_map = _find_map
PMF.map = property(_map)

所以现在我们的PMF类有一个map 属性，它将通过Powell优化找到或从之前的优化中加载。一旦我们有了MAP，我们就可以将其用作MCMC采样器的起点。我们需要一个采样函数来绘制MCMC样本，以近似PMF模型的后验分布。

# Draw MCMC samples.
def _draw_samples(self, **kwargs):
    # kwargs.setdefault("chains", 1)
    with self.model:
        self.trace = pm.sample(**kwargs)


# Update our class with the sampling infrastructure.
PMF.draw_samples = _draw_samples

我们可以像为MAP所做的那样定义某种默认的trace属性，但这意味着可能会使用对于nsamples和cores来说无意义的值。最好将其保留为对draw_samples的非可选调用。最后，我们需要一个函数来使用我们推断出的\(U\)和\(V\)值进行预测。对于用户\(i\)和电影\(j\)，预测是通过从\(\mathcal{N}(U_i V_j^T, \alpha)\)中抽取生成的。为了从采样器生成预测，我们为每个采样的\(U\)和\(V\)生成一个\(R\)矩阵，然后通过对\(K\)个样本进行平均来组合这些矩阵。

\[ P(R_{ij}^* \given R, \alpha, \alpha_U, \alpha_V) \approx \frac{1}{K} \sum_{k=1}^K \mathcal{N}(U_i V_j^T, \alpha) \]

我们希望在为了诊断目的对它们进行平均之前检查各个\(R\)矩阵。因此，我们将在评估期间编写用于平均部分的代码。下面的函数简单地根据给定的\(U\)和\(V\)以及PMF对象中存储的固定\(\alpha\)绘制一个\(R\)矩阵。

def _predict(self, U, V):
    """Estimate R from the given values of U and V."""
    R = np.dot(U, V.T)
    sample_R = rng.normal(R, self.std)
    # bound ratings
    low, high = self.bounds
    sample_R[sample_R < low] = low
    sample_R[sample_R > high] = high
    return sample_R


PMF.predict = _predict

最后一点需要注意的是：该模型中的点积通常使用逻辑函数 \(g(x) = 1/(1 + exp(-x))\) 进行约束，将预测值限制在范围 [0, 1] 内。为了便于这种限制，评分也使用 \(t(x) = (x + min) / range\) 映射到范围 [0, 1]。PMF 的作者还引入了一个约束版本，该版本在评分较少的用户上表现更好 [Salakhutdinov 和 Mnih, 2008]。这两种模型通常都是对这里介绍的基本模型的改进。然而，出于时间和空间的考虑，这里不会实现这些模型。

评估#

指标#

为了了解我们的模型有多有效，我们需要能够评估它们。我们将根据均方根误差（RMSE）进行评估，其形式如下：

\[ RMSE = \sqrt{ \frac{ \sum_{i=1}^N \sum_{j=1}^M I_{ij} (R_{ij} - R_{ij}^*)^2 } { \sum_{i=1}^N \sum_{j=1}^M I_{ij} } } \]

在这种情况下，RMSE 可以被认为是我们的预测与实际用户偏好之间的标准差。

# Define our evaluation function.
def rmse(test_data, predicted):
    """Calculate root mean squared error.
    Ignoring missing values in the test data.
    """
    I = ~np.isnan(test_data)  # indicator for missing values
    N = I.sum()  # number of non-missing values
    sqerror = abs(test_data - predicted) ** 2  # squared error array
    mse = sqerror[I].sum() / N  # mean squared error
    return np.sqrt(mse)  # RMSE

训练数据与测试数据#

接下来我们需要做的是将数据分为训练集和测试集。矩阵分解技术使用转导学习而非归纳学习。因此，我们通过从完整的\(N \times M\)数据矩阵中随机抽取单元格来生成测试集。在原始数据矩阵的副本中，被选为测试样本的值被替换为nan值，以生成训练集。由于我们将生成随机分割，我们还需要将生成的训练/测试集写出来。这将使我们能够复制结果。我们希望能够识别出哪个分割是哪个，因此我们将对选定的测试索引进行哈希处理，并使用它来保存数据。

# Define a function for splitting train/test data.
def split_train_test(data, percent_test=0.1):
    """Split the data into train/test sets.
    :param int percent_test: Percentage of data to use for testing. Default 10.
    """
    n, m = data.shape  # # users, # movies
    N = n * m  # # cells in matrix

    # Prepare train/test ndarrays.
    train = data.copy()
    test = np.ones(data.shape) * np.nan

    # Draw random sample of training data to use for testing.
    tosample = np.where(~np.isnan(train))  # ignore nan values in data
    idx_pairs = list(zip(tosample[0], tosample[1]))  # tuples of row/col index pairs

    test_size = int(len(idx_pairs) * percent_test)  # use 10% of data as test set
    train_size = len(idx_pairs) - test_size  # and remainder for training

    indices = np.arange(len(idx_pairs))  # indices of index pairs
    sample = rng.choice(indices, replace=False, size=test_size)

    # Transfer random sample from train set to test set.
    for idx in sample:
        idx_pair = idx_pairs[idx]
        test[idx_pair] = train[idx_pair]  # transfer to test set
        train[idx_pair] = np.nan  # remove from train set

    # Verify everything worked properly
    assert train_size == N - np.isnan(train).sum()
    assert test_size == N - np.isnan(test).sum()

    # Return train set and test set
    return train, test


train, test = split_train_test(dense_data)

结果#

# Let's see the results:
baselines = {}
for name in baseline_methods:
    Method = baseline_methods[name]
    method = Method(train)
    baselines[name] = method.rmse(test)
    print("{} RMSE:\t{:.5f}".format(method, baselines[name]))

Uniform Random Baseline RMSE:	1.68490
Global Mean Baseline RMSE:	1.11492
Mean Of Means Baseline RMSE:	1.00750

正如预期的那样：均匀随机基线是目前为止最差的，全局平均基线次之，而平均值方法是我们最好的基线。现在让我们看看PMF的表现如何。

# We use a fixed precision for the likelihood.
# This reflects uncertainty in the dot product.
# We choose 2 in the footsteps Salakhutdinov
# Mnihof.
ALPHA = 2

# The dimensionality D; the number of latent factors.
# We can adjust this higher to try to capture more subtle
# characteristics of each movie. However, the higher it is,
# the more expensive our inference procedures will be.
# Specifically, we have D(N + M) latent variables. For our
# Movielens dataset, this means we have D(2625), so for 5
# dimensions, we are sampling 13125 latent variables.
DIM = 10


pmf = PMF(train, DIM, ALPHA, std=0.05)

使用MAP进行预测#

# Find MAP for PMF.
pmf.find_map();

非常好。我们首先要做的是确保我们获得的MAP估计是合理的。我们可以通过计算从\(U\)和\(V\)的MAP值得到的预测评分的RMSE来做到这一点。首先，我们定义一个函数，用于从\(U\)和\(V\)生成预测评分\(R\)。我们通过将所有低于1的值设置为1，并将所有高于5的值设置为5，来确保实际评分的界限得到执行。最后，我们计算训练集和测试集的RMSE。我们预计测试集的RMSE会更高。两者之间的差异给出了我们对过拟合程度的某种概念。总是会存在一些差异，但如果在训练集上RMSE非常低，而在测试集上RMSE很高，这无疑是过拟合的明确迹象。

def eval_map(pmf_model, train, test):
    U = pmf_model.map["U"]
    V = pmf_model.map["V"]
    # Make predictions and calculate RMSE on train & test sets.
    predictions = pmf_model.predict(U, V)
    train_rmse = rmse(train, predictions)
    test_rmse = rmse(test, predictions)
    overfit = test_rmse - train_rmse

    # Print report.
    print("PMF MAP training RMSE: %.5f" % train_rmse)
    print("PMF MAP testing RMSE:  %.5f" % test_rmse)
    print("Train/test difference: %.5f" % overfit)

    return test_rmse


# Add eval function to PMF class.
PMF.eval_map = eval_map

# Evaluate PMF MAP estimates.
pmf_map_rmse = pmf.eval_map(train, test)
pmf_improvement = baselines["mom"] - pmf_map_rmse
print("PMF MAP Improvement:   %.5f" % pmf_improvement)

PMF MAP training RMSE: 2.63259
PMF MAP testing RMSE:  2.63524
Train/test difference: 0.00265
PMF MAP Improvement:   -1.62774

我们实际上看到在MAP估计和均值性能之间性能有所下降。我们还在训练集和测试集的RMSE值之间存在相当大的差异。这表明我们从数据中计算出的\(\alpha_U\)和\(\alpha_V\)的点估计在控制模型复杂性方面做得并不好。

让我们看看是否可以通过使用MCMC采样来近似我们的后验分布，从而改进我们的估计。我们将抽取500个样本，其中500个为调整样本。

使用MCMC进行预测#

# Draw MCMC samples.
pmf.draw_samples(draws=500, tune=500)

诊断和后验预测检查#

下一步是检查我们应该丢弃多少样本作为燃烧期。通常，我们会使用轨迹图来大致了解采样变量开始收敛的位置。在这种情况下，我们有多维样本，因此我们需要找到一种方法来近似它们。一种方法由Salakhutdinov 和 Mnih [2008]提出。我们可以计算每一步的\(U\)和\(V\)的Frobenius范数，并监控这些范数以观察其收敛情况。这基本上可以让我们了解潜在变量的平均幅度何时趋于稳定。下面显示了\(U\)和\(V\)的Frobenius范数的方程。我们将使用numpy的linalg包来计算这些范数。

\[ \|U\|_{Fro}^2 = \sqrt{\sum_{i=1}^N \sum_{d=1}^D |U_{id}|^2}, \hspace{40pt} \|V\|_{Fro}^2 = \sqrt{\sum_{j=1}^M \sum_{d=1}^D |V_{jd}|^2} \]

def _norms(pmf_model):
    """Return norms of latent variables at each step in the
    sample trace. These can be used to monitor convergence
    of the sampler.
    """

    norms = dict()
    norms["U"] = xr.apply_ufunc(
        np.linalg.norm,
        pmf_model.trace.posterior["U"],
        input_core_dims=[["users", "latent_factors"]],
        kwargs={"ord": "fro", "axis": (-2, -1)},
    )
    norms["V"] = xr.apply_ufunc(
        np.linalg.norm,
        pmf_model.trace.posterior["V"],
        input_core_dims=[["movies", "latent_factors"]],
        kwargs={"ord": "fro", "axis": (-2, -1)},
    )

    return xr.Dataset(norms)


def _traceplot(pmf_model):
    """Plot Frobenius norms of U and V as a function of sample #."""
    fig, ax = plt.subplots(2, 2, figsize=(12, 7))
    az.plot_trace(pmf_model.norms(), axes=ax)
    ax[0][1].set_title(label=r"$\|U\|_{Fro}^2$ at Each Sample", fontsize=10)
    ax[1][1].set_title(label=r"$\|V\|_{Fro}^2$ at Each Sample", fontsize=10)
    ax[1][1].set_xlabel("Sample Number", fontsize=10)


PMF.norms = _norms
PMF.traceplot = _traceplot

pmf.traceplot()

../_images/5b2551c6f4199dc0aca603ee728f30c661836b3959edf8ead2848973039ac1c4.png

看起来在默认调整后，\(U\) 和 \(V\) 达到了收敛。在测试收敛性时，我们还希望看到我们所寻找的特定统计量的收敛，因为后验的不同特征可能会以不同的速率收敛。让我们也绘制一个RSME的轨迹图。我们将为训练集和测试集计算RMSE，尽管收敛性仅由训练集上的RMSE指示。此外，让我们在训练/测试集上计算一个运行中的RMSE，以查看随着我们继续采样，整体性能是提高还是降低。

请注意，这里我们只从一个链中进行采样，这使得像\(\hat{R}\)这样的收敛统计量变得不可能（我们仍然可以计算split-rhat，但目的不同）。不从多个链中采样的原因是PMF可能没有唯一解。因此，在没有约束的情况下，解决方案充其量是对称的，最坏的情况下在任何旋转下都是相同的，无论如何都容易受到标签切换的影响。事实上，如果我们从多个链中采样，我们会看到较大的\(\hat{R}\)，表明采样器正在探索参数空间不同部分的不同的解决方案。

def _running_rmse(pmf_model, test_data, train_data, plot=True):
    """Calculate RMSE for each step of the trace to monitor convergence."""
    results = {"per-step-train": [], "running-train": [], "per-step-test": [], "running-test": []}
    R = np.zeros(test_data.shape)
    for cnt in pmf.trace.posterior.draw.values:
        U = pmf_model.trace.posterior["U"].sel(chain=0, draw=cnt)
        V = pmf_model.trace.posterior["V"].sel(chain=0, draw=cnt)
        sample_R = pmf_model.predict(U, V)
        R += sample_R
        running_R = R / (cnt + 1)
        results["per-step-train"].append(rmse(train_data, sample_R))
        results["running-train"].append(rmse(train_data, running_R))
        results["per-step-test"].append(rmse(test_data, sample_R))
        results["running-test"].append(rmse(test_data, running_R))

    results = pd.DataFrame(results)

    if plot:
        results.plot(
            kind="line",
            grid=False,
            figsize=(15, 7),
            title="Per-step and Running RMSE From Posterior Predictive",
        )

    # Return the final predictions, and the RMSE calculations
    return running_R, results


PMF.running_rmse = _running_rmse

predicted, results = pmf.running_rmse(test, train)

../_images/b64ff1f05a93e6953efa11de6aede55acb876803396c57c5f73447a2a6f5c6fa.png

# And our final RMSE?
final_test_rmse = results["running-test"].values[-1]
final_train_rmse = results["running-train"].values[-1]
print("Posterior predictive train RMSE: %.5f" % final_train_rmse)
print("Posterior predictive test RMSE:  %.5f" % final_test_rmse)
print("Train/test difference:           %.5f" % (final_test_rmse - final_train_rmse))
print("Improvement from MAP:            %.5f" % (pmf_map_rmse - final_test_rmse))
print("Improvement from Mean of Means:  %.5f" % (baselines["mom"] - final_test_rmse))

Posterior predictive train RMSE: 0.78127
Posterior predictive test RMSE:  0.90711
Train/test difference:           0.12584
Improvement from MAP:            1.72813
Improvement from Mean of Means:  0.10039

我们这里有一些有趣的结果。正如预期的那样，我们的MCMC采样器在训练集上提供了较低的误差。然而，这似乎是以过度拟合数据为代价的。这导致测试RMSE相对于MAP有所下降，尽管它仍然比我们最好的基线要好得多。那么为什么会这样呢？回想一下，我们为我们的精度参数\(\alpha_U\)和\(\alpha_V\)使用了点估计，并且我们选择了一个固定的精度\(\alpha\)。这样做很可能会以某种方式限制我们的后验，使其偏向于训练数据。实际上，用户评分和电影评分的方差不太可能等于我们使用的样本方差的平均值。此外，最合理的观测精度\(\alpha\)也可能不同。

结果摘要#

让我们总结一下我们的结果。

size = 100  # RMSE doesn't really change after 100th sample anyway.
all_results = pd.DataFrame(
    {
        "uniform random": np.repeat(baselines["ur"], size),
        "global means": np.repeat(baselines["gm"], size),
        "mean of means": np.repeat(baselines["mom"], size),
        "PMF MAP": np.repeat(pmf_map_rmse, size),
        "PMF MCMC": results["running-test"][:size],
    }
)
fig, ax = plt.subplots(figsize=(10, 5))
all_results.plot(kind="line", grid=False, ax=ax, title="RMSE for all methods")
ax.set_xlabel("Number of Samples")
ax.set_ylabel("RMSE");

../_images/936040c4347099c41577bbefec0ee8858c56046c30924aa7d4a61ff5c52453b0.png

概述#

我们着手预测用户对未观看电影的偏好。首先，我们讨论了协同过滤的用户-用户和物品-物品邻域方法背后的直观概念。然后，我们将这些直观概念形式化。在深入理解问题背景后，我们继续探索Movielens数据集的一个子集。在发现一些普遍模式后，我们定义了三种基线方法：均匀随机、全局均值和均值的均值。为了超越我们的基线方法，我们使用pymc实现了概率矩阵分解（PMF）的基本版本。

我们的结果表明，在预测任务中，均值的均值方法是我们的最佳基线。正如预期的那样，我们能够通过使用Powell优化的PMF MAP估计获得显著的RMSE减少。我们展示了一种使用采样变量的Frobenius范数来监控具有高维度采样空间的MCMC采样器收敛性的方法。使用这种方法的轨迹图似乎表明我们的采样器收敛到了后验分布。使用这个后验的结果显示，尝试通过MCMC采样来改进MAP估计实际上过度拟合了训练数据并增加了测试RMSE。这很可能是由于通过固定精度参数\(\alpha\)、\(\alpha_U\)和\(\alpha_V\)限制了后验分布。

作为对此分析的后续，实现PMF的逻辑和约束版本也将是有趣的。我们预计这两种模型都将优于基本的PMF模型。我们还可以实现PMF的完全贝叶斯版本（BPMF）[Salakhutdinov和Mnih，2008]，它在模型参数上放置超先验，以自动学习\(U\)和\(V\)的理想均值和精度参数。这可能会解决我们在本次分析中遇到的问题。我们预计BPMF将通过学习更合适的超参数和参数来改进此处生成的MAP估计。对于pymc中的基本（但可用的！）BPMF实现，请参见此gist。

如果你已经做到了这一步，那么恭喜你！你现在对如何构建一个基本的推荐系统有了一些了解。这些相同的想法和方法可以用于许多不同的推荐任务。物品可以是电影、产品、广告、课程，甚至是其他人。任何时候，只要你能够构建一个包含用户偏好的用户-物品矩阵，你就可以使用这些类型的协同过滤算法来预测缺失的值。如果你想了解更多关于推荐系统的知识，第一个参考资料是一个很好的起点。

作者#

本分析中讨论的模型由Ruslan Salakhutdinov和Andriy Mnih开发。代码和支持文本是Mack Sweeney的原创作品，由Colin Carroll和Rob Zinkov进行了修改，以适应MovieLens数据集。

参考资料#

[1] (1,2)

Yehuda Koren, Robert Bell, 和 Chris Volinsky. 推荐系统的矩阵分解技术。计算机，42(8):30–37, 2009. doi:10.1109/MC.2009.263.

[2]

F. Maxwell Harper 和 Joseph A. Konstan。MovieLens 数据集：历史和背景。2016年1月。URL: https://doi.org/10.1145/2827872.

[3] (1,2)

Andriy Mnih 和 Russ R Salakhutdinov. 概率矩阵分解。在 J. Platt, D. Koller, Y. Singer, 和 S. Roweis 编辑的 Advances in Neural Information Processing Systems 中，第 20 卷。Curran Associates, Inc., 2008. URL: https://proceedings.neurips.cc/paper/2007/file/d7322ed717dedf1eb4e6e52a37ea7bcd-Paper.pdf.

[4]

史蒂文·J·诺兰和杰弗里·E·辛顿。通过软权重共享简化神经网络。神经计算，4(4):473–493, 1992。

[5] (1,2,3)

Ruslan Salakhutdinov 和 Andriy Mnih。使用马尔可夫链蒙特卡罗的贝叶斯概率矩阵分解。在 第25届国际机器学习会议论文集，第25卷，880–887页。2008年。

[6]

肯·戈德堡、特蕾莎·罗德和克里斯·珀金斯。Eigentaste：一种常数时间的协同过滤算法。信息检索，4:133–151，2001年。

水印#

%load_ext watermark
%watermark -n -u -v -iv -w -p pytensor,aeppl,xarray

Last updated: Mon Jun 13 2022

Python implementation: CPython
Python version       : 3.9.7
IPython version      : 8.4.0

pytensor: 2.5.1
aeppl : 0.0.27
xarray: 2022.3.0

pandas    : 1.4.2
pymc      : 4.0.0b6
arviz     : 0.12.1
scipy     : 1.7.3
matplotlib: 3.5.2
logging   : 0.5.1.2
pytensor    : 2.5.1
numpy     : 1.22.4
xarray    : 2022.3.0

Watermark: 2.3.1

许可证声明#

本示例库中的所有笔记本均在MIT许可证下提供，该许可证允许修改和重新分发，前提是保留版权和许可证声明。

引用 PyMC 示例#

要引用此笔记本，请使用Zenodo为pymc-examples仓库提供的DOI。

重要

许多笔记本是从其他来源改编的：博客、书籍……在这种情况下，您应该引用原始来源。

同时记得引用代码中使用的相关库。

这是一个BibTeX的引用模板：

@incollection{citekey,
  author    = "<notebook authors, see above>",
  title     = "<notebook title>",
  editor    = "PyMC Team",
  booktitle = "PyMC examples",
  doi       = "10.5281/zenodo.5654871"
}

渲染后可能看起来像：

分类

标签