四、序列最小優化算法（smo算法）
1、smo算法基本思想
支持向量機的學習問題可以形式化為求解凸二次規划問題。 這樣的凸二次規划問題具有全局最優解， 並且有許多最優化算法可以用於這一問題的求解。 但是當訓練樣本容量很大時， 這些算法往往變得非常低效， 以致無法使用。 所以，如何高效地實現支持向量機學習就成為一個重要的問題。
目前人們已提出許多快速實現算法。其中最具代表的就是序列最小最優化算法（sequential minimaloptimization，smo）。
smo算法要解如下凸二次規划的對偶問題：
算法是一種啟發式算法， 其基本思路是：如果所有變量的解都滿足此最優化問題的KKT條件，那么這個最優化問題的解就得到了。因為KKT條件是該最優化問題的充分必要條件。否則,選擇兩個變量，固定其他變量，針對這兩個變量構建一個二次規划問題。這個二次規划問題關於這兩個變量的解應該更接近原始二次規划問題的解，因為這會使得原始二次規划問題的目標函數值變得更小。重要的是，這時子問題可以通過解析方法求解，這樣就可以大大提高整個算法的計算速度。子問題有兩個變量,一個是違反KKT條件最嚴重的那一個，另一個由約束條件自動確定。如此，smo算法將原問題不斷分解為子問題並對子問題求解，進而達到求解原問題的目的。具體如下：
選擇兩個變量 alpha1 , alpha2，其它變量 alphai , i = 3 , 4 , .. m是固定的。
於是，smo的最優化問題的子問題可以寫成：
2、SMO算法目標函數的優化
為了求解上面含有這兩個變量的目標優化問題，我們首先分析約束條件，所有的 alpha1 , alpha2 都要滿足約束條件，然后在約束條件下求最小。
3、SMO算法代碼實現
（1）數據集介紹
數據集共包括400個樣本，有年齡跟工資兩個特征，標簽為是否購買房產，1表示購買，0表示未購買。我們取300個樣本作為訓練集，100個樣本作為測試集。部分數據如下：
（2）獲取數據代碼
#獲取數據 import pandas as pd dataset = pd.read_csv('./step3/Social_Network_Ads.csv') x = dataset.iloc[:, [2, 3]].values x = x.astype(float) y = dataset.iloc[:, 4].values #將0替換為-1 for i in range(len(y)): if y[i]==0: y[i]=-1 from sklearn.model_selection import train_test_split train_data, test_data, train_label, test_label = train_test_split(x, y, test_size = 0.25, random_state = 61) #特征標准化 from sklearn.preprocessing import StandardScaler sc = StandardScaler() train_data = sc.fit_transform(train_data) test_data = sc.transform(test_data)
（3）整體代碼
#encoding=utf8 import numpy as np class smo: def __init__(self, max_iter=100, kernel='linear'): ''' input:max_iter(int):最大訓練輪數 kernel(str):核函數，等於'linear'表示線性，等於'poly'表示多項式 ''' self.max_iter = max_iter self._kernel = kernel #初始化模型 def init_args(self, features, labels): self.m, self.n = features.shape self.X = features self.Y = labels self.b = 0.0 # 將Ei保存在一個列表里 self.alpha = np.ones(self.m) self.E = [self._E(i) for i in range(self.m)] # 錯誤懲罰參數 self.C = 1.0 #kkt條件 def _KKT(self, i): y_g = self._g(i)*self.Y[i] if self.alpha[i] == 0: return y_g >= 1 elif 0 < self.alpha[i] < self.C: return y_g == 1 else: return y_g <= 1 # g(x)預測值，輸入xi（X[i]） def _g(self, i): r = self.b for j in range(self.m): r += self.alpha[j]*self.Y[j]*self.kernel(self.X[i], self.X[j]) return r # 核函數,多項式添加二次項即可 def kernel(self, x1, x2): if self._kernel == 'linear': return sum([x1[k]*x2[k] for k in range(self.n)]) elif self._kernel == 'poly': return (sum([x1[k]*x2[k] for k in range(self.n)]) + 1)**2 return 0 # E（x）為g(x)對輸入x的預測值和y的差 def _E(self, i): return self._g(i) - self.Y[i] #初始alpha def _init_alpha(self): # 外層循環首先遍歷所有滿足0<a<C的樣本點，檢驗是否滿足KKT index_list = [i for i in range(self.m) if 0 < self.alpha[i] < self.C] # 否則遍歷整個訓練集 non_satisfy_list = [i for i in range(self.m) if i not in index_list] index_list.extend(non_satisfy_list) for i in index_list: if self._KKT(i): continue E1 = self.E[i] # 如果E2是+，選擇最小的；如果E2是負的，選擇最大的 if E1 >= 0: j = min(range(self.m), key=lambda x: self.E[x]) else: j = max(range(self.m), key=lambda x: self.E[x]) return i, j #選擇alpha參數 def _compare(self, _alpha, L, H): if _alpha > H: return H elif _alpha < L: return L else: return _alpha #訓練 def fit(self, features, labels): ''' input:features(ndarray):特征 label(ndarray):標簽 ''' self.init_args(features, labels) for t in range(self.max_iter): i1, i2 = self._init_alpha() # 邊界 if self.Y[i1] == self.Y[i2]: L = max(0, self.alpha[i1]+self.alpha[i2]-self.C) H = min(self.C, self.alpha[i1]+self.alpha[i2]) else: L = max(0, self.alpha[i2]-self.alpha[i1]) H = min(self.C, self.C+self.alpha[i2]-self.alpha[i1]) E1 = self.E[i1] E2 = self.E[i2] # eta=K11+K22-2K12 eta = self.kernel(self.X[i1], self.X[i1]) + self.kernel(self.X[i2], self.X[i2]) - 2*self.kernel(self.X[i1], self.X[i2]) if eta <= 0: continue alpha2_new_unc = self.alpha[i2] + self.Y[i2] * (E2 - E1) / eta alpha2_new = self._compare(alpha2_new_unc, L, H) alpha1_new = self.alpha[i1] + self.Y[i1] * self.Y[i2] * (self.alpha[i2] - alpha2_new) b1_new = -E1 - self.Y[i1] * self.kernel(self.X[i1], self.X[i1]) * (alpha1_new-self.alpha[i1]) - self.Y[i2] * self.kernel(self.X[i2], self.X[i1]) * (alpha2_new-self.alpha[i2])+ self.b b2_new = -E2 - self.Y[i1] * self.kernel(self.X[i1], self.X[i2]) * (alpha1_new-self.alpha[i1]) - self.Y[i2] * self.kernel(self.X[i2], self.X[i2]) * (alpha2_new-self.alpha[i2])+ self.b if 0 < alpha1_new < self.C: b_new = b1_new elif 0 < alpha2_new < self.C: b_new = b2_new else: # 選擇中點 b_new = (b1_new + b2_new) / 2 # 更新參數 self.alpha[i1] = alpha1_new self.alpha[i2] = alpha2_new self.b = b_new self.E[i1] = self._E(i1) self.E[i2] = self._E(i2) def predict(self, data): ''' input:data(ndarray):單個樣本 output:預測為正樣本返回+1，負樣本返回-1 ''' r = self.b for i in range(self.m): r += self.alpha[i] * self.Y[i] * self.kernel(data, self.X[i]) return 1 if r > 0 else -1