Python · SVM（四）· SMO 演算法

02-02

（這裡是本章會用到的 GitHub 地址）

（這篇東西我真是覺得又臭又長 ┑(￣Д ￣)┍）

SMO 演算法概述

SMO 是由 Platt 在 1998 年提出的、針對軟間隔最大化 SVM 對偶問題求解的一個演算法，其基本思想很簡單：在每一步優化中，挑選出諸多參數（ $alpha_k(k=1,2,...,N)$ ）中的兩個參數（ $alpha_i$ 、 $alpha_j$ ）作為「真正的參數」，其餘參數都視為常數，從而問題就變成了類似於二次方程求最大值的問題，從而我們就能求出解析解

具體而言，SMO 要解決的是如下對偶問題：

$max_alpha L(alpha)=-frac12sum_{i=1}^Nsum_{j=1}^Nalpha_ialpha_jy_iy_j(x_icdot x_j) + sum_{i=1}^Nalpha_i$

使得對 $i=1,...,N$ 、都有

$sum_{i=1}^Nalpha_iy_i=0$ 、 $0lealpha_ile C$

其大致求解步驟則可以概括如下：

選出 $alpha_1,alpha_2,...,alpha_N$ 中「最不好的」兩個參數 $alpha_i$ 、 $alpha_j$
只把 $alpha_i$ 、 $alpha_j$ 視為參數並把其餘的 $alpha_k$ 視為常數，於是最大化 $L(alpha)$ 就變成了以 $alpha_i$ 、 $alpha_j$ 為參數的二次規劃問題，從而可以直接對其進行求解；但是，注意到 $alpha_i$ 、 $alpha_j$ 需滿足 $sum_{i=1}^Nalpha_iy_i=0$ 和 $0lealpha_i$ 、 $alpha_jle C$ ，所以求完解後需要檢查是否滿足約束；如不滿足，則進行調整

KKT 條件

先來看如何選取參數。在 SMO 演算法中，我們是依次選取參數的：

選出違反 KKT 條件最嚴重的樣本點、以其對應的參數作為第一個參數
第二個參數的選取有一種比較繁複且高效的方法，但對於一個樸素的實現而言、第二個參數即使隨機選取也無不可

這裡就有了一個叫 KKT 條件的東西，其詳細的陳列會放在文末，這裡就僅簡要的說明一下。具體而言，對於已有的模型 $f(x)=wcdot x + b$ 來說， $alpha_i$ 及其對應樣本 $(x_i,y_i)$ 的 KKT 條件為：

$begin{align}nalpha_i=0&Leftrightarrow y_if(x_i)>1 n0<alpha_i<C&Leftrightarrow y_if(x_i)=1 nalpha_i = C&Leftrightarrow y_if(x_i)<1nend{align}$

注意我們之前提過樣本到超平面的函數間隔為 $y(wcdot x+b)$ ，所以上述 KKT 條件可以直觀地敘述為：

$alpha_i=0Leftrightarrow$ 樣本離間隔超平面比較遠
$0<alpha_i<CLeftrightarrow$ 樣本落在間隔超平面上
$alpha_i=CLeftrightarrow$ 樣本在間隔超平面以內

【注意：這裡的間隔超平面即為滿足方程 $y(wcdot x+b)=1$ 的平面；由於 $y$ 可以取正負一兩個值，所以間隔超平面會有兩個—— $wcdot x+b=1$ 和 $wcdot x+b=-1$ 。而分類超平面則是滿足 $wcdot x+b=0$ 的平面，需要將它和間隔超平面加以區分】

可以以一張圖來直觀理解這裡提到的諸多概念：

（畫得有點亂，見諒……）

圖中外面有個黑圓圈的其實就是傳說中的「支持向量」，其定義會在文末給出

那麼我們到底應該如何刻畫「違反 KKT 條件」這麼個東西呢？從直觀上來說，我們可以有這麼一種簡單有效的定義：

計算三份「差異向量」 $c^{(k)}=(c^{(k)}_1,c^{(k)}_2,...,c^{(k)}_N)^T (k=1,2,3)$ ，其中第 $k$ 份對應於三個 KKT 條件中的第 $k$ 個，且 $c^{(k)}_i=y_if(x_i)-1 (i=1,2,...,N)$
針對不同的 KKT 條件，將 $c^{(k)}$ 的某些位置 $c^{(k)}_i$ 置為 0。具體而言：

對第一個 KKT 條件 $alpha_i=0Leftrightarrow y_if(x_i)>1Leftrightarrow c^{(1)}_i>0$ 而言，滿足以下兩種情況的 $c^{(1)}_i$ 將應該置為 0：

$alpha_i>0$ 且 $c^{(1)}_ile0$
$alpha_i=0$ 且 $c^{(1)}_i>0$

對第二個 KKT 條件 $0<alpha_i<CLeftrightarrow y_if(x_i)=1Leftrightarrow c^{(2)}_i=0$ 而言則是：

（ $alpha_i=0$ 或 $alpha_i=C$ ）且 $c^{(2)}_ine0$
$0<alpha_i<C$ 且 $c^{(2)}_i=0$

對第三個 KKT 條件 $alpha_i = CLeftrightarrow y_if(x_i)<1Leftrightarrow c^{(3)}_i<0$ 亦同理：

$alpha_i<C$ 且 $c^{(3)}_ige0$
$alpha_i=C$ 且 $c^{(3)}_i<0$

最後則可以簡單的將三份差異向量的平方相加來作為「損失」，從而直接選出使該損失最大的 $alpha_i$ 作為 SMO 的第一個參數即可。具體而言：
$i=argmax_ileft{ c^{(1)^2}_i + c^{(2)^2}_i + c^{(3)^2}_i | i=1,2,...,Nright}$

得益於 Numpy 強大的 Fancy Indexing，上述置 0 的實現非常好寫（？？？）：

# 得到 alpha > 0 和 alpha < C 的 maskncon1 = alpha > 0ncon2 = alpha < Cn# 算出「差異向量」並拷貝成三份nerr1 = y * y_pred - 1nerr2 = err1.copy()nerr3 = err1.copy()n# 依次根據三個 KKT 條件，將差異向量的某些位置設為 0n# 不難看出為了直觀、我做了不少重複的運算，所以這一步是可以優化的nerr1[(con1 & (err1 <= 0)) | (~con1 & (err1 > 0))] = 0nerr2[((~con1 | ~con2) & (err2 != 0)) | ((con1 & con2) & (err2 == 0))] = 0nerr3[(con2 & (err3 >= 0)) | (~con2 & (err3 < 0))] = 0n# 算出平方和並取出使得「損失」最大的 idxnerr = err1 ** 2 + err2 ** 2 + err3 ** 2nidx = np.argmax(err)n

第二個參數則可以簡單地隨機選取，雖然這不是特別好，但效果已然不錯，而且不僅實現起來更簡便、運行起來也更快（其實就是我太懶）（喂）。具體代碼如下：

idx = np.random.randint(len(self._y))n# 這裡的 idx1 是第一個參數對應的 idxnwhile idx == idx1:n idx = np.random.randint(len(self._y))nreturn idxn

至於 SMO 演算法的第二步，正如前文所說，它的本質就是一個帶約束的二次規劃，雖然求解過程可能會比較折騰，但其實難度不大。具體步驟會放在文末，這裡就暫時按下

SMO 的效果

仍是先看看螺旋線數據集上的訓練過程：

略顯糾結，不過還是不錯的

接下來看看蘑菇數據集上的表現；單就這個數據集而言，我們實現的樸素 SVM 和 sklearn 中的 SVM 表現幾乎是一致的（在使用 RBF 核時），比較具有代表性的訓練曲線則如下圖所示：

也算是符合 SMO 這種每次只取兩個參數進行更新的訓練方法的直觀

相關數學理論

1）KKT 條件的詳細陳列

注意到原始問題為

$min_{w,b,xi}L(w,b,xi)=frac12|w|^2+Csum_{i=1}^Nxi_i$ ，使得 $xi^*ge0$ 、 $y_i(wcdot x_i+b)ge1-xi_i$ （不妨稱這兩個約束為原始約束）

所以其拉格朗日運算元法對應的式子為

$L=L(w,b,xi,alpha,beta)=frac12|w|^2+Csum_{i=1}^Nxi_i-sum_{i=1}^Nalpha_i[y_i(wcdot x_i+b)-1+xi_i]-sum_{i=1}^Nbeta_ixi_i$

於是 KKT 條件的其中四個約束即為（不妨設最優解為 $w^*$ 、 $b^*$ 、 $xi^*$ 、 $alpha^*$ 和 $beta^*$ ）：

$alpha_i^*ge0,beta_i^*ge0$ （這是拉格朗日乘子法自身的要求）
$xi_i^*ge0$ 、 $y_i(w^*cdot x_i+b^*)-1+xi_i^*ge0$ （此即原始約束）
$alpha_i^*[y_i(w^*cdot x_i+b^*)-1+xi_i^*]=0$ （換句話說， $alpha_i^*$ 和 $y_i(w^*cdot x_i+b)-1+xi_i^*$ 中必有一個為 0）

該等式有著很好的直觀：設想它們同時不為 0，則必有 $y_i(w^*cdot x_i+b)-1+xi_i^*>0$ （注意原始約束）、從而 $alpha_i^*[y_i(w^*cdot x_i+b^*)-1+xi_i^*]ge0$ ，等號當且僅當 $alpha_i=0$ 時取得。然而由於 $alpha_i^*ne0$ ，所以若將 $alpha_i$ 取為 0、則上述 $L$ 將會變大。換句話說，將參數 $alpha_i$ 取為 0 將會使得目標函數比參數取 $alpha_i^*$ 時的目標函數要大，這與 $alpha_i^*$ 的最優性矛盾

$beta_i^*xi_i^*=0$ （換句話說， $beta_i^*$ 和 $xi_i^*$ 中必有一個為 0，理由同上）

從而原始問題轉為 $min_{w,b}max_alpha L$ ；而對偶問題的實質，其實就是將原始問題 $min_{w,b}max_alpha L$ 轉為 $max_alphamin_{w,b} L$ 。在求解 $nabla_wL=nabla_bL=nabla_xi L=0$ 後，可以得到如下對偶問題：

$max_alpha L(alpha)=-frac12sum_{i=1}^Nsum_{j=1}^Nalpha_ialpha_jy_iy_j(x_icdot x_j) + sum_{i=1}^Nalpha_i$ ，使得對 $i=1,...,N$ 、都有 $sum_{i=1}^Nalpha_iy_i=0$ 、 $0lealpha_ile C$

（雖然這些在 Python · SVM（二）· LinearSVM 中介紹過，不過為了連貫性，這裡還是再介紹一遍）

於是，最優解自然需要滿足這麼個條件：

$nabla_wL(w^*,b^*,xi^*,alpha^*,beta^*)=nabla_bL(w^*,b^*,xi^*,alpha^*,beta^*)=nabla_xi L(w^*,b^*,xi^*,alpha^*,beta^*)=0$

這個條件即是最後一個 KKT 條件

2）何謂「支持向量」

為方便說明，這裡再次放出上文給出過的圖：

圖中帶黑圈的樣本點即是支持向量，數學上來說的話，就是 $alpha_i>0$ 對應的樣本點即是支持向量。從圖中不難看出，支持向量從直觀上來說，就是比較難分的樣本點

此外，支持向量之所以稱之為「支持」向量，是因為在理想情況下，僅利用支持向量訓練出來的模型和利用所有樣本訓練出來的模型是一致的。這從直觀上是好理解的，粗略的證明則可以利用其定義來完成：非支持向量的樣本對應著 $alpha_i=0$ ，亦即它對最終模型—— $f(x)=sum_{i=1}^Nalpha_iy_i(x_icdot x)+b$ 沒有絲毫貢獻，所以可以直接扔掉

3）帶約束的二次規劃求解方法

不妨設我們選取出來的兩個參數就是 $alpha_1$ 和 $alpha_2$ ，那麼問題的關鍵就在於如何把 $alpha_1$ 和 $alpha_2$ 相關的東西抽取出來並把其它東西扔掉

注意到我們的對偶問題為

$max_alpha L(alpha)=-frac12sum_{i=1}^Nsum_{j=1}^Nalpha_ialpha_jy_iy_j(x_icdot x_j) + sum_{i=1}^Nalpha_i$ ，使得對 $i=1,...,N$ 、都有 $sum_{i=1}^Nalpha_iy_i=0$ 、 $0lealpha_ile C$

且我們在 Python · SVM（一）· 感知機的最後介紹過 Gram 矩陣：

$G=(x_icdot x_j)_{Ntimes N}$

所以 $L$ 就可以改寫為 $L(alpha)=-frac12sum_{i=1}^Nsum_{j=1}^Nalpha_ialpha_jy_iy_jG_{ij}+sum_{i=1}^Nalpha_i$

把和 $alpha_1$ 、 $alpha_2$ 無關的東西扔掉之後， $L$ 可以化簡為：

$L(alpha)=-frac12(G_{11}alpha_1^2+2y_1y_2G_{12}alpha_1alpha_2+G_{22}alpha_2^2)-left(y_1alpha_1sum_{i=3}^Ny_ialpha_iG_{i1}+y_2alpha_2sum_{i=3}^Ny_ialpha_iG_{i2}right)+(alpha_1+alpha_2)$

約束條件則可以化簡為對 $i=1$ 和 $i=2$ ，都有 $y_1alpha_1+y_2alpha_2=-sum_{i=3}^Ny_ialpha_i=c$ 、 $0lealpha_ile C$ ，其中 $c$ 是某個常數

而帶約束的二次規劃求解過程也算簡單：只需先求出無約束下的最優解，然後根據約束「裁剪」該最優解即可

無約束下的求解過程其實就是求偏導並令其為 0。以 $alpha_1$ 為例，注意到

$y_1alpha_1+y_2alpha_2=cRightarrowalpha_2=frac c{y_2}-frac{y_1}{y_2}alpha_1$

令 $c^*=frac c{y_2}, s=y_1y_2$ ，則 $c^*$ 亦是常數，且由於 $y_1$ 、 $y_2$ 都只能取正負 1，故不難發現 $frac{y_2}{y_1}=frac{y_1}{y_2}=s$ ，從而 $alpha_2=c^*-salpha_1Rightarrowfrac{partialalpha_2}{partialalpha_1}=-s$

於是

$begin{align}nfrac{partial L}{partialalpha_1}=&-G_{11}alpha_1-y_1y_2G_{12}(alpha_2+alpha_1frac{partialalpha_2}{partialalpha_1})-G_{22}alpha_2frac{partialalpha_2}{partialalpha_1} n&-y_1sum_{i=3}^Ny_ialpha_iG_{i1}-y_2frac{partialalpha_2}{partialalpha_1}sum_{i=3}^Ny_ialpha_iG_{i2}+1 n=&-G_{11}alpha_1-sG_{12}(c^*-salpha_1-alpha_1cdot s)-G_{22}(c^*-salpha_1)cdot(-s) n&-y_1sum_{i=3}^Ny_ialpha_iG_{i1}+sy_2sum_{i=3}^Ny_ialpha_iG_{i2}+left(frac{partialalpha_2}{partialalpha_1}+1right)nend{align}$

考慮到 $s^2=1$ 、 $sy_2=y_1$ 、Gram 矩陣是對稱陣、且模型在第 $k$ 個樣本 $x_k$ 處的輸出為 $f(x_k)=sum_{i=1}^Nalpha_iy_i(x_icdot x_k)+b=sum_{i=1}^Nalpha_iy_iG_{ik}+b$ ，從而可知

$begin{align}nfrac{partial L}{partialalpha_1}=&-G_{11}alpha_1-sG_{12}c^*+2G_{12}alpha_1+sG_{22}c^*-G_{22}alpha_1 n&-y_1[f(x_1)-y_1alpha_1G_{11}-y_2alpha_2G_{21}] n&+y_1[f(x_2)-y_1alpha_1G_{12}-y_2alpha_2G_{22}] +(1-s)nend{align}$

令 $v_i=(f(x_i)-b)-y_1alpha_1G_{1i}-y_2alpha_2G_{2i} (i=1,2)$ ，則

$frac{partial L}{partialalpha_1}=-(G_{11}-2G_{12}+G_{22})alpha_1-sc^*(G_{12}-G_{22})-y_1(v_1-v_2)+(1-s)$

於是

$begin{align}nfrac{partial L}{partialalpha_1}=0Rightarrowalpha_1&=-frac{sc^*(G_{12}-G_{22})+y_1(v_1-v_2)-(1-s)}{G_{11}-2G_{12}+G_{22}} n&=-frac{y_1[y_2c^*(G_{12}-G_{22})+(v_1-v_2)-(y_1-y_2)]}{G_{11}-2G_{12}+G_{22}}nend{align}$

注意到 $c^*=salpha_1+alpha_2$ ，從而

$y_2c^*(G_{12}-G_{22})=y_2(salpha_1+alpha_2)(G_{12}-G_{22})=(y_1alpha_1+y_2alpha_2)(G_{12}-G_{22})$

令 $dG=G_{11}-2G_{12}+G_{22}$ 、 $e_i=f(x_i)-y_i (i=1,2)$ ，則

$y_2c^*(G_{12}-G_{22})+(v_1-v_2)-(y_2+y_1)=... =e_1 - e_2 - y_1alpha_1dG$

從而

$alpha_1^{new,raw}=alpha_1^{old}-frac{y_1(e_1-e_2)}{dG}$

接下來就要對其進行裁剪了。注意到我們的約束為

$0lealpha_ile C$ 、 $alpha_1y_1+alpha_2y_2$ 為常數

所以我們需要分情況討論 $alpha_1$ 的下、上界

當 $y_1,y_2$ 異號（ $y_1y_2=-1$ ）時，可知 $alpha_1-alpha_2$ 為常數、亦即
$alpha_1^{new}-alpha_2^{new}=alpha_1^{old}-alpha_2^{old}Rightarrowalpha_2^{new}=alpha_1^{new}-(alpha_1^{old}-alpha_2^{old})$
結合 $0lealpha_2le C$ ，可知：

$alpha_1^{new}$ 不應小於 $alpha_1^{old}-alpha_2^{old}$ ，否則 $alpha_2$ 將小於 0
$alpha_1^{new}$ 不應大於 $C+alpha_1^{old}-alpha_2^{old}$ ，否則 $alpha_2$ 將大於 C

當 $y_1,y_2$ 同號（ $y_1y_2=1$ ）時，可知 $alpha_1+alpha_2$ 為常數、亦即
$alpha_1^{new}+alpha_2^{new}=alpha_1^{old}+alpha_2^{old}Rightarrowalpha_2^{new}=(alpha_1^{old}+alpha_2^{old}) - alpha_1^{new}$
結合 $0lealpha_2le C$ ，可知：

$alpha_1^{new}$ 不應小於 $alpha_1^{old}+alpha_2^{old}-C$ ，否則 $alpha_2$ 將大於 C
$alpha_1^{new}$ 不應大於 $alpha_1^{old}+alpha_2^{old}$ ，否則 $alpha_2$ 將小於 0

綜上可知

$alpha_1^{new}$ 的下界為 $U=left{begin{aligned}nmax{0,alpha_1^{old}-alpha_2^{old}} &y_1y_2=-1 nmax{0,alpha_1^{old}+alpha_2^{old}-C} &y_1y_2=1nend{aligned}nright.$
$alpha_1^{new}$ 的上界為 $V=left{begin{aligned}nmin{C,C+alpha_1^{old}-alpha_2^{old}} &y_1y_2=-1 nmax{C,alpha_1^{old}+alpha_2^{old}} &y_1y_2=1nend{aligned}nright.$

那麼直接做一個 clip 即可得到更新後的 $alpha_1$ ：

alpha1_new = np.clip(alpha1_new_raw, u, v)n

注意由於我們要保持 $alpha_1y_1+alpha_2y_2$ 為常數，所以（注意 $frac{y_1}{y_2}=y_1y_2$ ）

$begin{align}nalpha_2^{new}&=frac1{y_2}(alpha_1^{old}y_1+alpha_2^{old}y_2-alpha_1^{new}y_1) n&=alpha_2^{old}+y_1y_2(alpha_1^{old}-alpha_1^{new})nend{align}$

綜上所述，我們就完成了一次參數的更新，之後就不斷地更新直至滿足停機條件即可

此外，我在 Python · SVM（三）· 核方法這篇文章中提到過，對 SVM 的對偶形式應用核方法會非常自然。表現在 SMO 上的話就是，我們可以通過簡單粗暴地將核矩陣 $K$ 代替 Gram 矩陣 $G$ 來完成核方法的應用。直觀地說，我們只需將上面所有出現過的 $G$ 都換成 $K$ 就行了

至此，SVM 演算法的介紹就大致告一段落了。我們從感知機出發，依次介紹了「極大梯度法」、MBGD（Batch 梯度下降）法、核方法和 SMO 演算法；雖然都有點淺嘗輒止的味道，但覆蓋的東西……大概還是挺多的

希望觀眾老爺們能夠喜歡~