1. 什么是FM？
FM即Factor Machine，因子分解机。
2. 为什么需要FM？
1.特征组合是许多机器学习和建模过程中遇到的问题。如果直接建模特征，很可能会忽略特征与特征之间的相关信息。因此，模型的效果可以通过构建新的交叉特征组合来提高。
2.高维稀疏矩阵是实际工程中常见的问题，直接导致计算量过大，特征权值更新缓慢。想象一下1万*每列100表有8种元素，经过one-hot独热编码后，将产生1万个*800表。因此，表中每行元素只有100个值为1，700个值为0。
而FM其优点是处理这两个问题。一是特征组合，通过两个特征组合，引入交叉特征，提高模型得分；二是高维灾难，通过引入隐向量（参数矩阵分解），完成特征参数估计。
3. FM用在哪？
我们已经知道了FM它可以解决特征组合和高维稀疏矩阵的问题。在实际业务场景中，电子商务、豆瓣等推荐系统的场景是应用最广泛的领域。例如，小王只在豆瓣上浏览了20部电影，豆瓣上有2万部电影。如果构建一个基于小王的电影矩阵，毫无疑问，将有19980个元素为0。这样的问题可以通过FM来解决。
4. FM长什么样？
在展示FM在算法之前，让我们回顾一下最常见的线性表达：
其中w0对于初始权值，或理解为偏置项，wi为每个特征xi相应的权值。这种线性表达式只描述了每个特征与输出之间的关系。
FM表达式如下，可以观察到新的交叉特征和相应的权了新的交叉特征和相应的权重。
5. FM展开交叉项
5.1 寻找交叉项
FM表达式解决方案的核心在于对交叉项的解决方案。以下是许多人用来解决交叉项的扩展，对于第一次接触FM算法的人可能会有疑问，不知道公式是怎么展开的，然后作者会手动推导。
三个变量(特征)x1 x2 x3.每个特征的隐变量分别为v1=(1 2 3)、v2=(4 5 6)、v3=(1 2 1)，即：
由交叉项组成的权矩阵W是对称矩阵，它被设置为对称矩阵向量乘以向量转移所取代。
那么W=VVT，即
所以：
事实上，我们应该考虑的交叉项应该是排除自己组合的项目，即对于x1x1、x2x2、x3x3不认为是交叉项，所以真正的交叉项是x1x2、x1x3、x2x1、x2x3、x3x1、x3x2。
去重后，交叉项就是x1x2、x1x3、x2x3.这也是公式中1/2出现的原因。
5.2 交叉项权值转换
对交叉项有了基本了解后，下面将进行公式的分解，还是以n=3为例，
所以：
wij可记作
或
，这取决于vi是1*3 还是3*1 向量。
5.3 展开式交叉项
以上例子是三个特征的交叉项推导，因此具有n个特征，FM交叉项公式可推广为：
我们也可以进一步分解：
所以FM算法的交叉项最终可以扩展为：
5.4隐向量v是embedding vector?
假设训练数据集dataMatrix的shape为(20000，9)，取其中一行数据作为一条样本i，那么样本i 的shape假设隐向量为(1，9)vi的shape(9，8)(注:8为自定义值，代表embedding vector的长度)
所以5.三节中的交叉项可以表示为：
sum((inter_1)^2 - (inter_2)^2)/2
其中：
inter_1 =i*v shape为(1，8)
inter_2 =np.multiply(i)*np.multiply(v) shape为(1，8)
样本可见i 计算交叉项后，得到向量shape为(1，8)的inter_1和inter_2。
由于维度较低，这个计算过程可以被认为是交叉项中的样本i 进行了embedding vector转换。
因此，我们需要纠正以前的理解：
我们嘴里的隐向量vi它实际上是一个向量组，其形状是(输入特征One-hot自定义长度)；
隐向量vi不代表embedding vector，而是输入embedding vector向量组，也可以理解为权矩阵；
由输入i*vi得到的向量才是真正的embedding vector。
可以结合第7节点的代码来理解。
6. 权值求解
梯度通过求损函数对特征(输入项)的导数计算入项)的导数计算梯度，从而更新权值。将m为样本数，θ为权值。
如果是回归问题，损失函数通常是均方误差(MSE):
因此，回归问题的损失函数对权值的梯度（导数）为：
如果是二分类问题，损失函数通常是logit loss：
其中，
表示阶跃函数Sigmoid。
因此，权值的梯度（导数）为：
常数项、一次项、交叉项的导数分别为：
7. FM算法的Python实现
FM算法的Python实现流程图如下：
要注意以下四点：
1. 初始化参数包括偏置项权值w0.一次项权值w和交叉项辅助向量的初始化；
2. 定义FM算法；
3. 定义损失函数梯度；
4. 利用梯度下降更新参数。
以下代码片段是上述四点的描述，其中loss不是二分类的损失loss，而是分类loss梯度的一部分：
loss = self.sigmoid(classLabels[x] * p[0, 0]) -1
事实上，二分类的损失loss梯度可表示为：
gradient = (self.sigmoid(classLabels[x] * p[0, 0]) -1)*classLabels[x]*p_derivative
其中 p_derivative 代表常数项、一次项、交叉项的导数(详见本文第6节)。
FM算法代码片段
1 #初始化参数
2 w = zeros((n, 1)) #其中n是特征的个数
3 w_0 =0.4 v = normalvariate(0, 0.2) *ones((n, k))5 for it in range(self.iter): #迭代次数
6 #优化每个样本
7 for x inrange(m):8 #这里注意一个数学知识:通常会有相应点积的地方sum，通常没有相应位置积的地方，详见矩阵操作规则，本处计算逻辑如下：http://blog.csdn.net/google19890102/article/details/45532745
9 #xi·vi,xi与vi的矩阵点积
10 inter_1 = dataMatrix[x] *v11 #xi与xi的对应位置乘积 与 xi^2与vi^2相应位置的乘积 的点积
12 inter_2 = multiply(dataMatrix[x], dataMatrix[x]) * multiply(v, v) #multiply相应元素相乘
13 #完成交叉项，xi*vi*xi*vi - xi^2*vi^2
14 interaction = sum(multiply(inter_1, inter_1) - inter_2) / 2.15 #计算预测输出
16 p = w_0 dataMatrix[x] * w interaction17 print('classLabels[x]:',classLabels[x])18 print(预测输出p:', p)19 #计算sigmoid(y*pred_y)-准确地说不loss，原作者这边理解的有问题，只是作为更新w的中间参数，这边算出来的是越大越好，而下面却用了梯度下降而不是梯度上升的算法在
20 loss = self.sigmoid(classLabels[x] * p[0, 0]) - 1
21 if loss >= -1:22 loss_res = 正方向'
23 else:24 loss_res = '反方向'
25 #更新参数
26 w_0 = w_0 - self.alpha * loss *classLabels[x]27 for i inrange(n):28 if dataMatrix[x, i] !=0:29 w[i, 0] = w[i, 0] - self.alpha * loss * classLabels[x] *dataMatrix[x, i]30 for j inrange(k):31 v[i, j] = v[i, j] - self.alpha * loss * classLabels[x] *(32 dataMatrix[x, i] * inter_1[0, j] - v[i, j] * dataMatrix[x, i] * dataMatrix[x, i])
FM算法完整实现
1 #-*- coding: utf-8 -*-
2
3 from __future__ importdivision4 from math importexp5 from numpy import *
6 from random import normalvariate #正态分布
7 from sklearn importpreprocessing8 importnumpy as np9
10 '''
11 data : 数据的路径12 feature_potenital : 潜在分解维度数13 alpha ： 学习速率14 iter ： 迭代次数15 _w,_w_0,_v ： 拆分子矩阵的weight16 with_col : 是否带有columns_name17 first_col : 首列有价值的feature的index18 '''
19
20
21 classfm(object):22 def __init__(self):23 self.data =None24 self.feature_potential =None25 self.alpha =None26 self.iter =None27 self._w =None28 self._w_0 =None29 self.v =None30 self.with_col =None31 self.first_col =None32
33 defmin_max(self, data):34 self.data =data35 min_max_scaler =preprocessing.MinMaxScaler()36 returnmin_max_scaler.fit_transform(self.data)37
38 def loadDataSet(self, data, with_col=True, first_col=2):39 #我就是闲的蛋疼，明明pd.read_table()可以直接度，非要搞这样的，显得代码很长，小数据下完全可以直接读嘛，唉～
40 self.first_col =first_col41 dataMat =[]42 labelMat =[]43 fr =open(data)44 self.with_col =with_col45 ifself.with_col:46 N =047 for line infr.readlines():48 #N=1时干掉列表名
49 if N >0:50 currLine =line.strip().split()51 lineArr =[]52 featureNum =len(currLine)53 for i inrange(self.first_col, featureNum):54 lineArr.append(float(currLine[i]))55 dataMat.append(lineArr)56 labelMat.append(float(currLine[1]) * 2 - 1)57 N = N + 1
58 else:59 for line infr.readlines():60 currLine =line.strip().split()61 lineArr =[]62 featureNum =len(currLine)63 for i in range(2, featureNum):64 lineArr.append(float(currLine[i]))65 dataMat.append(lineArr)66 labelMat.append(float(currLine[1]) * 2 - 1)67 returnmat(self.min_max(dataMat)), labelMat68
69 defsigmoid(self, inx):70 #return 1.0/(1+exp(min(max(-inx,-10),10)))
71 return 1.0 / (1 + exp(-inx))72
73 #得到对应的特征weight的矩阵
74 def fit(self, data, feature_potential=8, alpha=0.01, iter=100):75 #alpha是学习速率
76 self.alpha =alpha77 self.feature_potential =feature_potential78 self.iter =iter79 #dataMatrix用的是mat, classLabels是列表
80 dataMatrix, classLabels =self.loadDataSet(data)81 print('dataMatrix:',dataMatrix.shape)82 print('classLabels:',classLabels)83 k =self.feature_potential84 m, n =shape(dataMatrix)85 #初始化参数
86 w = zeros((n, 1)) #其中n是特征的个数
87 w_0 =0.88 v = normalvariate(0, 0.2) *ones((n, k))89 for it in range(self.iter): #迭代次数
90 #对每一个样本，优化
91 for x inrange(m):92 #这边注意一个数学知识：对应点积的地方通常会有sum，对应位置积的地方通常都没有，详细参见矩阵运算规则，本处计算逻辑在：http://blog.csdn.net/google19890102/article/details/45532745
93 #xi·vi,xi与vi的矩阵点积
94 inter_1 = dataMatrix[x] *v95 #xi与xi的对应位置乘积 与 xi^2与vi^2对应位置的乘积 的点积
96 inter_2 = multiply(dataMatrix[x], dataMatrix[x]) * multiply(v, v) #multiply对应元素相乘
97 #完成交叉项,xi*vi*xi*vi - xi^2*vi^2
98 interaction = sum(multiply(inter_1, inter_1) - inter_2) / 2.99 #计算预测的输出
100 p = w_0 + dataMatrix[x] * w +interaction101 print('classLabels[x]:',classLabels[x])102 print('预测的输出p:', p)103 #计算sigmoid(y*pred_y)-1
104 loss = self.sigmoid(classLabels[x] * p[0, 0]) - 1
105 if loss >= -1:106 loss_res = '正方向'
107 else:108 loss_res = '反方向'
109 #更新参数
110 w_0 = w_0 - self.alpha * loss *classLabels[x]111 for i inrange(n):112 if dataMatrix[x, i] !=0:113 w[i, 0] = w[i, 0] - self.alpha * loss * classLabels[x] *dataMatrix[x, i]114 for j inrange(k):115 v[i, j] = v[i, j] - self.alpha * loss * classLabels[x] *(116 dataMatrix[x, i] * inter_1[0, j] - v[i, j] * dataMatrix[x, i] *dataMatrix[x, i])117 print('the no %s times, the loss arrach %s' %(it, loss_res))118 self._w_0, self._w, self._v =w_0, w, v119
120 defpredict(self, X):121 if (self._w_0 == None) or (self._w == None).any() or (self._v ==None).any():122 raise NotFittedError("Estimator not fitted, call `fit` first")123 #类型检查
124 ifisinstance(X, np.ndarray):125 pass
126 else:127 try:128 X =np.array(X)129 except:130 raise TypeError("numpy.ndarray required for X")131 w_0 =self._w_0132 w =self._w133 v =self._v134 m, n =shape(X)135 result =[]136 for x inrange(m):137 inter_1 = mat(X[x]) *v138 inter_2 = mat(multiply(X[x], X[x])) * multiply(v, v) #multiply对应元素相乘
139 #完成交叉项
140 interaction = sum(multiply(inter_1, inter_1) - inter_2) / 2.141 p = w_0 + X[x] * w + interaction #计算预测的输出
142 pre =self.sigmoid(p[0, 0])143 result.append(pre)144 returnresult145
146 defgetAccuracy(self, data):147 dataMatrix, classLabels =self.loadDataSet(data)148 w_0 =self._w_0149 w =self._w150 v =self._v151 m, n =shape(dataMatrix)152 allItem =0153 error =0154 result =[]155 for x inrange(m):156 allItem += 1
157 inter_1 = dataMatrix[x] *v158 inter_2 = multiply(dataMatrix[x], dataMatrix[x]) * multiply(v, v) #multiply对应元素相乘
159 #完成交叉项
160 interaction = sum(multiply(inter_1, inter_1) - inter_2) / 2.161 p = w_0 + dataMatrix[x] * w + interaction #计算预测的输出
162 pre =self.sigmoid(p[0, 0])163 result.append(pre)164 if pre < 0.5 and classLabels[x] == 1.0:165 error += 1
166 elif pre >= 0.5 and classLabels[x] == -1.0:167 error += 1
168 else:169 continue
170 #print(result)
171 value = 1 - float(error) /allItem172 returnvalue173
174
175 classNotFittedError(Exception):176 """
177 Exception class to raise if estimator is used before fitting178 """
179 pass
180
181
182 if __name__ == '__main__':183 fm()
=