推荐系统（三）Factorization Machines（FM）

推荐系统（三）Factorization Machines（FM）

推荐系统系列博客：
推荐系统（一）推荐系统整体概览
推荐系统（二）GBDT+LR模型

按照发表年份，这篇博客应该在GBDT+LR之前写的，但因为FM相比较GBDT+LR的内容稍微多些，所以就后写了这篇博客。言归正传，FM是推荐系统领域大佬rendle于2010年发表在ICDM上的论文，是一篇非常非常有影响力的论文，启发了此后10年学术界大量的工作，直接的改进就有引入神经网络的NFM，引入attention的AFM等（关于NFM和AFM这两个能不能像FM一样在工业界一样work，这个…，不过话又说回来，学术界的价值就在于提供idea，大规模商用另说了）。在详细介绍FM之前，先用一句话概括下FM：隐式向量特征交叉。这么一看FM在那个时候已经孕育了embedding的思想。接下来会介绍LR、SVM、MF，这三是FM的灵感源泉，这样看起来条理也比较清晰。

1. LR
先来回顾下LR：
$$\hat{y}(x)=w_0+\sum _{i=1}^n{w}_i{x}_i\text{\hspace{1em}(1)}$$
那么，如果我们想做特征的二阶交叉，很自然的可以把LR改造为如下的函数：
$$\hat{y}(x)=w_0+\sum _{i=1}^n{w}_i{x}_i+\sum _{i=1}^n\sum _{j=i+1}^n{w}_{ij}{x}_i{x}_j\text{\hspace{1em}(2)}$$
但是公式2明显有两个缺点：1. 时间段复杂度是O(N^2)，2.依赖于$x_i{x}_j$特征对的共现，如果这对特征在训练集中没有出现，那么这个参数学习不到。

2. SVM
我在很久很久以前的博客里写过一篇最基本最浅显的SVM（现在回头来看着实有些浅显，没有讲解到核函数的本质），参见：支持向量机SVM，支持向量机干的事就是既然低维空间下找不到一个超平面来划分两类样本，那么我映射下，把低维空间映射成高维空间来寻得一个超平面好了。所以，假如我们用$\varphi (x)$表示将$x$映射后的特征向量，支持向量机模型为：
$$\hat{y}(x)=w^T\varphi (x)+b\text{\hspace{1em}(3)}$$
一顿对偶猛如虎后，发现要算$\varphi (x_i{)}^T\varphi (x_j)$，也就是说要算$x_i$与$x_j$映射之后的内积，那这个复杂度有点略高啊（先映射成高维，再算内积），所以得想办法降低复杂度，那么核函数就出场了，假如我们有个核函数：
$$k(x_i,x_j)=<\varphi (x_i),\varphi (x_j)>=\varphi (x_i{)}^T\varphi (x_j)\text{\hspace{1em}(4)}$$
能够使得$x_i$与$x_j$在映射后的特征空间的内积等于他们在原始样本空间中通过$k(.,.)$函数计算的结果。

2.1 线性核函数
对于线性核$k(x_i,x_j)=1+<x_i,x_j>$，映射函数$\varphi (x)=(1,x_1,x_2,...,x_n)$，那么linear SVM模型为：
$$\hat{y}(x)=w_0+\sum _{i=1}^n{w}_i{x}_i\text{\hspace{1em}(5)}$$
乍一看仿佛和LR没啥区别啊，但实际上还是有区别的，因为loss不同和SVM需要支持向量等。

2.2 多项式核函数
多项式核 $k(x_i,x_j)=(<x_i,x_j>+1{)}^2$，映射函数 $\varphi (x)=(1,\sqrt{2}{x}_1,\sqrt{2}{x}_2,...,\sqrt{2}{x}_n,x_1^2,x_2^2,...,x_n^2,\sqrt{2}{x}_1{x}_2,\sqrt{2}{x}_1{x}_3,...,\sqrt{2}{x}_1{x}_n,\sqrt{2}{x}_2{x}_3...,\sqrt{2}{x}_{n-1}{x}_n)$，因此多项式核的支持向量机模型为：
$$\hat{y}(x)=w_0+\sqrt{2}\sum _{i=1}^n{w}_i{x}_i+\sum _{i=1}^n{w}_{ii}^2{x}_{ii}^2+\sqrt{2}\sum _{i=1}^n\sum _{j=i+1}^n{w}_{ij}^2{x}_i{x}_j\text{\hspace{1em}(6)}$$
我们看到二项式核函数的SVM和公式(2)中存在一样的问题，就是交叉项的参数是独立的，这会使得如果这个交叉特征值没有在样本里出现，这个参数是无法学到的。

所以，我们总结下，目前在模型层面做交叉特征的难点主要有以下两个方面：

• 交叉特征的参数独立，强依赖于在样本中的共现信息，如果交叉特征值没有出现，那么参数无法学习。
• 时间复杂度问题，直接做二阶交叉，时间复杂度为O(N^2)，复杂度过高。

3. FM
那么FM则解决了上面两个问题，直接上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\text{\hspace{1em}(7)}$$
其中$<\cdot ,\cdot >$为两个k维的向量的点积，即数量积。$v_i$表示第$i$个特征的向量。
$$<v_i,v_j>=\sum _{f=1}^k{v}_{i,f}\cdot v_{j,f}\text{\hspace{1em}(8)}$$
因此，公式（7）完整的为：
$$\hat{y}(x)=w_0+\sum _{i=1}^n{w}_i{x}_i+\sum _{i=1}^n\sum _{j=i+1}^n\sum _{f=1}^k{v}_{i,f}\cdot v_{j,f}{x}_i{x}_j\text{\hspace{1em}(9)}$$
因此，公式(7)这里实际可分为三部分，第一部分一个偏置单元 $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\text{\hspace{1em}(10)}$$
FM这里巧妙的把公式(2)中的独立参数$w_{ij}$分解成了$<v_i,v_j>$，实际上是通过学习每一个特征对应的隐向量（现在大家熟知的embedding向量），这样就不再依赖于交叉特征$x_i{x}_j$的共现信息，因为即使$x_i{x}_j$没有共现，假如$x_i{x}_k$有共现，那么$x_i$对应的隐向量$v_i$依然能够得到训练。

那么问题来了，这种方法是FM独创的吗，答案是：NO。这种思想来源于一个古老且有效的方法矩阵分解MF（matrix factorization），在推荐系统里，每个用户对每个物品的评分，可以构建出一个user-item矩阵，而矩阵分解的核心思想是用一个用户embedding矩阵和一个物品embedding矩阵的乘积来近似这个大矩阵，这两个embedding矩阵是可训练学习的。上一张图来形象的表示矩阵分解（图片来源：Introduction to Matrix Factorization - Collaborative filtering with Python 12）

到这里，可以看到FM解决了我们前面抛出的两个问题中的第一个问题（交叉特征参数独立，依赖于交叉特征的共现），下面来看看FM如何解决第二个问题（时间复杂度问题），再来看看公式（7）中的二阶部分，时间复杂度为$O(N^2)$，FM把这部分做了个公式推导，把时间复杂度降到了$O(KN)$，下面来看看FM的推到过程：
$$\begin{array}{rl}\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-\sum _{i=1}^n<v_i,v_i>x_i{x}_i]\\ & =\frac{1}{2}(\sum _{i=1}^n\sum _{j=1}^n\sum _{f=1}^k{v}_{i,f}\cdot v_{j,f}{x}_i{x}_j-\sum _{i=1}^n\sum _{f=1}^k{v}_{i,f}\cdot v_{i,f}{x}_i{x}_i)\\ & =\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}^2{x}_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}^2{x}_i^2)\end{array}\text{\hspace{1em}(11)}$$
关于上面这个公式，第一步到第二步，大家想象一个矩阵，行和列都是$x_1,....,x_n$，第一步为矩阵的上三角，所以等于全矩阵减去对角线，再折半，这样就比较好理解了。

4. FM实现

1.作者实现版本
rendle大佬用C++自己coding了个FM，并且开源了，地址：libFM，大佬就是大佬，不仅学术能力了解，coding能力也不逞多让，不得不服。下面来学习下大佬的代码，这里只列核心代码，讲解参见我的注释

double fm_model::predict(sparse_row<FM_FLOAT>& x, DVector<double> &sum, DVector<double> &sum_sqr) {
  double result = 0;
  // 公式(7)中的第一项，w0
  if (k0) {
    result += w0;
  }
  // 第二项，一阶部分w_i*x_i
  if (k1) {
    for (uint i = 0; i < x.size; i++) {
      assert(x.data[i].id < num_attribute); // num_attribute --> feature num
      result += w(x.data[i].id) * x.data[i].value;
    }
  }
  // 第三项，二阶部分
  for (int f = 0; f < num_factor; f++) {   // num_factor --> embedding size
    sum(f) = 0;
    sum_sqr(f) = 0;
    for (uint i = 0; i < x.size; i++) {
      double d = v(f,x.data[i].id) * x.data[i].value;
      sum(f) += d;
      sum_sqr(f) += d*d;
    }
    result += 0.5 * (sum(f)*sum(f) - sum_sqr(f));
  }
  return result;
}

// sgd，梯度部分
void fm_SGD(fm_model* fm, const double& learn_rate, sparse_row<DATA_FLOAT> &x, const double multiplier, DVector<double> &sum) {
  if (fm->k0) {
    double& w0 = fm->w0;
    w0 -= learn_rate * (multiplier + fm->reg0 * w0);
  }
  if (fm->k1) {
    for (uint i = 0; i < x.size; i++) {
      double& w = fm->w(x.data[i].id);
      w -= learn_rate * (multiplier * x.data[i].value + fm->regw * w);
    }
  }
  for (int f = 0; f < fm->num_factor; f++) {
    for (uint i = 0; i < x.size; i++) {
      double& v = fm->v(f,x.data[i].id);
      double grad = sum(f) * x.data[i].value - v * x.data[i].value * x.data[i].value;
      v -= learn_rate * (multiplier * grad + fm->regv * v);
    }
  }
}

2.paddlepaddle版本
百度官方开源了paddlepaddle版本的FM实现，文档介绍的相对详细点，参见：https://github.com/PaddlePaddle/PaddleRec/tree/release/2.1.0/models/rank/fm
因为使用的数据集包含离散特征和连续值特征两类，一种办法是连续特征离散化，这样所有特征都当做离散特征处理，处理起来比较简单，代码比较统一。另一种就是构建网络时离散值特征和连续值特征分开考虑，paddle版本FM采用了后者。因此，对应对于连续值特征直接w*x即可，而离散特征的embedding向量也比较容易获取，直接根据index查表取出来即可，这里要注意的是：特征分量$x_i$与$x_j$的交叉项系数就等于$x_i$对应的隐向量与$x_j$对应的隐向量的内积（关于这一点，可能有人比较迷惑，解释下，离散特征onehot后只有一个分量是1，所以直接取embedding向量没问题）

def forward(self, sparse_inputs, dense_inputs):
        # -------------------- first order term  --------------------
        sparse_inputs_concat = paddle.concat(sparse_inputs, axis=1)
        sparse_emb_one = self.embedding_one(sparse_inputs_concat)

        dense_emb_one = paddle.multiply(dense_inputs, self.dense_w_one)
        dense_emb_one = paddle.unsqueeze(dense_emb_one, axis=2)

        y_first_order = paddle.sum(sparse_emb_one, 1) + paddle.sum(
            dense_emb_one, 1)

        # -------------------- second order term  --------------------
        sparse_embeddings = self.embedding(sparse_inputs_concat)
        dense_inputs_re = paddle.unsqueeze(dense_inputs, axis=2)
        dense_embeddings = paddle.multiply(dense_inputs_re, self.dense_w)
        feat_embeddings = paddle.concat([sparse_embeddings, dense_embeddings],
                                        1)

        # sum_square part
        summed_features_emb = paddle.sum(feat_embeddings,
                                         1)  # None * embedding_size
        summed_features_emb_square = paddle.square(
            summed_features_emb)  # None * embedding_size

        # square_sum part
        squared_features_emb = paddle.square(
            feat_embeddings)  # None * num_field * embedding_size
        squared_sum_features_emb = paddle.sum(squared_features_emb,
                                              1)  # None * embedding_size

        y_second_order = 0.5 * paddle.sum(
            summed_features_emb_square - squared_sum_features_emb,
            1,
            keepdim=True)  # None * 1

        return y_first_order, y_second_order

3.tensorflow版本
网上实现的版本比较杂，挑了一个实现比较好的tensorflow版本DeepFM，参见：https://github.com/ChenglongChen/tensorflow-DeepFM，可以看看FM部分的实现。

# model
            self.embeddings = tf.nn.embedding_lookup(self.weights["feature_embeddings"],
                                                             self.feat_index)  # None * F * K
            feat_value = tf.reshape(self.feat_value, shape=[-1, self.field_size, 1])
            self.embeddings = tf.multiply(self.embeddings, feat_value)

            # ---------- first order term ----------
            self.y_first_order = tf.nn.embedding_lookup(self.weights["feature_bias"], self.feat_index) # None * F * 1
            self.y_first_order = tf.reduce_sum(tf.multiply(self.y_first_order, feat_value), 2)  # None * F
            self.y_first_order = tf.nn.dropout(self.y_first_order, self.dropout_keep_fm[0]) # None * F

            # ---------- second order term ---------------
            # sum_square part
            self.summed_features_emb = tf.reduce_sum(self.embeddings, 1)  # None * K
            self.summed_features_emb_square = tf.square(self.summed_features_emb)  # None * K

            # square_sum part
            self.squared_features_emb = tf.square(self.embeddings)
            self.squared_sum_features_emb = tf.reduce_sum(self.squared_features_emb, 1)  # None * K

            # second order
            self.y_second_order = 0.5 * tf.subtract(self.summed_features_emb_square, self.squared_sum_features_emb)  # None * K
            self.y_second_order = tf.nn.dropout(self.y_second_order, self.dropout_keep_fm[1])  # None * K

参考

[1]:Factorization Machines