這幾天我的主要任務是對論文《Towards Long-term Fairness in Recommendation》[1] 中所描述的算法進行編程實現，然后測試該算法的效果並記錄。以下分模型算法細節實現、數據集、模型評估准則、測試結果記錄四個部分來描述我的工作。
模型算法細節實現
由論文描述可知，論文算法最核心的部分就是以下帶約束優化問題的求解：
該問題被論文作者稱為約束策略優化問題。論文中提到該問題具體的求解算法可參考另一篇論文《Constrained policy optimization》[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條件理論)。論文利用強對偶性來求得該問題的解析解。以下是具體推導。我們先定義拉格朗日函數：
根據拉格朗日對偶性，原始問題的對偶問題是極大極小問題：
先求\(\underset{x}{\min}L(x, \lambda, \nu)\)，我們令\(\nabla_{x}L(x, \lambda, \nu)=0\)，解得
將\(x^*\)帶入拉格朗日函數，則問題進一步轉化為
我們將該式進行變量替換，令\(q=g^T H^{-1} g, r = g^TH^{-1}b, s=b^TH^{-1} b\),可進一步化簡為
我們又由\(\frac{\partial L(\lambda, \nu)}{\partial \nu} = - \frac{1}{2\lambda}(2r+2\nu s) + c\)，我們由\(\mathbb{R}_{+}\)上的單變量凸二次函數優化理論可以得到
進一步，我們將問題轉化為
注意，這里\(\Lambda_{a} = \{\lambda | \lambda c - r > 0, \lambda \geqslant 0\}\)，\(\Lambda_{b} = \{\lambda | \lambda c - r \leqslant 0, \lambda \geqslant 0\}\)。
於是，綜上所述，當至少有一個可行解時，最優解\(x^{*}\)滿足
這里\(\lambda^*\)和\(\nu^*\)定義如下
這里 \(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):
這里\(y_{t} = r_{t} + \gamma_{r}V_{w}(s_{t+1})\)。此外，為了提高准確率，我們針對約束策略優化引入了單獨的代價函數
評委\(V_{\phi}(s)\)，並按照與式\((12)\)相似的方式更新：
這里\(y_{t} = c_{t} + \gamma_{c}V_{\phi}(s_{t+1})\)。這樣，我們訓練部分的核心算法就可以描述如下面算法所示。
注意，其中的經驗回放(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)\)來計算
這里\(\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.