推荐算法之FM（Factorization Machines）算法

推荐系统算法大致发展流程如下，主要包括线性模型、树模型和神经网络模型。随着 DNN 在图像、语音、NLP 等领域取得突破，人们渐渐意识到 DNN 在特征表示上的天然优势。相继提出了使用 CNN 或 RNN 来做 CTR 预估的模型。但是，CNN 模型的缺点是：偏向于学习相邻特征的组合特征。 RNN 模型的缺点是：比较适用于有序列 (时序) 关系的数据。

特征组合

普通的线性模型，都是将各个特征独立考虑的，并没有考虑到特征与特征之间的相互关系。但实际上，大量的特征之间是有关联的。最简单的以电商为例，一般女性用户看化妆品服装之类的广告比较多，而男性更青睐各种球类装备。那很明显，女性这个特征与化妆品类服装类商品有很大的关联性，男性这个特征与球类装备的关联性更为密切。如果我们能将这些有关联的特征找出来，显然是很有意义的。

一般的线性模型为：

从上面的式子很容易看出，一般的线性模型压根没有考虑特征间的关联。为了表述特征间的相关性，我们采用多项式模型。在多项式模型中，特征

与

的组合用

表示。为了简单起见，我们只考虑二阶交叉的情况，具体的模型如下：

式中，n表示样本的特征数量,

表示第i个特征，与线性模型相比，FM的模型就多了后面特征组合的部分。

FM求解

从公式可以看出，组合特征的参数一共有 n(n−1)/2个，任意两个参数都是独立的。然而，在数据稀疏性普遍存在的实际应用场景中，二次项参数的训练是很困难的。其原因是，每个参数

的训练需要大量

和

都非零的样本；由于样本数据本来就比较稀疏，满足

和

都非零”的样本将会非常少。训练样本的不足，很容易导致参数

不准确，最终将严重影响模型的性能。

如何解决二次项参数的训练问题呢？矩阵分解提供了一种解决思路。在model-based的协同过滤中，一个rating矩阵可以分解为user矩阵和item矩阵，每个user和item都可以采用一个隐向量表示。比如在下图中的例子中，我们把每个user表示成一个二维向量，同时把每个item表示成一个二维向量，两个向量的点积就是矩阵中user对item的打分。

隐向量的长度为 k（k<<n），包含 k 个描述特征的因子，二次项的参数数量减少为 kn个，远少于多项式模型的参数数量。

公式2是一个通用的拟合方程，可以采用不同的损失函数用于解决回归、二元分类等问题，比如可以采用MSE（Mean Square Error）损失函数来求解回归问题，也可以采用Hinge/Cross-Entropy损失来求解分类问题。当然，在进行二元分类时，FM的输出需要经过sigmoid变换，这与Logistic回归是一样的。

直观上看，FM的复杂度是 O(kn^2^)。但是，通过公式(3)的等式，FM的二次项可以化简，其复杂度可以优化到 O(kn)。由此可见，FM可以在线性时间对新样本作出预测。

看一下FM的训练复杂度，利用SGD（Stochastic Gradient Descent）训练模型。模型各个参数的梯度如下：

因此，FM参数训练的复杂度也是O(kn)。综上可知，FM可以在线性时间训练和预测，是一种非常高效的模型。

FM思想及实现

FM的全称是Factorization Machines，就是因子分解机的意思，为什么叫因子分解呢，就是因为他对传统的线性回归模型加了一个因子交叉项，你可以理解为把每一个特征和其他特征相乘后求和。

其实就是把W分解了，分解成两个矩阵相乘，这样的话其实就剩kn个参数需要我们计算了，降低了算法复杂度。可以发现参数矩阵w是一个*实对称矩阵，可以使用矩阵分解的方法分解，通过引入辅助向量又称为隐向量V来表示。


可以使用常见的梯度下降法对参数进行求解，为了对参数进行梯度下降更新，需要计算模型各参数的梯度表达式。

案例演示：用python实现FM算法，数据场景为二分类问题

'''
FM(因子分解机)模型算法：稀疏数据下的特征二阶组合问题（个性化特征）
1、应用矩阵分解思想，引入隐向量构造FM模型方程
2、目标函数（损失函数复合FM模型方程）的最优问题：链式求偏导
3、SGD优化目标函数
'''
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import accuracy_score
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
'''二分类输出非线性映射'''
def sigmoid(x):
    return 1 / (1 + np.exp(-x))
'''计算logit损失函数：L = log(1 + e**(y_hat * y))'''
def logit(y, y_hat):
    return np.log(1 + np.exp(-y * y_hat))
'''计算logit损失函数的外层偏导(不含y_hat的一阶偏导)'''
def df_logit(y, y_hat):
    return sigmoid(-y * y_hat) * (-y)
'''FM的模型方程：LR线性组合 + 交叉项组合 = 1阶特征组合 + 2阶特征组合'''
def FM(Xi, w0, W, V):
    # 样本Xi的特征分量xi和xj的2阶交叉项组合系数wij = xi和xj对应的隐向量Vi和Vj的内积
    # 向量形式：Wij:= <Vi, Vj> * Xi * Xj
    interaction = np.sum((Xi.dot(V)) ** 2 - (Xi ** 2).dot(V ** 2))  # 二值硬核匹配->向量软匹配
    y_hat = w0 + Xi.dot(W) + interaction / 2  # FM预测函数
    return y_hat[0]
'''SGD更新FM模型的参数列表：[w0, W, V]'''
def FM_SGD(X, y, k=2, alpha=0.01, iter=50):
    m, n = np.shape(X)
    w0, W = 0, np.zeros((n, 1))  # 初始化wo=R、W=(n, 1)
    V = np.random.normal(loc=0, scale=1, size=(n, k))  # 初始化隐向量矩阵V=(n, k)~N(0, 1)，其中Vj是第j维特征的隐向量
    all_FM_params = []  # FM模型的参数列表：[w0, W, V]
    for it in range(iter):
        total_loss = 0  # 当前迭代模型的损失值
        for i in range(m):  # 遍历训练集
            y_hat = FM(Xi=X[i], w0=w0, W=W, V=V)  # FM的模型方程
            total_loss += logit(y=y[i], y_hat=y_hat)  # 计算logit损失函数值
            dloss = df_logit(y=y[i], y_hat=y_hat)  # 计算logit损失函数的外层偏导
            dloss_w0 = dloss * 1  # l(y, y_hat)中y_hat展开w0，求关于w0的内层偏导
            w0 = w0 - alpha * dloss_w0  # 梯度下降更新w0
            for j in range(n):  # 遍历n维向量X[i]
                if X[i, j] != 0:
                    dloss_Wj = dloss * X[i, j]  # l(y, y_hat)中y_hat展开y_hat，求关于W[j]的内层偏导
                    W[j] = W[j] - alpha * dloss_Wj  # 梯度下降更新W[j]
                    for f in range(k):  # 遍历k维隐向量Vj
                        # l(y, y_hat)中y_hat展开V[j, f]，求关于V[j, f]的内层偏导
                        dloss_Vjf = dloss * (X[i, j] * (X[i].dot(V[:, f])) - V[j, f] * X[i, j] ** 2)
                        V[j, f] = V[j, f] - alpha * dloss_Vjf  # 梯度下降更新V[j, f]
        print('FM第{}次迭代，当前损失值为：{:.4f}'.format(it + 1, total_loss / m))
        all_FM_params.append([w0, W, V])  # 保存当前迭代下FM的参数列表:[w0, W, V]
    return all_FM_params
'''FM模型预测测试集分类结果'''
def FM_predict(X, w0, W, V):
    predicts, threshold = [], 0.5  # sigmoid阈值设置
    for i in range(X.shape[0]):  # 遍历测试集
        y_hat = FM(Xi=X[i], w0=w0, W=W, V=V)  # FM的模型方程
        predicts.append(-1 if sigmoid(y_hat) < threshold else 1)  # 分类结果非线性映射
    return np.array(predicts)
'''FM在不同迭代次数下的参数列表中，训练集的损失值和测试集的准确率变化'''
def draw_research(all_FM_params, X_train, y_train, X_test, y_test):
    all_total_loss, all_total_accuracy = [], []
    for w0, W, V in all_FM_params:
        total_loss = 0
        for i in range(X_train.shape[0]):
            total_loss += logit(y=y_train[i], y_hat=FM(Xi=X_train[i], w0=w0, W=W, V=V))
        all_total_loss.append(total_loss / X_train.shape[0])
        all_total_accuracy.append(accuracy_score(y_test, FM_predict(X=X_test, w0=w0, W=W, V=V)))
    plt.plot(np.arange(len(all_FM_params)), all_total_loss, color='#FF4040', label='训练集的损失值')
    plt.plot(np.arange(len(all_FM_params)), all_total_accuracy, color='#4876FF', label='测试集的准确率')
    plt.xlabel('SGD迭代次数')
    plt.title('FM模型:二阶互异特征组合')
    plt.legend()
    plt.show()
if __name__ == '__main__':
    np.random.seed(123)
    df = pd.read_csv(r'D:\\FM-master\\data\\xg.csv')
    df['Class'] = df['Class'].map({0: -1, 1: 1})  # 标签列从[0, 1]离散到[-1, 1]
    X_train, X_test, y_train, y_test = train_test_split(df.iloc[:, :-1].values, df.iloc[:, -1].values, test_size=0.3, random_state=123)
    X_train = MinMaxScaler().fit_transform(X_train)  # 归一化训练集，返回[0, 1]区间
    X_test = MinMaxScaler().fit_transform(X_test)  # 归一化测试集，返回[0, 1]区间
    '''*****************FM预测模型*****************'''
    all_FM_params = FM_SGD(X=X_train, y=y_train, k=2, alpha=0.01, iter=45)  # SGD更新FM模型的参数列表：[w0, W, V]
    w0, W, V = all_FM_params[-1]  # FM模型的参数列表
    predicts = FM_predict(X=X_test, w0=w0, W=W, V=V)  # FM模型预测测试集分类结果 80.52%  80.09%
    print('FM在测试集的分类准确率为: {:.2%}'.format(accuracy_score(y_test, predicts)))
    # draw_research(all_FM_params=all_FM_params, X_train=X_train, y_train=y_train, X_test=X_test, y_test=y_test)

import pandas as pd
import numpy as np

pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', None)

rnames = ['user_id', 'movie_id', 'rating', 'timestamp']
df = pd.read_csv('D:\\u.data', sep='\t', header=None, names=rnames, engine='python')
#构造2分类数据集

df['rating']=df['rating'].map(lambda x: -1 if x>=3 else 1) #1,2是label=1  3，4,5是label=0

#one-hot encoder
from sklearn.preprocessing import OneHotEncoder
columns=['user_id', 'movie_id']

for i in columns:
    get_dummy_feature=pd.get_dummies(df[i])
    df=pd.concat([df, get_dummy_feature],axis=1)
    df=df.drop(i, axis=1)

df=df.drop(['timestamp'], axis=1)
#这些特征可以进一步挖掘。这里都不要了，只保留one-hot特征

from sklearn.model_selection import train_test_split

X=df.drop('rating', axis=1)
Y=df['rating']

X_train,X_val,Y_train,Y_val=train_test_split(X, Y, test_size=0.3, random_state=123)

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

def logit(y, y_hat): #对每一个样本计算损失
    if y_hat == 'nan':
        return 0
    else:
        return np.log(1 + np.exp(-y * y_hat))

def df_logit(y, y_hat):
    return sigmoid(-y * y_hat) * (-y)

from sklearn.base import BaseEstimator, ClassifierMixin
from collections import Counter


class FactorizationMachine(BaseEstimator):
    def __init__(self, k=5, learning_rate=0.01, iternum=2):
        self.w0 = None
        self.W = None
        self.V = None
        self.k = k
        self.alpha = learning_rate
        self.iternum = iternum

    def _FM(self, Xi):
        interaction = np.sum((Xi.dot(self.V)) ** 2 - (Xi ** 2).dot(self.V ** 2))
        y_hat = self.w0 + Xi.dot(self.W) + interaction / 2
        return y_hat[0]

    def _FM_SGD(self, X, y):
        m, n = np.shape(X)
        # 初始化参数
        self.w0 = 0
        self.W = np.random.uniform(size=(n, 1))
        self.V = np.random.uniform(size=(n, self.k))  # Vj是第j个特征的隐向量  Vjf是第j个特征的隐向量表示中的第f维

        for it in range(self.iternum):
            total_loss = 0
            for i in range(m):  # 遍历训练集
                y_hat = self._FM(Xi=X[i])  # X[i]是第i个样本  X[i,j]是第i个样本的第j个特征

                total_loss += logit(y=y[i], y_hat=y_hat)  # 计算logit损失函数值
                dloss = df_logit(y=y[i], y_hat=y_hat)  # 计算logit损失函数的外层偏导

                dloss_w0 = dloss * 1  # 公式中的w0求导，计算复杂度O(1)
                self.w0 = self.w0 - self.alpha * dloss_w0

                for j in range(n):
                    if X[i, j] != 0:
                        dloss_Wj = dloss * X[i, j]  # 公式中的wi求导，计算复杂度O(n)
                        self.W[j] = self.W[j] - self.alpha * dloss_Wj
                        for f in range(self.k):  # 公式中的vif求导，计算复杂度O(kn)
                            dloss_Vjf = dloss * (X[i, j] * (X[i].dot(self.V[:, f])) - self.V[j, f] * X[i, j] ** 2)
                            self.V[j, f] = self.V[j, f] - self.alpha * dloss_Vjf

            print('iter={}, loss={:.4f}'.format(it+1, total_loss / m))

        return self

    def _FM_predict(self, X):
        predicts, threshold = [], 0.5  # sigmoid阈值设置
        for i in range(X.shape[0]):  # 遍历测试集
            y_hat = self._FM(Xi=X[i])  # FM的模型方程
            predicts.append(-1 if sigmoid(y_hat) < threshold else 1)
        return np.array(predicts)

    def fit(self, X, y):
        if isinstance(X, pd.DataFrame):
            X = np.array(X)
            y = np.array(y)

        return self._FM_SGD(X, y)

    def predict(self, X):
        if isinstance(X, pd.DataFrame):
            X = np.array(X)

        return self._FM_predict(X)

    def predict_proba(self, X):
        pass

from sklearn.metrics import roc_auc_score, confusion_matrix

model=FactorizationMachine(k=10, learning_rate=0.001, iternum=2)
model.fit(X_train, Y_train)

y_pred=model.predict(X_train)

print('训练集roc: {:.2%}'.format(roc_auc_score(Y_train.values, y_pred)))
print('混淆矩阵: \n',confusion_matrix(Y_train.values, y_pred))

y_true=model.predict(X_val)

print('验证集roc: {:.2%}'.format(roc_auc_score(Y_val.values, y_true)))
print('混淆矩阵: \n',confusion_matrix(Y_val.values, y_true))

from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import recall_score
from sklearn.metrics import precision_score
from sklearn.metrics import accuracy_score

X_val = MinMaxScaler().fit_transform(X_val)#归一化测试集，返回[0,1]区间

val_predicts = model._FM_predict(X_val)
print('FM测试集的分类准确率为: {:.2%}'.format(accuracy_score(Y_val,val_predicts)))
print("FM测试集均方误差mse：{:.2%}".format(mean_squared_error(Y_val,val_predicts)))
print("FM测试集召回率recall：{:.2%}".format(recall_score(Y_val,val_predicts)))
print("FM测试集的精度precision：{:.2%}".format(precision_score(Y_val,val_predicts)))