Factorization Machine 详解（三）—— FM与高阶特征训练

Factorization Machine 详解（三）—— FM与高阶特征训练

之前的文章回顾：

• FM基础概念的介绍与训练方法
• FM开源库LibFM简介与使用

可以看到，FM主要用于推测推荐系统中不同特征之间的关联，尤其适用于稀疏数据的推荐。然而，Rendle等提出的基础FM模型一般只能推测2-way的特征关系，而对于更高阶的特征关系则缺乏更为有效的训练手段。另外，在传统的FM模型中，目标函数是非凸的，因此很难进行优化，尤其是当涉及高阶特征关系的挖掘时，目标函数的优化更为复杂。
对于高阶特征的训练，目前有以下两个思路：

• 使用更高阶的Kernel函数对目标函数进行建模
• 使用深度神经网络对输入数据进行处理

目前，在第一种解决方案中，我们需要对目标函数进行建模：
$${\hat{y}}_{\mathcal{K}}(\mathbf{x};\mathbf{\lambda },\mathbf{P}):=\sum _{s=1}^k{\lambda }_s\mathcal{K}\left({\mathbf{p}}_s,\mathbf{x}\right)$$
总的loss function：
$$D_{\mathcal{K}}(\mathbf{\lambda },\mathbf{P}):=\sum _{i=1}^n\ell \left(y_i,{\hat{y}}_{\mathcal{K}}\left({\mathbf{x}}_i;\mathbf{\lambda },\mathbf{P}\right)\right)$$
其中函数$\ell$为凸的loss function。从中我们可以看到，其使用的目标函数是一个Kernel，即为从$\mathbf{P}$映射到输入向量$\mathbf{x}$空间的Kernel。在正式阐述高阶FM的训练方式之前，先说明一下上述参数的具体意义：

• ${\hat{y}}_{\mathcal{K}}$表示的是基于Kernel $\mathcal{K}$的目标函数。
• $\mathbf{\lambda }={\left[{\lambda }_1,\dots ,{\lambda }_k\right]}^{\mathrm{T}}\in {\mathbb{R}}^k$ 是超参，具体意义是决定不同维度Kernel的权重。在FM中${\lambda }_s=1(s=1,2,...,k)$。
• $\mathbf{P}=\left[{\mathbf{p}}_1,\dots ,{\mathbf{p}}_k\right]\in {\mathbb{R}}^{d\times k}$表示的是已经分解的向量（矩阵）
• $\mathcal{K}$表示的是Kernel函数。目前研究的Kernel包括ANOVA Kernel以及同质的多项式Kernel

ANOVA Kernel

2016年，Mathieu Blondel等[1]解决了Polynomial Network(PN) 以及Factorization Machine(FM)之间的关系，并且证明了提出了coordinate descent(CD)算法，以训练三阶的FM。关于ANOVA Kernel，定义如下：
$${\mathcal{A}}^m(\mathbf{p},\mathbf{x}):=\sum _{j_m>\cdots >j_1}{p}_{j_1}{x}_{j_1}\dots p_{j_m}{x}_{j_m}$$
其中${\mathcal{A}}^0(\mathbf{p},\mathbf{x}):=1$以及${\mathcal{A}}^1(\mathbf{p},\mathbf{x}):=⟨\mathbf{p},\mathbf{x}⟩$。可以看到ANOVA Kernel使用的是多种distinct 的特征（也就是没有替代的特征组合）。ANOVA Kernel具有多线性的性质。通过展开ANOVA Kernel的定义式，我们可以使用数学归纳法证明，对于$\mathbf{p},\mathbf{x}\in {\mathbb{R}}^d$，以及$j\in [d]$和$1\le m\le d$，有如下的递归多线性（multi-linearity）关系：
$${\mathcal{A}}^m(\mathbf{p},\mathbf{x})={\mathcal{A}}^m\left({\mathbf{p}}_{\neg j},{\mathbf{x}}_{\neg j}\right)+p_j{x}_j{\mathcal{A}}^{m-1}\left({\mathbf{p}}_{\neg j},{\mathbf{x}}_{\neg j}\right)$$
其中$\neg j$表示的是向量中排除掉第$j$个元素的内容，因此当我们对其中第$j$个坐标$p_j$计算导数时，上式是一个线性函数。因此，我们可以理解多线性的含义为对某个特定的方向$j$，都满足既凸又凹的线性关系。那么，根据Rendle的FM定义，我们可以改写FM的目标函数如下：
$${\hat{y}}_{FM}(\mathbf{x};\mathbf{w},\mathbf{P})=⟨\mathbf{w},\mathbf{x}⟩+{\hat{y}}_{{\mathcal{A}}^2}(\mathbf{x};1,\mathbf{P})$$
因此，我们可以使用坐标梯度法（Coordinate Descent, CD）来训练。算法如下：

可以看到，内循环遍历了各个坐标方向，并分别使用梯度下降方法进行训练。这样显著提升了训练的效果。

Homogeneous Polynomial Kernel

以上讨论了目标函数使用ANOVA Kernel的情况。当$\mathcal{K}$使用同质的多项式Kernel之时，其定义如下：对于$\mathbf{p}={\left[p_1,\dots ,p_d\right]}^{\mathrm{T}}$以及$\mathbf{x}={\left[x_1,\dots ,x_d\right]}^{\mathrm{T}}$
$${\mathcal{H}}^m(\mathbf{p},\mathbf{x})=\sum _{j_1=1}^d\dots \sum _{j_m=1}^d{p}_{j_1}{x}_{j_1}\dots p_{j_m}{x}_{j_n}=⟨p,x{⟩}^m$$
可以看到，${\mathcal{H}}^m$与${\mathcal{A}}^m$不同的是前者考虑了所有不同特征之间的组合。因此，其并不具有多线性的性质。其预测函数的表达式如下：
$${\hat{y}}_{PN}(\mathbf{x};\mathbf{w},\mathbf{\lambda },\mathbf{P})=⟨\mathbf{w},\mathbf{x}⟩+{\hat{y}}_{{\mathcal{H}}^2}(\mathbf{x};\mathbf{\lambda },\mathbf{P})$$

Lifted Approach

正是因为$\mathcal{H}$没有多线性，因此作者提供了一种Lifted的方法。具体而言，Lifted方法包含以下两个部分:

1. 首先，将整个Kernel转化为低阶的张量。
2. 将低阶张量对称化。具体方法为特征值分解（这种方法在解决Convex FM[2]时被广泛使用）

具体的Tensor转换可以通过引理证明得到，最后ANOVA Kernel和Homogeneous Polynomial Kernel的低阶张量变换结果为：
$$\mathcal{W}=\sum _{s=1}^k{\lambda }_s{p}_s^{\otimes m}$$
并且令$\lambda ={\left[{\lambda }_1,\dots ,{\lambda }_k\right]}^{\mathrm{T}}$以及$P={\left[p_1,\dots ,p_k\right]}^{\mathrm{T}}$，得到：
$$\begin{array}{rl}⟨\mathcal{W},{\mathbf{x}}^{\otimes m}⟩& ={\hat{y}}_{{\mathcal{H}}^m}(\mathbf{x};\mathbf{\lambda },\mathbf{P})\\ {⟨\mathcal{W},{\mathbf{x}}^{\otimes m}⟩}_{>}& ={\hat{y}}_{{\mathcal{A}}^m}(\mathbf{x};\mathbf{\lambda },\mathbf{P})\end{array}$$
其loss function为:
$$\begin{array}{rl}{L}_{{\mathcal{H}}^m}(\mathcal{W})& :=\sum _{i=1}^n\ell \left(y_i,⟨\mathcal{W},x_i^{\otimes m}⟩\right)\\ L_{{\mathcal{A}}^m}(\mathcal{W})& :=\sum _{i=1}^n\ell \left(y_i,{⟨\mathcal{W},x_i^{\otimes m}⟩}_{>}\right)\end{array}$$
可以看到，其通过张量化的操作将复杂的非凸目标函数问题转化为了凸优化问题。

HOFM

上文的工作从另一个思路总结了FM的解决思路，也为更高阶FM提供了一种真正有用的工具——ANOVA Kernel。正是由于此Kernel所具有的多线性，且具有递归特性，在[3]中，定义了如下目标函数：
$${\hat{y}}_{\mathrm{H}\mathrm{O}\mathrm{F}\mathrm{M}}(\mathbf{x})=⟨\mathbf{w},\mathbf{x}⟩+\sum _{s=1}^{k_2}{\mathcal{A}}^2\left({\mathbf{p}}_s^{(2)},\mathbf{x}\right)+\cdots +\sum _{s=1}^{k_m}{\mathcal{A}}^m\left({\mathbf{p}}_s^{(m)},\mathbf{x}\right)$$
那么当$m$比较大时，为了估计${\mathcal{A}}^m$，我们可以简写为：$a_{j,t}:={\mathcal{A}}^t\left(p_{1:j},x_{1:j}\right)$，那么我们可以定义状态转移方程：
$$a_{j,t}=a_{j-1,t}+p_j{x}_j{a}_{j-1,t-1}\forall d\ge j\ge t\ge 1$$
那么我们就可以使用动态规划的方法来解决这个问题。通过SGD优化，我们可以解决高阶的特征训练的问题。利用chain rule, 高阶梯度的计算如下：
$${\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}\forall d-1\ge j\ge t\ge 1$$
使用下述算法，可以实现最基本的HOFM的训练。

当然，根据[1]的叙述，并考虑到DP table的更新的同步性，可以得到坐标下降算法的递归公式如下：
$${\mathcal{A}}^m(\mathbf{p},\mathbf{x})=\frac{1}{m}\sum _{t=1}^m(-1{)}^{t+1}{\mathcal{A}}^{m-t}(\mathbf{p},\mathbf{x}){\mathcal{D}}^t(\mathbf{p},\mathbf{x})$$
其中${\mathcal{D}}^t(\mathbf{p},\mathbf{x}):={\sum }_{j=1}^d{\left(p_j{x}_j\right)}^t$，可以得到对每个维度的梯度为：
$$\frac{\partial {\mathcal{A}}^m(\mathbf{p},\mathbf{x})}{\partial p_j}=\frac{1}{m}\sum _{t=1}^m(-1{)}^{t+1}\left[\frac{\partial {\mathcal{A}}^{m-t}(\mathbf{p},\mathbf{x})}{\partial p_j}{\mathcal{D}}^t(\mathbf{p},\mathbf{x})+{\mathcal{A}}^{m-t}(\mathbf{p},\mathbf{x})\frac{\partial {\mathcal{D}}^t(\mathbf{p},\mathbf{x})}{\partial p_j}\right]$$
这样，考虑到Loss function，可以得到每一步的更新条件为：
$$p_{j,s}\leftarrow p_{j,s}-{\eta }_{j,s}^{-1}\frac{\partial F(\mathbf{P})}{\partial p_{j,s}}$$
其中：
$${\eta }_{j,s}:=\frac{\mu }{n}\sum _{i=1}^n{\left(\frac{\partial {\mathcal{A}}^m\left({\mathbf{p}}_s,{\mathbf{x}}_i\right)}{\partial p_{j,s}}\right)}^2+\beta$$
$$\frac{\partial F(\mathbf{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({\mathbf{p}}_s,{\mathbf{x}}_i\right)}{\partial p_{j,s}}+\beta p_{j,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.