简介

本文要介绍的是S.Rendle在2010年提出的FM（Factorization Machine）模型，此模型的提出是为了解决在数据极其稀疏的情况下的特征组合问题。FM模型跟SVM模型类似，都是一个通用的预测器，但是FM模型可以在数据极其稀疏的情况下估计可靠的模型参数。FM模型对变量之间的嵌套交互进行建模（类似多项式核函数SVM），但是却是用因子分解参数化的方式，而SVM中用的是稠密参数化的方式，这使得FM相比SVM的参数少了很多，更加容易计算，并且无需存储额外的训练数据（比如SVM中的支持向量）。

稀疏数据下的预测

通用的预测任务就是要估计一个函数：
$\hat y: \mathbb R^n → \mathbb T$ 该函数将一个 $n$ 维的实值特征向量 $x \in \mathbb R^n$ 映射到一个目标域 $T$ 。例如，对于回归 $T = \mathbb R$ ，对于分类问题 $T = \{+, -\}$ 。在有监督学习场景中，通常还有一个带标签的训练数据集：
$D = \{ (x^{(1)},y^{(1)}), (x^{(2)},y^{(2)}),...,(x^{(N)},y^{(N)}) \}$ 其中 $x^{(i) \in \mathbb R}$ 表示输入数据，对应样本的特征向量， $y^{(i)}$ 对应标签， $N$ 为样本的数量。
FM模型要处理的数据 $x$ 是高度稀疏的。举个例子，向量 $x$ 中的元素 $x_i$ 几乎大多数都是0。令 $m(x)$ 为向量 $x$ 中所有非零元素的总和， $\overline m_D$ 是所有的向量 $x \in D$ 中非零元素 $m(x)$ 的平均值。现实世界中的数据中存在着巨大的稀疏性，即( $\overline m_D <<n$ )。比如文本分析，推荐系统等。
下面以电影评分系统为例，给出一个关于高度稀疏数据的实例。

例1 考虑电影评分中的数据，每一条都包含了哪个用户 $u \in U$ 在什么时候 $t \in \mathbb R$ 对哪部电影 $i \in I$ 打了多少分 $r \in \{1,2,3,4,5\}$ 这样的信息。假定用户集 $U$ 和电影集 $I$ 如下：

设观察到的数据集

S

为：

利用观测数据集

S

来进行预测任务的一个实例是：估计一个函数

\hat y

来预测某个用户在某个时刻对某部电影的打分行为。根据这些数据，我们可以创建一个关于特征向量和标签的表，如下：

图1

图1中的标签部分比较简单，就是把用户对当前电影的评分数据当做标签。比如

Alice

对电影TI的评分是5，因此第一行最后一列的数据为

y^{(1)}=5

。特征向量

x

包含5个部分，分别对应图1中5个不同颜色的矩形框。

蓝色矩形框表示为当前用户信息，其维度为 $|U|$ ，该部分的分量中，当前用户所在的位置为1，其余为0，即是one-hot向量的表示方法。
红色矩形框表示当前评分电影信息，其维度为 $|I|$ 。当前被评分的电影所在的位置为1，其余为0，也是one-hot向量的表示方式。
黄色矩形框代表当前评分用户评分过的所有电影信息，其维度为 $I$ 。该部分的分量中，被当前用户评分过的所有电影（设个数为n）的位置为 $1/n$ ，其余位置为0。比如 $Alice$ 一共评价了3部电影TI，NH和SW，因此这3部电影所在的位置均为1/3=0.3。
绿色矩形框代表评分日期信息，其维度为1，表示当前用户评分的时间。表示方式是将2009年1月作为基数1，以后每增加1个月就加1，如2009年5月可以换算为5。
绿色矩形框代表当前评分用户最近评分过的一部电影的信息，其维度为 $I$ ，在该部分分量中，若当前用户评价当前电影之前还评价过其他电影，则当前用户评价的上一部电影的位置取1，其他为0。若当前用户评价当前电影之前没有评价过其他电影，则所有分量取0。

在本例中，特征向量 $x$ 的维度为 $|U|+|I|+|I|+1+|I| = |U| + 3|I| + 1$ 。在真实的电影评分系统中，用户数量 $|U|$ 和电影数目 $|I|$ 都非常大，而每个用户参与评分的电影数目则相对很小，于是，可想而知，每一条记录对应的特征向量会是多么的稀疏。

模型方程

1.线性回归模型

介绍FM模型方程之前，先回顾一下线性回归方程。对于一个给定的特征向量 $x = (x_1,x_1,...,x_n)^T$ ，线性回归的建模方程为：

其中 $w_0$ 和 $w=(w_1,w_2,...,w_n)^T$ 为模型参数。其中 $w_i$ 分别对应不同特征分量的权重，即代替了不同特征分量的重要程度。线性回归的优点在于可解释性强，易扩展，效率高，因而在CTR领域还有着较为广泛的运用。不过从上式中也可以很容易地看出它的缺陷，每个特征分量 $x_i$ 和 $x_j(i \ne j)$ 之间是相互独立的，即 $\hat y(x)$ 中仅仅考虑单个特征分量，而没有考虑特征分量之间的相互关系。

2.二阶多项式回归模型

在实际应用中，组合特征是很有意义的。比如在新闻推荐系统中, 喜欢军事新闻的也很可能喜欢国际新闻, 喜欢时尚新闻的也很可能喜欢娱乐新闻。因此，我们把特征的两两组合加到线性回归模型中，就可以得到二阶多项式回归模型：

其中 $n$ 代表样本特征维度，截距 $w_0 \in \mathbb R$ ， $w = {w_1, w_2,...,w_n} \in \mathbb R^n, w_{ij} \in \mathbb R^n$ 为模型参数。这样一来，就可以将任意两个互异的特征分量之间的关系也考虑进来了。
但是从上式中可以发现，交叉项的系数 $w_{ij}$ 一共有 $C_n^2 = \frac {n(n-1)}{2}$ 个，并且彼此之间是互相独立的，在数据高度稀疏的情况下，这会导致交叉项系数学习困难。原因如下：
我们知道，回归模型的参数 $w$ 的学习结果就是从训练样本中计算充分统计量（凡是符合指数簇分布的模型都具有此性质），而在这里交叉项的每一个参数 $w_{ij}$ 的学习过程都需要大量的 $x_i$ 、 $x_j$ 同时非零的训练样本数据。由于样本数据本来就很稀疏，能够满足" $x_i \ne 0 \ and \ x_j \ne 0$ "的样本数据就更少。训练样本不充分，学到的参数 $w_{ij}$ 就不是充分统计的结果，导致参数 $w_{ij}$ 不准确，而这会严重影响模型预测的结果和稳定性。

为什么参数 $w_{ij}$ 的更新依赖于满足" $x_i \ne 0 \ and \ x_j \ne 0$ "条件的的样本数据？
其实只需要自己动手计算 $\hat y(x)$ 对 $w_{ij}$ 的偏导数就能知道为什么了，留给读者自己思考。

下面再举个例子实际说明上面描述的问题。
仍以上文的电影评分数据为例，在观测集 $S$ 中， $Alice$ 没有对电影 $Stark Trek$ 的评分记录，如果要直接估计 $Alice(A)$ 和 $Stark Trek(ST)$ 之间，或者说特征分类 $x_A$ 和 $x_{ST}$ 之间的相互关系，显然会得到系数 $w_{A,ST}=0$ 。即对于观察样本数据中未出现过交互的特征分量，不能对相应的参数进行估计。在高度稀疏的数据场景中，由于数据量的不足，样本中出现未交互的特征分类是很常见的。

3.FM模型

为了解决这个问题，我们对系数 $w_{ij}$ 做文章，将其表示成其他形式。为此，针对每个维度的特征分量 $x_i$ ，引入辅助向量：
$v_i = (v_{i1}, v_{i2},...,v_{ik})^T, i = 1,2,...,n$ 其中 $k \in \mathbb N^{+}$ 为超参数，并将 $w_{ij}$ 改写成：
$\hat w_{i,j} = \Bbb v_i^T \Bbb v := \sum_{l=1}^k v_{il}v_{jl}$ 于是，函数 $\hat y(x)$ 可以被写成：

从数学的角度来看，这样做的好处在哪里呢？对于交叉项系数，二阶多项式回归模型使用的是参数

w_{ij}

，而后者用的是

\hat w_{ij}

，为了更好地看清它们之间的关系，引入几个矩阵：

将 $\{v_i\}_{i=1}^n$ 按行拼成的长方阵 $V$
交互矩阵 $W$
交互矩阵 $\hat W$

由此可见，二阶多项式回归模型和改进后的模型分别采用了交互矩阵 $W$ 和 $\hat W$ 的非对角线元素作为 $x_ix_j$ 的系数。由于 $\hat W = VV^T$ 对应一种矩阵分解。因此，我们把这种改进二阶多项式模型的方法叫做因子分解机(Factorization Machines)方法。
读者可能要问了， $VV^T$ 的表达能力怎么样呢？即任意给定一个交互矩阵 $\hat W$ ，能否找到相应的(分解)矩阵 $V$ 呢？答案是肯定的，这里需要用到一个结论“当 $k$ 足够大时，对于任意对称正定的实矩阵 $\hat W \in \mathbb R^{n \times n}$ ，均存在实矩阵 $V \in \mathbb R^{n \times k}$ ,使得 $\hat W = VV^T$ 成立”。

将每个交叉项系数 $w_{ij}$ 用隐向量的内积<v_i,v_j>表示，是FM模型的核心思想。具体地说，FM为每个特征学习了一个隐权重向量，在特征交叉时，使用两个特征隐向量的内积作为交叉特征的权重。
下面给出论文中的FM方程的形式：

其中

v_i

代表第

i

个特征的隐向量，

<.,.>

代表的是两个长度为

k

的向量的内积，计算公式为：

隐向量的长度

k

是一个超参数

(k \in \mathbb N^+, k<<n)

，隐向量

v_i = (v_{i,1}, v_{i,2},...,v_{i,k})

的含义是用

k

个描述特征的因子来描述第

i

维特征。这样做之后，交叉项的参数由原来的

\frac {n(n-1)}{2}

降低为

nk

个，远少于二阶多项式模型中的参数数量。
此外，参数因子化表示后，使得 $x_hx_i$ 的参数与 $x_ix_j$ 的参数不再相互独立。这样我们就可以在样本稀疏的情况下相对合理地估计FM交叉项的参数。

FM复杂度分析

先列出FM的方程：

其中需要估计的参数包括：

1+n+nk

个，其中

w_0

是整体的偏移量，

w

是对特征向量的各个分量的强度进行建模，

V

是对特征向量中任意两个分量之间的关系进行建模。那么根据公式来看，其计算复杂度为：

其中第一个花括号对应

\sum_{i=1}^{n} w_{i} x_{i}

的加法和乘法操作数，第二个花括号对应

\sum_{i=1}^{n-1} \sum_{j=i+1}^{n}\left(\mathbf{v}_{i}^{\top} \mathbf{v}_{j}\right) x_{i} x_{j}

的加法和乘法操作数，最后那个2代表整体的两次加法。
这样指数级的复杂度显然是不能够运用到大规模的实际应用中的，不过好在通过对上式改写，可以将计算复杂度降低为线性的

O(kn)

，论文给出了具体过程:

分析一下改进后的时间复杂度，为：

注意，在高度稀疏数据场景中，特征向量

x

中绝大部分分量为0，即

m(x)

很小，而上式中关于

i

的求和只需要对非零元素进行。于是，计算复杂度变成

\mathcal{O} \left(k m( \mathbf{ x }) \right) \

了。

这里需要讲一下这个公式的第（1）步到第（2）步的过程，其他的推导由于篇幅原因略过，如果读者想深入了解，请参考https://zhuanlan.zhihu.com/p/58160982。
先计算一下：

将结果写成矩阵的形式：

其中红色的上三角部分就是我们要求的。
设该上三角为 $A$ ，则可推导如下：

其实也可以用排列组合的角度来思考这个问题， $\sum_{i=1}^{n} \sum_{j=i+1}^{n}<v_{i}, v_{j}>x_{i} x_{j}$ 一共有 $C_n^2 = \frac {n(n-1)}{2} = \frac {1}{2} n^2 - \frac {1}{2} n$ 种组合方式，而
$\sum_{i=1}^{n} \sum_{j=1}^{n}\left\langle\mathbf{v}_{i} \mathbf{v}_{j}\right\rangle x_{i} x_{j} \$
一共有 $n^2$ 种组合方式。
$\sum_{i=1}^{n} \left \langle \mathbf{v}_{i}, \mathbf{v}_{i} \right \rangle x_{i} x_{i}$
一共有 $n$ 种组合方式，故从排列组合的方式来看，上式也是成立的。

如果用随机梯度下降的方式来学习模型参数。那么，模型各个参数的梯度如下：

其中，

v_{j,f}

是隐向量

v_j

的第

f

个元素。

最优化问题

FM的优化目标是整体损失函数：

其中，对于回归问题，loss取最小平方误差：

对于二分类问题，loss取logit函数：

于是，FM的最优化问题变成了：

我们通常还需要加上

L2

正则来防止过拟合，因此最优化问题变成了：

其中

\lambda_0

表示参数

\theta

的正则化系数。

代码实践

虽然FM的公式看起来很复杂，但是代码实现起来其实很简单，模型部分代码如下：

import torch
import torch.nn as nn
from BaseModel.basemodel import BaseModel

class FM(BaseModel):
    def __init__(self, config):
        super(FM, self).__init__(config)
        # 特征的个数
        self.p = 13 #config['num_features']
        # 隐向量的维度
        self.k = config['latent_dim']
        # FM的线性部分，即∑WiXi
        self.linear = nn.Linear(self.p, 1, bias=True)
        # 隐向量的大小为nxk,即为每个特征学习一个维度为k的隐向量
        self.v = nn.Parameter(torch.randn(self.k, self.p))
        # nn.init.uniform_(self.v, -0.1, 0.1)

    def forward(self, x):
        x = x[:, :13]
        # 线性部分
        linear_part = self.linear(x)
        # 矩阵相乘 (batch*p) * (p*k)
        inter_part1 = torch.mm(x, self.v.t())
        # 矩阵相乘 (batch*p)^2 * (p*k)^2
        inter_part2 = torch.mm(torch.pow(x, 2), torch.pow(self.v, 2).t())
        output = linear_part + 0.5 * torch.sum(torch.pow(inter_part1, 2) - inter_part2)
        return output

使用criteo数据集来做训练，但是对于类别数据要先转成one-hot向量编码的形式，criteo数据集的一条记录包含39个特征，其中前13个是连续特征，后26个是类别特征。13个连续特征进行归一化处理之后保持不动，将26个类别特征先转换成one-hot向量形式，然从concatenate起来，组成一个高维的稀疏向量输入到FM模型中去训练。采用SGD优化器，MSE损失函数来训练。

在训练的时候，发现训练一两次，loss就变成了nan。这是由于loss数值太大，已经无法使用float或者double来表示了。主要可能是数据本身有问题，学习率设置过大等。我自己降低了学习率后，模型可以优化训练。但是不保证代码本身没有问题，仅供参考。

测试部分代码：

import torch
from FM.trainer import Trainer
from FM.network import FM
from Utils.criteo_loader import getTestData, getTrainData
import torch.utils.data as Data

fm_config = \
{
    'latent_dim': 10,
    'num_dense_features': 13, # for criteo data set
    'num_epoch': 10,
    'batch_size': 64,
    'lr': 1e-6,
    'l2_regularization': 1e-3,
    'device_id': 0,
    'use_cuda': False,
    'train_file': '../Data/criteo/processed_data/train_set.csv',
    'fea_file': '../Data/criteo/processed_data/fea_col.npy',
    'validate_file': '../Data/criteo/processed_data/val_set.csv',
    'test_file': '../Data/criteo/processed_data/test_set.csv',
    'model_name': '../TrainedModels/fm.model'
}

def toOneHot(x, MaxList):
    res = []
    for i in range(len(x)):
        t = torch.zeros(MaxList[i])
        t[int(x[i])] = 1
        res.append(t)
    return torch.cat(res, -1)

if __name__ == "__main__":
    ####################################################################################
    # FM 模型
    ####################################################################################
    training_data, training_label, dense_features_col, sparse_features_col = getTrainData(fm_config['train_file'], fm_config['fea_file'])
    p = sum(sparse_features_col) + fm_config['num_dense_features']

    rows, cols = training_data.shape
    train_x = torch.zeros((rows, p))
    for row in range(rows):
        dense = torch.Tensor(training_data[row][:fm_config['num_dense_features']])
        sparse = training_data[row][fm_config['num_dense_features']:]
        sparse = toOneHot(sparse, sparse_features_col)
        train_x[row] = torch.cat((dense, sparse),0)

    train_dataset = Data.TensorDataset(train_x.float().clone().detach().requires_grad_(True), torch.tensor(training_label).float())
    test_data = getTestData(fm_config['test_file'])
    test_dataset = Data.TensorDataset(torch.tensor(test_data).float())

    fm = FM(fm_config, p)

    ####################################################################################
    # 模型训练阶段
    ####################################################################################
    # # 实例化模型训练器
    trainer = Trainer(model=fm, config=fm_config)
    # 训练
    trainer.train(train_dataset)
    # 保存模型
    trainer.save()

完整代码见https://github.com/HeartbreakSurvivor/RsAlgorithms/tree/main/FM。

参考

S. Rendle. Factorization machines. In ICDM, 2010.
《深度学习推荐系统》-- 王喆
https://www.cnblogs.com/pinard/p/6370127.html
http://www.52caml.com/head_first_ml/ml-chapter9-factorization-family/
https://www.jianshu.com/p/25f0139637b7
http://stackbox.cn/2018-12-factorization-machine/
https://zhuanlan.zhihu.com/p/58160982
https://zhuanlan.zhihu.com/p/72586223
https://zhuanlan.zhihu.com/p/91151350
https://www.cnblogs.com/wyhluckdog/p/12168087.html
http://shomy.top/2018/12/31/factorization-machine/
https://github.com/1JasonZhang/FM-by-pytorch/blob/master/fm_model.py

推荐系统之FM（因子分解机）模型原理以及代码实践