Fast Greedy MAP for DPP 论文精读

论文引用: Chen, Laming, Guoxin Zhang, and Eric Zhou. “Fast greedy map inference for determinantal point process to improve recommendation diversity.” Advances in Neural Information Processing Systems. 2018.

内容推荐系统的设计宗旨是帮助用户从近乎无限的内容中找到自己喜欢的，为了这个目标，推荐系统最主要的两个任务就是探索与利用。利用的含义是当推荐系统有足够的依据推测用户当下可能喜欢什么的时候，将相关的内容推荐给用户；探索的含义是当推荐系统不知道用户是否喜欢某个品类的内容时，推荐少量这些品类的精选内容给用户，试探用户的反馈。利用的任务我们通常交给排序算法去解决；而探索的任务则困难得多，或者说风险高得多，它需要借助排序之外的手段去实现，例如基于用户实时反馈进行快速试探（冷启动），或者保持一定比例的流量专门做试探，等等。

无论使用哪种方式，这些试探的内容还是需要和排序的结果进行融合后，以一定的顺序展示给用户。由于这两类任务是不可比的，在进行融合时很容易就引入一些人工规则。人工规则设计的好其实也是一项技术活，更简单的办法是用算法来将这些结果进行融合，既能最大化用户的喜好，又能引入尽可能多样的内容。本文主要介绍基于 DPP 的多样性算法。

DPP 简介

DPP 最早是在量子物理学中引入 $^{[1]}$，用于估计费米子系统在热平衡状态下的分布，它能够精确的描述粒子间的斥力。在内容推荐领域，我们可以把斥力与物品间的相似度进行类比：对于一个推荐列表，我们希望列表内的物品不是那么雷同，换句话说，我们希望相似的物品能够相互排斥，这就很符合 DPP 的建模场景，因此 DPP 经常被用于做物品的多样性调整。

DPP 用于推荐系统中来提升推荐列表多样性时，最大的优势在于，它能够直接估计整个列表的得分，而不是两两之间进行比较，因此它对所有物品间的相似关系考虑的更加全面。它的精确定义可以参考 wiki $^{[2]}$。以内容推荐场景来说，给定一个离散的内容集合 $Z=\{1,2,\cdots,m\}$ 及其关联的一个内核矩阵 $\boldsymbol{L}\in \mathbb{R}^{m\times m}$，一个 DPP $\mathcal{P}$ 是对该集合所有子集的一个概率度量，并且每个子集 $Y\subseteq Z$ 的概率与该子集投影在内核矩阵中的主子矩阵的行列式 $^{[3]} $ 成正比：

$\mathcal{P}(Y)\propto \text{det}(\boldsymbol{L}_Y)$

内核矩阵马上会介绍如何构造，这里只需要知道它是半正定的 (PSD, positive semidefinite)，因此 $\det(\boldsymbol{L}_Y) \ge 0$。举个例子，给定 $Z={1,2}$， $\boldsymbol{L}=\left[\begin{matrix}2,1\\1,1\end{matrix}\right]$ ，则 DPP $\mathcal{P}$ 对 $Z$ 所有子集的概率度量如下：

$\begin{matrix} \begin{align} (未归一化)\\ \mathcal{P}(\varnothing)&\propto \text{det}(\varnothing)=1 \\ \mathcal{P}(\{1\})&\propto \text{det}([2])=2\\ \mathcal{P}(\{2\})&\propto \text{det}([1])=1\\ \mathcal{P}(\{1,2\})&\propto \text{det}(\boldsymbol{L})=1 \end{align} \end{matrix} \Rightarrow \begin{matrix} (归一化后)\\ \mathcal{P}(\varnothing)&=0.2\\ \mathcal{P}(\{1\})&=0.4\\ \mathcal{P}(\{2\})&=0.2\\ \mathcal{P}(\{1,2\})&=0.2 \end{matrix}$

在这个例子里，我们想找到最大概率的集合就是 $\{1\}$ ，寻找这个集合的过程就是最大后验推断 (MAP, maximum a posteriori)。不幸的是，DPP 的 MAP 问题是 NP-hard 的；幸运的是，$\log(\det(\cdot))$ 是一个 submodular 函数 $^{[4]}$，而submodular 函数的贪心算法有 1/2 的最优近似保证。由此，文章基于 DPP 模型提出了一种 MAP 的贪心近似推断方法，将时间复杂度降到 $O(n^2m)$，其中，$m$ 表示候选集合内容数量，$n$ 表示待推荐的内容数量。本文只介绍单个短列表内的 DPP 模型推断，基于窗口的长列表的推断可以参考论文或者其他博客。

DPP 建模

前面简单介绍了 DPP 模型的概念和特点，接下来考虑将 DPP 模型应用到以下典型场景：从排序模型输出的列表 $Z$ 中选择 top $n$ 的内容 $Y$ 推荐给用户，既能保证列表里的内容都是用户感兴趣的，同时提升推荐列表的多样性。DPP 建模的核心在于内核矩阵的构建。

要保证用户感兴趣，也就是要保证推荐内容与用户的相关性，一般可以通过排序得分来度量，我们把列表 $Z$ 中内容的排序分记为 $\boldsymbol{r}=[r_1, r_2, \cdots, r_m]$ 。

而内容的多样性一般通过内容间是不是不相似来度量。相似度评估可以是隐式的，例如基于物品的 embedding；也可以是规则式的，比如只考察可读的几个维度：标签、分类、作者等，然后基于 Jaccard 相似度来计算。我们把列表 $Z$ 中内容两两间相似度记为 $\boldsymbol{S}$，其中，$\boldsymbol{S}_{i,j}$ 表示内容 $i,j$ 之间的相似度，相似度的取值需要缩放到区间 $[0,1]$。

那怎么生成内核矩阵呢？前面提到，内核矩阵需要是半正定的，因此可以被分解为 $\boldsymbol{L}=\boldsymbol{B}^\top \boldsymbol{B}$，其中，$\boldsymbol{B}\in\mathbb{R}^{m\times m}$ 的每一列 $\boldsymbol{b}_i$ 可以看成是相应内容的表示向量，这个表示向量既要考虑与用户的相关性，又能用于相似度计算，因此可以设计成 $\boldsymbol{b}_i=r_i \boldsymbol{f}_i$，$\boldsymbol{f}_i$ 表示内容 $i$ 的 embedding 向量，或者是元数据的 onehot 编码等。这样，

$\boldsymbol{L}_{i,j}=\langle\boldsymbol{b}_i,\boldsymbol{b}_j\rangle=\langle r_i\boldsymbol{f}_i,r_j\boldsymbol{f}_j\rangle=r_ir_j\langle\boldsymbol{f}_i,\boldsymbol{f}_j\rangle=r_ir_j\boldsymbol{S}_{i,j}$

写成矩阵的形式就是：

$\boldsymbol{L}=\text{diag}(\boldsymbol{r}) \cdot \boldsymbol{S}\cdot \text{diag}(\boldsymbol{r})$

为了能够贪心求解，我们对目标行列式取对数，得到 submodular 的目标，即找到最优的列表 $Y$ 来最大化：

$\begin{align} \log\mathcal{P}(Y) \propto \log\det(\boldsymbol{L}_{Y})&=\log\Big(\det\big(\text{diag}(\boldsymbol{r}_Y) \cdot \boldsymbol{S}_Y\cdot \text{diag}(\boldsymbol{r}_Y)\big)\Big) \\ &=\log\Big(\det\big(\text{diag}(\boldsymbol{r}_Y)\big) \cdot\det\big( \boldsymbol{S}_Y\big) \cdot\det\big( \text{diag}(\boldsymbol{r}_Y)\big)\Big) \\ &=2\sum_{i\in Y}\log(r_i) + \log \det\big( \boldsymbol{S}_Y\big) \end{align}$

这里矩阵或者向量的下标为 $Y$ 表示得到 $Y$ 中所有元素在矩阵/向量的索引得到的主子矩阵/向量。这样，我们就能看出来，最终推荐列表 $Y$ 的 DPP 概率得分实际上是 $\sum_{i\in R}\log(r_i)$ 与 $\log \det\big( \boldsymbol{S}_Y\big)$ 的线性组合。前者表示内容 $i$ 的排序得分，后者表示 $Y$ 中内容的差异得分。

为了能够调整排序分与相似分的权重，可以将上式进一步调整成下式：

$\begin{align} \log\mathcal{P}(Y) &\propto \theta \cdot \sum_{i\in Y}r_i + (1-\theta)\cdot\log \det\big( \boldsymbol{S}_Y\big) \\ &\propto\frac{\theta}{2(1-\theta)} \cdot \sum_{i\in Y}r_i + \log \det\big( S_R\big)+ \frac{\theta}{2(1-\theta)} \cdot \sum_{i\in Y}r_i\\ &\overset{\alpha=\frac{\theta}{2(1-\theta)}}{=} \sum_{i\in Y}\log\big(e^{\alpha r_i}) + \log \det\big( \boldsymbol{S}_Y\big)+ \sum_{i\in Y}\log\big(e^{\alpha r_i})\\ &=\log\Big(\det\big(\text{diag}(e^{\alpha\boldsymbol{r}_Y}) \cdot \boldsymbol{S}_Y\cdot \text{diag}(e^{\alpha\boldsymbol{r}_Y})\big)\Big) \end{align}$

再回到内核矩阵上来，我们得到了一个新的 $\boldsymbol{L}_{Y} = \text{diag}(e^{\alpha\boldsymbol{r}_Y}) \cdot \boldsymbol{S}_Y\cdot \text{diag}(e^{\alpha\boldsymbol{r}_Y})$，其中，$\alpha=\frac{\theta}{2(1-\theta)}$，$\theta\in [0,1)$，$\theta$ 越大，排序得分的影响越大，$\theta$ 越小，列表中的内容差异越大。

需要说明的是，这种构造方法只有在 $\boldsymbol{S}$ 是半正定的情况下，才能保证 $\boldsymbol{L}$ 是半正定的。

Fast Greedy MAP Inference

假设我们已经构造出内容推荐场景下的内核矩阵 $\boldsymbol{L}$，下面就需要贪心每次向最终推荐列表 $Y$ 中添加内容 $j$，使得我们的 submodular 的目标增益最大：

$j = \arg \max_{i\in Z-Y} \log \det(\boldsymbol{L}_{Y\cup\{i\}}) - \log \det(\boldsymbol{L}_{Y})$

首先，我们推导一下，要最大化目标增益，其实我们每次只需要找到一个内容 $i$，它在 $\boldsymbol{L} $ 主子矩阵 $\boldsymbol{L}_{Y\cup\{i\}}$ 的 Cholesky 分解矩阵的对角元素最大即可。

具体来讲，如果 $\boldsymbol{L}$ 是半正定的，则 $\boldsymbol{L}$ 的任意主子矩阵 $\boldsymbol{L}_Y$ 也是半正定的，因此可以进行 Cholesky 分解。但是由于半正定矩阵的 Cholesky 分解不唯一，我们姑且假设 $\boldsymbol{L}_Y$ 是正定的（后面会讨论异常如何处理），则它可以被分解为：

$\boldsymbol{L}_Y=\boldsymbol{V}\boldsymbol{V}^\top$

这里 $\boldsymbol{V}$ 是一个下三角矩阵，并且对角线上的元素都是正实数；$Y$ 是我们目前根据贪心规则已经从候选列表 $Z$ 选出的部分的推荐结果。对于内容 $i\in Z-Y$， $Y'=Y\cup \{i\}$ 在 $\boldsymbol{L} $ 上投影的主子矩阵 $\boldsymbol{L}_{Y’}$ 可以表示为：

$\boldsymbol{L}_{Y'} = \left [\begin{matrix} \boldsymbol{L}_{Y} & \boldsymbol{L}_{Y,i}\\ \boldsymbol{L}_{i,Y} & \boldsymbol{L}_{i,i} \end{matrix}\right ]= \left [\begin{matrix} \boldsymbol{V} & \boldsymbol{0}\\ \boldsymbol{c}_i & d_i \end{matrix}\right ] \left [\begin{matrix} \boldsymbol{V} & \boldsymbol{0}\\ \boldsymbol{c}_i & d_i \end{matrix}\right ]^\top= \left [\begin{matrix} \boldsymbol{V}\boldsymbol{V}^\top & \boldsymbol{V}\boldsymbol{c}_i^\top\\ \boldsymbol{c}_i\boldsymbol{V}^\top+d_i\boldsymbol{c}_i^\top & \boldsymbol{c}_i\boldsymbol{c}_i^\top+d_i^2 \end{matrix}\right ] \qquad(1)$

根据式 $(1)$ 可以得 $\det(\boldsymbol{L}_{Y'})$ 与 $\det(\boldsymbol{L}_{Y})$ 的关系：

$\begin{align}\det(\boldsymbol{L}_{Y'}) = \det\left(\left [\begin{matrix}\boldsymbol{V} & \boldsymbol{0}\\\boldsymbol{c}_i & d_i\end{matrix}\right ]\left [\begin{matrix}\boldsymbol{V} & \boldsymbol{0}\\\boldsymbol{c}_i & d_i\end{matrix}\right ]^\top \right)&=\det\left(\left [\begin{matrix}\boldsymbol{V} & \boldsymbol{0}\\\boldsymbol{c}_i & d_i\end{matrix}\right ]\right)\det\left(\left [\begin{matrix}\boldsymbol{V} & \boldsymbol{0}\\\boldsymbol{c}_i & d_i\end{matrix}\right ]^\top\right) \\&=\left(\det(\boldsymbol{V})\cdot d_i\right)^2 \\&=\det(\boldsymbol{V}\boldsymbol{V}^\top)\cdot d_i^2\\&=\det(\boldsymbol{L}_Y)\cdot d_i^2\end{align}$

因此，基于当前局部最优结果 $Y$，下一个被选中的内容 $j$ 满足：

$j = \arg \max_{i\in Z-Y} \log \det(\boldsymbol{L}_{Y\cup\{i\}}) - \log \det(\boldsymbol{L}_{Y}) = \arg\max_{i\in Z-Y} \log(d_i^2)$

也就是说，在给定当前推荐列表 $Y$ 的情况下，我们只需要计算所有剩下候选内容 $i\in Z-Y$ 对应的 $d_i$ 值即可，可以认为 $\log(d_i^2)$ 就是我们每次选中内容 $i$ 得到的边际增益，为了保证这个增益始终为正，需要保证 $d_i^2>1$ (当然，这个增益物理含义不明确，实际上小于 $0$ 问题也不大)。

观察式 $(1)$ 可得：

$d_i^2=\boldsymbol{L}_{i,i}-\parallel\boldsymbol{c}_i\parallel^2_2$

正常来讲，我们需要从矩阵 $\boldsymbol{V}$ 的左上角开始，按行或者按列依次算出所有的 $\boldsymbol{V}_{k,l}$，最后才能得出 $d_i$ $^{[6]}$，但是这样的复杂度太高了。文章提出了一种迭代的计算方法，在每次得到最优推荐结果的同时，记录所有的 $\boldsymbol{c}_i$，并在计算下一轮最优结果时更新这些 $\boldsymbol{c}_i$。

怎么更新 $\boldsymbol{c}_i$ 呢？假设我们在某一轮已经找到最优的内容 $j$，在找的过程中，已经计算并保存所有的 $\boldsymbol{c}_k, \ k\in Z$；当我们再找下一个最优内容 $i$ 时，由于 $Y$ 中的内容数量会增加一个，$\boldsymbol{c}_i$ 的维度也会相应的增加一维，记待更新的向量为 $\boldsymbol{c}_i’=[\boldsymbol{c}_i, e_i]$。还是观察式 $(1)$ 的右上角项，易知在本轮：

$\begin{align} \left [\begin{matrix} \boldsymbol{V} & \boldsymbol{0}\\ \boldsymbol{c}_j & d_j \end{matrix}\right ]\boldsymbol{c}_i'^\top = \boldsymbol{L}_{Y',i} = \left [\begin{matrix} \boldsymbol{L}_{Y,i}\\ \boldsymbol{L}_{j,i} \end{matrix}\right ] &\Rightarrow \langle[\boldsymbol{c}_j d_j], \boldsymbol{c}_i'\rangle=\boldsymbol{L}_{j,i} \\ &\Rightarrow \langle\boldsymbol{c}_j , \boldsymbol{c}_i\rangle+ d_j\cdot e_i=\boldsymbol{L}_{j,i}\\ &\Rightarrow e_i = \frac{\boldsymbol{L}_{j,i}-\langle\boldsymbol{c}_j , \boldsymbol{c}_i\rangle}{d_j} \end{align}$

由此可以进一步算出，

$d_i'^2=\boldsymbol{L}_{i,i}-\parallel\boldsymbol{c}_i'\parallel^2_2=\boldsymbol{L}_{i,i}-\parallel\boldsymbol{c}_i\parallel^2_2 -e_i^2=d_i^2-e_i^2$

初始时，$Y=\varnothing$， $d_i^2=\boldsymbol{L}_{i,i}$ ， $j=\arg \max_{i\in Z} d_i^2$ 。

算法整体流程

上面我们推导了内核矩阵如何构建，基于内核矩阵如何贪心的找到增益最高的推荐列表。最后总结一下整个算法流程（python 实现可参考 $[6]$）：

Fast Greedy MAP Inference for DPP

Input: candidates $Z=[0,1,\cdots m-1]$, similarity matrix of candidates $\boldsymbol{S}\in \mathbb{R}^{m\times m}$, ranking vector of candidates $\boldsymbol{r}=[r_0,r_1,\cdots,r_{m-1}]$ , recommend item count $n$
Initialize: $\boldsymbol{L} = \text{diag}(e^{\alpha\boldsymbol{r}}) \cdot \boldsymbol{S}\cdot \text{diag}(e^{\alpha\boldsymbol{r}})$, $\boldsymbol{C}=\boldsymbol{0}_{m\times n}$, $\boldsymbol{d}^2=\text{diag}(\boldsymbol{L})$, $j=\arg \max \boldsymbol{d}^2$, $Y=[j]$
for $k$ in range($0$, $n-1$):
$\boldsymbol{e}=(\boldsymbol{L}_{j,:}-\boldsymbol{C}_{:,:k}\cdot \boldsymbol{C}_{j,:k})/ \sqrt{\boldsymbol{d}^2_j}$
$\boldsymbol{C}_{:,k} = \boldsymbol{e}$
$\boldsymbol{d}^2_j = -\infty$
$j=\arg \max \boldsymbol{d}^2$
if $\boldsymbol{d}^2_j<\epsilon$:
break
$Y=[Y, j]$
return $Y$

前面提到，我们假设 $\boldsymbol{S}$ 和 $\boldsymbol{L} $ 是半正定的，但是在推导迭代公式的时候却假设 $\boldsymbol{L}_Y $ 是正定的，因此在实际计算的时候可能得到 $d_j=0$，导致第 4 步的数值不稳定，因此在第 8 步引入一个很小的 $\epsilon>0$，当 $d_j^2<\epsilon$ 时，就进入异常处理程序。这里写的是直接返回结果，实际也可以将列表用某种策略填充到 $n$ 个再返回。

什么情况下会发生这个问题呢？回顾前文，如果有两个内容 $i,j$ 的表征向量线性相关/元数据完全相同，这个时候 $\boldsymbol{b}_i,\boldsymbol{b}_j$ 是线性相关的，如果它们都被选到推荐列表 $Y$ 中时，就会导致 $\det(\boldsymbol{L}_Y )=0$，说明此时 $\boldsymbol{L}_Y $ 半正定，其 Cholesky 分解的矩阵 $\boldsymbol{V}$ 对角线元素为 $0$，即 $d_j=0$。

这种情况在现实中是可能出现的，可以考虑在计算相似矩阵 $\boldsymbol{S}$ 的时候，为每两个内容的相似度加上一个小的随机数（要对称的加），保证 $\boldsymbol{S}$ 是正定的即可，但是这种方式可能对推荐结果有一定的影响。当然，如果 $Z$ 足够多，而且里面内容足够多样，这种情况出现的概率就非常低了。