用Keras实现libFM

主要来自于kaggle比赛 Talking Data competition on Kaggle,构建一个模型来预测用户在移动端点击广告之后是否会下载应用app，主要数据字段如下:

ip: 点击广告行为的ip地址
app: app id
device: 用户手机的设备类型，比如（iPhone6 plus，iPhone7 等）
os：用户手机的系统版本号
channel：移动广告发布的channel id
click_time: 用户点击广告时间
attributed_time: 如果用户在点击广告之后下载了应用程序，则表示用户的下载时间
is_attributed: 预测变量，应用程序是否被下载

这份数据集中，ip，device和os组合一起基本上可以表示一名用户，但是某些情况下，存在一个ip地址被许多用户共享。在这次比赛中，这样的数据会增加问题的复杂性，为了简单化，忽略掉。另外channel数据似乎没啥作用，对其进行删除。赛方提供大约2亿条数据，目标是预测用户是否下载应用程序。

这让人想起推荐系统，它的目标是预测用户是否会购买产品(app就是这里的产品)。

推荐系统中最常用的方法是矩阵分解，主要计算用户和产品之间的表征矩阵信息，以便具有类似购买模式的用户的推荐相似的产品。一个非常简单的方法是为每个用户和每个产品计算一个k值的向量，K表示k个隐藏或潜在的特征。两向量间内积(点积)表示这两个向量的亲和力。

矩阵分解方法实际上相当简单。对于推荐系统，矩阵捕获过去的用户x产品交互(例如销售)信息。矩阵行代表用户，列代表产品，比如有M个用户和N个产品，则得到一个MxN矩阵C。其中元素C(i,j)表示用户i购买产品j的次数。C的大部分元素都是0。为了获得更好的准确性，通常对次数进行log化（当为0时，使用1进行代替），这样可以使用次数分布不偏向与大的数值。

表征向量可以通过对M×N矩阵C近似因子分解计算得到。我们的任务是找出两个矩阵,U (M x K)和P (K x N) (U – 用户, D – 产品)，使得 U x PT 近似等于矩阵C，其中U的每一行表示一个用户信息，P的每一列表征一个产品信息。如何得到U和P，一种方法是使用截断奇异值分解(tSVD)[scikit-learn实现了tSVD]。其他方法比如ALS也是可能的。

更多关于FM算法的详细信息可以参考文献:

Steffen Rendle (2010): Factorization Machines, in Proceedings of the 10th IEEE International Conference on Data Mining (ICDM 2010), Sydney, Australia [PDF](https://www.csie.ntu.edu.tw/~b97053/paper/Rendle2010FM.pdf) 此外，作者Rendle提供了算法的源码[libFM](http://www.libfm.org/)。我们可以直接使用libFM，但它需要不同于我使用的数据格式(熊猫或numpy)。它也是一个批处理命令行工具，不与我熟悉的Python环境集成。因此，我决定使用Keras作为一个深度学习模型来实现libFM方法。

下面我们FM模型代码进行说明，主要是：

k_latent = 2
embedding_reg = 0.0002
kernel_reg = 0.1
def get_embed(x_input, x_size, k_latent):
    if x_size > 0:  # category
        embed = Embedding(x_size, k_latent, input_length=1,embeddings_regularizer=l2(embedding_reg))(x_input)
        embed = Flatten()(embed)
    else:
        embed = Dense(k_latent, kernel_regularizer=l2(embedding_reg))(x_input)
    return embed

def build_model_1(X, f_size):
    dim_input = len(f_size)
    input_x = [Input(shape=(1,)) for i in range(dim_input)]
    biases = [get_embed(x, size, 1) for (x, size) in zip(input_x, f_size)]
    factors = [get_embed(x, size, k_latent) for (x, size) in zip(input_x, f_size)]
    s = Add()(factors)
    diffs = [Subtract()([s, x]) for x in factors]
    dots = [Dot(axes=1)([d, x]) for d, x in zip(diffs, factors)]
    x = Concatenate()(biases + dots)
    x = BatchNormalization()(x)
    output = Dense(1, activation='relu', kernel_regularizer=l2(kernel_reg))(x)
    model = Model(inputs=input_x, outputs=[output])
    opt = Adam(clipnorm=0.5)
    model.compile(optimizer=opt, loss='mean_squared_error')
    output_f = factors + biases
    model_features = Model(inputs=input_x, outputs=output_f)
    return model, model_features

接下来，分析每一部分，对于输入：

1	input_x = [Input(shape=(1,)) for i in range(dim_input)]

第一层是创造潜在的向量。这里我们对每一个类别变量创建两个潜在的向量，一个因子向量，和一个偏置向量：

1 2	biases = [get_embed(x, size, 1) for (x, size) in zip(input_x, f_size)] factors = [get_embed(x, size, k_latent) for (x, size) in zip(input_x, f_size)]

我们可以对所有的组合进行循环，计算两两因子相乘，假设我们有N个类别，就会得到N(N -1)/2的乘积。随着类别数量的增加，这可能导致内存问题。Rendle设计了一种巧妙的方法来线性地计算，我用一种稍微不同的方法来简化。让我来解释一下。我们想要计算:

$F =\sum_{i=1}^n \sum_{j = i+1}^n x_i x_j for all i, j such that i < j$

其中 $x_i$ 是第i类的因子向量。如果我们将

$S = \sum x_i$

得到

$F =\sum x_j (S - x_j) for all j$

这样我们可以看到变成了线性。

1
2
3

s = Add()(factors)
diffs = [Subtract()([s, x]) for x in factors]
dots = [Dot(axes=1)([d, x]) for d,x in zip(diffs, factors)]

添加bias：

1	x = Concatenate()(biases + dots)

我们几乎已经完成了。如果我们实现了原始的libFM，那么我们就应该使用x的和。当我们使用一个灵活的深度学习框架，我们可以使用一个mode来计算复杂的最后一层。而且可以添加batch normalization以得到更好的结果:

1 2	x = BatchNormalization()(x) output = Dense(1,activation='relu', kernel_regularizer=l2(kernel_reg))(x)

当模型搭建完成之后，进行编译。我们使用adam优化器，并增加clipnorm参数以避免出现梯度爆炸。

1
2
3

model = Model(inputs=input_x, outputs=[output])
opt = Adam(clipnorm=0.5)
model.compile(optimizer=opt, loss='mean_squared_error')

考虑到这个模型的目的是计算潜在的向量，我们需要通过训练模型来检索它们。为此，我们创建了一个二级模型，它与前面的模型共享嵌入层。

1 2	output_f = factors + biases model_features = Model(inputs=input_x, outputs=output_f)

接着训练模型。我使用一个相当大的batch大小2¹⁷。一般不建议使用这么大的值作为默认的batch大小。

n_epochs = 100
P = 17
batch_size = 2**P
earlystopper = EarlyStopping(patience=0, verbose=1)
model.fit(X_train,
          y_train,
          epochs=n_epochs,
          batch_size=batch_size,
          verbose=1,
          shuffle=True,
          validation_data=(X_train, y_train),
          sample_weight=w_train,
          callbacks=[earlystopper],
         )