FM算法解析及Python實現
阿新 • • 發佈：2018-09-08
line 類型檢查 lse pytho spl ont 核心 ict 作者 7. FM算法的Python實現
1. 什麽是FM？
FM即Factor Machine，因子分解機。
2. 為什麽需要FM？
1、特征組合是許多機器學習建模過程中遇到的問題，如果對特征直接建模，很有可能會忽略掉特征與特征之間的關聯信息，因此，可以通過構建新的交叉特征這一特征組合方式提高模型的效果。
2、高維的稀疏矩陣是實際工程中常見的問題，並直接會導致計算量過大，特征權值更新緩慢。試想一個10000*100的表，每一列都有8種元素，經過one-hot獨熱編碼之後，會產生一個10000*800的表。因此表中每行元素只有100個值為1，700個值為0。
而FM的優勢就在於對這兩方面問題的處理。首先是特征組合，通過對兩兩特征組合，引入交叉項特征，提高模型得分；其次是高維災難，通過引入輔助向量（對參數矩陣進行矩陣分解），完成對特征的參數估計。
3. FM用在哪？
我們已經知道了FM可以解決特征組合以及高維稀疏矩陣問題，而實際業務場景中，電商、豆瓣等推薦系統的場景是使用最廣的領域，打個比方，小王只在豆瓣上瀏覽過20部電影，而豆瓣上面有20000部電影，如果構建一個基於小王的電影矩陣，毫無疑問，裏面將有199980個元素全為0。而類似於這樣的問題就可以通過FM來解決。
4. FM長什麽樣？
在展示FM算法前，我們先回顧一下最常見的線性表達式：
其中w0 為初始權值，或者理解為偏置項，wi 為每個特征xi 對應的權值。可以看到，這種線性表達式只描述了每個特征與輸出的關系。
FM的表達式如下，可觀察到，只是在線性表達式後面加入了新的交叉項特征及對應的權值。
5. FM交叉項的展開
5.1 尋找交叉項
FM表達式的求解核心在於對交叉項的求解。下面是很多人用來求解交叉項的展開式，對於第一次接觸FM算法的人來說可能會有疑惑，不知道公式怎麽展開的，接下來筆者會手動推導一遍。
設有3個變量（特征）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 交叉項展開式
上面的例子是對3個特征做的交叉項推導，因此對具有n個特征，FM的交叉項公式就可推廣為：
我們還可以進一步分解：
所以FM算法的交叉項最終可展開為：
6. 權值求解
利用梯度下降法，通過求損失函數對特征（輸入項）的導數計算出梯度，從而更新權值。設m為樣本個數，θ為權值。
如果是回歸問題，損失函數一般是均方誤差（MSE）:
所以回歸問題的損失函數對權值的梯度（導數）為：
如果是二分類問題，損失函數一般是logit loss：
其中，
所以分類問題的損失函數對權值的梯度（導數）為：
相應的，對於常數項、一次項、交叉項的導數分別為:
7. FM算法的Python實現
FM算法的Python實現流程圖如下：
我們需要註意以下四點：
1. 初始化參數，包括對偏置項權值w0、一次項權值w以及交叉項輔助向v量的初始化；
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 in range(m): 8 # 這邊註意一個數學知識：對應點積的地方通常會有sum，對應位置積的地方通常都沒有，詳細參見矩陣運算規則，本處計算邏輯在：http://blog.csdn.net/google19890102/article/details/45532745 9 # xi·vi,xi與vi的矩陣點積 10 inter_1 = dataMatrix[x] * v 11 # 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 + interaction 17 print(‘classLabels[x]:‘,classLabels[x]) 18 print(‘預測的輸出p:‘, p) 19 # 計算sigmoid(y*pred_y)-1準確的說不是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 in range(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 in range(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__ import division 4 from math import exp 5 from numpy import * 6 from random import normalvariate # 正態分布 7 from sklearn import preprocessing 8 import numpy as np 9 10 ‘‘‘ 11 data : 數據的路徑 12 feature_potenital : 潛在分解維度數 13 alpha ： 學習速率 14 iter ： 叠代次數 15 _w,_w_0,_v ： 拆分子矩陣的weight 16 with_col : 是否帶有columns_name 17 first_col : 首列有價值的feature的index 18 ‘‘‘ 19 20 21 class fm(object): 22 def __init__(self): 23 self.data = None 24 self.feature_potential = None 25 self.alpha = None 26 self.iter = None 27 self._w = None 28 self._w_0 = None 29 self.v = None 30 self.with_col = None 31 self.first_col = None 32 33 def min_max(self, data): 34 self.data = data 35 min_max_scaler = preprocessing.MinMaxScaler() 36 return min_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_col 41 dataMat = [] 42 labelMat = [] 43 fr = open(data) 44 self.with_col = with_col 45 if self.with_col: 46 N = 0 47 for line in fr.readlines(): 48 # N=1時幹掉列表名 49 if N > 0: 50 currLine = line.strip().split() 51 lineArr = [] 52 featureNum = len(currLine) 53 for i in range(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 in fr.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 return mat(self.min_max(dataMat)), labelMat 68 69 def sigmoid(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 = alpha 77 self.feature_potential = feature_potential 78 self.iter = iter 79 # dataMatrix用的是mat, classLabels是列表 80 dataMatrix, classLabels = self.loadDataSet(data) 81 print(‘dataMatrix:‘,dataMatrix.shape) 82 print(‘classLabels:‘,classLabels) 83 k = self.feature_potential 84 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 in range(m): 92 # 這邊註意一個數學知識：對應點積的地方通常會有sum，對應位置積的地方通常都沒有，詳細參見矩陣運算規則，本處計算邏輯在：http://blog.csdn.net/google19890102/article/details/45532745 93 # xi·vi,xi與vi的矩陣點積 94 inter_1 = dataMatrix[x] * v 95 # 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 + interaction 101 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 in range(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 in range(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, v 119 120 def predict(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 if isinstance(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_0 132 w = self._w 133 v = self._v 134 m, n = shape(X) 135 result = [] 136 for x in range(m): 137 inter_1 = mat(X[x]) * v 138 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 return result 145 146 def getAccuracy(self, data): 147 dataMatrix, classLabels = self.loadDataSet(data) 148 w_0 = self._w_0 149 w = self._w 150 v = self._v 151 m, n = shape(dataMatrix) 152 allItem = 0 153 error = 0 154 result = [] 155 for x in range(m): 156 allItem += 1 157 inter_1 = dataMatrix[x] * v 158 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) / allItem 172 return value 173 174 175 class NotFittedError(Exception): 176 """ 177 Exception class to raise if estimator is used before fitting 178 """ 179 pass 180 181 182 if __name__ == ‘__main__‘: 183 fm()
FM算法解析及Python實現