一、引入FM

在传统的线性模型如LR中，每个特征都是独立的，如果需要考虑特征与特征直接的交互作用，可能需要人工对特征进行交叉组合；非线性SVM可以对特征进行kernel映射，但是在特征高度稀疏的情况下，并不能很好地进行学习。

因子分解机(FactorizationMachine,FM)是由SteffenRendle在2010年提出的一种基于矩阵分解的机器学习算法。算法的核心在于特征组合，以此来减少人工参与特征组合工作。 在计算广告和推荐系统中，CTR预估(click-throughrate)是非常重要的一个环节，判断一个物品是否进行推荐需要根据CTR预估的点击率排序决定。 下面以一个示例引入FM模型。假设一个广告分类的问题，根据用户和广告位相关的特征，预测用户是否点击了广告。

我们手动生成源数据如下：

import pandas as pd
train = [
{"country": "China", "type": "Game", "age": 19,"Clicked":1},
{"country": "USA", "type": "Game", "age": 33,"Clicked":1},
{"country": "UK", "type": "Music", "age": 55,"Clicked":1},
{"country": "China", "type": "Movie", "age": 20,"Clicked":0},
{"country": "China", "type": "Music", "age": 25,"Clicked":1},
]
df = pd.DataFrame(train)

“Clicked”是label，Age、Country、Type是特征。由于前后种特征都是categorical类型的，需要经过独热编码（One-HotEncoding）转换成数值型特征。

pd.get_dummies(df)

由上表可以看出，经过One-Hot编码之后，大部分样本数据特征是比较稀疏的。

上面的样例中，每个样本有7维特征，但平均仅有3维特征具有非零值。实际上，这种情况并不是此例独有的，在真实应用场景中这种情况普遍存在。例如，CTR/CVR预测时，用户的性别、职业、教育水平、品类偏好，商品的品类等，经过One-Hot编码转换后都会导致样本数据的稀疏性。特别是商品品类这种类型的特征，如商品的末级品类约有几百个，采用One-Hot编码生成几百个数值特征，但每个样本的这几百个特征，有且仅有一个是有效的（非零）。由此可见，数据稀疏性是实际问题中不可避免的挑战。 One-Hot编码的另一个特点就是导致特征空间大。

例如，商品品类有100维特征，一个categorical特征转换为100维数值特征，特征空间剧增。 同时通过观察大量的样本数据可以发现，某些特征经过关联之后，与label之间的相关性就会提高。例如，“USA”与“Thanksgiving”、 “China” 与 “ChineseNewYear” 这样的关联特征，对用户的点击有着正向的影响。

换句话说：来自“China”的用户很可能会在“ChineseNewYear”有大量的浏览、购买行为；而在“Thanksgiving”却不会有特别的消费行为。这种关联特征与label的正向相关性在实际问题中是普遍存在的，如<“化妆品”类商品，“女”性>，<“球类运动配件”的商品，“男”性>，<“电影票”的商品，“电影”>品类偏好等。因此，引入两个特征的组合是非常有意义的。

 

二、引入隐变量V

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

其中，n代表样本的特征数量，xi是第i个特征的值，w0、wi、wij是模型参数，只有当xi与xj都不为0时，交叉才有意义在数据稀疏的情况下，满足交叉项不为0的样本将非常少，当训练样本不足时，很容易导致参数wij训练不充分而不准确，最终影响模型的效果。 那么，交叉项参数的训练问题可以用矩阵分解来近似解决，有下面的公式。

对任意正定矩阵W，只要k足够大，就存在矩阵W，使得W=VVT。然而在数据稀疏的情况下，应该选择较小的k，因为没有足够的数据来估计wj。限制k的大小提高了模型更好的泛化能力。 为什么说可以提高模型的泛化能力呢？ 我们用原论文的例子看一下

假设我们要计算用户A与电影ST的交叉项，很显然，训练数据里没有这种情况，这样会导致=0，但是我们可以近似计算出。首先，用户B和C有相似的向量VB和Vc ，因为他们对SW的预测评分比较相似， 所以 和 是相似的。用户A和C有不同的向量，因为对TI和SW的预测评分完全不同。接下来，SW和ST的向量可能相似，因为用户B对这两部电影的评分也相似。最后可以看出，与 是相似的。

我们先用代码把我们文章开头的例子进行测算一下

#说明：DictVectorizer的处理对象是符号化(非数字化)的但是具有一定结构的特征数据，如字典等，将符号转成数字0/1表示
v = DictVectorizer()
X = v.fit_transform(train)
y = array([1,1,1,0,1])
fm = pylibfm.FM(num_factors=5, 
                num_iter=200, 
                #verbose=True, 
                task="classification", 
                initial_learning_rate=0.001, 
                learning_rate_schedule="optimal")
fm.fit(X,y) 

fm.predict(v.transform({"country": "USA", "type": "Movie", "age": 33}))

array([0.05668886])

对于预测数据，"country": "USA"和"type": "Movie"，是没有交集的，但是我们可以发现第一条数据和第四条数据中，type=Game 和type=Movie，一个点击一个不点击，说明这两个是相反的；而第一条数据和第二条数据说明"country"= "USA"和"country"= "China"相似，那么测试中的这个也不会点击。

fm.predict(v.transform({"country": "USA", "type": "Music", "age": 20}))

 array([0.98671174])

同样我们可以发现第一条数据和第五条数据中，type=Game 和type=Music，都进行了点击，说明这两个是相同的。而第一条数据和第二条数据说明"country"= "USA"和"country"= "China"相似，那么测试中的这个很有可能点击。

 

三、公式变换

通过公式变换，可以减少到线性复杂度，方法如下(利用了xy= [(x+y)^2 – x^2 – y^2]/2 的思路)：

FM可以应用于很多预测任务，比如回归、分类、排序等等。

1.回归Regression：y(x)直接作为预测值，损失函数可以采用least square error；

2.二值分类Binary Classification：y(x)需转化为二值标签，如0,1。损失函数可以采用hinge loss或logit loss；

3.排序Rank：x可能需要转化为pair-wise的形式如(Xa,Xb)，损失函数可以采用pairwise loss

 

四、损失函数

在二分类问题中使用logit loss作为优化的标准，即

注解：

可能有人会问这个损失函数为什么跟我们之前见的损失函数不太一样？

一般情况下，损失函数是

这种情况是Y∈{0,1}但是，这个函数在我们进行代码计算时会比较繁琐，所以我们进行转换

Y∈{-1,1}，这样就需要用上面的损失函数来进行梯度下降求解了

基于随机梯度的方式求解

其中

五、算法流程

1、初始化权重w0w1……wn和V

2、对每个样本进行计算：

 

六、python实战

我们分别用梯度下降迭代和pylibfm直接预测两种方法进行

我们先用梯度下降迭代进行演示

1>梯度下降法

我们用糖尿病数据集，先倒入训练集

#定义导入数据的函数
def loadTrainDataSet(data):
    global max_list
    global min_list
    dataMat = []
    labelMat = []
    
    fr = open(data)#打开文件
    
    for line in fr.readlines():
        currLine = line.strip().split(',')
        #lineArr = [1.0]
        lineArr = []
        
        for i in range(featureNum):
            #数据集有9列数据，前8列为变量，最后为标签
            lineArr.append(float(currLine[i]))

        dataMat.append(lineArr)
        #标签从【0，1】转换为【-1，1】
        labelMat.append(float(currLine[-1]) * 2 - 1)
    
    data_array = np.array(dataMat)
    #取出每列最大值最小值，进行归一化处理准备
    max_list = np.max(data_array,axis=0)
    min_list = np.min(data_array,axis=0)

    scalar_dataMat = []
    for row in dataMat:
        #归一化处理
        scalar_row = normalize(row,max_list,min_list)
        scalar_dataMat.append(scalar_row)
    return scalar_dataMat, labelMat

#梯度下降迭代
def stocGradAscent(dataMatrix, classLabels, k, iter):
    #dataMatrix用的是mat, classLabels是列表
    m, n = shape(dataMatrix)
    alpha = 0.01#学习率
    #初始化参数
    #w = random.randn(n, 1)#其中n是特征的个数
    w = zeros((n, 1))
    w_0 = 0.
    v = normalvariate(0, 0.2) * ones((n, k))
    
    for it in range(iter):
        #print (it)
        for x in range(m):#随机优化，对每一个样本而言的
            inter_1 = dataMatrix[x] * v
            inter_2 = multiply(dataMatrix[x], dataMatrix[x]) * multiply(v, v)#multiply对应元素相乘
            #完成交叉项
            interaction = sum(multiply(inter_1, inter_1) - inter_2) / 2.
            
            p = w_0 + dataMatrix[x] * w + interaction#计算预测的输出
            #print "y: ",p 
            loss = sigmoid(classLabels[x] * p[0, 0]) - 1
            #print "loss: ",loss
            
            #更新w_0
            w_0 = w_0 - alpha * loss * classLabels[x]
            
            for i in range(n):
                if dataMatrix[x, i] != 0:
                    #更新w_i
                    w[i, 0] = w[i, 0] - alpha * loss * classLabels[x] * dataMatrix[x, i]
                    #更新v
                    for j in range(k):
                        v[i, j] = v[i, j] - alpha * loss * classLabels[x] * (dataMatrix[x, i] * inter_1[0, j] - v[i, j] * dataMatrix[x, i] * dataMatrix[x, i])
        
    
    return w_0, w, v

最后进行准确率的测算

def getAccuracy(dataMatrix, classLabels, w_0, w, v):
    m, n = shape(dataMatrix)
    allItem = 0
    error = 0
    result = []
    for x in range(m):
        allItem += 1
        #带入参数w_0, w, v
        inter_1 = dataMatrix[x] * v
        inter_2 = multiply(dataMatrix[x], dataMatrix[x]) * multiply(v, v)#multiply对应元素相乘
        #完成交叉项
        interaction = sum(multiply(inter_1, inter_1) - inter_2) / 2.
        p = w_0 + dataMatrix[x] * w + interaction#计算预测的输出        
        pre = sigmoid(p[0, 0])
        
        result.append(pre) 
        if pre < 0.5 and classLabels[x] == 1.0:
            error += 1
        elif pre >= 0.5 and classLabels[x] == -1.0:
            error += 1
        else:
            continue            
    print (result)   
    return float(error) / allItem
        

最后我们就可以进行测试了

if __name__ == '__main__':
    dataTrain, labelTrain = loadTrainDataSet(trainData)
    dataTest, labelTest = loadTestDataSet(testData)
    w_0, w, v = stocGradAscent(mat(dataTrain), labelTrain, 20, 200)
    print ("测试准确性为：%f" % (1 - getAccuracy(mat(dataTest), labelTest, w_0, w, v))  )

最终预测的准确率为：

测试准确性为：0.791045

2、pylibfm进行预测

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from pyfm import pylibfm
from sklearn.feature_extraction import DictVectorizer
from scipy.sparse import csr_matrix
#压缩矩阵
dataTrain = csr_matrix(dataTrain)
fm= pylibfm.FM(num_factors=8, 
                num_iter=200, 
                verbose=True, 
                task="classification", 
                initial_learning_rate=0.001, 
                learning_rate_schedule="optimal")
fm1.fit(dataTrain, labelTrain)
dataTest = csr_matrix(dataTest)
y_preds = fm.predict(dataTest)
y_preds_label = y_preds > 0.5
from sklearn.metrics import log_loss,accuracy_score
print ("Validation log loss: %.4f" % log_loss(labelTest, y_preds))
print ("accuracy: %.4f" % accuracy_score(labelTest, y_preds_label))

最终输出结果为：

Validation log loss: 0.4851
accuracy: 0.7948

 

最后补充一点关于pylibfm的知识

安装 pip install git+https://github.com/coreylynch/pyFM

参数设置

num_factors : int, The dimensionality of the factorized 2-way interactions
num_iter : int, 迭代次数
learning_rate_schedule（可选参数） : string类型, 学习率迭代参数值:

initial_learning_rate : double, 初始时的学习率，默认等于 0.01
task : 模型任务类型，string
regression: Labels 为实数值.
classification: Labels 为1或0.
verbose : bool，是否需要打印当前迭代时的训练误差（training error）.
power_t : double, The exponent for inverse scaling learning rate [默认值 0.5].
t0 : double，learning_rate_schedule的初始值. 默认值 0.001.

 

完整代码及数据可以关注公众号：不懂乱问

回复FM获取

 

 

参考：

https://www.csie.ntu.edu.tw/~b97053/paper/Rendle2010FM.pdf

https://zhuanlan.zhihu.com/p/50426292

https://blog.csdn.net/google19890102/article/details/45532745

https://zhuanlan.zhihu.com/p/37963267

 https://www.csuldw.com/2019/02/08/2019-02-08-fm-algorithm-theory/