Factorization Machine 詳解（三）—— FM與高階特征訓練

之前的文章回顧：

• FM基礎概念的介紹與訓練方法
• FM開源庫LibFM簡介與使用

可以看到，FM主要用于推測推薦系統中不同特征之間的關聯，尤其适用于稀疏資料的推薦。然而，Rendle等提出的基礎FM模型一般隻能推測2-way的特征關系，而對于更高階的特征關系則缺乏更為有效的訓練手段。另外，在傳統的FM模型中，目标函數是非凸的，是以很難進行優化，尤其是當涉及高階特征關系的挖掘時，目标函數的優化更為複雜。

對于高階特征的訓練，目前有以下兩個思路：

• 使用更高階的Kernel函數對目标函數進行模組化
• 使用深度神經網絡對輸入資料進行處理

目前，在第一種解決方案中，我們需要對目标函數進行模組化：

y ^ K ( x ; λ , P ) : = ∑ s = 1 k λ s K ( p s , x ) \hat{y}_{\mathcal{K}}(\boldsymbol{x} ; \boldsymbol{\lambda}, \boldsymbol{P}) :=\sum_{s=1}^{k} \lambda_{s} \mathcal{K}\left(\boldsymbol{p}_{s}, \boldsymbol{x}\right) y^​K​(x;λ,P):=s=1∑k​λs​K(ps​,x)

總的loss function：

D K ( λ , P ) : = ∑ i = 1 n ℓ ( y i , y ^ K ( x i ; λ , P ) ) D_{\mathcal{K}}(\boldsymbol{\lambda}, \boldsymbol{P}) :=\sum_{i=1}^{n} \ell\left(y_{i}, \hat{y}_{\mathcal{K}}\left(\boldsymbol{x}_{i} ; \boldsymbol{\lambda}, \boldsymbol{P}\right)\right) DK​(λ,P):=i=1∑n​ℓ(yi​,y^​K​(xi​;λ,P))

其中函數 ℓ \ell ℓ為凸的loss function。從中我們可以看到，其使用的目标函數是一個Kernel，即為從 P \boldsymbol P P映射到輸入向量 x \boldsymbol x x空間的Kernel。在正式闡述高階FM的訓練方式之前，先說明一下上述參數的具體意義：

• y ^ K \hat y_{\mathcal K} y^​K​表示的是基于Kernel K \mathcal K K的目标函數。
• λ = [ λ 1 , … , λ k ] T ∈ R k \boldsymbol{\lambda}=\left[\lambda_{1}, \dots, \lambda_{k}\right]^{\mathrm{T}} \in \mathbb{R}^{k} λ=[λ1​,…,λk​]T∈Rk 是超參，具體意義是決定不同次元Kernel的權重。在FM中 λ s = 1 ( s = 1 , 2 , . . . , k ) \lambda_s = 1(s = 1, 2, ..., k) λs​=1(s=1,2,...,k)。
• P = [ p 1 , … , p k ] ∈ R d × k \boldsymbol{P}=\left[\boldsymbol{p}_{1}, \ldots, \boldsymbol{p}_{k}\right] \in \mathbb{R}^{d \times k} P=[p1​,…,pk​]∈Rd×k表示的是已經分解的向量（矩陣）
• K \mathcal K K表示的是Kernel函數。目前研究的Kernel包括ANOVA Kernel以及同質的多項式Kernel

ANOVA Kernel

2016年，Mathieu Blondel等[1]解決了Polynomial Network(PN) 以及Factorization Machine(FM)之間的關系，并且證明了提出了coordinate descent(CD)算法，以訓練三階的FM。關于ANOVA Kernel，定義如下：

A m ( p , x ) : = ∑ j m > ⋯ > j 1 p j 1 x j 1 … p j m x j m \mathcal{A}^{m}(\boldsymbol{p}, \boldsymbol{x}) :=\sum_{j_{m}>\cdots>j_{1}} p_{j_{1}} x_{j_{1}} \ldots p_{j_{m}} x_{j_{m}} Am(p,x):=jm​>⋯>j1​∑​pj1​​xj1​​…pjm​​xjm​​

其中 A 0 ( p , x ) : = 1 \mathcal{A}^{0}(\boldsymbol{p}, \boldsymbol{x}) :=1 A0(p,x):=1以及 A 1 ( p , x ) : = ⟨ p , x ⟩ \mathcal{A}^{1}(\boldsymbol{p}, \boldsymbol{x}) :=\langle\boldsymbol{p}, \boldsymbol{x}\rangle A1(p,x):=⟨p,x⟩。可以看到ANOVA Kernel使用的是多種distinct 的特征（也就是沒有替代的特征組合）。ANOVA Kernel具有多線性的性質。通過展開ANOVA Kernel的定義式，我們可以使用數學歸納法證明，對于 p , x ∈ R d \boldsymbol{p}, \boldsymbol{x} \in \mathbb{R}^{d} p,x∈Rd，以及 j ∈ [ d ] j \in[d] j∈[d]和 1 ≤ m ≤ d 1 \leq m \leq d 1≤m≤d，有如下的遞歸多線性（multi-linearity）關系：

A m ( p , x ) = A m ( p ¬ j , x ¬ j ) + p j x j A m − 1 ( p ¬ j , x ¬ j ) \mathcal{A}^{m}(\boldsymbol{p}, \boldsymbol{x})=\mathcal{A}^{m}\left(\boldsymbol{p}_{\neg j}, \boldsymbol{x}_{\neg j}\right)+p_{j} x_{j} \mathcal{A}^{m-1}\left(\boldsymbol{p}_{\neg j}, \boldsymbol{x}_{\neg j}\right) Am(p,x)=Am(p¬j​,x¬j​)+pj​xj​Am−1(p¬j​,x¬j​)

其中 ¬ j \neg j ¬j表示的是向量中排除掉第 j j j個元素的内容，是以當我們對其中第 j j j個坐标 p j p_j pj​計算導數時，上式是一個線性函數。是以，我們可以了解多線性的含義為對某個特定的方向 j j j，都滿足既凸又凹的線性關系。那麼，根據Rendle的FM定義，我們可以改寫FM的目标函數如下：

y ^ F M ( x ; w , P ) = ⟨ w , x ⟩ + y ^ A 2 ( x ; 1 , P ) \hat{y}_{F M}(\boldsymbol{x} ; \boldsymbol{w}, \boldsymbol{P})=\langle\boldsymbol{w}, \boldsymbol{x}\rangle+\hat{y}_{\mathcal{A}^{2}}(\boldsymbol{x} ; \mathbf{1}, \boldsymbol{P}) y^​FM​(x;w,P)=⟨w,x⟩+y^​A2​(x;1,P)

是以，我們可以使用坐标梯度法（Coordinate Descent, CD）來訓練。算法如下：

可以看到，内循環周遊了各個坐标方向，并分别使用梯度下降方法進行訓練。這樣顯著提升了訓練的效果。

Homogeneous Polynomial Kernel

以上讨論了目标函數使用ANOVA Kernel的情況。當 K \mathcal K K使用同質的多項式Kernel之時，其定義如下：對于 p = [ p 1 , … , p d ] T \boldsymbol{p}=\left[p_{1}, \dots, p_{d}\right]^{\mathrm{T}} p=[p1​,…,pd​]T以及 x = [ x 1 , … , x d ] T \boldsymbol{x}=\left[x_{1}, \dots, x_{d}\right]^{\mathrm{T}} x=[x1​,…,xd​]T

H m ( p , x ) = ∑ j 1 = 1 d … ∑ j m = 1 d p j 1 x j 1 … p j m x j n = ⟨ p , x ⟩ m \mathcal{H}^{m}(\boldsymbol{p}, \boldsymbol{x})=\sum_{j_{1}=1}^{d} \ldots \sum_{j_{m}=1}^{d} p_{j_{1}} x_{j_{1}} \ldots p_{j_{m}} x_{j_{n}} = \langle p, x\rangle^{m} Hm(p,x)=j1​=1∑d​…jm​=1∑d​pj1​​xj1​​…pjm​​xjn​​=⟨p,x⟩m

可以看到， H m \mathcal{H}^{m} Hm與 A m \mathcal{A}^{m} Am不同的是前者考慮了所有不同特征之間的組合。是以，其并不具有多線性的性質。其預測函數的表達式如下：

y ^ P N ( x ; w , λ , P ) = ⟨ w , x ⟩ + y ^ H 2 ( x ; λ , P ) \hat{y}_{P N}(\boldsymbol{x} ; \boldsymbol{w}, \boldsymbol{\lambda}, \boldsymbol{P})=\langle\boldsymbol{w}, \boldsymbol{x}\rangle+\hat{y}_{\mathcal{H}^{2}}(\boldsymbol{x} ; \boldsymbol{\lambda}, \boldsymbol{P}) y^​PN​(x;w,λ,P)=⟨w,x⟩+y^​H2​(x;λ,P)

Lifted Approach

正是因為 H \mathcal H H沒有多線性，是以作者提供了一種Lifted的方法。具體而言，Lifted方法包含以下兩個部分:

1. 首先，将整個Kernel轉化為低階的張量。
2. 将低階張量對稱化。具體方法為特征值分解（這種方法在解決Convex FM[2]時被廣泛使用）

具體的Tensor轉換可以通過引理證明得到，最後ANOVA Kernel和Homogeneous Polynomial Kernel的低階張量變換結果為：

W = ∑ s = 1 k λ s p s ⊗ m \mathcal{W}=\sum_{s=1}^{k} \lambda_{s} p_{s}^{\otimes m} W=s=1∑k​λs​ps⊗m​

并且令 λ = [ λ 1 , … , λ k ] T \lambda=\left[\lambda_{1}, \ldots, \lambda_{k}\right]^{\mathrm{T}} λ=[λ1​,…,λk​]T以及 P = [ p 1 , … , p k ] T P=\left[p_{1}, \ldots, p_{k}\right]^{\mathrm{T}} P=[p1​,…,pk​]T，得到：

⟨ W , x ⊗ m ⟩ = y ^ H m ( x ; λ , P ) ⟨ W , x ⊗ m ⟩ > = y ^ A m ( x ; λ , P ) \begin{aligned}\left\langle\mathcal{W}, \boldsymbol{x}^{\otimes m}\right\rangle &=\hat{y}_{\mathcal{H}^{m}}(\boldsymbol{x} ; \boldsymbol{\lambda}, \boldsymbol{P}) \\\left\langle\mathcal{W}, \boldsymbol{x}^{\otimes m}\right\rangle_{>} &=\hat{y}_{\mathcal{A}^{m}}(\boldsymbol{x} ; \boldsymbol{\lambda}, \boldsymbol{P}) \end{aligned} ⟨W,x⊗m⟩⟨W,x⊗m⟩>​​=y^​Hm​(x;λ,P)=y^​Am​(x;λ,P)​

其loss function為:

L H m ( W ) : = ∑ i = 1 n ℓ ( y i , ⟨ W , x i ⊗ m ⟩ ) L A m ( W ) : = ∑ i = 1 n ℓ ( y i , ⟨ W , x i ⊗ m ⟩ > ) \begin{aligned} L_{\mathcal{H}^{m}}(\mathcal{W}) & :=\sum_{i=1}^{n} \ell\left(y_{i},\left\langle\mathcal{W}, x_{i}^{\otimes m}\right\rangle\right) \\ L_{\mathcal{A}^{m}}(\mathcal{W}) & :=\sum_{i=1}^{n} \ell\left(y_{i},\left\langle\mathcal{W}, x_{i}^{\otimes m}\right\rangle_{>}\right) \end{aligned} LHm​(W)LAm​(W)​:=i=1∑n​ℓ(yi​,⟨W,xi⊗m​⟩):=i=1∑n​ℓ(yi​,⟨W,xi⊗m​⟩>​)​

可以看到，其通過張量化的操作将複雜的非凸目标函數問題轉化為了凸優化問題。

HOFM

上文的工作從另一個思路總結了FM的解決思路，也為更高階FM提供了一種真正有用的工具——ANOVA Kernel。正是由于此Kernel所具有的多線性，且具有遞歸特性，在[3]中，定義了如下目标函數：

y ^ H O F M ( x ) = ⟨ w , x ⟩ + ∑ s = 1 k 2 A 2 ( p s ( 2 ) , x ) + ⋯ + ∑ s = 1 k m A m ( p s ( m ) , x ) \hat{y}_{\mathrm{HOFM}}(\boldsymbol{x})=\langle\boldsymbol{w}, \boldsymbol{x}\rangle+\sum_{s=1}^{k_{2}} \mathcal{A}^{2}\left(\boldsymbol{p}_{s}^{(2)}, \boldsymbol{x}\right)+\cdots+\sum_{s=1}^{k_{m}} \mathcal{A}^{m}\left(\boldsymbol{p}_{s}^{(m)}, \boldsymbol{x}\right) y^​HOFM​(x)=⟨w,x⟩+s=1∑k2​​A2(ps(2)​,x)+⋯+s=1∑km​​Am(ps(m)​,x)

那麼當 m m m比較大時，為了估計 A m \mathcal A^m Am，我們可以簡寫為： a j , t : = A t ( p 1 : j , x 1 : j ) a_{j, t} :=\mathcal{A}^{t}\left(p_{1 : j}, x_{1 : j}\right) aj,t​:=At(p1:j​,x1:j​)，那麼我們可以定義狀态轉移方程：

a j , t = a j − 1 , t + p j x j a j − 1 , t − 1 ∀ d ≥ j ≥ t ≥ 1 a_{j, t}=a_{j-1, t}+p_{j} x_{j} a_{j-1, t-1} \quad \forall d \geq j \geq t \geq 1 aj,t​=aj−1,t​+pj​xj​aj−1,t−1​∀d≥j≥t≥1

那麼我們就可以使用動态規劃的方法來解決這個問題。通過SGD優化，我們可以解決高階的特征訓練的問題。利用chain rule, 高階梯度的計算如下：

a ~ j , t = ∂ a d , m ∂ a j + 1 , t ∂ a j + 1 , t ∂ a j , t + ∂ a d , m ∂ a j + 1 , t + 1 ∂ a j + 1 , t + 1 ∂ a j , t = a ~ j + 1 , t + p j + 1 x j + 1 a ~ j + 1 , t + 1 ∀ d − 1 ≥ j ≥ t ≥ 1 \tilde{a}_{j, t}=\frac{\partial a_{d, m}}{\partial a_{j+1, t}} \frac{\partial a_{j+1, t}}{\partial a_{j, t}}+\frac{\partial a_{d, m}}{\partial a_{j+1, t+1}} \frac{\partial a_{j+1, t+1}}{\partial a_{j, t}}=\tilde{a}_{j+1, t}+p_{j+1} x_{j+1} \tilde{a}_{j+1, t+1} \quad \forall d-1 \geq j \geq t \geq 1 a~j,t​=∂aj+1,t​∂ad,m​​∂aj,t​∂aj+1,t​​+∂aj+1,t+1​∂ad,m​​∂aj,t​∂aj+1,t+1​​=a~j+1,t​+pj+1​xj+1​a~j+1,t+1​∀d−1≥j≥t≥1

使用下述算法，可以實作最基本的HOFM的訓練。

當然，根據[1]的叙述，并考慮到DP table的更新的同步性，可以得到坐标下降算法的遞歸公式如下：

A m ( p , x ) = 1 m ∑ t = 1 m ( − 1 ) t + 1 A m − t ( p , x ) D t ( p , x ) \mathcal{A}^{m}(\boldsymbol{p}, \boldsymbol{x})=\frac{1}{m} \sum_{t=1}^{m}(-1)^{t+1} \mathcal{A}^{m-t}(\boldsymbol{p}, \boldsymbol{x}) \mathcal{D}^{t}(\boldsymbol{p}, \boldsymbol{x}) Am(p,x)=m1​t=1∑m​(−1)t+1Am−t(p,x)Dt(p,x)

其中 D t ( p , x ) : = ∑ j = 1 d ( p j x j ) t \mathcal{D}^{t}(\boldsymbol{p}, \boldsymbol{x}) :=\sum_{j=1}^{d}\left(p_{j} x_{j}\right)^{t} Dt(p,x):=∑j=1d​(pj​xj​)t，可以得到對每個次元的梯度為：

∂ A m ( p , x ) ∂ p j = 1 m ∑ t = 1 m ( − 1 ) t + 1 [ ∂ A m − t ( p , x ) ∂ p j D t ( p , x ) + A m − t ( p , x ) ∂ D t ( p , x ) ∂ p j ] \frac{\partial \mathcal{A}^{m}(\boldsymbol{p}, \boldsymbol{x})}{\partial p_{j}}=\frac{1}{m} \sum_{t=1}^{m}(-1)^{t+1}\left[\frac{\partial \mathcal{A}^{m-t}(\boldsymbol{p}, \boldsymbol{x})}{\partial p_{j}} \mathcal{D}^{t}(\boldsymbol{p}, \boldsymbol{x})+\mathcal{A}^{m-t}(\boldsymbol{p}, \boldsymbol{x}) \frac{\partial \mathcal{D}^{t}(\boldsymbol{p}, \boldsymbol{x})}{\partial p_{j}}\right] ∂pj​∂Am(p,x)​=m1​t=1∑m​(−1)t+1[∂pj​∂Am−t(p,x)​Dt(p,x)+Am−t(p,x)∂pj​∂Dt(p,x)​]

這樣，考慮到Loss function，可以得到每一步的更新條件為：

p j , s ← p j , s − η j , s − 1 ∂ F ( P ) ∂ p j , s p_{j, s} \leftarrow p_{j, s}-\eta_{j, s}^{-1} \frac{\partial F(\boldsymbol{P})}{\partial p_{j, s}} pj,s​←pj,s​−ηj,s−1​∂pj,s​∂F(P)​

其中：

η j , s : = μ n ∑ i = 1 n ( ∂ A m ( p s , x i ) ∂ p j , s ) 2 + β \eta_{j, s} :=\frac{\mu}{n} \sum_{i=1}^{n}\left(\frac{\partial \mathcal{A}^{m}\left(\boldsymbol{p}_{s}, \boldsymbol{x}_{i}\right)}{\partial p_{j, s}}\right)^{2}+\beta ηj,s​:=nμ​i=1∑n​(∂pj,s​∂Am(ps​,xi​)​)2+β

∂ F ( P ) ∂ p j , s = 1 n ∑ i = 1 n ℓ ′ ( y i , y ^ i ) ∂ A m ( p s , x i ) ∂ p j , s + β p j , s \frac{\partial F(\boldsymbol{P})}{\partial p_{j, s}}=\frac{1}{n} \sum_{i=1}^{n} \ell^{\prime}\left(y_{i}, \hat{y}_{i}\right) \frac{\partial \mathcal{A}^{m}\left(\boldsymbol{p}_{s}, \boldsymbol{x}_{i}\right)}{\partial p_{j, s}}+\beta p_{j, s} ∂pj,s​∂F(P)​=n1​i=1∑n​ℓ′(yi​,y^​i​)∂pj,s​∂Am(ps​,xi​)​+βpj,s​

高階FM的實作可以使用tffm，它是基于Tensorflow的。關于HOFM的實作和代碼解析将在以後詳細說明。

Reference

[1] Blondel M, Ishihata M, Fujino A, et al. Polynomial networks and factorization machines: New insights and efficient training algorithms[J]. arXiv preprint arXiv:1607.08810, 2016.

[2] Blondel M, Fujino A, Ueda N. Convex factorization machines[C]//Joint European Conference on Machine Learning and Knowledge Discovery in Databases. Springer, Cham, 2015: 19-35.

[3] Blondel M, Fujino A, Ueda N, et al. Higher-order factorization machines[C]//Advances in Neural Information Processing Systems. 2016: 3351-3359.