因子分解机（Factorization Machine, FM）及其在推荐系统中的应用

Mike123
62 阅读7分钟

实验和完整代码

完整代码实现和jupyter运行：github.com/Myolive-Lin…

引言

在推荐系统和计算广告领域，准确预测用户行为是核心任务之一。因子分解机（Factorization Machines，简称FM）作为一种高效的预测模型，自2010年由Steffen Rendle提出以来，在处理稀疏数据和特征交互方面表现出色，广泛应用于点击率（CTR）预测和评分预测等任务

一、因子分解机的基本原理

模型公式

本质，FM引入隐向量的做法，与矩阵分解用隐向量代表用户和物品的做法异曲同工。可以说FM是将矩阵分解隐向量的思想进行了进一步扩展，从单纯的用户、物品隐向量扩展到所有特征上来

1.1 FM模型预测公式

FM的预测公式包含线性部分和二阶交互部分：

y^(x)=w0+∑i=1nwixi+∑i=1n∑j=i+1n(vi⋅vj)xixj\hat{y}(x) = w_0 + \sum_{i=1}^{n} w_i x_i + \sum_{i=1}^{n} \sum_{j=i+1}^{n} (v_i \cdot v_j) x_i x_j

其中：

  • w0w_0：全局偏置项
  • wiw_i：第 ii 个特征的线性权重
  • viv_i：第 ii 个特征的隐向量
  • xix_i：第 ii 个特征的输入值
  • vi⋅vjv_i \cdot v_j：第 ii 和第 jj 特征的隐向量内积

其中二阶段交互部分还可以进行优化从由O(kn2)O(kn^2) 到 O(kn)O(kn)

原本这个其实就是部分上三角，由于上下三角都是一样，用全部的一半减对角部分进行优化

∑i=1n∑j=i+1n(vi⋅vj)xixj=12∑i=1n∑j=1n(vi⋅vj)xixj−12∑i=1n(vi2)xi2=12∑i=1n∑j=1n∑i=1k(vi,f⋅vj,f)xixj−12∑i=1n∑i=1k(vi,f2)xi2=12∑i=1k(∑i=1n∑j=1n(vif⋅vjf)xixj−∑i=1n(vif)2xi2)=12∑i=1k((∑i=1nvifxi)(∑j=1nvjfxj)−∑i=1n(vif)2xi2)=12∑i=1k((∑i=1nvifxi)2−∑i=1n(vif)2xi2)\begin{aligned} \sum_{i=1}^{n} \sum_{j=i+1}^{n} (v_i \cdot v_j) x_i x_j &= \frac{1}{2} \sum_{i=1}^{n} \sum_{j=1}^{n} (v_i \cdot v_j) x_i x_j - \frac{1}{2} \sum_{i=1}^{n} (v_i^2) x_i^2\\ &= \frac{1}{2} \sum_{i=1}^{n} \sum_{j=1}^{n}\sum_{i=1}^{k} (v_{i,f} \cdot v_{j,f}) x_i x_j - \frac{1}{2} \sum_{i=1}^{n} \sum_{i=1}^{k}(v_{i,f}^2) x_i^2\\ &= \frac{1}{2} \sum_{i=1}^{k} \left( \sum_{i=1}^{n} \sum_{j=1}^{n} (v_{if} \cdot v_{jf}) x_i x_j - \sum_{i=1}^{n} (v_{if})^2 x_i^2 \right)\\ &= \frac{1}{2} \sum_{i=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)\\ &= \frac{1}{2} \sum_{i=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) \end{aligned}

1.2 优化后预测函数

y^(x)=w0+∑i=1nwixi+12∑i=1k(∑i=1nvifxi)2−∑i=1n(vif)2xi2)\hat{y}(x) = w_0 + \sum_{i=1}^{n} w_i x_i + \frac{1}{2} \sum_{i=1}^{k} ( \sum_{i=1}^{n} v_{if} x_i )^2 - \sum_{i=1}^{n} (v_{if})^2 x_i^2)

1.3 损失函数

FM模型的损失函数通常是均方误差（MSE）损失函数，表示为：

L(θ)=12∑t=1T(yt−yt^)2L(\theta) = \frac{1}{2} \sum_{t=1}^{T} (y_t - \hat{y_t})^2

其中：

  • TT 是训练样本的数量
  • yty_t 是第 tt 个样本的真实值
  • yt^\hat{y_t} 是第 tt 个样本的预测值
  • θ\theta 是模型的所有参数（包括线性权重和隐向量）

1.4 FM模型的梯度公式

对于FM模型的参数（包括线性权重 ww 和隐向量 vv），我们需要计算损失函数的梯度，以便使用梯度下降法进行优化。我们可以根据损失函数来推导出每个参数的梯度。

偏置项 w0w_0 梯度：

∂L∂w0=yt^−yt\frac{\partial L}{\partial w_0} = \hat{y_t} - y_t

线性权重 wiw_i 梯度：

∂L∂wi=(yt^−yt)⋅xi\frac{\partial L}{\partial w_i} = (\hat{y_t} - y_t) \cdot x_i

隐向量 viv_i 梯度：

∂L∂vi=(yt^−yt)⋅(∑j=1nvjxj−vixi)⋅xi=(yt^−yt)⋅(xi∑j=1nvjxj−vixi2)\begin{aligned} \frac{\partial L}{\partial v_i} = (\hat{y_t} - y_t) \cdot \left( \sum_{j=1}^{n} v_j x_j - v_i x_i \right) \cdot x_i \\ = (\hat{y_t} - y_t) \cdot \left( x_i \sum_{j=1}^{n} v_j x_j - v_i x_i^2\right) \end{aligned}

二、代码和实现部分

numpy 实现

class FactorizationMachine:
    def __init__(self, n_features, k, learning_rate=0.01, n_iter=1000):
        """
        初始化因子分解机模型。
        参数：
        - n_features: 特征数量
        - k: 隐向量的维度
        - learning_rate: 学习率
        - n_iter: 迭代次数
        """

        self.n_features = n_features
        self.k = k
        self.learning_rate = learning_rate
        self.n_iter = n_iter

        # 初始化参数
        self.w0 = 0
        self.w = np.zeros(n_features) # 线性权重
        self.V = np.random.normal(scale = 1.0 / k, size = (n_features, k))# 隐向量矩阵


    def _predict_instance(self, x):
        """
        预测单个样本的输出。
        参数：
        - x: 输入特征向量
        返回：
        - 预测值
        """
        linear_terms = self.w0 + np.dot(self.w, x)
        interactions = 0.5 * np.sum(
            np.square(np.dot(x, self.V)) - np.dot(x ** 2, self.V ** 2)
        )
        return linear_terms + interactions
         


    def fit(self,X,y):
        """
        训练因子分解机模型。
        参数：
        - X: 特征矩阵，形状为(n_samples, n_features)
        - y: 目标向量，形状为(n_samples,)
        """

        pbar = tqdm(range(self.n_iter), desc = "Training FM", ncols = 100, unit = "iter")
        for _ in pbar:
            rows, cols = X.shape
            mse = 0
            for i in range(rows):
                prediction = self._predict_instance(X[i])
                yi = y[i]
                error = (prediction - yi).astype(np.float32)
                mse += error ** 2
                

                #更新参数
                self.w0 -= self.learning_rate * error
                self.w -= self.learning_rate * error *X[i]

                # 计算 V 的更新，使用矩阵运算来代替逐元素的运算
                x_square = X[i] ** 2  # (n_features,)
                x_V = np.dot(X[i], self.V)  # (1,k)
                grad_V = error * (
                    #这纬度一致其实是矩阵乘法
                    np.dot(X[i].reshape(-1, 1), x_V.reshape(1, -1)) #(n_features,k)
                    - self.V * x_square.reshape(-1, 1) #(n_features,k) 广播机制
                )
                self.V -= self.learning_rate * grad_V
            pbar.set_postfix({"MSE": mse})

            
    def predict(self,X):
        """
        预测样本的输出。
        参数：
        - X: 特征矩阵，形状为(n_samples, n_features)
        返回：
        - 预测值数组，形状为(n_samples,)
        """
        return np.array([self._predict_instance(xi) for xi in X])

测试部分可以在上面github链接获取

torch实现

class FM(torch.nn.Module):
    def __init__(self,n_features, k):
        """
        初始化因子分解机模型。
        参数：
        - n_features: 特征数量
        - k: 隐向量的维度
        """

        super(FM,self).__init__()
        self.n_features = n_features
        self.k = k

        #初始化模型参数
        self.w0 = nn.Parameter(torch.zeros(1))
        self.w = nn.Parameter(torch.zeros(n_features))
        self.V = nn.Parameter(torch.zeros(n_features,k) * 0.01)

    def forward(self, X):
        """
        定义前向传播。
        参数：
        - X: 输入特征矩阵，形状为 (batch_size, n_features)
        返回：
        - 预测值 (batch_size,)
        """

        #线性部分
        linear_terms = self.w0 + torch.matmul(X, self.w)

        # 二阶交互部分 (优化计算)
        interactions = 0.5 * torch.sum(
            torch.pow(torch.matmul(X,self.V),2) - - torch.matmul(torch.pow(X, 2), torch.pow(self.V, 2)), dim=1 
        )

        return linear_terms + interactions


def train_fm(model,optimizer, criterion, X, y, n_epochs):
    """
    训练FM模型。
    参数：
    - model: FM模型实例
    - optimizer: 优化器
    - criterion: 损失函数
    - X: 训练数据 (torch.Tensor)
    - y: 目标值 (torch.Tensor)
    - n_epochs: 训练轮数
    """

    model.train()
    pbar = tqdm(range(n_epochs), desc="Training FM", ncols=100, unit="epoch")
    
    for epoch in pbar:
        predictions = model(X)
        loss = criterion(predictions, y)

        #反向传播   
        optimizer.zero_grad() #清除梯度
        loss.backward()       #计算梯度
        optimizer.step()      #更新参数
        
        #显示损失
        pbar.set_postfix({'Loss':loss.item()})

def predict_fm(model, X):
    model.eval()
    with torch.no_grad():
        predictions = model(X)
    return predictions

测试代码

#数据准备
n_features = 10
k = 5
X = torch.randn(100, n_features)
y = torch.randn(100)

fm_model = FM(n_features = n_features, k = k)
criterion = nn.MSELoss() # 使用均方误差损失函数
optimizer = optim.Adam(fm_model.parameters(), lr = 0.01)

train_fm(fm_model, optimizer = optimizer, criterion = criterion, X=X, y=y, n_epochs= 100)

y_pred = predict_fm(fm_model,X)

res_torch = pd.DataFrame(X)
res_torch['y'] = y
res_torch['y_pred'] = y_pred

res_torch
IndexFeature 0Feature 1Feature 2Feature 3Feature 4Feature 5Feature 6Feature 7Feature 8Feature 9yy_pred
0tensor(1.9269)tensor(1.4873)tensor(0.9007)tensor(-2.1055)tensor(0.6784)tensor(-1.2345)tensor(-0.0431)tensor(-1.6047)tensor(-0.7521)tensor(1.6487)0.1010670.237655
1tensor(-0.3925)tensor(-1.4036)tensor(-0.7279)tensor(-0.5594)tensor(-0.7688)tensor(0.7624)tensor(1.6423)tensor(-0.1596)tensor(-0.4974)tensor(0.4396)-1.309492-0.170638
2tensor(-0.7581)tensor(1.0783)tensor(0.8008)tensor(1.6806)tensor(1.2791)tensor(1.2964)tensor(0.6105)tensor(1.3347)tensor(-0.2316)tensor(0.0418)-0.4103580.360829
3tensor(-0.2516)tensor(0.8599)tensor(-1.3847)tensor(-0.8712)tensor(-0.2234)tensor(1.7174)tensor(0.3189)tensor(-0.4245)tensor(0.3057)tensor(-0.7746)0.468094-0.080337
4tensor(-1.5576)tensor(0.9956)tensor(-0.8798)tensor(-0.6011)tensor(-1.2742)tensor(2.1228)tensor(-1.2347)tensor(-0.4879)tensor(-0.9138)tensor(-0.6581)-0.234628-0.057512
.......................................
95tensor(0.5692)tensor(-0.7911)tensor(-0.1990)tensor(-1.3616)tensor(-0.5194)tensor(0.0765)tensor(0.3401)tensor(1.4557)tensor(-0.3461)tensor(-0.2634)-0.014565-0.286557
96tensor(-0.4477)tensor(-0.7288)tensor(-0.1607)tensor(-0.3206)tensor(-0.6308)tensor(-0.7888)tensor(1.3062)tensor(-0.9276)tensor(-0.2627)tensor(0.9315)-0.748688-0.001467
97tensor(-0.4593)tensor(-0.9419)tensor(-0.7089)tensor(2.1861)tensor(-0.6493)tensor(0.4521)tensor(0.8521)tensor(-1.6947)tensor(1.1806)tensor(-2.8929)-0.384403-0.337158
98tensor(-0.3876)tensor(-0.7124)tensor(-1.6171)tensor(-0.3590)tensor(-0.4137)tensor(-0.5285)tensor(-0.5082)tensor(1.1478)tensor(0.2401)tensor(-0.7907)0.143502-0.181565
99tensor(-0.8132)tensor(0.1415)tensor(1.5630)tensor(0.5983)tensor(-0.5407)tensor(0.7842)tensor(0.6991)tensor(-0.0363)tensor(-0.3798)tensor(-0.8535)-0.081268-0.101379

三、FM 的特点与优势

  1. 灵活性与高效性
    FM 能够自动建模二阶特征交互，无需显式生成交互特征，避免了特征工程的复杂性。

  2. 适应高维稀疏数据
    FM 的参数数量与特征维度呈线性关系，因此在处理高维稀疏数据时具有显著优势。

  3. 易扩展性
    FM 可与其他模型结合，如深度神经网络（DNN），形成如 DeepFM 的深度交互模型，以提升预测性能。

  4. 良好的解释性
    FM 的线性部分与特征交互部分能够独立建模，因此在解释模型输出时较为直观。

推荐系统的核心挑战

推荐系统的核心任务是预测用户对物品的兴趣程度。由于用户行为数据通常具有以下特点，传统模型往往面临挑战：

  • 高维稀疏性：用户和物品的交互矩阵通常是稀疏的。
  • 非线性特征交互：用户兴趣往往由多个特征交互决定，例如用户年龄与物品类型的组合。

FM 能够高效建模特征交互，同时避免显式构造组合特征带来的维度爆炸问题，因此成为在过去推荐系统的核心方法之一。

应用场景

  • 个性化推荐 FM 在个性化推荐中表现优异。通过对用户特征（如年龄、性别、偏好）和物品特征（如类型、价格、评分）之间的交互建模，FM 能够预测用户对物品的评分或点击率。

  • 点击率（CTR）预测 在广告推荐领域，FM 被广泛用于 CTR 预测任务。CTR 预测的目标是估计用户点击广告的概率，其关键在于建模用户与广告特征的复杂交互关系。FM 能够高效捕捉这些关系，从而显著提升预测精度。

  • 冷启动问题 对于新用户或新物品，FM 能够通过隐向量分解，结合其他侧信息（如人口统计学特征或商品属性），有效缓解冷启动问题。

四、FM 与其他推荐算法的比较

算法优点局限性
矩阵分解适用于评分预测，参数少无法直接建模特征交互
因子分解机（FM）自动建模特征交互，适应高维稀疏数据仅建模二阶特征交互
深度因子分解机（DeepFM）能够建模高阶特征交互，预测能力更强计算复杂度更高，训练时间较长

Reference

王喆 《深度学习推荐系统》
FM因子分解机的原理、公式推导、Python实现和应用

avatar
文章
阅读
粉丝