原创 | 想做推荐算法？先把FM模型搞懂再说

2020-10-13 14:07:02

在上一篇文章当中我们剖析了Facebook的论文GBDT+LR，虽然这篇paper在业内广受好评，但是毕竟GBDT已经是有些老旧的模型了。今天我们要介绍一个业内使用得更多的模型，它诞生于2010年，原作者是Steffen Rendle。虽然诞生得更早，但是它的活力更强，并且衍生出了多种版本。我们今天剖析的就是这篇2010年经典的原版论文。

说到推荐、广告的算法模型，几乎很难绕开FM，它是一个非常强的模型。理论简单、推导严谨、实现容易，并且效果不俗。即使是目前仍然在各大厂商当中发挥用场，在一些小厂当中甚至是主力模型。我们初看之前也许还有疑惑，但是相信当大家看完之后，必然会有全新的认识。

同样，这是一篇长文，希望大家耐心读到后。

简介

FM的全称是Factorization Machines，翻译过来就是因式分解机，这么翻译很拗口，不得其神也不得其形，所以我们一般都不翻译简称为FM模型。

论文开篇就讲了FM模型的种种好处，比如和SVM对比，它能够在稀疏的特征当中仍然拥有很好的表现，并且FM的效率非常高，可以在线性时间内获得结果。并且不像非线性的SVM（核函数），FM并不需要将特征转换成对偶形式，并且模型的参数可以直接训练，也不用借助支持向量或者是其他的方法。

除了SVM之外，FM模型和其他的因式分解模型比如SVD、PITF、FPMC比较起来都有非常明显的优势。这些模型往往只针对一个特定的场景或者是特定的数据集，并且它们的训练以及优化方案都是根据场景定制化的。FM模型不需要定制化就可以实现同样好的效果，可以非常简易地应用在各个预测问题当中。

这段摘要其实没什么太多的内涵，主要就是吹嘘了一通FM。不过作者所言不虚，和他列举的这些模型比较起来，FM的确更加出色。

总结一下，FM模型的优点有如下几点：

FM模型的参数支持非常稀疏的特征，而SVM等模型不行
FM的时间复杂度为，并且可以直接优化原问题的参数，而不需要依靠支持向量或者是转化成对偶问题解决
FM是通用的模型，可以适用于任何实数特征的场景，其他的模型不行

paper在一开始的时候就表明了FM模型的优点，给足了好处，之后才是对FM模型的深入分析，让大家明白这些优点是如何得到的，以及它的工作原理。

稀疏场景

在有监督学习的场景当中，常见的任务就是给定一个向量x，要求预测一个目标T。如果T是一个实数，那么这就是回归模型，如果T是一个类别的常量，就是一个分类模型。

这些都是机器学习的基础知识，相信大家都了解，然而对于线上排序的功能来说，我们重要的是给商品排序，而不是分类。常规来说我们可以将impression和click看成是两个类别进行分类预测，也可以直接用回归模型预测商品的点击率，根据点击率进行排序。这两种其实本质上是一样的，都是针对商品本身进行学习。还有一种方法是更加侧重商品之间的排序，我们的模型学习的是商品之间的优劣关系。

举个例子，比如我们用向量表示i商品的特征向量，表示j商品的特征向量，那么我们可以将这两个向量合并一起作为输入，进行分类预测。分类的结果表示的是i商品在j商品之前还是反之。这样我们可以通过多次预测，直接得到商品之间的排序关系，而不是根据一个分数进行排序。这样可以在只有正样本的情况下进行训练。这种方法直接训练的模型商品的优劣，业内叫做learning to rank。

然而不论是哪一种做法，都有一个问题绕不开就是特征的稀疏。举个很简单的例子，比如我们对商品的类目进行one-hot处理，在电商场景下商品的类目动辄上万个，那么one-hot之后的结果就是一个长度上万的数组，并且这个数组当中只有一位为1，其他均为0。当然这只是一个例子，除此之外还有许多许多的特征有可能是非常稀疏的。

我们用表示x向量当中非零的元素的个数，用表示所有x样本平均的非零元素的个数，在大多数真实数据场景下，我们都可以得到，这里的n是特征的维数。

真实场景例子

paper当中举了一个真实场景的例子，我们的问题是预测用户对于一部电影的评分。我们先来观察一下下图，下图是样本矩阵。

很明显上图的左边一大块是特征，右边的Target y表示的预测结果，也就是用户可能对电影做出的评价。这里一共有[1, 2, 3, 4, 5]这5种可能，也就是说这是一个多分类的问题。

接着我们再来看特征，特征一共也有5个部分，其中蓝色的部分表示的用户的one-hot。那么这个数组的长度应该是用户的数量，只有代表当前用户的那一维为1，其他均为0。

红色部分表示电影，同样是电影的one-hot，和用户的one-hot是一样的逻辑。代表当前电影的那一维度为1，其他为0。

之后是黄色的部分，表示的之前用户对于电影的评分行为，维度同样是电影的数量。凡是用户评分过的电影分数大于0，没有评分的等于0。得分的计算逻辑是1除以用户评论过的电影数量。比如行当中，个用户评价过前三部电影，所以前三部电影每一部分到了的分数。

绿色的部分只有1维，表示的是用户评论电影的时间。计算方法是将记录当中早的日期作为基数（这里是2009年1月），之后每过一个月，增加1。比如2009年5就可以折算成5。

后棕色的部分表示的是用户近一次评论的电影，这同样是一个one-hot的数组，它的维度和电影的数量是一致的。

我们假设用户的数量是U，电影的数量是M，那么后得到的整个特征的维度数应该是U+3M+1。即使是小众一些的电影评分网站，用户数也至少是以上百万起的，再加上电影的数量，这显然是一个非常庞大的数字。而在这么庞大的维度当中只有少数的一些维度是有值的，其余均为0。

对于这样稀疏的特征矩阵而言，一般的模型是很难保证效果的。为什么FM模型可以置身其外呢？下面，我们就来看看FM模型的原理。

FM模型原理

在我们介绍FM模型的方程之前，我们先来回顾一下线性回归的表达式，其实非常简单，只有一行：

也就是说，W是这样一个n+1维的向量，X是一个n x m的特征矩阵。这里的n是特征的维数，m是样本的数量。所以我们也经常把它写成，不管怎么写，形式都是一样的。

这个式子称为线性回归的原因也很简单，因为它就是一个简单的线性方程，只有一次项。但是一次项有时候效果不好，尤其是在特别稀疏的场景当中，刻画能力不够。我们做特征的时候经常会把两项特征组合起来做成新的组合特征，由于我们这样操作引入了新的特征，找到了新的特征组合，所以能够挖掘出之前无法得到的信息，因此模型也会有更好的效果。

当我们把特征进行二项组合之后，会得到这样的式子：

这里的和分别表示两个不同的特征取值，对于n维的特征来说，这样的组合应该一共有种，这个值是，也就意味着我们需要同样数量的权重参数。但是有一个小问题，我们前面已经说过了，由于特征可能非常稀疏，导致n非常大，比如上百万，这里两两特征组合的特征量级大约是n的平方，那么因此带来的参数数量就是一个天文数字。想想看对于一个上百亿甚至更多参数空间的模型来说，我们需要多少训练样本才可以保证完全收敛？这是不可想象的。

既然如此，那么针对性的优化就是必不可少的。FM为人称道的也正是这一点，既引入了特征交叉，又解决了复杂度以及模型参数的问题，真的可以说是非常棒了。

FM解决这个问题的方法非常简单，它不再是简单地为交叉之后的特征对设置参数，而是设置了一种计算特征参数的方法。FM模型引入了新的矩阵V，矩阵V是一个n x k的二维矩阵。这里的k是我们设置的参数，一般不会很大，比如16、32之类。对于特征每一个维度i，我们都可以找到一个，它表示一个长度为k的向量。

于是我们可以用和来计算得出上式当中的：

也就是说我们用向量的内积来计算得到了就交叉特征的系数，相比于原先量级的参数而言，我们将参数的量级降低到了n x k。由于k是一个常数值，所以可以看成我们的参数数量是。通过这种方式我们把参数的数量降低了一个数量级。

有了V矩阵之后，刚才的公式可以改写成：

我们很容易可以知道，当k这个参数足够大的时候，我们可以保证成立。所以我们引入的参数矩阵V可以看成是对W矩阵做了一个因子分解，这也是FM得名的由来。当然在实际的应用场景当中，我们并不需要设置非常大的K，因为特征矩阵往往非常稀疏，我们可能没有足够多的样本来训练这么大量的参数，并且限制K也可以一定程度上提升FM模型的泛化能力。

此外这样做还有一个好处就是有利于模型训练，因为对于有些稀疏的特征组合来说，我们所有的样本当中可能都是空的。比如在刚才的例子当中用户A和电影B的组合，可能用户A在电影B上就没有过任何行为，那么这个数据就是空的，我们也不可能训练出任何参数来。但是引入了V之后，虽然这两项缺失，但是我们仍然可以得到一个参数。因为我们针对用户A和电影B训练出了一个向量参数，我们用这两个向量参数点乘，就得到了这个交叉特征的系数。

这种做法有两种理解方式，一种就是刚才说的，我们对于一些稀疏的组合也可以很好地计算出系数。另外一种理解方式是这其实也是一种embedding的方式，将某一个id映射成向量。所以业内也有使用FM来做embedding的。

复杂度优化

接下来我们来看另外一个很精髓的内容，就是关于FM模型的复杂度优化。我们观察一下刚才上面的式子，不难发现，目前对于预测一条样本的计算复杂度为。

n我们前文说过了是一个比较大的数字，那么显然级的计算也是我们不能接受的。所以对它进行优化也是必须的，并且这里的优化非常简单，可以直接通过数学公式的变形推导得到。

这个是它的原式，对于这个式子来说，前面两项的复杂度是，我们可以先忽略，重点来看后一项。我们要做的就是通过数学公式的变形来对这一项进行化简：

简单来解释一下这个推导过程，行我想大家应该都能看懂，第二行也很好理解，其实就是把向量内积展开。第二行到第三行的转化也不难理解，这里有三个，我们提取出的是里面的，因为是有限项求和，我们调换顺序并不会影响结果。提取出了公因式之后，得到的结果是两个平方项。我们观察一下可以发现，这两个平方项的计算复杂度都是，再加上外面一层的复杂度，整体的复杂度是。

这样我们就完成了FM模型预测的优化。

模型训练

我们再来看下模型训练的过程，我们写出变形之后的原式：

针对FM模型我们一样可以使用梯度下降算法来进行优化，既然要使用梯度下降，那么我们就需要写出模型当中所有参数的偏导。

我们可以把参数分成三个部分，个部分自然是，。第二个部分是，对于其中每一个来说它的系数是。由于样本是固定的，我们要把样本的特征值看成是常数。所以。后一个部分就是看起来复杂很多，但其实偏导也不难求，我们首先明确这里要优化的参数是，所以我们可以忽略外层的一层，直接对括号内的进行求导，得出的结果是

我们把这三种综合一下，就得到了：

其中和i是独立的，所以它是可以提前算好的，这样一来对于所有参数项，我们都可以在的时间内计算出它们的梯度。由于我们所有参数的量级是，意味着我们可以在时间内对所有的参数完成梯度下降的更新。

拓展到d维

我们刚才介绍的FM模型专注的是二维特征交叉的情况，如果我们愿意也可以将它拓展到d维特征交叉的情况。这样的话我们的特征可以拟合任意维特征交叉的相互关系。

我们仿照刚才的公式，可以写出FM模型推广到d维的方程：

前面两项都很好理解，我们着重来看第三项。第三项当中包含了从2维到d维交叉特征的情况，我们以d=3为例，那么这一项当中应该包含二维的交叉项以及三维的交叉项，应该是这样的：

这个式子整体上和之前的形式是一样的，我们不难分析出它的复杂度是。当d=2的时候，我们通过一系列变形将它的复杂度优化到了，而当d>2的时候，没有很好的优化方法，而且三重特征的交叉往往没有意义，并且会过于稀疏，所以我们一般情况下只会使用d = 2的情况。

代码实现

上面这段大家了解一下知道有这么回事就可以了，实际上意义不大。后的时候paper还比较了FM和其他一些经典的模型的效果，比如SVD、SVM等，也没太多价值，因为现在在推荐领域已经几乎没有人直接使用这些模型了。后我贴一下我用pytorch和TensorFlow两个框架分别实现的FM模型，给大家做一个参考。

Pytorch

import torch
from torch import nn

ndim = len(feature_names)
k = 4

class FM(nn.Module):
    def __init__(self, dim, k):
        super(FM, self).__init__()
        self.dim = dim
        self.k = k
        self.w = nn.Linear(self.dim, 1, bias=True)
        # 初始化V矩阵
        self.v = nn.Parameter(torch.rand(self.dim, self.k) / 100)
        
    def forward(self, x):
        linear = self.w(x)
        # 二次项
        quadradic = 0.5 * torch.sum(torch.pow(torch.mm(x, self.v), 2) - torch.mm(torch.pow(x, 2), torch.pow(self.v, 2)))
        # 套一层sigmoid转成分类模型，也可以不加，就是回归模型
        return torch.sigmoid(linear + quadradic)
    
    
fm = FM(ndim, k)
loss_fn = nn.BCELoss()
optimizer = torch.optim.SGD(fm.parameters(), lr=0.005, weight_decay=0.001)
iteration = 
epochs = 10

for epoch in range(epochs):
    fm.train()
    for X, y in data_iter:
        output = fm(X)
        l = loss_fn(output.squeeze(dim=1), y)
        optimizer.zero_grad()
        l.backward()
        optimizer.step()
        iteration += 1        
        if iteration % 200 == 199:
            with torch.no_grad():
                fm.eval()
                output = fm(X_eva_tensor)
                l = loss_fn(output.squeeze(dim=1), y_eva_tensor)
                acc = ((torch.round(output).long() == y_eva_tensor.view(-1, 1).long()).sum().float().item()) / 1024
                print('Epoch: {}, iteration: {}, loss: {}, acc: {}'.format(epoch, iteration, l.item(), acc))
            fm.train()

TensorFlow

import tensorflow as tf
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import accuracy_score, roc_auc_score
import numpy as np


class FactorizationMachine:

    def __init__(self, n_dim=1, k=4, learning_rate=0.05, epochs=8):
        self._learning_rate = learning_rate
        self._n_dim = n_dim
        self._k = k
        self._epochs = epochs
        self.sess = tf.Session()
        self.x_input = tf.placeholder(shape=[None, self._n_dim], dtype=tf.float32)
        self.y_input = tf.placeholder(shape=[None, 1], dtype=tf.float32)
        # 初始化W和V
        self.w = tf.Variable(tf.truncated_normal(shape=[self._n_dim, 1], dtype=tf.float32))
        self.V = tf.Variable(tf.truncated_normal(shape=[self._n_dim, self._k], dtype=tf.float32))
        self.b = tf.Variable(tf.truncated_normal(shape=[1, 1]))
        self.linear = tf.add(self.b, tf.matmul(self.x_input, self.w))
        self.quadratic = 1/2 * tf.reduce_sum(tf.square(tf.matmul(self.x_input, self.V)) - tf.matmul(tf.square(self.x_input), tf.square(self.V)), axis=1, keepdims=True)
        self.y_out = self.linear + self.quadratic
        self.y_pred = tf.round(tf.sigmoid(self.y_out))
        self.loss = tf.reduce_mean(tf.nn.sigmoid_cross_entropy_with_logits(logits=self.y_out, labels=self.y_input))
        self.train_op = tf.train.GradientDescentOptimizer(self._learning_rate).minimize(self.loss)
        self.accuracy = tf.reduce_mean(tf.cast(tf.equal(self.y_pred, self.y_input), tf.float32))
        init = tf.global_variables_initializer()
        self.sess.run(init)

    def train(self, X, Y, iterations=1000, batch_size=16, validation_size=0.1):
        x_train, x_test, y_train, y_test = train_test_split(X, Y, test_size=validation_size)
        for epoch in range(self._epochs):
            for i in range(iterations):
                rand_idx = np.random.choice(x_train.shape[], size=batch_size)
                rand_x = x_train[rand_idx]
                rand_y = y_train[rand_idx]
                self.sess.run(self.train_op, feed_dict={self.x_input: rand_x, self.y_input: rand_y})
                if i % 100 == 99:
                    loss = self.sess.run(self.loss, feed_dict={self.x_input: x_test, self.y_input: y_test})
                    acc = self.sess.run(self.accuracy, feed_dict={self.x_input: x_test, self.y_input: y_test})
                    print('epoch = {}, iteration ={}, loss = {}, accuracy ={}'.format(epoch, i, loss, acc))

    def predict(self, x):
        return self.sess.run(self.y_pred, feed_dict={self.x_input: x})