在上一篇文章當中我們剖析了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 = 0
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[0], 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})

註：TensorFlow代碼使用的1.0版本

FM的paper到這裡我們就剖析完了，也給出了代碼實現。但是這並沒有結束，後續關於FM又有了好幾個變種的更新版本。像是AFM、FFM、PFM等等。關於這些paper，將會在後續給大家一一介紹。

今天的文章就到這裡，衷心祝願大家每天都有所收穫。如果還喜歡今天的內容的話，請來一個三連支持吧~（點贊、關注、轉發）

原文鏈接，求個關注