large{<--------收藏別忘記點贊-------->} ag{^_^}

前言

上一篇文章講了LR+FTRL演算法原理以及工程化實現。在實際的項目開發中，常常使用的是LR+組合特徵+FTRL的方式進行建模。這種方式需要人工組合特徵，非常考驗經驗，而且存在線下組合的有效特徵線上不一定有效、當前有效的特徵未來也不一定有效，所以逐漸被其它的可以自動組合特徵的模型取代。業界常用的兩種組合特徵的方式是：FM系列以及Tree系列。本文主要講解FM系列裡面的基本款：Factor Machine模型。

特徵組合

我們採用多項式模型表述組合特徵的作用。在多項式模型中，特徵 x_i? 與 x_j? 的組合用 x_ix_j? 表示。為了簡單起見，我們討論二階多項式模型。具體的模型表達式如下：

y=w_{0} + sum_{i=1}^{N}w_{i}x_{i} + sum_{i=1}^{N}sum_{j=i+1}^{N}w_{i,j}x_{i}x_{j} ag1

公式1的前兩個屬於標準的線性模型，第三個部分對應的就是特徵組合。傳統的多項式模型中，係數 w_{i,j}? 是相互獨立的，需要組合特徵 x_{i} 和 x_{j}? 出現足夠置信的情況下才能學習出有效的值，而實際的推薦、廣告系統數據中訓練數據往往非常的稀疏，因此簡單的多項式組合特徵獲得的模型往往失效，比如高階kernel的SVM模型。

那麼如何解決二次項參數的訓練問題？矩陣分解提供了一種解決思路。在model-based的協同過濾中，一個rating矩陣可以分解為user矩陣和item矩陣，每個user和item都可以用一個隱向量表示。如下圖所示，我們把每一個user表示成一個二維向量，同時把每個item表示成一個二維向量，兩個向量的乘積就是就是矩陣中user對item的打分。

因為 w_{i,j}->W? 是一個NxN的實對稱矩陣，根據實對稱矩陣的性質，可以知道W有N個線性無關的特徵向量，並且這些特徵向量都可以正則化而得到一組模為1的向量。故而實對稱矩陣可以分解為:

W=Q*Lambda*Q ag{2}

其中Q?為正交矩陣， Lambda? 為實對角矩陣。所以W可以分解為 W=V^{T}V? ，其中V的第j列( v_{j}? )對應的便是第j維特徵 x_{j}? 的隱向量。換句話說，特徵 x_{i}? 與 x_{j}? 的交叉項係數 w_{i,j}=<v_{i},v_{j}>

因此FM的組合特徵的係數更新為一下的形式：

hat{y}(x) = w_{0} + sum_{i=1}^{N}w_{i}x_{i} + sum_{i=1}^{N}sum_{j=i+1}^{N}<v_{i},v_{j}>x_{i}x_{j} ag3

其中： v_{i} = (v_{i1},v_{i2},v_{i3},v_{i4},v_{i5},....v_{ik})^T , <v_{i}, v_{j}> 代表隱向量的點積。隱向量的長度為k(k<?<n，包含k個描述特徵的因子。參數因子化使得 x_{i}x_{j}? 的參數不再是相互獨立的，比如 x_{i}x_{j}? 和 x_{i}x_{h}? 的係數分別是 <v_{i},v_{j}>? 以及 <v_{i},v_{h}>? ，他們之間有共同項 v_{i}? 。也就是說，所有包含 x_{i}? 的非零組合特徵（存在某個 j eq{i}? ，使得 x_{i,j} eq{0}? ）的樣本都可以用來學習隱向量 v_{i}? ，這很大程度上避免了數據稀疏性造成的影響。而在多項式模型中， w_{hi} 和 w_{ij} 相互獨立的。

舉個例子，如上圖所示，假設我們想預測組合特徵A與ST對應的target的值，很明顯上面並沒有A與ST同時為1的數據，如果 w_{a,st}? 是獨立的，那麼 w_{a,st}=0? 。但是對於FM而言， v_{a}, v_{st}? 可以從其它的組合（比如 <v_{a}, v_{ti}>?, <v_{b},v_{st}>? )中學習出來，因此 w_{a,st} eq{0}? 。

所以FM能很好的解決高維稀疏數據的特徵組合問題。

參數優化

公式2是一個通用的擬合方程，使用不同的損失函數可以用於解決分類、回歸、排序問題，比如使用MSE求解回歸問題，使用合頁損失函數/cross-entropy損失函數來解決分類問題。

公式2的計算複雜度，從第一眼上看是 kn^{2}? ，但是可以通過以下的方式優化到 kn? ：

egin{align*} &sum_{i=1}^{n}sum_{j=i+1}^{n}<v_{i},v_{j}>x_{i}x_{j}\ &=frac{1}{2}sum_{i=1}^{n}sum_{j=1}^{n}<v_{i},v_{j}>x_{i}x_{j}-frac{1}{2}sum_{i=1}^{n}x_{i}x_{i}\ &=frac{1}{2}(sum_{i=1}^{n}sum_{j=1}^{n}sum_{f=1}^{k}<v_{i,f},v_{j,f}>x_{i}x_{j}-sum_{i=1}^{n}sum_{f=1}^{k}v_{i,f}^2x_{i}^2)\ &=frac{1}{2}sum_{f=1}^{k}((sum_{i=1}^{n}v_{i,f}x_{i})(sum_{j=1}^{n}v_{j,f}x_{j})-(sum_{i=1}^{n}v_{i,f}^2x_{i}^2)))\ &=frac{1}{2}sum_{f=1}^{k}((sum_{i=1}^{n}v_{i,f}x_{i})^2-(sum_{i=1}^{n}v_{i,f}^2x_{i}^2)) end{align*} ag4

公式4第一步推理是根據點積的對稱性，把雙重求和想像成實對稱矩陣的上三角（不包括對角線）元素求和，那麼可以得到上三角矩陣元素的和等於整個矩陣元素求和減去矩陣的trace，然後除以2的值。

可以看到FM的組合特徵係數的計算複雜度已經降低為 kn? 。其中k一般設為2或者4，通過抑制FM的表達能力，來提高模型的泛化能力。當K固定以後，FM的計算時間複雜度直隨著特徵的數量線性增加。結合公式3和公式4得到FM模型的參數方程為:

egin{align*} hat{y}(x) &= w_{0} + sum_{i=1}^{N}w_{i}x_{i} + frac{1}{2}sum_{f=1}^{k}((sum_{i=1}^{n}v_{i,f}x_{i})^2-(sum_{i=1}^{n}v_{i,f}^2x_{i}^2))\ end{align*} ag{5}

公式5的參數空間為 heta={w_{0},w_{1},w_{2}...,v_{1,1},....v_{n,k}}? 。以二分類為例，在公式5的基礎上再外面在套一個sigmod函數即可：

egin{align*} P(Y=1|X) &=pi(x)\ &= frac{exp(hat{y(x| heta)})}{1+exp(hat{y(x| heta)})} \ end{align*} ag{6}

egin{align*} P(Y=1|X) &=1-pi(x)\ &= frac{1}{1+exp(hat{y(x| heta)})} \ end{align*} ag{7}

似然函數是：

prod_{i=1}^{n}[pi(x_i)]^{y_{i}}*[1-pi(x_{i})]^{1-y_{i}} ag{8}

對數似然函數是:

egin{align*} L( heta) =& sum_{i=1}^{N}[y_{i}*log(pi(x_{i}))+(1-y_{i})*log(1-pi(x_{i})]))\ =&sum_{i=1}^{N}[y_{i}logfrac{pi(x_{i})}{1-pi_(x_{i})}+log(1-pi(x_{i}))]\ =&sum_{i=1}^{N}[y_{i}*(hat{y(x| heta)})+log(1+exp(hat{y(x| heta)}))] ag{9} end{align*}

這樣FM參數優化問題就變成了以對數似然函數為目標的最優化問題，也就是對公式9求最大值時的參數 heta ；實際應用中，通常是使用最小化損失函數的概念來尋找最佳參數空間，所以對公式10取負號得到邏輯回歸模型的損失函數：

egin{align*} l( heta) =&-sum_{i=1}^{N}[y_{i}*(hat{y(x| heta)})-log(1+exp(hat{y(x| heta)}))]\ end{align*} ag{10}

也就是我們常說的交叉熵損失函數，是一個高階可導連續凸函數。根據凸函數優化理論，可以使用梯度下降法、牛頓法、trust-region等方法進行優化。 l( heta)? 對應的 模型參數的偏導是： frac{partial}{partial heta}=left{ egin{array}{fm} &1 & ext{if $ heta$ is $w_{0}$}\ &x_{i} & ext{if $ heta$ is $w_{i}$}\ &x_{i}sum_{j=1}^{j=n}v_{j,f}x_{j}-x_{i,f}x_{i}^2 & ext{if $ heta$ is $v_{i,f}$}\ end{array} ag{11} ight.

這裡我們使用更加高效的FTRL優化演算法。簡單回顧一下FTRL演算法工程化偽代碼：

egin{align*} hline &Large{ ext{Algorithm1 Per-Cordinate FTRL-Proximal with $L_{1}$ and}}\ &Large{ ext{$L_2$ Regularization for Logistic Regression}}\ hline &large{ ext{# With per-cordinate learning rates of Eq.(2)}}\ &largeold{Input:} ext{parameters $alpha$,$eta$,$lambda_{1}$,$lambda_{2}$}\ & ext{$forall{i}in{1,.....,d}$,initialize $z_{i}=0$ and $n_{i}=0$}\ &old{for} ext{ t=1 $old{to}$ T } old{do}\ &spacespacespacespace ext{Recieve feature vector $x_{t}$ and let $I={i|x_{i} eq{0}}$}\ &spacespacespacespace ext{For $iin{I}$ compute}\ &spacespacespacespacespacespacespacespace{w_{t,i}=left{ egin{array}{fm} &0 & ext{if $|z_{i}|$ < $lambda_{1}$}\ &-(frac{eta+sqrt{n_{i}}}{alpha}+lambda_{2})^{-1}(z_{i}-sgn(z_{i})lambda_{1}) & ext{otherwise.}\ end{array} ight.}\ &spacespacespacespace ext{Predict $P_t = sigma(x_{t}*w)$ using the $w_{t,i}$ compute above} \ &spacespacespacespace ext{Observe label $y_{t}in{0,1}$}\ &spacespacespacespace ext{$old{for}$ all $iin{I}$ do}\ &spacespacespacespacespacespacespacespace ext{$g_{i}=(p_{t}-y_{t})x_{i}$ #gradient of loss w.r.t. $w_{i}$}\ &spacespacespacespacespacespacespacespace ext{$sigma_{i}=frac{1}{alpha}(sqrt{n_{i}+g_{i}^2}-sqrt{n_{i}})$ # equal $frac{1}{eta_{t,i}}-frac{1}{eta_{t-1,i}}$}\ &spacespacespacespacespacespacespacespace ext{$z_{i}leftarrow{z_{i}}+g_{i}-sigma_{i}omega_{t,i}$}\ &spacespacespacespacespacespacespacespace ext{$n_{i}leftarrow{n_{i}}+g_{i}^2$}\ &spacespacespacespace ext{$old{endspace{for}}$}\ & ext{$old{endspace{for}}$}\ hline end{align*}

其實FTRL是一個online learning的框架，能解決的問題絕不僅僅是LR，已經成了一個通用的優化運算元，比如TensorFlow的optimizer中都包含了FTRL。我們只要把截圖中的偽代碼修改， p_{t}? 的計算改為 y(x| heta)? ，對於每一輪的特徵向量 x? 的每一維非0特徵 x_{i}? ，都要相應的更新模型參數 w_{0},w_{1},w_{2}...,v_{1,1},....v_{n,k}? ，更新公式不變和截圖一致，梯度的計算即為損失函數對每個參數的偏導，前面已經給出。 sigma,z,n? 的更新公式不變。偽代碼如下(labtex代碼有點長，知乎無法顯示，只好省略一部分偽代碼)：

egin{align*} hline &Large ext{Algorithm: FM with FTRL}\ hline &large ext{Input: Parameters $alpha^{w},alpha^{v},eta^{w},eta^{v},lambda_{1}^{w},lambda_2^{w},sigma$}\ & ext{Init:$w_{0}=0;n_{0}^{w}=0,z_{0}^{w}=0;$}\ & ext{Init:$forall{i},forall{f},w_{i}=0;n_{i}^{w}=0;z_{i}^{w}=0;v_{i,f}sim{N(0,sigma)};n_{i,f}^v=0;z_{i,f}^v=0$;}\ & ext{for t=1 to T, do}\ &spacespacespacespace ext{Recive feature vector x and let $I={i|x_{i} eq{0}}$}\ &spacespacespacespacespacespacespacespacespacespacespacespace ext{$w_{0}=left{ egin{array}{fm} &0 & ext{if $|z_{0}^w|$ < $lambda_{1}^w$}\ &-(frac{eta^w+sqrt{n_{0}^w}}{alpha^w}+lambda_{2}^w)^{-1}(z_{0}^w-sgn(z_{0}^w)lambda_{1}^w) & ext{otherwise.}\ end{array} ight.$}\ &spacespacespacespace ext{for $iin{I}$ compute:}\ &spacespacespacespacespacespacespacespacespacespacespacespace ext{$w_{i}=left{ egin{array}{fm} &0 & ext{if $|z_{i}^w|$ < $lambda_{i}^w$}\ &-(frac{eta^w+sqrt{n_{i}^w}}{alpha^w}+lambda_{2}^w)^{-1}(z_{i}^w-sgn(z_{i}^w)lambda_{1}^w) & ext{otherwise.}\ end{array} ight.$}\ &spacespacespacespacespacespacespacespace ext{for $f=1space{to}space{k}$ compute:}\ &spacespacespacespacespacespacespacespacespacespacespacespace ext{$v_{i,f}=left{ egin{array}{fm} &0 & ext{if $|z_{i,f}^v|$ < $lambda_{1}^v$}\ &-(frac{eta^v+sqrt{n_{i,f}^v}}{alpha^v}+lambda_{2}^v)^{-1}(z_{i,f}^v-sgn(z_{i,f}^v)lambda_{1}^v) & ext{otherwise.}\ end{array} ight.$}\ &spacespacespacespacespacespacespacespace ext{end for}\ &spacespacespacespace ext{end for}\ &spacespacespacespace ext{Compute $hat{y}(x| heta)$}\ &spacespacespacespace ext{Observe label $yin{-1,1}$}\ &spacespacespacespace ext{Compute $g_{0}^{w}$ #對應$ heta$為$w_{0}$}\ &spacespacespacespace ext{$ heta_{0}^{w}=frac{1}{alpha^w}(sqrt{n_{0}^{w}+(g_{0}^{w})^2}-sqrt{n_{0}^{w}})$}\ &spacespacespacespace ext{$z_{0}^{w}leftarrow{z_{0}^w+g_{0}^w-sigma_{0}^womega_{0}^w}$}\ &spacespacespacespace ext{$n_{0}^{w}leftarrow{n_{0}^w+(g_{0}^w)^2}$}\ &spacespacespacespace ext{for $iin{I}, do$}\ &spacespacespacespacespacespacespacespace ext{compute $g_{i}^{w}$}\ &spacespacespacespacespacespacespacespace ext{$sigma_{i}^{w}=frac{1}{alpha^{w}}(sqrt{n_{i}^{w}+(g_{i}^{w})^2}-sqrt{n_{i}^{w}})$}\ &spacespacespacespacespacespacespacespace ext{$z_{i}^{w}leftarrow{z_{i}^w+g_{i}^w-sigma_{i}^womega_{i}^w}$}\ &spacespacespacespacespacespacespacespace ext{$n_{i}^{w}leftarrow{n_{i}^w+(g_{i}^w)^2}$}\ &spacespacespacespacespacespacespacespace ext{for $f=1$ to k do}\ &spacespacespacespacespacespacespacespace ext{……}\ &spacespacespacespacespacespacespacespace ext{end for}\ &spacespacespacespace ext{end for}\ & ext{end for}\ hline end{align*}

工程實現

通過上面的講述，理解了FM+FTRL的基本原理，按照之前的LR+FTRL的流程，也能給出對應的FM+FTRL的代碼實現。這裡主要給出的是基於TensorFlow的FM+FTRL代碼實現，其中FTRL的部分使用的是TensorFlow的通用運算元：

# -*- coding: utf-8 -*-

"""
this is a fm_ftrl model with structured tensorflow coding style, and support online feature encoding
"""

import functools
import tensorflow as tf
import numpy as np
import os
import pandas as pd

def doublewrap(function):
"""
A decorator decorator, allowing to use the decorator to be used without
parentheses if no arguments are provided. All arguments must be optional.
"""

@functools.wraps(function)
def decorator(*args, **kwargs):
if len(args) == 1 and len(kwargs) == 0 and callable(args[0]):
return function(args[0])
else:
return lambda wrapee: function(wrapee, *args, **kwargs)

return decorator

@doublewrap
def define_scope(function, scope=None, *args, **kwargs):
"""
A decorator for functions that define TensorFlow operations. The wrapped
function will only be executed once. Subsequent calls to it will directly
return the result so that operations are added to the graph only once.
The operations added by the function live within a tf.variable_scope(). If
this decorator is used with arguments, they will be forwarded to the
variable scope. The scope name defaults to the name of the wrapped
function.
"""
attribute = _cache_ + function.__name__
name = scope or function.__name__

@property
@functools.wraps(function)
def decorator(self):
if not hasattr(self, attribute):
with tf.variable_scope(name, *args, **kwargs):
setattr(self, attribute, function(self))
return getattr(self, attribute)

return decorator

class FM_FTRL:
def __init__(self, x, y, p, k):
"""

:param x: input x
:param y: label
:param p: number of columns
:param k: dim of v for FM pair interaction vector
"""
self.x = x
self.y = y
self.p = p
self.k = k
self.predict
self.optimize
self.w0
self.W
self.V
self.norm
self.error
self.loss

@define_scope
def predict(self):
"""
this function used to predict data
:return:
"""
x = self.x
self.w0 = tf.Variable(tf.zeros([1]))
self.W = tf.Variable(tf.zeros([self.p]))
self.V = tf.Variable(tf.random_normal([self.k, self.p], stddev=0.01))
liner_terms = tf.add(self.w0,
tf.reduce_sum(tf.multiply(self.W, x), 1, keepdims=True)
)
pair_terms = tf.multiply(0.5,
tf.reduce_sum(
tf.subtract(
tf.pow(tf.matmul(x, tf.transpose(self.V)), 2),
tf.matmul(tf.pow(x, 2), tf.transpose(tf.pow(self.V, 2)))
)
))
predict = tf.add(liner_terms, pair_terms)
return predict

@define_scope
def norm(self):
"""

:return:
"""
lambda_w = tf.constant(0.001, name="lambda_w")
lambda_v = tf.constant(0.001, name="lambda_v")
l2_norm = tf.reduce_sum(
tf.add(
tf.multiply(lambda_w, tf.pow(self.W, 2)),
tf.multiply(lambda_v, tf.pow(self.V, 2))
)
)
return l2_norm

@define_scope
def error(self):
y = self.y
y_hat = self.predict
error = tf.reduce_mean(
tf.square(
tf.subtract(y, y_hat)
)
)
return error

@define_scope
def loss(self):
loss = tf.add(self.error, self.norm)
return loss

@define_scope
def optimize(self, mode="ftrl"):
if mode == ftrl:
opt = tf.train.FtrlOptimizer(learning_rate=0.1).minimize(self.loss)
else:
opt = tf.train.AdamOptimizer(learning_rate=0.001).minimize(self.loss)
return opt

def hash_java(key):
"""
hash equal to jaha hash funtion,which hash valus > 0, this is very import for engineer and ont-hot encode
:param key:
:return:
"""
h = 0
for c in key:
h = ((h * 37) + ord(c)) & 0xFFFFFFFF
return h

def main():
"""

:return:
"""
epochs = 20
batch_size = 1000

D = 3000
p = 2
k = 2

cols = [user, item, rating, timestamp]
use_cols = [user, item, rating]
features = [user, item]

data_dir = os.path.abspath(f"{os.path.abspath(os.path.dirname(os.path.realpath(__file__)))}/../../Data/fm/ml-100k")

x = tf.placeholder(float, shape=[None, D])
y = tf.placeholder(float, shape=[None, 1])
model = FM_FTRL(x=x, y=y, p=D, k=k)
sess = tf.Session()
sess.run(tf.global_variables_initializer())
num_lines = sum(1 for l in open(f{data_dir}/ua.base)) - 1
print(f"total train lines number is {num_lines}")
for epoch in range(epochs):
total_bacth = 0
avg_cost = 0.
# create random data based on random index
index_random = np.random.permutation(num_lines)

for row_index in range(0, index_random.shape[0], batch_size):
skip_rows = np.concatenate([index_random[:row_index], index_random[row_index+batch_size:]])
row = pd.read_csv(f{data_dir}/ua.base, delimiter= , names=cols,
usecols=[user, item, rating],
skiprows=skip_rows)
total_bacth += 1
bY = row[rating].values.reshape(-1, 1)
bX = np.zeros([D])
for f in features:
hash_index = hash_java(str(row[f]) + f) % D
if hash_index < 0:
raise ValueError("index for one-hot should be bigger than 0")
bX[hash_index] = 1
bX = bX.reshape(-1, D)
mse, loss_val, w, v, _ = sess.run([model.error, model.loss, model.W, model.V, model.optimize],
feed_dict={x: bX, y: bY})
avg_cost += loss_val
# Display logs per epoch step
if (epoch + 1) % 1 == 0:
print(f"total batch is {total_bacth}")
print(f"Epoch:{epoch + 1}, cost={avg_cost/total_bacth}")
print(MSE: , mse)
print(Learnt weights:, w, w.shape)
print(Learnt factors:, v, v.shape)
# print(f"auc value is {tf.summary.scalar(AUC, auc)}")

errors = []
test = pd.read_csv(f{data_dir}/ua.test, delimiter= , names=cols, usecols=[user, item, rating],
chunksize=batch_size)
for row in test:
bY = row[rating].values.reshape(-1, 1)
bX = np.zeros([D])
for f in features:
hash_index = hash_java(str(row[f]) + "_" + f) % D
bX[hash_index] = 1
bX = bX.reshape(-1, D)
errors.append(sess.run(model.error, feed_dict={x: bX, y: bY}))

RMSE = np.sqrt(np.array(errors).mean())
print(RMSE)
sess.close()

if __name__ == __main__:
main()

數據地址： MovieLens 100K Dataset

參考文獻：

  1. Factorization Machines (Steffen Rendle)
  2. Factorization Machines with libFM (STEFFEN RENDLE, University of Konstanz)
  3. 美團技術團隊
  4. Introductory Guide – Factorization Machines & their application on huge datasets (with codes in Python)
  5. 推薦系統遇上深度學習(一)--FM模型理論和實踐
  6. FM, FTRL, Softmax
  7. LR+FTRL演算法原理以及工程化實現

推薦閱讀：

相關文章