支持向量機(SVM)推導以及代碼實現

SVM 是一個非常優雅的算法，具有完善的數學理論，雖然如今工業界用到的不多，但還是決定花點時間去寫篇文章整理一下。

1. 支持向量

1.1 線性可分

首先我們先來了解下什么是線性可分。

在二維空間上，兩類點被一條直線完全分開叫做線性可分。

嚴格的數學定義是：

$D_0$ 和 $D_1$ 是 n 維歐氏空間中的兩個點集。如果存在 n 維向量 w 和實數 b，使得所有屬於 $D_0$ 的點 $x_i$ 都有 $wx_i + b > 0$ ，而對於所有屬於 $D_1$ 的點 $x_j$ 則有 $wx_j + b < 0$ ，則我們稱 $D_0$ 和 $D_1$ 線性可分。

1.2 最大間隔超平面

從二維擴展到多維空間中時，將 $D_0$ 和 $D_1$ 完全正確地划分開的 $wx+b=0$ 就成了一個超平面。

為了使這個超平面更具魯棒性，我們會去找最佳超平面，以最大間隔把兩類樣本分開的超平面，也稱之為最大間隔超平面。

兩類樣本分別分割在該超平面的兩側；
兩側距離超平面最近的樣本點到超平面的距離被最大化了。

1.3 支持向量

樣本中距離超平面最近的一些點，這些點叫做支持向量。

1.4 SVM 最優化問題

SVM 想要的就是找到各類樣本點到超平面的距離最遠，也就是找到最大間隔超平面。任意超平面可以用下面這個線性方程來描述：

$w^Tx+b=0 \\$

二維空間點 $(x,y)$ 到直線 $Ax+By+C=0$ 的距離公式是：

$\frac{|Ax+By+C|}{\sqrt{A^2+B^2}} \\$

擴展到 n 維空間后，點 $x=(x_1,x_2…x_n)$ 到直線 $w^Tx+b=0$ 的距離為：

$\frac{|w^Tx+b|}{||w||} \\$

其中 $||w||=\sqrt{w_1^2+…w_n^2}$ 。

如圖所示，根據支持向量的定義我們知道，支持向量到超平面的距離為 d，其他點到超平面的距離大於 d。

於是我們有這樣的一個公式：

$\left\{ \begin{aligned} \frac{w^Tx+b}{||w||} &\geq d \quad y=1 \\ \frac{w^Tx+b}{||w||} &\leq -d \quad y=-1 \end{aligned} \right. \\$

稍作轉化可以得到：

$\left\{ \begin{aligned} \frac{w^Tx+b}{||w||d} &\geq 1 \quad y=1 \\ \frac{w^Tx+b}{||w||d} &\leq -1 \quad y=-1 \end{aligned} \right. \\$

$||w|| d$ 是正數，我們暫且令它為 1（之所以令它等於 1，是為了方便推導和優化，且這樣做對目標函數的優化沒有影響），故：

$\left\{ \begin{aligned} w^Tx+b &\geq 1 \quad y=1 \\ w^Tx+b &\leq -1 \quad y=-1 \end{aligned} \right. \\$

將兩個方程合並，我們可以簡寫為：

$y(w^Tx+b) \geq 1 \\$

至此我們就可以得到最大間隔超平面的上下兩個超平面：

每個支持向量到超平面的距離可以寫為：

$d=\frac{|w^Tx+b|}{||w||} \\$

由上述 $y(w^Tx+b) > 1 > 0$ 可以得到 $y(w^Tx+b) = |w^Tx+b|$ ，所以我們得到：

$d = \frac{y(w^Tx+b)}{||w||} \\$

最大化這個距離：

$\max 2* \frac{y(w^Tx+b)}{||w||} \\$

這里乘上 2 倍也是為了后面推導，對目標函數沒有影響。剛剛我們得到支持向量 $y(w^Tx+b) = 1$ ，所以我們得到：

$\max \frac{2}{||w||} \\$

再做一個轉換：

$\min \frac{1}{2}||w|| \\$

為了方便計算（去除 $||w||$ 的根號），我們有：

$\min \frac{1}{2}||w||^2\\$

所以得到的最優化問題是：

$\min \frac{1}{2} ||w||^2 \ s.t. \quad y_i（w^Tx_i+b）\geq 1 \\$

2. 對偶問題

2.1 拉格朗日乘數法

2.1.1 等式約束優化問題

本科高等數學學的拉格朗日程數法是等式約束優化問題：

$\min f(x_{1} ,x_{2} ,...,x_{n} ) \\ s.t. \quad h_{k} (x_{1} ,x_{2} ,...,x_{n} )=0 \quad k =1,2,...,l\\$

我們令 $L(x,\lambda ) = f(x) + \sum\limits_{k = 1}^l \lambda _k h_k (x)$ ，函數 $L(x,y)$ 稱為 Lagrange 函數，參數 $\lambda$ 稱為 Lagrange 乘子沒有非負要求。

利用必要條件找到可能的極值點：

$\left\{ \begin{aligned} \frac{\partial L}{\partial x_i} = 0 \quad i=1,2,...,n \\ \frac{\partial L}{\partial \lambda_k} = 0 \quad k=1,2,...,l \end{aligned} \right. \\$

具體是否為極值點需根據問題本身的具體情況檢驗。這個方程組稱為等式約束的極值必要條件。

等式約束下的 Lagrange 乘數法引入了 $l$ 個 Lagrange 乘子，我們將 $x_{i}$ 與 $\lambda_{k}$ 一視同仁，把 $\lambda_{k}$ 也看作優化變量，共有 $(n+l)$ 個優化變量。

2.1.2 不等式約束優化問題

而我們現在面對的是不等式優化問題，針對這種情況其主要思想是將不等式約束條件轉變為等式約束條件，引入松弛變量，將松弛變量也是為優化變量。

以我們的例子為例：

$min f(w) = min\frac{1}{2} ||w||^2 \\ s.t. \quad g_i(w) = 1 - y_i(w^Tx_i+b)\leq 0 \\$

我們引入松弛變量 $a_i^2$ 得到 $h_i(w,a_i) = g_i(w) + a_i^2 = 0$ 。這里加平方主要為了不再引入新的約束條件，如果只引入 $a_i$ 那我們必須要保證 $a_i \geq 0$ 才能保證 $h_i(w,a_i) = 0$ ，這不符合我們的意願。

由此我們將不等式約束轉化為了等式約束，並得到 Lagrange 函數：

$\begin{aligned} L(w,\lambda,a) &= {f(w)} + \sum\limits_{i = 1}^n \lambda_i h_i (w) \\ &= {f(w)} + \sum\limits_{i = 1}^n \lambda_i [g_i(w) + a_i^2] \quad \lambda_i \geq 0 \end{aligned} \\$

由等式約束優化問題極值的必要條件對其求解，聯立方程：

$\left\{ \begin{aligned} \frac{\partial L}{\partial w_i} &= \frac{\partial f}{\partial w_i} + \sum\limits_{i=1}^{n} \lambda_i \frac{\partial g_i}{\partial w_i}= 0, \\ \frac{\partial L}{\partial a_i} &= 2 \lambda_i a_i = 0, \\ \frac{\partial L}{\partial \lambda_i} &= g_i(w) + a_i^2 = 0, \\ \lambda_i &\geq 0 \end{aligned} \right. \\$

（為什么取 $\lambda_i \geq 0$ ，可以通過幾何性質來解釋，有興趣的同學可以查下 KKT 的證明）。

針對 $\lambda_i a_i = 0$ 我們有兩種情況：

情形一： $\lambda_i = 0, a_i \neq 0$

由於 $\lambda_i=0$ ，因此約束條件 $g_i(w)$ 不起作用，且 $g_i(w)<0$

情形二： $\lambda_i \neq 0, a_i = 0$

此時 $g_i(w)=0$ 且 $\lambda_i>0$ ，可以理解為約束條件 $g_i(w)$ 起作用了，且 $g_i(w)=0$

綜合可得： $\lambda_ig_i(w)=0$ ，且在約束條件起作用時 $\lambda_i>0,g_i(w)=0$ ；約束不起作用時 $\lambda_i = 0,g_i(w) < 0$

由此方程組轉換為：

$\left\{ \begin{aligned} \frac{\partial L}{\partial w_i} &= \frac{\partial f}{\partial w_i} + \sum\limits_{j=1}^{n} \lambda_j \frac{\partial g_j}{\partial w_i}= 0, \\ \lambda_ig_i(w) &= 0, \\ g_i(w)&\leq 0 \\ \lambda_i &\geq 0 \end{aligned} \right. \\$

以上便是不等式約束優化優化問題的 KKT(Karush-Kuhn-Tucker) 條件， $\lambda_i$ 稱為 KKT 乘子。

這個式子告訴了我們什么事情呢？

直觀來講就是，支持向量 $g_i(w)=0$ ，所以 $\lambda_i > 0$ 即可。而其他向量 $g_i(w)<0, \lambda_i=0$ 。

我們原本問題時要求： $min \frac{1}{2} ||w||^2$ ，即求 $minL(w,\lambda,a)$

$\begin{aligned} L(w,\lambda,a) &= {f(w)} + \sum\limits_{i = 1}^n \lambda_i [g_i(w) + a_i^2] \quad \\ &= {f(w)} + \sum\limits_{i = 1}^n \lambda_i g_i(w) + \sum\limits_{i = 1}^n \lambda_i a_i^2 \end{aligned} \\$

由於 $\sum\limits_{i = 1}^n \lambda_i a_i^2 \geq 0$ ，故我們將問題轉換為： $minL(w,\lambda)$ ：

$L(w,\lambda)={f(w)} + \sum\limits_{i = 1}^n \lambda_i g_i(w) \\$

假設找到了最佳參數是的目標函數取得了最小值 p。即 $\frac{1}{2} ||w||^2 =p$ 。而根據 $\lambda_{i} \geq 0$ ，可知 $\sum\limits_{i = 1}^n \lambda_i g_i(w) \leq 0$ ，因此 $L(w,\lambda) \leq p$ ，為了找到最優的參數 ${\lambda}$ ，使得 $L(w,\lambda)$ 接近 p，故問題轉換為出 $\max\limits_{\lambda}L(w,\lambda)$ 。

故我們的最優化問題轉換為：

$\min\limits_w \max\limits_{\lambda} L(w,\lambda) \\ s.t. \quad \lambda_i \geq 0 \\$

出了上面的理解方式，我們還可以有另一種理解方式：由於 $\lambda_i \geq 0$ ，

$\max\limits_{\lambda} L(w, \lambda) = \left\{ \begin{aligned} \infty \quad g_i(w) \geq 0 \\ \frac{1}{2} {||w||^2} \quad g_i(w) \leq 0 \end{aligned} \right. \\$

所以 $\min(\infty, \frac{1}{2} {||w||^2}) = \frac{1}{2} {||w||^2}$ ，所以轉化后的式子和原來的式子也是一樣的。

2.2 強對偶性

對偶問題其實就是將：

$\min\limits_w \max\limits_{\lambda} L(w,\lambda) \\ s.t. \quad \lambda_i \geq 0 \\$

變成了：

$\max\limits_{\lambda} \min\limits_w L(w,\lambda) \\ s.t. \quad \lambda_i \geq 0 \\$

假設有個函數 $f$ 我們有：

$\min\max f \geq \max\min f \\$

也就是說，最大的里面挑出來的最小的也要比最小的里面挑出來的最大的要大。這關系實際上就是弱對偶關系，而強對偶關系是當等號成立時，即：

$\min\max f = \max\min f \\$

如果 $f$ 是凸優化問題，強對偶性成立。而我們之前求的 KKT 條件是強對偶性的充要條件。

3. SVM 優化

我們已知 SVM 優化的主問題是：

$\min\limits_{w} \frac{1}{2} ||w||^2 \\ s.t.\quad g_i(w,b) = 1 - y_i（w^Tx_i+b）\leq 0, \quad i=1,2,...,n \\$

那么求解線性可分的 SVM 的步驟為：

步驟 1：

構造拉格朗日函數：

$\min\limits_{w,b}\max\limits_{\lambda} L(w,b,\lambda)= \frac{1}{2}{||w||}^2 + \sum\limits_{i = 1}^n \lambda_i [1-y_i(w^Tx_i+b)] \\ s.t. \quad \lambda_i \geq 0 \\$

步驟 2：

利用強對偶性轉化：

$\max\limits_{\lambda}\min\limits_{w,b} L(w,b,\lambda) \\$

現對參數 w 和 b 求偏導數：

$\begin{aligned} \frac{\partial L}{\partial w} &= w - \sum_{i=1}^{n}\lambda_ix_iy_i = 0 \\ \frac{\partial L}{\partial b} &= \sum_{i=1}^{n}\lambda_iy_i = 0 \\ \end{aligned} \\$

得到：

$\begin{aligned} \sum_{i=1}^{n}\lambda_ix_iy_i &= w\\ \sum_{i=1}^{n}\lambda_iy_i &= 0 \\ \end{aligned} \\$

我們將這個結果帶回到函數中可得：

$\begin{aligned} L(w,b,\lambda) &= \frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j) + \sum_{i = 1}^n \lambda_i - \sum_{i = 1}^n \lambda_i y_i(\sum_{j = 1}^n \lambda_j y_j (x_i \cdot x_j) + b) \\ &= \frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j) + \sum_{i = 1}^n \lambda_i - \sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j)-\sum_{i = 1}^n \lambda_i y_i b \\ &= \sum_{j=1}^{n}\lambda_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j) \end{aligned} \\$

也就是說：

$\min\limits_{w,b}L(w,b,\lambda) = \sum_{j=1}^{n}\lambda_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j) \\$

步驟 3：

由步驟 2 得：

$\max\limits_{\lambda} [\sum_{j=1}^{n}\lambda_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j)] \\ s.t. \quad \sum_{i=1}^{n}\lambda_iy_i = 0 \quad \lambda_i \geq 0 \\$

我們可以看出來這是一個二次規划問題，問題規模正比於訓練樣本數，我們常用 SMO(Sequential Minimal Optimization) 算法求解。

SMO(Sequential Minimal Optimization)，序列最小優化算法，其核心思想非常簡單：每次只優化一個參數，其他參數先固定住，僅求當前這個優化參數的極值。我們來看一下 SMO 算法在 SVM 中的應用。

我們剛說了 SMO 算法每次只優化一個參數，但我們的優化目標有約束條件： $\sum\limits_{i=1}^{n}\lambda_iy_i = 0$ ，沒法一次只變動一個參數。所以我們選擇了一次選擇兩個參數。具體步驟為：

選擇兩個需要更新的參數 $\lambda_i$ 和 $\lambda_j$ ，固定其他參數。於是我們有以下約束：

這樣約束就變成了：

$\lambda_i y_i+\lambda_j y_j = c \quad \lambda_i \geq 0,\lambda_j \geq 0 \\$

其中 $c=-\sum\limits_{k \ne i,j}\lambda_ky_k$ ，由此可以得出 $\lambda_j=\frac{c-\lambda_iy_i}{y_j}$ ，也就是說我們可以用 $\lambda_i$ 的表達式代替 $\lambda_{j}$ 。這樣就相當於把目標問題轉化成了僅有一個約束條件的最優化問題，僅有的約束是 $\lambda_i \geq 0$ 。

2. 對於僅有一個約束條件的最優化問題，我們完全可以在 $\lambda_{i}$ 上對優化目標求偏導，令導數為零，從而求出變量值 $\lambda_{i_{new}}$ ，然后根據 $\lambda_{i_{new}}$ 求出 $\lambda_{j_{new}}$ 。

3. 多次迭代直至收斂。

通過 SMO 求得最優解 $\lambda^*$ 。

步驟 4 ：

我們求偏導數時得到：

$w = \sum_{i=1}^m \lambda_i y_i x_i \\$

由上式可求得 w。

我們知道所有 $\lambda_i > 0$ 對應的點都是支持向量，我們可以隨便找個支持向量，然后帶入： $y_s(wx_s+b) = 1$ ，求出 b 即可，

兩邊同乘 $y_s$ ，得 $y_s^2(wx_s+b) = y_s$

因為 $y_s^2=1$ ，所以： $b=y_s - wx_s$

為了更具魯棒性，我們可以求得支持向量的均值：

$b = \frac{1}{|S|}\sum_{s \in S} (y_s -wx_s) \\$

步驟 5： w 和 b 都求出來了，我們就能構造出最大分割超平面： $w^Tx+b=0$

分類決策函數： $f(x)=sign(w^Tx+b)$

其中 $sign( \cdot )$ 為階躍函數：

$sign(x) = \left\{ \begin{aligned} -1 \quad x<0 \\ 0 \quad x=0 \\ 1 \quad x>0 \end{aligned} \right. \\$

將新樣本點導入到決策函數中既可得到樣本的分類。

4. 軟間隔

4.1 解決問題

在實際應用中，完全線性可分的樣本是很少的，如果遇到了不能夠完全線性可分的樣本，我們應該怎么辦？比如下面這個：

於是我們就有了軟間隔，相比於硬間隔的苛刻條件，我們允許個別樣本點出現在間隔帶里面，比如：

我們允許部分樣本點不滿足約束條件：

$1-y_i(w^Tx_i + b) \leq 0 \\$

為了度量這個間隔軟到何種程度，我們為每個樣本引入一個松弛變量 $\xi_{i}$ ，令 $\xi_{i} \geq 0$ ，且 $1 - y_i(w^Tx_i + b)-\xi_i \leq 0$ 。對應如下圖所示：

4.2 優化目標及求解

增加軟間隔后我們的優化目標變成了：

$\min\limits_{w} \frac{1}{2} ||w||^2 + C\sum_{i=1}^{m}\xi_i \\ s.t.\quad g_i(w,b) = 1 - y_i(w^Tx_i+b) - \xi_i\leq 0, \quad \xi_i \geq 0, \quad i=1,2,...,n \\$

其中 C 是一個大於 0 的常數，可以理解為錯誤樣本的懲罰程度，若 C 為無窮大， $\xi_{i}$ 必然無窮小，如此一來線性 SVM 就又變成了線性可分 SVM；當 C 為有限值的時候，才會允許部分樣本不遵循約束條件。

接下來我們將針對新的優化目標求解最優化問題：

步驟 1：

構造拉格朗日函數：

$\min\limits_{w,b,\xi}\max\limits_{\lambda, \mu} L(w,b,\xi,\lambda,\mu)= \frac{1}{2}{||w||}^2 + C\sum_{i=1}^{m}\xi_i+ \sum\limits_{i = 1}^n \lambda_i [1-\xi_i-y_i(w^Tx_i+b)] - \sum_{i=1}^{n}\mu_i\xi_i \\ s.t. \quad \lambda_i \geq 0 \quad \mu_i \geq 0 \\$

其中 $\lambda_{i}$ 和 $\mu_{i}$ 是拉格朗日乘子，w、b 和 $\xi_{i}$ 是主問題參數。

根據強對偶性，將對偶問題轉換為：

$\max\limits_{\lambda, \mu}\min\limits_{w,b,\xi} L(w,b,\xi,\lambda,\mu) \\$

步驟 2：

分別對主問題參數w、b 和 $\xi_{i}$ 求偏導數，並令偏導數為 0，得出如下關系：

$w = \sum_{i=1}^{m}\lambda_i y_i x_i \\ 0 = \sum_{i=1}^{m}\lambda_i y_i \\ C = \lambda_i + \mu_i \\$

將這些關系帶入拉格朗日函數中，得到：

$\min\limits_{w,b,\xi}L(w,b,\xi,\lambda,\mu) = \sum_{j=1}^{n}\lambda_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (x_i \cdot x_j) \\$

最小化結果只有 $\lambda$ 而沒有 $\mu$ ，所以現在只需要最大化 $\lambda$ 就好：

我們可以看到這個和硬間隔的一樣，只是多了個約束條件。

然后我們利用 SMO 算法求解得到拉格朗日乘子 $\lambda^{*}$ 。

步驟 3 ：

$w = \sum_{i=1}^m \lambda_i y_i x_i \\ b = \frac{1}{|S|}\sum_{s \in S} (y_s -wx_s) \\$

然后我們通過上面兩個式子求出 w 和 b，最終求得超平面 $w^Tx+b=0$ ，

這邊要注意一個問題，在間隔內的那部分樣本點是不是支持向量？

我們可以由求參數 w 的那個式子可看出，只要 $\lambda_{i} > 0$ 的點都能夠影響我們的超平面，因此都是支持向量。

5. 核函數

5.1 線性不可分

我們剛剛討論的硬間隔和軟間隔都是在說樣本的完全線性可分或者大部分樣本點的線性可分。

但我們可能會碰到的一種情況是樣本點不是線性可分的，比如：

這種情況的解決方法就是：將二維線性不可分樣本映射到高維空間中，讓樣本點在高維空間線性可分，比如：

對於在有限維度向量空間中線性不可分的樣本，我們將其映射到更高維度的向量空間里，再通過間隔最大化的方式，學習得到支持向量機，就是非線性 SVM。

我們用 x 表示原來的樣本點，用 $\phi(x)$ 表示 x 映射到特征新的特征空間后到新向量。那么分割超平面可以表示為： $f(x)=w \phi(x)+b$ 。

對於非線性 SVM 的對偶問題就變成了：

$\min\limits_{\lambda} [\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\lambda_i \lambda_j y_i y_j (\phi(x_i) \cdot \phi(x_j))-\sum_{j=1}^{n}\lambda_i] \\ s.t. \quad \sum_{i=1}^{n}\lambda_iy_i = 0, \quad \lambda_i \geq 0, \quad C-\lambda_i-\mu_i=0 \\$

可以看到與線性 SVM 唯一的不同就是：之前的 $(x_i \cdot x_j)$ 變成了 $(\phi(x_i) \cdot \phi(x_j))$ 。

5.2 核函數的作用

我們不禁有個疑問：只是做個內積運算，為什么要有核函數的呢？

這是因為低維空間映射到高維空間后維度可能會很大，如果將全部樣本的點乘全部計算好，這樣的計算量太大了。

但如果我們有這樣的一核函數 $k(x,y) = (\phi(x),\phi(y))$ ， $x_i$ 與 $x_j$ 在特征空間的內積等於它們在原始樣本空間中通過函數 $k( x, y)$ 計算的結果，我們就不需要計算高維甚至無窮維空間的內積了。

舉個例子：假設我們有一個多項式核函數：

$k(x,y)=(x \cdot y + 1)^2 \\$

帶進樣本點的后：

$k(x,y) = (\sum_{i=1}^n(x_i \cdot y_i) + 1)^2 \\$

而它的展開項是：

$\sum_{i=1}^nx_i^2y_i^2+\sum_{i=2}^n\sum_{j=1}^{i-1}(\sqrt2x_ix_j)(\sqrt2y_iy_j)+\sum_{i=1}{n}(\sqrt2x_i)(\sqrt2y_i)+1 \\$

如果沒有核函數，我們則需要把向量映射成：

$x^{'} = (x_1^2,...,x_n^2,...\sqrt2x_1,...,\sqrt2x_n,1) \\$

然后在進行內積計算，才能與多項式核函數達到相同的效果。

可見核函數的引入一方面減少了我們計算量，另一方面也減少了我們存儲數據的內存使用量。

5.3 常見核函數

我們常用核函數有：

線性核函數

$k(x_i,x_j) = x_i^Tx_j \\$

多項式核函數

$k(x_i,x_j) = (x_i^Tx_j)^d\\$

高斯核函數

$k(x_i,x_j) = exp(-\frac{||x_i-x_j||}{2\delta^2}) \\$

這三個常用的核函數中只有高斯核函數是需要調參的。

6. 優缺點

6.1 優點

有嚴格的數學理論支持，可解釋性強，不依靠統計方法，從而簡化了通常的分類和回歸問題；
能找出對任務至關重要的關鍵樣本（即：支持向量）；
采用核技巧之后，可以處理非線性分類/回歸任務；
最終決策函數只由少數的支持向量所確定，計算的復雜性取決於支持向量的數目，而不是樣本空間的維數，這在某種意義上避免了“維數災難”。

6.2 缺點

訓練時間長。當采用 SMO 算法時，由於每次都需要挑選一對參數，因此時間復雜度為 $O(N^2)$ ，其中 N 為訓練樣本的數量；
當采用核技巧時，如果需要存儲核矩陣，則空間復雜度為 $O(N^2)$ ；
模型預測時，預測時間與支持向量的個數成正比。當支持向量的數量較大時，預測計算復雜度較高。

因此支持向量機目前只適合小批量樣本的任務，無法適應百萬甚至上億樣本的任務。

import random
from numpy import *


# SMO算法相關輔助中的輔助函數
# 1 解析文本數據函數，提取每個樣本的特征組成向量，添加到數據矩陣
def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    with open(fileName) as f:
        for line in f.readlines():
            lineArr = line.strip().split('\t')
            dataMat.append([float(lineArr[0]), float(lineArr[1])])
            labelMat.append(float(lineArr[2]))
    return dataMat, labelMat


# 2 在樣本集中采取隨機選擇的方法選取第二個不等於第一個alphai的
# 優化向量alphaj
def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j


# 3 約束范圍L<=alphaj<=H內的更新后的alphaj值
def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


# @dataMat    ：數據列表
# @classLabels：標簽列表
# @C          ：權衡因子（增加松弛因子而在目標優化函數中引入了懲罰項）
# @toler      ：容錯率
# @maxIter    ：最大迭代次數
def smoSimple(dataMat, classLabels, C, toler, maxIter):
    # 將列表形式轉為矩陣或向量形式
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    # 初始化b=0，獲取矩陣行列
    b = 0
    m, n = shape(dataMatrix)
    # 新建一個m行1列的向量
    alphas = mat(zeros((m, 1)))
    # 迭代次數為0
    iter = 0
    while (iter < maxIter):
        # 改變的alpha對數
        alphaPairsChanged = 0
        # 遍歷樣本集中樣本
        for i in range(m):
            # 計算支持向量機算法的預測值
            fXi = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b
            # 計算預測值與實際值的誤差
            Ei = fXi - float(labelMat[i])
            # 如果不滿足KKT條件，即labelMat[i]*fXi<1(labelMat[i]*fXi-1<-toler)
            # and alpha<C 或者labelMat[i]*fXi>1(labelMat[i]*fXi-1>toler)and alpha>0
            if ((labelMat[i] * Ei < -toler) and (alpha < C)) or ((labelMat[i] * Ei > toler) and (alpha[i] > 0)):
                # 隨機選擇第二個變量alphaj
                j = selectJrand(i, m)
                # 計算第二個變量對應數據的預測值
                fXj = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :]).T) + b
                # 計算與測試與實際值的差值
                Ej = fXj - float(label[j])
                # 記錄alphai和alphaj的原始值，便於后續的比較
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                # 如何兩個alpha對應樣本的標簽不相同
                if (labelMat[i] != labelMat[j]):
                    # 求出相應的上下邊界
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                if L == H: print("L==H");continue
                # 根據公式計算未經剪輯的alphaj
                # ------------------------------------------
                eta = 2.0 * dataMatrix[i, :] * dataMatrix[j, :].T - ataMatrix[i, :] * dataMatrix[i, :].T - dataMatrix[j,:] * dataMatrix[j, :].T
                # 如果eta>=0,跳出本次循環
                if eta >= 0: print("eta>=0");continue
                alphas[j] -= labelMat[j] * (Ei - Ej) / eta
                alphas[j] = clipAlpha(alphas[j], H, L)
                # ------------------------------------------
                # 如果改變后的alphaj值變化不大，跳出本次循環
                if (abs(alphas[j] - alphaJold) < 0.00001): print("j not moving enough");continue
                # 否則，計算相應的alphai值
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])
                # 再分別計算兩個alpha情況下對於的b值
                b1 = b - Ei - labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMat[i, :].T - labelMat[j] * (alphas[j] - alphaJold) * dataMatrix[i, :] * dataMatrix[j, :].T
                b2 = b - Ej - labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMatrix[j, :].T - labelMat[j] * (alphas[j] - alphaJold) * dataMatrix[j, :] * dataMatrix[j, :].T
                # 如果0<alphai<C,那么b=b1
                if (0 < alphas[i]) and (C > alphas[i]):
                    b = b1
                # 否則如果0<alphai<C,那么b=b1
                elif (0 < alphas[j]) and (C > alphas[j]):
                    b = b2
                # 否則，alphai，alphaj=0或C
                else:
                    b = (b1 + b2) / 2.0
                # 如果走到此步，表面改變了一對alpha值
                alphaPairsChanged += 1
                print("iter: %d i:%d,paird changed %d" % (iter, i, alphaPairsChanged))
        # 最后判斷是否有改變的alpha對，沒有就進行下一次迭代
        if (alphaPairsChanged == 0):
            iter += 1
        # 否則，迭代次數置0，繼續循環
        else:
            iter = 0
        print("iteration number: %d" % iter)
    # 返回最后的b值和alpha向量
    return b, alphas

#啟發式SMO算法的支持函數
#新建一個類的收據結構，保存當前重要的值
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        self.eCache = mat(zeros((self.m, 2)))


# 格式化計算誤差的函數，方便多次調用
def calcEk(oS, k):
    fXk = float(multiply(oS.alphas, oS.labelMat).T * (oS.X * oS.X[k, :].T)) + oS.b
    Ek = fXk - float(oS.labelMat[k])
    return Ek


# 修改選擇第二個變量alphaj的方法
def selectJ(i, oS, Ei):
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    # 將誤差矩陣每一行第一列置1，以此確定出誤差不為0
    # 的樣本
    oS.eCache[i] = [1, Ei]
    # 獲取緩存中Ei不為0的樣本對應的alpha列表
    validEcacheList = nonzero(oS.Cache[:, 0].A)[0]
    # 在誤差不為0的列表中找出使abs(Ei-Ej)最大的alphaj
    if (len(validEcacheList) > 0):
        for k in validEcacheList:
            if k == i: continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k;
                maxDeltaE = deltaE;
                Ej = Ek
        return maxK, Ej
    else:
        # 否則，就從樣本集中隨機選取alphaj
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej


# 更新誤差矩陣
def updateEk(oS, k):
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]


        # 內循環尋找alphaj

    def innerL(i, oS):
        # 計算誤差
        Ei = calcEk(oS, i)
        # 違背kkt條件
        if (((oS.labelMat[i] * Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i] * Ei > oS.tol) and (oS.alphas[i] > 0))):
            j, Ej = selectJ(i, oS, Ei)
            alphaIold = alphas[i].copy()
            alphaJold = alphas[j].copy()
            # 計算上下界
            if (oS.labelMat[i] != oS.labelMat[j]):
                L = max(0, oS.alphas[j] - oS.alphas[i])
                H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
            else:
                L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
                H = min(oS.C, oS.alphas[j] + oS.alphas[i])
            if L == H: print("L==H");return 0
            # 計算兩個alpha值
            eta = 2.0 * oS.X[i, :] * oS.X[j, :].T - oS.X[i, :] * oS.X[i, :].T -  oS.X[j, :] * oS.X[j, :].T
            if eta >= 0: print("eta>=0");return 0
            oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
            oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
            updateEk(oS, j)
            if (abs(oS.alphas[j] - alphaJold) < 0.00001):
                print("j not moving enough")
                return 0
            oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
            updateEk(oS, i)
            # 在這兩個alpha值情況下，計算對應的b值
            # 注，非線性可分情況，將所有內積項替換為核函數K[i,j]
            b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[i, :].T - oS.labelMat[j] * (oS.alphas[j] - alphaJold) * oS.X[i, :] * oS.X[j, :].T
            b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[j, :].T - oS.labelMat[j] * (oS.alphas[j] - alphaJold) * oS.X[j, :] * oS.X[j, :].T
            if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
                oS.b = b1
            elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]):
                oS.b = b2
            else:
                oS.b = (b1 + b2) / 2.0
            # 如果有alpha對更新
            return 1
            # 否則返回0
        else: return 0

# SMO外循環代碼
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
    # 保存關鍵數據
    oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler)
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    # 選取第一個變量alpha的三種情況，從間隔邊界上選取或者整個數據集
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        # 沒有alpha更新對
        if entireSet:
            for i in range(oS.m):
                alphaPairsChanged += innerL(i, oS)
                print("fullSet,iter: %d i:%d,pairs changed %d" %(iter, i, alphaPairsChanged))
            iter += 1
        else:
            # 統計alphas向量中滿足0<alpha<C的alpha列表
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i, oS)
                print("non-bound,iter: %d i:%d,pairs changed %d"% (iter, i, alphaPairsChanged))
            iter += 1
        if entireSet:
            entireSet = False
            # 如果本次循環沒有改變的alpha對，將entireSet置為true，
            # 下個循環仍遍歷數據集
        elif (alphaPairsChanged == 0):
            entireSet = True
            print("iteration number: %d" % iter)
    return oS.b, oS.alphas

# 求出了alpha值和對應的b值，就可以求出對應的w值，以及分類函數值
def predict(alphas, dataArr, classLabels):
    X = mat(dataArr)
    labelMat = mat(classLabels)
    m, n = shape(X)
    w = zeros((n, 1))
    for i in range(m):
        w += multiply(alphas[i] * labelMat[i], X[i, :].T)
    result = dataArr[0] * mat(ws) + b
    return sign(result)