推薦系統系列（一）：FM理論與實踐

本文轉載自查看原文 2019-10-03 15:17 417 機器學習/ 推薦系統/ 深度學習/ FM/ 機器學習&深度學習

背景

在推薦領域CTR（click-through rate）預估任務中，最常用到的baseline模型就是LR（Logistic Regression）。對數據進行特征工程，構造出大量單特征，編碼之后送入模型。這種線性模型的優勢在於，運算速度快可解釋性強，在特征挖掘完備且訓練數據充分的前提下能夠達到一定精度。但這種模型的缺點也是較為明顯的：

模型並未考慮到特征之間的關系 $y=w_0+\sum_{i=1}^{n}w_ix_i$ 。在實踐經驗中，對特征進行交叉組合往往能夠更好地提升模型效果。
對於多取值的categorical特征進行one-hot編碼，具有高度稀疏性，帶來維度災難問題。

FM（Factorization Machine）模型就是針對在特征組合過程中遇到的上述問題而提出的一種高效的解決方案[1]。由於FM優越的性能表現，后續出現了一系列FM變種模型，從淺層模型到深度推薦模型中都有FM的影子。

分析

1. FM定義

FM以特征組合進行切入點，在公式定義中引入特征交叉項，彌補了一般線性模型未考慮特征間關系的缺憾。公式如下（FM模型可拓展到高階，但為簡化且不失一般性，這里只討論二階交叉）[1]：

\[y = w_0+\sum_{i=1}^nw_ix_i+\sum_{i=1}^{n-1}\sum_{j=i+1}^nw_{ij}x_ix_j \tag{1} \]

與一般線性模型相比，公式（1）僅多了一個二階交叉項，模型參數多了 $\frac{n(n+1)}{2}$ 個。雖然這種顯式交叉的方式能夠刻畫特征間關系，但是對公式求解帶來困難。

因為大量特征進行one-hot表示之后具有高度稀疏性的問題，所以公式（1）中的 $x_ix_j$ 同樣會產生大量的0值。參數學習不充分，直接導致$w_{ij}$ 無法通過訓練得到。（解釋：令$x_ix_j=X$，則$\frac{\partial{y}}{\partial{w_{ij}}}=X$ ，又因 $X=0$，所以$w_{ij}^{new}=w_{ij}^{old}+{\alpha}X=w_{ij}^{old}$ ，梯度為0參數無法更新。）

導致這種情況出現的根源在於：特征過於稀疏。我們期望的是找到一種方法，使得 $w_{ij}$ 的求解不受特征稀疏性的影響。

2. 公式改寫

為了克服上述困難，需要對FM公式進行改寫，使得求解更加順利。受 矩陣分解 的啟發，對於每一個特征 $x_i$ 引入輔助向量（隱向量）$V_i=(v_{i1},v_{i2},\cdots,v_{ik})$，然后利用$V_iV_j^T$ 對$w_{ij}$ 進行求解。即，做如下假設： $w_{ij} \approx V_iV_j^T$ 。

引入隱向量的好處是：

二階項的參數量由原來的 $\frac{n(n-1)}{2}$ 降為 $kn$ 。
原先參數之間並無關聯關系，但是現在通過隱向量可以建立關系。如，之前 $w_{ij}$ 與 $w_{ik}$ 無關，但是現在 $w_{ij}=\langle V_i,V_j\rangle ,w_{ik}=\langle V_i,V_k\rangle $ 兩者有共同的 $V_i$ ，也就是說，所有包含 $x_ix_j$ 的非零組合特征（存在某個 $j\neq i$ ，使得 $x_ix_j\neq 0$ ）的樣本都可以用來學習隱向量 $V_i$ ，這很大程度上避免了數據稀疏性造成的影響。[2]

現在可以將公式（1）進行改寫：

\[y = w_0+\sum_{i=1}^nw_ix_i+\sum_{i=1}^{n-1}\sum_{j=i+1}^n\langle V_i,V_j\rangle x_ix_j \tag{2} \]

重心轉移到如何求解公式（2）后面的二階項。

預備知識：

首先了解對稱矩陣上三角求和，設矩陣為 $M$：

\[M=\left( \begin{matrix} m_{11} & m_{12} & \cdots & m_{1n} \\ m_{21} & m_{22} & \cdots & m_{1n} \\ \vdots & \vdots & \ddots & \vdots \\ m_{n1} & m_{n2} & \cdots & m_{nn} \\ \end{matrix} \right)_{n*n} \]

其中，$m_{ij}=m_{ji}$ 。

令上三角元素和為 $A$ ，即 $\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}m_{ij}=A$ 。那么，$M$ 的所有元素之和等於 $2*A+tr(M)$ ，$tr(M)$為矩陣的跡。

\[\sum_{i=1}^n\sum_{j=1}^nm_{ij}=2*\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}m_{ij} + \sum_{i=1}^{n}m_{ii} \]

可得，

\[A=\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}m_{ij}=\frac{1}{2}*\left\{\sum_{i=1}^n\sum_{j=1}^nm_{ij}-\sum_{i=1}^{n}m_{ii}\right\} \]

正式改寫：

有了上述預備知識，可以對公式（2）的二階項進行推導：

\[\begin{align} & \sum_{i=1}^{n-1}\sum_{j=i+1}^n\langle V_i,V_j\rangle x_ix_j \notag \\ ={} & \frac{1}{2}*\left\{\sum_{i=1}^{n}\sum_{j=1}^{n}\langle V_i,V_j\rangle x_ix_j-\sum_{i=1}^{n}\langle V_i,V_i\rangle x_ix_i\right\} \notag \\ ={} & \frac{1}{2}*\left\{\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{f=1}^{k}v_{if}v_{jf}x_ix_j-\sum_{i=1}^{n}\sum_{f=1}^{k}v_{if}v_{if}x_ix_i\right\} \notag \\ ={} & \frac{1}{2}*\sum_{f=1}^{k}\left\{\sum_{i=1}^{n}\sum_{j=1}^{n}v_{if}x_iv_{jf}x_j-\sum_{i=1}^{n}v_{if}^{2}x_{i}^2\right\} \notag \\ ={} & \frac{1}{2}*\sum_{f=1}^{k}\left\{\left(\sum_{i=1}^{n}v_{if}x_i\right)\left(\sum_{j=1}^{n}v_{jf}x_j\right)-\sum_{i=1}^{n}v_{if}^{2}x_{i}^2\right\} \notag \\ ={} & \frac{1}{2}*\sum_{f=1}^{k}\left\{\left(\sum_{i=1}^{n}v_{if}x_i\right)^{2}-\sum_{i=1}^{n}v_{if}^{2}x_{i}^2\right\} \notag \\ \end{align}\tag{3} \]

結合（2）（3），可以得到：

\[\begin{align} y ={} & w_0+\sum_{i=1}^nw_ix_i+\sum_{i=1}^{n-1}\sum_{j=i+1}^n\langle V_i,V_j\rangle x_ix_j \notag \\ ={} & w_0+\sum_{i=1}^nw_ix_i+\frac{1}{2}*\sum_{f=1}^{k}\left\{\left(\sum_{i=1}^{n}v_{if}x_i\right)^{2}-\sum_{i=1}^{n}v_{if}^{2}x_{i}^2\right\} \notag \\ \end{align} \tag{4} \]

至此，我們得到了想要的模型表達式。

為什么要將公式（2）改寫為公式（4），是因為在改寫之前，計算 $y$ 的復雜度為 $O(kn^2)$ ，改寫后的計算復雜度為 $O(kn)$ ，提高模型推斷速度。

3. FM求解

到目前為止已經得到了FM的模型表示（4），如何對模型參數求解呢？可以使用常見的梯度下降法對參數進行求解，為了對參數進行梯度下降更新，需要計算模型各參數的梯度表達式：

當參數為 $w_0$ 時，$\frac{\partial{y}}{\partial{w_0}}=1$ 。

當參數為 $w_i$ 時，$\frac{\partial{y}}{\partial{w_i}}=x_i$ 。

當參數為 $v_{if}$ 時，只需要關注模型高階項，當計算參數 $v_{if}$ 的梯度時，其余無關參數可看做常數。

\[\begin{align} \frac{\partial{y}}{\partial{v_{if}}} ={} & \partial{\frac{1}{2}\left\{\left(\sum_{i=1}^{n}v_{if}x_i\right)^{2}-\sum_{i=1}^{n}v_{if}^{2}x_{i}^2\right\}}/\partial{v_{if}} \notag \\ ={} & \frac{1}{2}* \left\{ \frac{ \partial{ \left\{ \sum_{i=1}^{n}v_{if}x_i \right\}^2 } }{\partial{v_{if}}} - \frac{ \partial{ \left\{ \sum_{i=1}^{n}v_{if}^{2}x_{i}^2 \right\} } }{\partial{v_{if}}} \right\} \notag \\ \end{align} \tag{5} \]

其中：

\[\frac{ \partial{ \left\{ \sum_{i=1}^{n}v_{if}^{2}x_{i}^2 \right\} } }{\partial{v_{if}}} = 2x_{i}^2v_{if} \tag{6} \]

令 $\lambda=\sum_{i=1}^{n}v_{if}x_i$ ，則：

\[\begin{align} \frac{ \partial{ \left\{ \sum_{i=1}^{n}v_{if}x_i \right\}^2 } }{\partial{v_{if}}} ={} & \frac{\partial{\lambda^2}}{\partial{v_{if}}} \notag \\ ={} & \frac{\partial{\lambda^2}}{\partial{\lambda}} \frac{\partial{\lambda}}{\partial{v_{if}}} \notag \\ ={} & 2\lambda*\frac{\partial{\sum_{i=1}^{n}v_{if}x_i}}{\partial{v_{if}}} \notag \\ ={} & 2\lambda*x_i \notag \\ ={} & 2*x_i*\sum_{j=1}^{n}v_{jf}x_j \notag \\ \end{align} \tag{7} \]

結合公式（5~7），可得：

\[\frac{\partial{y}}{\partial{v_{if}}} = x_i\sum_{j=1}^{n}v_{jf}x_j-x_{i}^2v_{if} \tag{8} \]

綜上，最終模型各參數的梯度表達式如下：

\[\begin{equation} \frac{\partial{y}}{\partial{\theta}} = \begin{cases} 1, & \text{if } \theta \text{ is } w_0; \\ x_i, & \text{if } \theta \text{ is } w_i; \\ x_i\sum_{j=1}^{n}v_{jf}x_j-x_{i}^2v_{if}, & \text{if } \theta \text{ is } v_{if}. \end{cases} \end{equation} \notag \]

4. 性能分析

由第2小節可知，FM進行推斷的時間復雜度為 $O(kn)$ 。分析訓練的復雜度，依據參數的梯度表達式，$\sum_{j=1}^{n}v_{jf}x_{j}$ 與 $i$ 無關，在參數更新時可以首先將所有的 $\sum_{j=1}^{n}v_{jf}x_{j}$ 計算出來，復雜度為 $O(kn)$ ，后續更新所有參數的時間復雜度均為 $O(1)$ ，參數量為 $1+n+kn$ ，所以最終訓練的時間復雜度同樣為 $O(kn)$ ，其中 $n$ 為特征數，$k$ 為隱向量維數。

FM訓練與預測的時間復雜度均為 $O(kn)$ ，是一種十分高效的模型。

5. 優缺點

優點 [1]：

In total, the advantages of our proposed FM are:

FMs allow parameter estimation under very sparse data where SVMs fail.

FMs have linear complexity, can be optimized in the primal and do not rely on support vectors like SVMs. We show that FMs scale to large datasets like Netflix with 100 millions of training instances.

FMs are a general predictor that can work with any real valued feature vector. In contrast to this, other state-of- the-art factorization models work only on very restricted input data. We will show that just by defining the feature vectors of the input data, FMs can mimic state-of-the-art models like biased MF, SVD++, PITF or FPMC.

缺點：

每個特征只引入了一個隱向量，不同類型特征之間交叉沒有區分性。FFM模型正是以這一點作為切入進行改進。

實驗

FM既可以應用在回歸任務，也可以應用在分類任務中。如，在二分類任務中只需在公式（2）最外層套上 $sigmoid$ 函數即可，上述解析都是基於回歸任務來進行推導的。

關於模型最終的損失函數同樣可以有多種形式，如回歸任務可以使用 $MSE$ ，分類任務可以使用 $Cross Entropy$ 等。

1. 代碼演示

雖然知道可以通過引入輔助向量進行計算，但是輔助向量是如何與特征 $x_i$ 建立聯系的，換句話說，如何通過 $x_i$ 得到輔助向量 $V_i$ ？在使用神經網絡實現FM的過程中，將 $x_i$ 的 $embedding$ 作為輔助向量，最終得到的 $embedding$ 向量組也可以看作是對應特征的低維稠密表征，可以應用到其他下游任務中。

1.1 回歸任務

本文使用了 $MovieLens 100K Dataset$ [3] 作為實驗輸入，特征組分別為用戶編號、電影編號，用戶對電影的歷史評分作為 $Label$ 。

具體代碼實現如下：

# -*- coding:utf-8 -*-
import pandas as pd
import numpy as np
from scipy.sparse import csr
from itertools import count
from collections import defaultdict
import tensorflow as tf


def vectorize_dic(dic, label2index=None, hold_num=None):
  
    if label2index == None:
        d = count(0)
        label2index = defaultdict(lambda: next(d))  # 數值映射表

    sample_num = len(list(dic.values())[0])  # 樣本數
    feat_num = len(list(dic.keys()))  # 特征數
    total_value_num = sample_num * feat_num

    col_ix = np.empty(total_value_num, dtype=int)

    i = 0
    for k, lis in dic.items():
        col_ix[i::feat_num] = [label2index[str(k) + str(el)] for el in lis]
        i += 1

    row_ix = np.repeat(np.arange(sample_num), feat_num)
    data = np.ones(total_value_num)

    if hold_num is None:
        hold_num = len(label2index)

    left_data_index = np.where(col_ix < hold_num)  # 為了剔除不在train set中出現的test set數據

    return csr.csr_matrix(
        (data[left_data_index], (row_ix[left_data_index], col_ix[left_data_index])),
        shape=(sample_num, hold_num)), label2index

def batcher(X_, y_, batch_size=-1):

    assert X_.shape[0] == len(y_)

    n_samples = X_.shape[0]
    if batch_size == -1:
        batch_size = n_samples
    if batch_size < 1:
        raise ValueError('Parameter batch_size={} is unsupported'.format(batch_size))

    for i in range(0, n_samples, batch_size):
        upper_bound = min(i + batch_size, n_samples)
        ret_x = X_[i:upper_bound]
        ret_y = y_[i:upper_bound]
        yield(ret_x, ret_y)

def load_dataset():
    cols = ['user', 'item', 'rating', 'timestamp']
    train = pd.read_csv('data/ua.base', delimiter='\t', names=cols)
    test = pd.read_csv('data/ua.test', delimiter='\t', names=cols)

    x_train, label2index = vectorize_dic({'users': train.user.values, 'items': train.item.values})
    x_test, label2index = vectorize_dic({'users': test.user.values, 'items': test.item.values}, label2index, x_train.shape[1])

    y_train = train.rating.values
    y_test = test.rating.values

    x_train = x_train.todense()
    x_test = x_test.todense()

    return x_train, x_test, y_train, y_test

x_train, x_test, y_train, y_test = load_dataset()

print("x_train shape: ", x_train.shape)
print("x_test shape: ", x_test.shape)
print("y_train shape: ", y_train.shape)
print("y_test shape: ", y_test.shape)

vec_dim = 10
batch_size = 1000
epochs = 10
learning_rate = 0.001
sample_num, feat_num = x_train.shape

x = tf.placeholder(tf.float32, shape=[None, feat_num], name="input_x")
y = tf.placeholder(tf.float32, shape=[None,1], name="ground_truth")

w0 = tf.get_variable(name="bias", shape=(1), dtype=tf.float32)
W = tf.get_variable(name="linear_w", shape=(feat_num), dtype=tf.float32)
V = tf.get_variable(name="interaction_w", shape=(feat_num, vec_dim), dtype=tf.float32)

linear_part = w0 + tf.reduce_sum(tf.multiply(x, W), axis=1, keep_dims=True)
interaction_part = 0.5 * tf.reduce_sum(tf.square(tf.matmul(x, V)) - tf.matmul(tf.square(x), tf.square(V)), axis=1, keep_dims=True)
y_hat = linear_part + interaction_part
loss = tf.reduce_mean(tf.square(y - y_hat))
train_op = tf.train.AdamOptimizer(learning_rate).minimize(loss)

with tf.Session() as sess:
    sess.run(tf.global_variables_initializer())
    for e in range(epochs):
        step = 0
        print("epoch:{}".format(e))
        for batch_x, batch_y in batcher(x_train, y_train, batch_size):
            sess.run(train_op, feed_dict={x:batch_x, y:batch_y.reshape(-1, 1)})
            step += 1
            if step % 10 == 0:
                for val_x, val_y in batcher(x_test, y_test):
                    train_loss = sess.run(loss, feed_dict={x:batch_x, y:batch_y.reshape(-1, 1)})
                    val_loss = sess.run(loss, feed_dict={x:val_x, y:val_y.reshape(-1, 1)})
                    print("batch train_mse={}, val_mse={}".format(train_loss, val_loss))

    for val_x, val_y in batcher(x_test, y_test):
        val_loss = sess.run(loss, feed_dict={x: val_x, y: val_y.reshape(-1, 1)})
        print("test set rmse = {}".format(np.sqrt(val_loss)))

實驗結果：

epoch:0
batch train_mse=19.54930305480957, val_mse=19.687997817993164
batch train_mse=16.957233428955078, val_mse=19.531404495239258
batch train_mse=18.544944763183594, val_mse=19.376962661743164
batch train_mse=18.870519638061523, val_mse=19.222412109375
batch train_mse=18.769777297973633, val_mse=19.070764541625977
batch train_mse=19.383392333984375, val_mse=18.915040969848633
batch train_mse=17.26403045654297, val_mse=18.75937843322754
batch train_mse=17.652183532714844, val_mse=18.6033935546875
batch train_mse=18.331804275512695, val_mse=18.447608947753906
......
epoch:9
batch train_mse=1.394300103187561, val_mse=1.4516444206237793
batch train_mse=1.2031371593475342, val_mse=1.4285767078399658
batch train_mse=1.1761484146118164, val_mse=1.4077649116516113
batch train_mse=1.134848952293396, val_mse=1.3872103691101074
batch train_mse=1.2191411256790161, val_mse=1.3692644834518433
batch train_mse=1.572729468345642, val_mse=1.3509554862976074
batch train_mse=1.3323310613632202, val_mse=1.3339732885360718
batch train_mse=1.1601723432540894, val_mse=1.3183823823928833
batch train_mse=1.2751621007919312, val_mse=1.3023829460144043
test set rmse = 1.1405380964279175

1.2 分類任務

使用更全的 $MovieLens 100K Dataset$ 特征，將評分大於3分的樣本作為正類，其他為負類，構造二分類任務。核心代碼如下:

class FM(object):
    def __init__(self, vec_dim, feat_num, lr, lamda):
        self.vec_dim = vec_dim
        self.feat_num = feat_num
        self.lr = lr
        self.lamda = lamda

        self._build_graph()

    def _build_graph(self):
        self.add_input()
        self.inference()

    def add_input(self):
        self.x = tf.placeholder(tf.float32, shape=[None, self.feat_num], name='input_x')
        self.y = tf.placeholder(tf.float32, shape=[None], name='input_y')

    def inference(self):
        with tf.variable_scope('linear_part'):
            w0 = tf.get_variable(name='bias', shape=[1], dtype=tf.float32)
            self.W = tf.get_variable(name='linear_w', shape=[self.feat_num], dtype=tf.float32)
            self.linear_part = w0 + tf.reduce_sum(tf.multiply(self.x, self.W), axis=1)
        with tf.variable_scope('interaction_part'):
            self.V = tf.get_variable(name='interaction_w', shape=[self.feat_num, self.vec_dim], dtype=tf.float32)
            self.interaction_part = 0.5 * tf.reduce_sum(
                tf.square(tf.matmul(self.x, self.V)) - tf.matmul(tf.square(self.x), tf.square(self.V)),
                axis=1
            )
        self.y_logits = self.linear_part + self.interaction_part
        self.y_hat = tf.nn.sigmoid(self.y_logits)
        self.pred_label = tf.cast(self.y_hat > 0.5, tf.int32)
        self.loss = -tf.reduce_mean(self.y*tf.log(self.y_hat+1e-8) + (1-self.y)*tf.log(1-self.y_hat+1e-8))
        self.reg_loss = self.lamda*(tf.reduce_mean(tf.nn.l2_loss(self.W)) + tf.reduce_mean(tf.nn.l2_loss(self.V)))
        self.total_loss = self.loss + self.reg_loss

        self.train_op = tf.train.AdamOptimizer(self.lr).minimize(self.total_loss)

實驗結果：

Iter:  59400, Train acc: 0.7812, Val acc: 0.6867, Val auc: 0.7285, Val loss: 0.614005, Flag: 
Iter:  59600, Train acc: 0.8125, Val acc:  0.684, Val auc: 0.7294, Val loss: 0.615628, Flag: *
Iter:  59800, Train acc:  0.875, Val acc: 0.6665, Val auc: 0.7282, Val loss: 0.625017, Flag: 
Iter:  60000, Train acc: 0.9375, Val acc: 0.6767, Val auc: 0.7282, Val loss: 0.617686, Flag: 
Iter:  60200, Train acc:   0.75, Val acc: 0.6815, Val auc: 0.7277, Val loss: 0.614763, Flag: 
Iter:  60400, Train acc: 0.9062, Val acc:  0.681, Val auc: 0.7283, Val loss: 0.614414, Flag: 
Iter:  60600, Train acc: 0.6875, Val acc: 0.6853, Val auc: 0.7291, Val loss: 0.621548, Flag: 
Iter:  60800, Train acc:  0.625, Val acc:  0.679, Val auc: 0.7288, Val loss: 0.617327, Flag: 
Iter:  61000, Train acc: 0.7812, Val acc: 0.6835, Val auc: 0.7293, Val loss: 0.616952, Flag: 
Iter:  61200, Train acc: 0.8125, Val acc:  0.686, Val auc: 0.7292, Val loss: 0.614379, Flag: 
Iter:  61400, Train acc: 0.6562, Val acc:  0.688, Val auc: 0.7284, Val loss: 0.613859, Flag: 
Iter:  61600, Train acc: 0.6875, Val acc: 0.6725, Val auc: 0.7279, Val loss: 0.618824, Flag: 
No optimization for a long time, auto-stopping...
====== let's test =====
Test acc: 0.6833, Test auc: 0.7369

2. 注意事項

雖然FM可以應用於任意數值類型的數據上，但是需要注意對輸入特征數值進行預處理。優先進行特征歸一化，其次再進行樣本歸一化。[2]
FM不僅可以用於rank階段，同時可以用於向量召回：好文推薦

reference

[1] Rendle, S. (2010, December). Factorization machines. In 2010 IEEE International Conference on Data Mining (pp. 995-1000). IEEE.
[2] https://tech.meituan.com/2016/03/03/deep-understanding-of-ffm-principles-and-practices.html
[3] https://grouplens.org/datasets/movielens/100k/
[4] https://www.jianshu.com/p/152ae633fb00
[5] http://nowave.it/factorization-machines-with-tensorflow.html
[6] https://blog.csdn.net/xxiaobaib/article/details/92801244
[7] [https://github.com/wyl6/Recommender-Systems-Samples/blob/master/RecSys Traditional/MF/FM/fm_tensorflow.py](https://github.com/wyl6/Recommender-Systems-Samples/blob/master/RecSys Traditional/MF/FM/fm_tensorflow.py)
[8] https://zhuanlan.zhihu.com/p/58160982

知識分享

個人知乎專欄：https://zhuanlan.zhihu.com/c_1164954275573858304

歡迎關注微信公眾號：SOTA Lab
專注知識分享，不定期更新計算機、金融類文章

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。