線性支持向量機
SVM（Support Vector Machine）是數據挖掘中常用的分類算法。事實上，在2012年深度學習算法提出之前，SVM一直被認為是機器學習領域近幾十年來最成功的算法。SVM的難度相較於其他分類算法要高一些兒，特別是涉及到的數學知識比較多，不過不要害怕，本文盡量把每個知識點講清楚，力求更多的人能看懂SVM分類算法。現在，請深呼吸幾次，給大腦充足的氧氣，我們即將一起去探究SVM的點點滴滴。
我們根據數據集的分布情況，可將數據集分成兩種：線性可分與線性不可分。所謂的線性可分，是指我們可以找到一個線性方程（2維情況下是一條直線，3維情況下是一個平面），將數據分隔開來。下圖展示了2維數據集線性可分與線性不可分的情況。
如圖所示，在左圖中可以找到一條直線，將兩類數據分隔開，而右圖則找不到這樣一條直線，因此我們稱其是線性不可分的。
本文將着重討論SVM在線性可分數據集上的應用，對於線性不可分的數據集，留待在后續的文章中詳細討論。
SVM的模型
對於給定的訓練數據集，我們可能有多條不同的線性方程可以將數據分隔開，如上圖所示。SVM分類算法的目標，就是找到最優划分的那條線性方程。因此，線性方程的數學模型，即是SVM的分類模型，其可以表示為：
其中，X是數據集的特征向量，W是權重向量，b是偏置項。
我們說到，要尋找最優划分的線性方程，那么，該如何度量線性方程划分結果的優劣呢？這里我們引入一個概念：邊緣（Margin）
有時候，我們也會稱線性方程為超平面（hyperplane，因為3維空間里的線性方程就是指一個平面，而對於高維空間，則是超平面）。對於一個超平面，我們做一個這樣的操作：將超平面向上平移，直到藍色類的第一個點落在超平面上時停止（上方虛線）；將超平面向下平移，直到紅色類的第一個點落在超平面上時停止（下方虛線）。兩條虛線間的距離即為該超平面的邊緣。
那么這個距離得怎么算呢？我們知道，將超平面平移時，權重向量W是不會改變的，僅僅是改變偏置項b。（平移時直線的斜率不會變），因此，兩條虛線的方程可以這樣表示：(我們把偏置項的變化量移到了方程右邊)
我們不妨用（X，y）表示一條數據記錄，其中X是記錄的特征向量，y是記錄的類標號，我們把藍色類的類標號記為1，把紅色類的類標號記為-1，則有：
為了之后計算方便，我們可以調整W和b的取值，對方程做縮放變換，使得k=1，k'=-1，則上式可改寫為：
我們可以將上面兩個不等式概括為更加簡潔的形式：
假設超平面 W * X + b = 0上有p1，p2兩點，則有 W * Xp1 + bp1 = 0，W * Xp2 + bp2 = 0，兩式相減有 W * ( Xp1 - Xp2 ) = 0。其中 Xp1 - Xp2 是該超平面上的向量，說明權重向量W和該超平面垂直。
假設上下虛線上分別有p3，p4兩點，則對於p3點，有 W * Xp3 + bp3 = 1；對於p4點，有 W * Xp4 + bp4 = -1，兩式相減有 W * ( Xp3 - Xp4 ) = 2。Xp3 - Xp4 是p3，p4兩點間的向量，其在W方向上的投影長度即是兩條虛線間的距離d，所以有：
即：
其中，||W||是向量W的模長。至此，我們已經知道了如何計算一個超平面邊緣的大小，SVM的目標，就是找到邊緣最大的那個超平面，我們稱之為最大邊緣超平面（Maximum Margin Hyperplane）。也就是說，我們要優化距離函數d，使得取得最大值，此時對應的W0，就是MMH的權重向量。
為什么要尋找邊緣最大的超平面呢？因為具有較大邊緣的超平面具有更小的泛化誤差（也就是將應用在其它樣本上的誤差），邊緣較小的超平面，對訓練數據的過分擬合程度會更高，從而在未知樣本上的泛化能力較差。
目標函數
通過上面的討論，我們知道了需要最大化邊緣，最大化邊緣等價於最小化下面的函數：
像 f(W) 這樣我們要優化的目標，在機器學習領域，有一個專門的術語：目標函數
同時，不要忘了還有式1的不等式約束。總結來說，SVM的學習任務為以下被約束的優化問題：
可以看到，目標函數 f(W) 是二次的，而不等式約束在參數W和b上面是線性的，這在數學上屬於一個凸優化問題，對於凸優化問題，我們可以通過標准的拉格朗日乘子（Lagrange Multiplier）方法求解。下面簡要介紹一下求解這個優化問題的主要思想。
首先，必須改寫目標函數，將施加在解上的約束包括進來，新目標函數稱為該優化問題的拉格朗日函數：
其中，N為訓練數據集的大小，λi稱為拉格朗日乘子。拉格朗日函數的第一項為原函數，第二項則捕獲了不等式約束。為了理解改寫原目標函數的必要性，考慮式2，容易證明當 W=0（零向量）時函數取得最小值，然而，這樣的解違背了約束條件，因為b此時沒有解。事實上，如果 W，b 的解違反了約束，即 yi ( W Xi + b ) - 1 < 0 ，假定λi≥ 0,則任何不可行的解僅僅是增加了拉格朗日函數的值，當我們求得函數的最小值時，即可保證此時的解W0，b0是滿足約束條件的。
為了最小化拉格朗日函數，需要對L求關於W和b的偏導，並令它們等於0：
因為拉格朗日乘子λi是未知的，所以現在我們仍然求解不了W和b的值。
如果我們能將式1中的不等式約束轉換為等式約束，則我們可以從該等式約束得到的N的方程，加上式5,6（共N+2個方程），從而求得W，b，λi（共N+2個參數）的可行解。幸運的是，這種轉換在數學上是可行的，稱為Karuch-Kuhn-Tucher（KKT）轉換。式1的不等式約束等價於以下的等式約束：
乍一看，拉格朗日乘子的數量好像和訓練數據集的規模一樣大。不過我們冷靜想一下，絕大部分的點，滿足 y ( W * X + b ) - 1 > 0，因為它們都在虛線的一側，並沒有落在虛線上，因此，為了滿足式6的等式約束，這些點對應的拉格朗日乘子必須等於0，我們所需要求的，只是那些落在虛線上的點對應的λ。這些點被稱為支持向量，因為式5中W的取值只受那些λ≠0的點的影響，這也是這個算法被稱為支持向量機的原因。
現在，我們已經過濾掉了大量為0的參數，但是對於式4的優化問題求解仍然是一項十分棘手的任務，因為方程中同時涉及了W，b，λ三類參數。我們可以將拉格朗日函數轉化成對應的對偶函數（duel function）來進一步簡化該問題。在對偶函數中，參數只有拉格朗日乘子λ。
對偶定理陳述如下：
如果原問題有最優解，則其對偶問題也有最優解，並且相應的最優解是相同的。
我們可以將式5,6代入到式4中求得原拉格朗日函數的對偶拉格朗日函數。
首先，將式4展開：
代入式5,6：
因為，λi=0時，對結果沒有影響，所以我們可以把數據的范圍從整個訓練數據集縮小到支持向量的集合，用S表示支持向量的集合。
此時有約束條件如下：
至此，通過層層的轉化，我們已經找到了最終的目標函數：對偶拉格朗日函數LD。可以看到，LD僅涉及到拉格朗日乘子和支持向量，而原拉格朗日函數L除了拉格朗日乘子外還涉及超平面的參數W，b。奇妙的是，這兩個優化問題的解是等價的。同時，需要注意的是，LD中的二次項前面的符號是負的，這說明原來涉及L的最小化問題已經變成了涉及LD的最大化問題。
SMO算法
首次，對堅持到這里的同學們點個贊，經過各種數學公式的層層洗禮之后還未放棄實屬不易。告訴大家一個好消息，SMO算法是求解SVM分類器的最后一步了，現在，大家可以站起來活動活動筋骨，再次深吸幾口氣，讓大腦休息一會兒再一鼓作氣攻克最后一步吧！
那么，如何求解式7呢？這在數學上屬於二次規划問題，可以使用通用的二次規划算法來求解，然而，該問題的復雜程度正比於訓練數據集的規模，這會在實際的任務中造成很大的開銷，為了避開這個障礙，人們通過利用問題本身的特性，提出了很多高效的算法，SMO（Sequential Minimal Optimization）就是其中一個著名的代表。
SMO的基本思想是選取一對變量λi，λj，並固定其他參數，由於存在式8的約束，所以λi可以用λj和其他固定的拉格朗日乘子導出。這樣LD就變成了一個只關於單變量λi的二次規划問題，僅有的約束是λi > 0。我們可以將參數全都初始化為0，然后不斷執行如下步驟直至收斂。
1.選取一對需要更新的變量λi，λj
2.固定除λi，λj外的參數，通過式7更新λi，λj
我們假定更新前的λ記為λold，更新后的λ記為λnew 。則式8的約束可以改寫為：(因為得保證式8一直成立)
為了求式7的導數時容易理解，我們讓 i=1，j=2，則有：(其中Constant是和λ1，λ2無關的常數項，並不影響后面導數的計算)
將式9代入式10並令其關於λ1的導數為0可以得到未截斷（unclipped）的最優解：(其中E是誤差項，E = W * X + b - y)
對於未截斷的最優解λ1unclipped，通過二次規划的約束λ1≥0，λ2≥0便可截斷λ1unclipped，得到λ1new的值。然后通過式9計算出λ2new，至此，我們便完成了一組拉格朗日乘子的更新。重復以上計算過程，直到各拉格朗日乘子不再變化。
當計算出支持向量對應的拉格朗日乘子后，我們就可以通過式5計算出權重向量W的值了。
說了這么多，一直在說W如何計算，那么偏置項b呢？我們可以通過兩條虛線的方程：W*X+b=±1 來計算b的值，可以知道，兩條虛線的偏置項b計算出來的結果是不一樣的，我們不妨記為b1，b2。常用的做法是，取它們的平均值作為b的最終取值，也即是：b=(b1+b2)/2。
到這里，線性SVM算是基本講完了，只剩下一個小問題：模型訓練誤差為0，也即是說模型對訓練數據是完全擬合的。但有時候，可能有一定訓練誤差的模型效果要更好一些。在接下來的文章中，我們將一起討論如何訓練允許一定訓練誤差的SVM。好了，今天就先到這吧，快去讓大腦放松一下。
代碼（Python3.7）：
線性SVM的實現類代碼：
1 import numpy as np 2 3 4 class LinearSVM(object): 5 """ 6 線性SVM的實現類 7 """ 8 def __init__(self, dataset_size, vector_size): 9 """ 10 初始化函數 11 :param dataset_size:數據集規模，用於初始化 ‘拉格朗日乘子’ 的數量 12 :param vector_size:特征向量的維度，用於初始化 ‘權重向量’ 的維度 13 """ 14 self.__multipliers = np.zeros(dataset_size, np.float_) 15 self.weight_vec = np.zeros(vector_size, np.float_) 16 self.bias = 0 17 18 def train(self, dataset, iteration_num): 19 """ 20 訓練函數 21 :param dataset:數據集，每條數據的形式為（X，y），其中X是特征向量，y是類標號 22 :param iteration_num: 23 """ 24 dataset = np.array(dataset) 25 for k in range(iteration_num): 26 self.__update(dataset, k) 27 28 def __update(self, dataset, k): 29 """ 30 更新函數 31 :param dataset: 32 :param k: 33 """ 34 for i in range(dataset.__len__() // 2): 35 j = (dataset.__len__() // 2 + i + k) % dataset.__len__() 36 record_i = dataset[i] 37 record_j = dataset[j] 38 self.__sequential_minimal_optimization(dataset, record_i, record_j, i, j) 39 self.__update_weight_vec(dataset) 40 self.__update_bias(dataset) 41 42 def __sequential_minimal_optimization(self, dataset, record_i, record_j, i, j): 43 """ 44 SMO函數，每次選取兩條記錄，更新對應的‘拉格朗日乘子’ 45 :param dataset: 46 :param record_i:記錄i 47 :param record_j:記錄j 48 :param i: 49 :param j: 50 """ 51 label_i = record_i[-1] 52 vector_i = np.array(record_i[0]) 53 label_j = record_j[-1] 54 vector_j = np.array(record_j[0]) 55 56 # 計算出截斷前的記錄i的‘拉格朗日乘子’unclipped_i 57 error_i = np.dot(self.weight_vec, vector_i) + self.bias - label_i 58 error_j = np.dot(self.weight_vec, vector_j) + self.bias - label_j 59 eta = np.dot(vector_i - vector_j, vector_i - vector_j) 60 unclipped_i = self.__multipliers[i] + label_i * (error_j - error_i) / eta 61 62 # 截斷記錄i的`拉格朗日乘子`並計算記錄j的`拉格朗日乘子` 63 constant = -self.__calculate_constant(dataset, i, j) 64 multiplier = self.__quadratic_programming(unclipped_i, label_i, label_j, i, j) 65 if multiplier >= 0: 66 self.__multipliers[i] = multiplier 67 self.__multipliers[j] = (constant - multiplier * label_i) * label_j 68 69 def __update_bias(self, dataset): 70 """ 71 計算偏置項bias，使用平均值作為最終結果 72 :param dataset: 73 """ 74 sum_bias = 0 75 count = 0 76 for k in range(self.__multipliers.__len__()): 77 if self.__multipliers[k] != 0: 78 label = dataset[k][-1] 79 vector = np.array(dataset[k][0]) 80 sum_bias += 1 / label - np.dot(self.weight_vec, vector) 81 count += 1 82 if count == 0: 83 self.bias = 0 84 else: 85 self.bias = sum_bias / count 86 87 def __update_weight_vec(self, dataset): 88 """ 89 計算權重向量 90 :param dataset: 91 """ 92 weight_vector = np.zeros(dataset[0][0].__len__()) 93 for k in range(dataset.__len__()): 94 label = dataset[k][-1] 95 vector = np.array(dataset[k][0]) 96 weight_vector += self.__multipliers[k] * label * vector 97 self.weight_vec = weight_vector 98 99 def __calculate_constant(self, dataset, i, j): 100 label_i = dataset[i][-1] 101 label_j = dataset[j][-1] 102 dataset[i][-1] = 0 103 dataset[j][-1] = 0 104 sum_constant = 0 105 for k in range(dataset.__len__()): 106 label = dataset[k][-1] 107 sum_constant += self.__multipliers[k] * label 108 dataset[i][-1] = label_i 109 dataset[j][-1] = label_j 110 return sum_constant 111 112 def __quadratic_programming(self, unclipped_i, label_i, label_j, i, j): 113 """ 114 二次規划，截斷`拉格朗日乘子` 115 :param unclipped_i: 116 :param label_i: 117 :param label_j: 118 :return: 119 """ 120 multiplier = -1 121 if label_i * label_j == 1: 122 boundary = self.__multipliers[i] + self.__multipliers[j] 123 if boundary >= 0: 124 if unclipped_i <= 0: 125 multiplier = 0 126 elif unclipped_i < boundary: 127 multiplier = unclipped_i 128 else: 129 multiplier = boundary 130 else: 131 boundary = max(0, self.__multipliers[i] - self.__multipliers[j]) 132 if unclipped_i <= boundary: 133 multiplier = boundary 134 else: 135 multiplier = unclipped_i 136 return multiplier 137 138 def predict(self, vector): 139 result = np.dot(self.weight_vec, np.array(vector)) + self.bias 140 if result >= 0: 141 return 1 142 else: 143 return -1 144 145 def __str__(self): 146 return "multipliers:" + self.__multipliers.__str__() + '\n' + \ 147 "weight_vector:" + self.weight_vec.__str__() + '\n' + \ 148 "bias:" + self.bias.__str__()
數據集和測試代碼：
1 from SVM import LinearSVM 2 import numpy as np 3 from matplotlib import pyplot as plt 4 5 dataset = [ 6 [[0.3858, 0.4687], 1], 7 [[0.4871, 0.611], -1], 8 [[0.9218, 0.4103], -1], 9 [[0.7382, 0.8936], -1], 10 [[0.1763, 0.0579], 1], 11 [[0.4057, 0.3529], 1], 12 [[0.9355, 0.8132], -1], 13 [[0.2146, 0.0099], 1] 14 ] 15 16 linearSVM = LinearSVM.LinearSVM(dataset.__len__(), dataset[0][0].__len__()) 17 linearSVM.train(dataset, 100) 18 print(linearSVM) 19 20 for record in dataset: 21 vector = record[0] 22 label = record[-1] 23 if label == 1: 24 plt.plot(vector[0], vector[1], 'r-o') 25 else: 26 plt.plot(vector[0], vector[1], 'g-o') 27 28 predict = linearSVM.predict(vector) 29 print(record.__str__() + predict.__str__() + '\n') 30 31 x1 = np.linspace(0, 1, 50) 32 x2 = (-linearSVM.bias - linearSVM.weight_vec[0] * x1) / linearSVM.weight_vec[1] 33 plt.plot(x1, x2) 34 plt.show()