推荐系统公平性论文阅读(五)

这几天我的主要任务是对论文《Towards Long-term Fairness in Recommendation》[1] 中所描述的算法进行编程实现,然后测试该算法的效果并记录。以下分模型算法细节实现、数据集、模型评估准则、测试结果记录四个部分来描述我的工作。

模型算法细节实现

由论文描述可知,论文算法最核心的部分就是以下带约束优化问题的求解:

\[\begin{matrix}
\theta_{k+1} = \underset{\theta}{\arg\max}g^{T}(\theta-\theta_{k}) \\
s.t. \quad c+b^{T}(\theta – \theta_{k}) \leqslant 0 \\
\frac{1}{2}(\theta – \theta_{k})^{T}H(\theta-\theta_{k}) \leqslant \delta
\end{matrix}
\tag{1}
\]

该问题被论文作者称为约束策略优化问题。论文中提到该问题具体的求解算法可参考另一篇论文《Constrained policy optimization》[2]。我查阅了另一篇论文,该论文中详细地论证了形如以下问题的求解方法:

\[\begin{matrix}
p^{*} = \underset{x}{\min} g^{T}x \\
s.t. \quad b^{T}x + c \leqslant 0 \\
x^{T}Hx \leqslant \delta
\end{matrix}
\tag{2}
\]

这里\(g, b, x \in \mathbb{R}^{n}, c, \delta \in \mathbb{R}, \delta > 0, H \in \mathbb{S}^{n}, \text{且} H \succ 0\)
该问题是线性目标函数+线性与二次约束的最优化问题,是一个典型的凸优化问题。当至少有一个严格的可行解时,强对偶性得到满足(根据Slater’s条件理论)。论文利用强对偶性来求得该问题的解析解。以下是具体推导。我们先定义拉格朗日函数:

\[L(x, \lambda, \nu) = g^{T}x + \frac{\lambda}{2}(x^THx-\delta) + \nu(b^Tx+c)
\tag{3}
\]

根据拉格朗日对偶性,原始问题的对偶问题是极大极小问题:

\[\underset{\nu \geqslant 0}{\underset{\lambda \geqslant 0}{\max}} \underset{x}{\min}L(x, \lambda, \nu)
\tag{4}
\]

先求\(\underset{x}{\min}L(x, \lambda, \nu)\),我们令\(\nabla_{x}L(x, \lambda, \nu)=0\),解得

\[ x^* = – \frac{1}{\lambda}H^{-1}(g+\nu b)
\tag{5}
\]

\(x^*\)带入拉格朗日函数,则问题进一步转化为

\[ \underset{\nu \geqslant 0}{\underset{\lambda \geqslant 0}{\max}}
– \frac{1}{2\lambda}(g+\nu b)^TH^{-1}(g+\nu b) + (\nu c – \frac{1}{2}\lambda \delta)
\tag{6}
\]

我们将该式进行变量替换,令\(q=g^T H^{-1} g, r = g^TH^{-1}b, s=b^TH^{-1} b\),可进一步化简为

\[ \underset{\nu \geqslant 0}{\underset{\lambda \geqslant 0}{\max}}
– \frac{1}{2\lambda}(q + 2\nu r + \nu^2 s) + (\nu c – \frac{1}{2} \lambda \delta)
\tag{7}
\]

我们又由\(\frac{\partial L(\lambda, \nu)}{\partial \nu} = – \frac{1}{2\lambda}(2r+2\nu s) + c\),我们由\(\mathbb{R}_{+}\)上的单变量凸二次函数优化理论可以得到

\[ \nu = (\frac{\lambda c- r}{s})_{+}
\tag{8}
\]

进一步,我们将问题转化为

\[\underset{\lambda \geqslant 0}{\max}
\left\{
\begin{matrix}
\frac{1}{2 \lambda}(\frac{r^2}{s}-q) + \frac{\lambda}{2}(\frac{c^2}{s}-\delta) – \frac{rc}{s} \quad \text{如果} \lambda \in \Lambda_{a} \\
-\frac{1}{2}(\frac{q}{\lambda} + \lambda \delta) \quad \text{如果} \lambda \in \Lambda_{b}
\end{matrix}
\right.
\tag{9}
\]

注意,这里\(\Lambda_{a} = \{\lambda | \lambda c – r > 0, \lambda \geqslant 0\}\)\(\Lambda_{b} = \{\lambda | \lambda c – r \leqslant 0, \lambda \geqslant 0\}\)

于是,综上所述,当至少有一个可行解时,最优解\(x^{*}\)满足

\[ x^* = – \frac{1}{\lambda^*}H^{-1}(g+v^*b)
\tag{10}
\]

这里\(\lambda^*\)\(\nu^*\)定义如下

\[ \begin{matrix}
\nu^* = ( \frac{\lambda^*c – r}{s})_{+} \\
\lambda^* = \underset{\lambda \geqslant 0}{\arg\max}
\left\{
\begin{matrix}
f_{a}(\lambda) = \frac{1}{2\lambda}(\frac{r^2}{s}-q) + \frac{\lambda}{2}(\frac{c^2}{s} – \delta) – \frac{rc}{s} \text{如果} \lambda c – r >0 \\
f_{b}(\lambda) = -\frac{1}{2}(\frac{q}{\lambda}+\lambda \delta) \text{其他情况}
\end{matrix}
\right.
\end{matrix}
\tag{11}
\]

这里 \(q=g^TH^{-1}g, r=g^TH^{-1}b\)\(s=b^TH^{-1}b\)
接下来是编程细节部分。首先我们需要先计算\(q\)\(r\), \(s\)的值。而这需要高效地计算\(H^{-1}g, H^{-1}b\),而对Hessian矩阵求逆的时间复杂度是很高的,我们通过求解方程组\(Hx_{1}=g,Hx_{2}=b\)来实现。但Hessian矩阵很可能是稀疏的,这进一步加大了
我们的求解难度,因此我们采用共轭梯度法来完成方程组的求解,具体实现代码如下:

def cg_solver(Avp_fun, b, device, max_iter=10):
    device = device
    x = torch.zeros_like(b).to(device)
    r = b.clone()
    p = b.clone()
    for i in range(max_iter):
        Avp = Avp_fun(p, retain_graph=True)
        alpha = torch.matmul(r, r) / torch.matmul(p, Avp)
        x += alpha * p
        if i == max_iter - 1:
            return x
        r_new = r - alpha * Avp
        beta = torch.matmul(r_new, r_new) / torch.matmul(r, r
        r = r_new
        p = r + beta * p

注意,这里我们按照共轭梯度法求解方程\(Ax=b\)时,中途有一个步骤需要计算\(A\)和搜索方向向量\(p\)的乘积\(Ap\),这里我们采用专用的函数实现,我们将该函数定义为一个函数闭包(这里用到了Python的语言特性),在这里闭包内层定义的变量对外界完全隔离,具体实现如下:

def get_Hvp_fun(functional_output, inputs, damping_coef=0.0):
    inputs = list(inputs)
    grad_f = flat_grad(functional_output, inputs, create_graph=True)
    def Hvp_fun(v, retain_graph=True):
        gvp = torch.matmul(grad_f, v)
        Hvp = flat_grad(gvp, inputs, retain_graph=retain_graph)
        Hvp += damping_coef * v
        return Hvp
    return Hvp_fun 

这样,我们高效地计算出了\(H^{-1}g\)\(H^{-1}b\),从而得到\(q,r,s\)的值。又由原始论文的定义,可知\(c=J_{C}(\pi_{k})-d\)。这样,我们算法所需要的基础变量q, r, s, c已经得到,若问题是可解的,即如果\(c^2/s-\delta >0\)\(c>0\)时,我们可以使用上面推导的解析公式\((11)\)来计算对偶变量\(\lambda^*\)\(\nu^*\),然后依照解析公式\((10)\)计算出\(x^*\)。其中计算对偶变量\(\lambda^*\)\(\nu^*\)的函数实现如下:

def calc_dual_vars(self, q, r, s, c):
    if c < 0.0 and c ** 2 / s - 2 * self.max_kl > 0.0:
        lam = torch.sqrt(q / (2 * self.max_kl))
        nu = 0.0
        return lam, nu
    A = q - r ** 2 / s
    B = 2 * self.max_kl - c ** 2 / s
    lam_mid = r / c
    # lam_a*
    lam_a = torch.sqrt(A / B)
    # lam_b*
    lam_b = torch.sqrt(q / (2 * self.max_kl))
    f_mid = -0.5 * (q / lam_mid + 2 * lam_mid * self.max_kl)
    f_a = -torch.sqrt(A * B) - r * c / s
    f_b = -torch.sqrt(2 * q * self.max_kl)
    if lam_mid > 0:
        if c < 0:
            if lam_a > lam_mid:
                lam_a = lam_mid
                f_a = f_mid
            if lam_b < lam_mid:
                lam_b = lam_mid
                f_b = f_mid
        else:
            if lam_a < lam_mid:
                lam_a = lam_mid
                f_a = f_mid
            if lam_b > lam_mid:
                lam_b = lam_mid
                f_b = f_mid
    else:
        if c < 0:
            lam = lam_b
        else:
            lam = lam_a
    # lam*
    lam = lam_a if f_a >= f_b else lam_b
    # v*
    nu = max(0.0, (lam * c - r) / s)
    return lam, nu

至此,我们已经完成了论文算法的核心部分,即策略函数\(\pi\)的参数\(\theta\)如何做到服从约束地更新(尽管我们还没有完整实现其训练算法)。
接下来我们描述更完整的框架,即强化学习的模型的策略函数\(\pi\)究竟怎样训练。论文中完整的训练算法采用演员-评委模式,其示意图如下图所示:
演员-评委训练方法示意图
其中,演员中的策略函数\(\pi_{\theta}\)、状态价值函数\(V_{\omega}^{\pi}(s_{t})\)
代价函数\(V_{\phi}^{\pi}(s_{t})\)都由神经网络实现。其中策略函数、状态价值函数和代价函数的参数分别为\(\theta\)\(\omega\)\(\phi\)
我们构建状态价值函数评委\(V_{\omega}(s_t)\)来近似真实的状态价值函数\(V_{w}^{\pi}(s_{t})\),然后被用于优化演员。评委网络按照时间差分学习来优化,其中需要最小化均方误差(MSE):

\[ L(\omega) = \sum_{t}(y_{t} – V_{\omega}(s_{t}))^2
\tag{12}
\]

这里\(y_{t} = r_{t} + \gamma_{r}V_{w}(s_{t+1})\)。此外,为了提高准确率,我们针对约束策略优化引入了单独的代价函数
评委\(V_{\phi}(s)\),并按照与式\((12)\)相似的方式更新:

\[ L(\phi) = \sum_{t}(y_{t}-V_{\phi}(s_{t}))^2
\tag{13}
\]

这里\(y_{t} = c_{t} + \gamma_{c}V_{\phi}(s_{t+1})\)。这样,我们训练部分的核心算法就可以描述如下面算法所示。

fcpo算法伪代码

注意,其中的经验回放(Replay)数据集存储了总共\(T\)个时间步的轨迹(Trajectory),我们从中采样得到状态转换五元组\((s_{t}, a_{t}, r_{t}, c_{t}, s_{t+1})\)。而经验回放数据集我们先对数据进行一遍处理才能得到。
这样,我们的主体算法流程已描述完毕,接下来需要对模型进行测试。我们采用还是和在上一次复现论文中所做的一样,采用Movielens-1M数据集中的用户交易数据(包括用户id,物品id, 用户评分,时间戳)来测试论文提出的FCPO模型。

数据集

我们先按照时间戳的顺序对用户交易数据排序,然后将数据按照4:1的比例划分为训练集和测试集,其中我们将训练集中每个用户被推送的最后一个物品另外划分为验证集。Movielens-1M数据集的一些基础统计量如下表所示。我们基于人气(比如根据物品的曝光次数)将数据划分为\(G_{0}\)\(G_{1}\)两个群组。具体地,曝光次数位居前20%的物品属于受欢迎的群组\(G_{0}\),后80%的物品属于长尾群组\(G_{1}\)

用户数 物品数 动作数/用户数 动作数/物品数 动作数 密度
6040 3706 166 270 1000209 4.468

此外,对于该基于强化学习的推荐系统,每个用户在训练时的初始状态是在训练集中最先点击的5个物品,在测试时的初始状态是在训练集中
最后点击的5个物品。出于简便,我们这里测试时让RL智能体每次向用户推荐一个物品,不过实际上该算法的推荐列表长度可调整。

模型评估准则

我们采用一些常见的准则,比如召回率、F1-Score等来评估我们的推荐系统的表现。除了基于模型精准程度
的准则,我们同时也增设了两个公平性的度量。对于物品曝光给特定用户,我们度量基尼系数;对于物品曝光给特定群组,我们度量流行比率。
基尼系数度量频率分布(比如物品的曝光次数)之间的不公平程度,这可以被视为个体层面的公平性度量。给定来自各物品曝光记录,
\(\mathcal{M}=[g_{1}, g_{2}, …, g_{|\mathcal{I}|}]\),基尼系数可以按照式\((14)\)来计算

\[ Gini\quad Index(\mathcal{G}) = \frac{1}{2 |\mathcal{I}|^2 \bar{g}}\sum_{i=1}^{|\mathcal{I}|}\sum_{j=1}^{|\mathcal{I}|}|g_{i}-g_{j}|
\tag{14}
\]

这里\(\bar{g}\)代表所有物品曝光次数的平均值。流行比率指的是在推荐列表中流行物品的比例,可以被看做是人气层面的公平性度量。
这两个公平性的度量相比原始推荐系统都更加公平。

测试结果

接下来是我们在Movielens-1M上的测试结果,测试结果如下表所示。这里FCPO-1,FCPO-2和FCPO-3对应的公平性约束的参数分别为\(\alpha^{‘}=1, \alpha^{‘}=0.8\)\(\alpha^{‘}=0.4\)
且我们假定推荐列表的长度K=5。

方法 召回率(%)\(\uparrow\) F1(%) \(\uparrow\) 基尼系数(%) \(\downarrow\) 关注度排名(%) \(\uparrow\)
FCPO-1 2.033 2.668 99.81 99.28
FCPO-2 1.520 2.015 99.47 99.28
FCPO-3 0.998 1.449 99.47 99.28

参考文献

  • [1] Ge Y, Liu S, Gao R, et al. Towards Long-term Fairness in Recommendation[C]//Proceedings of the 14th ACM International Conference on Web Search and Data Mining. 2021: 445-453.
  • [2] Achiam J, Held D, Tamar A, et al. Constrained policy optimization[C]//International Conference on Machine Learning. PMLR, 2017: 22-31.