在上一篇文章當中我們剖析了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，將會在後續給大家一一介紹。
今天的文章就到這裡，衷心祝願大家每天都有所收穫。如果還喜歡今天的內容的話，請來一個三連支援吧~（點贊、關注、轉發）
