《機器學習_07_01_svm_硬間隔支援向量機與SMO》
阿新 • • 發佈：2020-05-21
### 一.簡介
支援向量機(svm)的想法與前面介紹的感知機模型類似，找一個超平面將正負樣本分開，但svm的想法要更深入了一步，它要求正負樣本中離超平面最近的點的距離要儘可能的大，所以svm模型建模可以分為兩個子問題：
（1）分的對：怎麼能讓超平面將正負樣本分的開；
（2）分的好：怎麼能讓距離超平面最近的點的距離儘可能的大。
**對於第一個子問題**：將樣本分開，與感知機模型一樣，我們也可以定義模型目標函式為：
$$
f(x)=sign(w^Tx+b)
$$
所以對每對樣本$(x,y)$，只要滿足$y\cdot (w^Tx+b)>0$，即表示模型將樣本正確分開了
**對於第二個子問題**：怎麼能讓離超平面最近的點的距離儘可能的大，對於這個問題，又可以拆解為兩個小問題：
（1）怎麼度量距離？
（2）距離超平面最近的點如何定義？
距離的度量很簡單，可以使用高中時代就知道的點到面的距離公式：
$$
d=\frac{|w^Tx+b|}{||w||}
$$
距離超平面最近的點，我們可以強制定義它為滿足$|w^Tx+b|=1$的點（注意，正負樣本都要滿足），為什麼可以這樣定義呢？我們可以反過來看，一個訓練好的模型可以滿足：（1）要使得正負樣本距離超平面最近的點的距離都儘可能大，那麼這個距離必然要相等，（2）引數$w,b$可以等比例的變化，而不會影響到模型自身，所以$|w^Tx+b|=1$自然也可以滿足，所以這時最近的點的距離可以表示為：
$$
d^*=\frac{1}{||w||}
$$
同時第一個子問題的條件要調整為$y\cdot(w^Tx+b)\geq1$，而$\max d^*$可以等價的表示為$\min \frac{1}{2}w^Tw$，所以svm模型的求解可以表述為如下優化問題：
$$
\min_{w,b} \frac{1}{2}w^Tw \\
s.t.y_i(w^Tx_i+b)\geq 1,i=1,2,...,N
$$
### 二.原優化問題的對偶問題
對於上面優化問題的求解往往轉化為對其對偶問題的求解，首先，構造其拉格朗日函式：
$$
L(w,b,\alpha)=\frac{1}{2}w^Tw+\sum_{i=1}^N \alpha_i(1-y_i(w^Tx_i+b)),\alpha=[\alpha_1,\alpha_2,...,\alpha_N]
$$
這時，原優化問題（設為$P$）就等價於：
$$
\min_{w,b}\max_{\alpha}L(w,b,\alpha)\\
s.t.\alpha_i\geq 0,i=1,2,...,N
$$
這裡簡單說明一下為什麼等價，首先看裡面$\max$那一層
$$\max_{\alpha}L(w,b,\alpha)\\
s.t.\alpha_i\geq 0,i=1,2,...,N$$
對每個樣本都有約束條件$1-y_i(w^Tx_i+b)$，如果滿足約束，即$\leq 0$，必有$\alpha_i(1-y_i(w^Tx_i+b))=0$，如果不滿足，必有$\alpha_i(1-y_i(w^Tx_i+b))\rightarrow 正無窮$，所以，（1）如果所有樣本均滿足約束條件(即$w,b$在可行域內時)，原問題與上面的$\min\max$問題等價，（2）如果有任意一個樣本不滿足約束，這時上面$\max$問題的函式取值為正無窮，外層再對其求$\min$會約束其只能在可行域內求最小值，所以兩問題是等價的，簡單手繪演示一下（兩個問題的最優解都是紅點標記）：
假設對於問題$P$我們求得了最優解$w^*,b^*,\alpha^*$，則必有$L(w^*,b^*,\alpha^*)=L(w^*,b^*,0)$，所以有：
$$
\sum_{i=1}^N\alpha_i^*(1-y_i({w^*}^Tx_i+b^*))=0(條件1)
$$
而最優解自然也滿足原始的約束條件，即：
$$
1-y_i({w^*}^Tx_i+b)\leq0,i=1,2,...,N(條件2)\\
\alpha_i^*\geq0,i=1,2,...,N(條件3)\\
$$
由條件1，2，3，我們可以得出更強地約束條件：
$$
\alpha_i^*(1-y_i({w^*}^Tx_i+b^*))=0,i=1,2,...,N(條件4)
$$
證明也很簡單，由條件2,3可以知道，$\forall i,\alpha_i^*(1-y_i({w^*}^Tx_i+b^*))\leq0$都成立，要使條件1成立，則只能$\alpha_i^*(1-y_i({w^*}^Tx_i+b^*))=0,i=1,2,...,N$。
進一步的，可以推匯出這樣的關係：
$$
\forall \alpha_i^*>0\Rightarrow 1-y_i({w^*}^Tx_i+b^*)=0(關係1)\\
\forall 1-y_i({w^*}^Tx_i+b^*)<0\Rightarrow \alpha_i^*=0(關係2)
$$
所以條件4有個很形象的稱呼：**互補鬆弛條件**，而對於滿足關係1的樣本，也有個稱呼，叫**支援向量**
好的，我們繼續看svm的對偶問題（設為$Q$）的定義：
$$
\max_{\alpha}\min_{w,b}L(w,b,\alpha)\\
s.t.\alpha_i\geq 0,i=1,2,...,N
$$
很幸運，svm的對偶問題$\max\min$與原問題$\min\max$等價（等價是指兩個問題的最優值、最優解$w,b,\alpha$均相等，**具體證明需要用到原問題為凸以及slater條件，可以參看《凸優化》**），先看裡層的$\min_{w,b} L(w,b,\alpha)，$由於$L(w,b,\alpha)$是關於$w,b$的凸函式，所以對偶問題的最優解必然滿足：$L(w,b,\alpha)$關於$w,b$的偏導為0，即：
$$
w=\sum_{i=1}^N\alpha_iy_ix_i(條件5)\\
0=\sum_{i=1}^N\alpha_iy_i(條件6)
$$
消去$w,b$，可得對偶問題關於$\alpha$的表示式：
$$
\max_{\alpha} \sum_{i=1}^N\alpha_i-\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_jx_i^Tx_j\\
s.t.\sum_{i=1}^N\alpha_iy_i=0,\\
\alpha_i\geq0,i=1,2,...,N
$$
顯然，等價於如下優化問題（設為$Q^*$）：
$$
\min_{\alpha} \frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_jx_i^Tx_j-\sum_{i=1}^N\alpha_i\\
s.t.\sum_{i=1}^N\alpha_iy_i=0,\\
\alpha_i\geq0,i=1,2,...,N
$$
該問題是關於$\alpha$的凸二次規劃(QP)問題，可以通過一些優化計算包(比如cvxopt)直接求解最優的$\alpha^*$，再由條件5，可知：
$$
w^*=\sum_{i=1}^N\alpha_i^*y_ix_i
$$
而關於$b^*$，我們可以巧妙求解：找一個樣本點$(x_i,y_i)$，它滿足對應的$\alpha_i^*>0$（即支援向量），利用關係1，可知$1-y_i({w^*}^Tx_i+b^*)=0$，所以：$b^*=y_i-{w^*}^Tx_i$
這裡，條件2,3,4,5,6即是**KKT條件**，而且對於該優化問題，**KKT條件**還是最優解的充分條件（**證明部分...可以參考《凸優化》**），即滿足KKT條件的解就是最優解。
### 三.SMO求解對偶問題最優解
關於對偶問題($Q^*$)可以使用軟體包暴力求解，而且一定能得到最優解，但它的複雜度有點高：（1）變數數與樣本數相同，每個變數$\alpha_i$對應樣本$(x_i,y_i)$；（2）約束條件數也與樣本數相同；而序列最小最優化化(sequential minimal optimization,SMO)演算法是求解SVM對偶問題的一種啟發式演算法，它的思路是：**每次只選擇一個變數優化，而固定住其他變數**，比如選擇$\alpha_1$進行優化，而固定住$\alpha_i,i=2,3,...,N$，但由於我們的問題中有一個約束：$\sum_i^N\alpha_iy_i=0$，需要另外選擇一個$\alpha_2$來配合$\alpha_1$做改變，當兩者中任何一個變數確定後，另外一個也就隨之確定了，比如確定$\alpha_2$後：
$$
\alpha_1=-y_i\sum_{i=2}^N\alpha_iy_i(關係3)
$$
**選擇兩個變數後，如果優化？**
我們在選擇好兩個變數後，如何進行優化呢？比如選擇的$\alpha_1,\alpha_2$，由於剩餘的$\alpha_3,\alpha_4,...,\alpha_N$都視作常量，在$Q^*$中可以忽略，重新整理一下此時的$Q^*$：
$$
\min_{\alpha_1,\alpha_2}\frac{1}{2}\alpha_1^2 x_1^Tx_1+\frac{1}{2}\alpha_2^2x_2^Tx_2+\alpha_1\alpha_2y_1y_2x_1^Tx_2+\frac{1}{2}\alpha_1y_1x_1^T\sum_{i=3}^N\alpha_iy_ix_i+\frac{1}{2}\alpha_2y_2x_2^T\sum_{i=3}^N\alpha_iy_ix_i-\alpha_1-\alpha_2\\
s.t.\alpha_1y_1+\alpha_2y_2=-\sum_{i=3}^N\alpha_iy_i=\eta\\
\alpha_1\geq0,\alpha_2\geq0
$$
這裡求解其實就很簡單了，將關係3帶入，消掉$\alpha_1$後可以發現，優化的目標函式其實是關於$\alpha_2$的二次函式（且開口朝上）：
$$
\min_{\alpha_2}\frac{1}{2}(x_1-x_2)^T(x_1-x_2)\alpha_2^2+(-y_2\eta x_1^Tx_1+y_1\eta x_1^Tx_2+\frac{1}{2}y_2x_2^T\gamma-\frac{1}{2}y_2x_1^T\gamma-1+y_1y_2)\alpha_2\\
s.t.\alpha_2\geq0,y_1(\eta-\alpha_2y_2)\geq0
$$
這裡，$\eta=-\sum_{i=3}^N\alpha_iy_i,\gamma=\sum_{i=3}^N\alpha_iy_ix_i$
所以該問題無約束的最優解為：
$$
\alpha_2^{unc}=-\frac{-y_2\eta x_1^Tx_1+y_1\eta x_1^Tx_2+\frac{1}{2}y_2x_2^T\gamma-\frac{1}{2}y_2x_1^T\gamma-1+y_1y_2}{(x_1-x_2)^T(x_1-x_2)}(公式1)
$$
接下來，我們對上面的表示式做一些優化，大家注意每次迭代時，$\gamma,\eta$都有大量的重複計算（每次僅修改了$\alpha$的兩個變數，剩餘部分其實無需重複計算），而且對於$\alpha_1,\alpha_2$的更新也沒有有效利用它上一階段的取值（記作$\alpha_1^{old},\alpha_2^{old}$）：
我們記：
$$
g(x)=\sum_{i=1}^N\alpha_iy_ix_i^Tx+b
$$
記：
$$
E_i=g(x_i)-y_i
$$
這裡$g(x)$表示模型對$x$的預測值，$E_i$表示預測值與真實值之差，於是我們有：
$$
x_1^T\gamma=g(x_1)-\alpha_1^{old}y_1x_1^Tx_1-\alpha_2^{old}y_2x_2^Tx_1-b^{old}\\
x_2^T\gamma=g(x_2)-\alpha_1^{old}y_1x_1^Tx_2-\alpha_2^{old}y_2x_2^Tx_2-b^{old}
$$
另外：
$$
\eta=\alpha_1^{old}y_1+\alpha_2^{old}y_2
$$
帶入公式1，可得：
$$
\alpha_2^{unc}=\alpha_2^{old}+\frac{y_2(E_1^{old}-E_2^{old})}{\beta}
$$
這裡$\beta=(x_1-x_2)^T(x_1-x_2)$，到這一步，可以發現計算量大大降低，因為$E_1^{old},E_2^{old}$可先快取到記憶體中，但別忘了$\alpha_2$還有約束條件$\alpha_2\geq0,y_1(\eta-\alpha_2y_2)\geq0$，所以需要進一步對它的最優解分情況討論：
當$y_1y_2=1$時，
$$
\alpha_2^{new}=\left\{\begin{matrix}
0 & \alpha_2^{unc}<0\\
\alpha_2^{unc} & 0\leq\alpha_2^{unc}\leq \alpha_1^{old}+\alpha_2^{old}\\
\alpha_1^{old}+\alpha_2^{old} & \alpha_2^{unc}>\alpha_1^{old}+\alpha_2^{old}
\end{matrix}\right.
$$
當$y_1y_2=-1$時，
$$
\alpha_2^{new}=\left\{\begin{matrix}
\alpha_2^{unc} & \alpha_2^{unc}\geq max\{0,\alpha_2^{old}-\alpha_1^{old}\}\\
max\{0,\alpha_2^{old}-\alpha_1^{old}\} & \alpha_2^{unc}< max\{0, \alpha_2^{old}-\alpha_1^{old}\}
\end{matrix}\right.
$$
到這兒，我們可以發現，SMO演算法可以極大的方便$Q^*$的求解，而且是以解析解方式，得到$\alpha_2^{new}$後，由於$\alpha_1^{new}y_1+\alpha_2^{new}y_2=\alpha_1^{old}y_1+\alpha_2^{old}y_2$，可得到$\alpha_1^{new}$的更新公式：
$$
\alpha_1^{new}=\alpha_1^{old}+y_1y_2(\alpha_2^{old}-\alpha_2^{new})
$$
最後，得到$w$的更新公式：
$$
w^{new}=w^{old}+(\alpha_1^{new}-\alpha_1^{old})y_1x_1+(\alpha_2^{new}-\alpha_2^{old})y_2x_2
$$
**對$b$和$E$的更新**
而對於$b$的更新同樣藉助於$\alpha_1,\alpha_2$更新，在更新後，傾向於$\alpha_1^{new}>0,\alpha_2^{new}>0$，還記得前面的互補鬆弛條件吧（條件4），即對於$\alpha_i>0$的情況，必然要有$1-y_i(w^Tx_i+b)=0$成立，即$w^Tx_i+b=y_i$，所以對$(x_1,y_1),(x_2,y_2)$有如下關係：
$$
{w^{new}}^Tx_1+b=y_1(關係4)\\
{w^{new}}^Tx_2+b=y_2(關係5)\\
$$
對關係4和關係5可以分別計算出$b_1^{new}=y_1-{w^{new}}^Tx_1,b_2^{new}=y_2-{w^{new}}^Tx_2$，對$b$的更新，可以取兩者的均值：
$$
b^{new}=\frac{b_1^{new}+b_2^{new}}{2}
$$
接下來，對於$E_1,E_2$的更新就很自然了：
$$
E_1^{new}={w^{new}}^Tx_1+b^{new}-y_1\\
E_2^{new}={w^{new}}^Tx_2+b^{new}-y_2
$$
那接下來還有一個問題，那就是$\alpha_1,\alpha_2$如何選擇的問題
**如何選擇兩個優化變數？**
這可以啟發式選擇，分為兩步：第一步是如何選擇$\alpha_1$，第二步是在選定$\alpha_1$時，如何選擇一個不錯的$\alpha_2$：
**$\alpha_1$的選擇**
選擇$\alpha_1$同感知機模型類似，選擇一個不滿足KKT條件的點$(x_i,y_i)$，即不滿足如下兩種情況之一的點：
$$
\left\{\begin{matrix}
\alpha_i=0\Leftrightarrow y_i(w^Tx_i+b)\geq1\\
\alpha_i>0\Leftrightarrow y_i(w^Tx_i+b)=1
\end{matrix}\right.
$$
**$\alpha_2$的選擇**
對$\alpha_2$的選擇傾向於選擇使其變化儘可能大的點，由前面的更新公式可知是使$\mid E_1^{old}-E_2^{old}\mid$最大的點，所以選擇的兩個點$(x_1,y_1)$和$(x_2,y_2)$會更傾向於選擇異類的點
### 四.程式碼實現
```python
import numpy as np
import matplotlib.pyplot as plt
import copy
import random
import os
os.chdir('../')
from ml_models import utils
%matplotlib inline
```
```python
#定義一個繪製決策邊界以及支援向量的函式（並放到utils中）
def plot_decision_function(X, y, clf, support_vectors=None):
plot_step = 0.02
x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, plot_step),
np.arange(y_min, y_max, plot_step))
Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
plt.contourf(xx, yy, Z, alpha=0.4)
plt.scatter(X[:, 0], X[:, 1], alpha=0.8, c=y, edgecolor='k')
# 繪製支援向量
if support_vectors is not None:
plt.scatter(X[support_vectors, 0], X[support_vectors, 1], s=80, c='none', alpha=0.7, edgecolor='red')
```
```python
"""
硬間隔支援向量機的smo實現,放到ml_models.svm模組
"""
class HardMarginSVM(object):
def __init__(self, epochs=100):
self.w = None
self.b = None
self.alpha = None
self.E = None
self.epochs = epochs
# 記錄支援向量
self.support_vectors = None
def init_params(self, X, y):
"""
:param X: (n_samples,n_features)
:param y: (n_samples,) y_i\in\{0,1\}
:return:
"""
n_samples, n_features = X.shape
self.w = np.zeros(n_features)
self.b = .0
self.alpha = np.zeros(n_samples)
self.E = np.zeros(n_samples)
# 初始化E
for i in range(0, n_samples):
self.E[i] = np.dot(self.w, X[i, :]) + self.b - y[i]
def _select_j(self, best_i):
"""
選擇j
:param best_i:
:return:
"""
valid_j_list = [i for i in range(0, len(self.alpha)) if self.alpha[i] > 0 and i != best_i]
best_j = -1
# 優先選擇使得|E_i-E_j|最大的j
if len(valid_j_list) > 0:
max_e = 0
for j in valid_j_list:
current_e = np.abs(self.E[best_i] - self.E[j])
if current_e > max_e:
best_j = j
max_e = current_e
else:
# 隨機選擇
l = list(range(len(self.alpha)))
seq = l[: best_i] + l[best_i + 1:]
best_j = random.choice(seq)
return best_j
def _meet_kkt(self, w, b, x_i, y_i, alpha_i):
"""
判斷是否滿足KKT條件
:param w:
:param b:
:param x_i:
:param y_i:
:return:
"""
if alpha_i < 1e-7:
return y_i * (np.dot(w, x_i) + b) >= 1
else:
return abs(y_i * (np.dot(w, x_i) + b) - 1) < 1e-7
def fit(self, X, y2, show_train_process=False):
"""
:param X:
:param y2:
:param show_train_process: 顯示訓練過程
:return:
"""
y = copy.deepcopy(y2)
y[y == 0] = -1
# 初始化引數
self.init_params(X, y)
for _ in range(0, self.epochs):
if_all_match_kkt = True
for i in range(0, len(self.alpha)):
x_i = X[i, :]
y_i = y[i]
alpha_i_old = self.alpha[i]
E_i_old = self.E[i]
# 外層迴圈：選擇違反KKT條件的點i
if not self._meet_kkt(self.w, self.b, x_i, y_i, alpha_i_old):
if_all_match_kkt = False
# 內層迴圈，選擇使|Ei-Ej|最大的點j
best_j = self._select_j(i)
alpha_j_old = self.alpha[best_j]
x_j = X[best_j, :]
y_j = y[best_j]
E_j_old = self.E[best_j]
# 進行更新
# 1.首先獲取無裁剪的最優alpha_2
eta = np.dot(x_i - x_j, x_i - x_j)
# 如果x_i和x_j很接近，則跳過
if eta < 1e-3:
continue
alpha_j_unc = alpha_j_old + y_j * (E_i_old - E_j_old) / eta
# 2.裁剪並得到new alpha_2
if y_i == y_j:
if alpha_j_unc < 0:
alpha_j_new = 0
elif 0 <= alpha_j_unc <= alpha_i_old + alpha_j_old:
alpha_j_new = alpha_j_unc
else:
alpha_j_new = alpha_i_old + alpha_j_old
else:
if alpha_j_unc < max(0, alpha_j_old - alpha_i_old):
alpha_j_new = max(0, alpha_j_old - alpha_i_old)
else:
alpha_j_new = alpha_j_unc
# 如果變化不夠大則跳過
if np.abs(alpha_j_new - alpha_j_old) < 1e-5:
continue
# 3.得到alpha_1_new
alpha_i_new = alpha_i_old + y_i * y_j * (alpha_j_old - alpha_j_new)
# 4.更新w
self.w = self.w + (alpha_i_new - alpha_i_old) * y_i * x_i + (alpha_j_new - alpha_j_old) * y_j * x_j
# 5.更新alpha_1,alpha_2
self.alpha[i] = alpha_i_new
self.alpha[best_j] = alpha_j_new
# 6.更新b
b_i_new = y_i - np.dot(self.w, x_i)
b_j_new = y_j - np.dot(self.w, x_j)
if alpha_i_new > 0:
self.b = b_i_new
elif alpha_j_new > 0:
self.b = b_j_new
else:
self.b = (b_i_new + b_j_new) / 2.0
# 7.更新E
for k in range(0, len(self.E)):
self.E[k] = np.dot(self.w, X[k, :]) + self.b - y[k]
# 顯示訓練過程
if show_train_process is True:
utils.plot_decision_function(X, y2, self, [i, best_j])
utils.plt.pause(0.1)
utils.plt.clf()
# 如果所有的點都滿足KKT條件，則中止
if if_all_match_kkt is True:
break
# 計算支援向量
self.support_vectors = np.where(self.alpha > 1e-3)[0]
# 利用所有的支援向量，更新b
self.b = np.mean([y[s] - np.dot(self.w, X[s, :]) for s in self.support_vectors.tolist()])
# 顯示最終結果
if show_train_process is True:
utils.plot_decision_function(X, y2, self, self.support_vectors)
utils.plt.show()
def get_params(self):
"""
輸出原始的係數
:return: w
"""
return self.w, self.b
def predict_proba(self, x):
"""
:param x:ndarray格式資料: m x n
:return: m x 1
"""
return utils.sigmoid(x.dot(self.w) + self.b)
def predict(self, x):
"""
:param x:ndarray格式資料: m x n
:return: m x 1
"""
proba = self.predict_proba(x)
return (proba >= 0.5).astype(int)
```
### 五.效果檢驗
```python
from sklearn.datasets import make_classification
# 生成分類資料
data, target = make_classification(n_samples=100, n_features=2, n_classes=2, n_informative=1, n_redundant=0,
n_repeated=0, n_clusters_per_class=1, class_sep=2.0)
```
```python
plt.scatter(data[:,0],dat,1],c=target)
```
假設對於問題$P$我們求得了最優解$w^*,b^*,\alpha^*$，則必有$L(w^*,b^*,\alpha^*)=L(w^*,b^*,0)$，所以有：
$$
\sum_{i=1}^N\alpha_i^*(1-y_i({w^*}^Tx_i+b^*))=0(條件1)
$$
而最優解自然也滿足原始的約束條件，即：
$$
1-y_i({w^*}^Tx_i+b)\leq0,i=1,2,...,N(條件2)\\
\alpha_i^*\geq0,i=1,2,...,N(條件3)\\
$$
由條件1，2，3，我們可以得出更強地約束條件：
$$
\alpha_i^*(1-y_i({w^*}^Tx_i+b^*))=0,i=1,2,...,N(條件4)
$$
證明也很簡單，由條件2,3可以知道，$\forall i,\alpha_i^*(1-y_i({w^*}^Tx_i+b^*))\leq0$都成立，要使條件1成立，則只能$\alpha_i^*(1-y_i({w^*}^Tx_i+b^*))=0,i=1,2,...,N$。
進一步的，可以推匯出這樣的關係：
$$
\forall \alpha_i^*>0\Rightarrow 1-y_i({w^*}^Tx_i+b^*)=0(關係1)\\
\forall 1-y_i({w^*}^Tx_i+b^*)<0\Rightarrow \alpha_i^*=0(關係2)
$$
所以條件4有個很形象的稱呼：**互補鬆弛條件**，而對於滿足關係1的樣本，也有個稱呼，叫**支援向量**
好的，我們繼續看svm的對偶問題（設為$Q$）的定義：
$$
\max_{\alpha}\min_{w,b}L(w,b,\alpha)\\
s.t.\alpha_i\geq 0,i=1,2,...,N
$$
很幸運，svm的對偶問題$\max\min$與原問題$\min\max$等價（等價是指兩個問題的最優值、最優解$w,b,\alpha$均相等，**具體證明需要用到原問題為凸以及slater條件，可以參看《凸優化》**），先看裡層的$\min_{w,b} L(w,b,\alpha)，$由於$L(w,b,\alpha)$是關於$w,b$的凸函式，所以對偶問題的最優解必然滿足：$L(w,b,\alpha)$關於$w,b$的偏導為0，即：
$$
w=\sum_{i=1}^N\alpha_iy_ix_i(條件5)\\
0=\sum_{i=1}^N\alpha_iy_i(條件6)
$$
消去$w,b$，可得對偶問題關於$\alpha$的表示式：
$$
\max_{\alpha} \sum_{i=1}^N\alpha_i-\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_jx_i^Tx_j\\
s.t.\sum_{i=1}^N\alpha_iy_i=0,\\
\alpha_i\geq0,i=1,2,...,N
$$
顯然，等價於如下優化問題（設為$Q^*$）：
$$
\min_{\alpha} \frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_jx_i^Tx_j-\sum_{i=1}^N\alpha_i\\
s.t.\sum_{i=1}^N\alpha_iy_i=0,\\
\alpha_i\geq0,i=1,2,...,N
$$
該問題是關於$\alpha$的凸二次規劃(QP)問題，可以通過一些優化計算包(比如cvxopt)直接求解最優的$\alpha^*$，再由條件5，可知：
$$
w^*=\sum_{i=1}^N\alpha_i^*y_ix_i
$$
而關於$b^*$，我們可以巧妙求解：找一個樣本點$(x_i,y_i)$，它滿足對應的$\alpha_i^*>0$（即支援向量），利用關係1，可知$1-y_i({w^*}^Tx_i+b^*)=0$，所以：$b^*=y_i-{w^*}^Tx_i$
這裡，條件2,3,4,5,6即是**KKT條件**，而且對於該優化問題，**KKT條件**還是最優解的充分條件（**證明部分...可以參考《凸優化》**），即滿足KKT條件的解就是最優解。
### 三.SMO求解對偶問題最優解
關於對偶問題($Q^*$)可以使用軟體包暴力求解，而且一定能得到最優解，但它的複雜度有點高：（1）變數數與樣本數相同，每個變數$\alpha_i$對應樣本$(x_i,y_i)$；（2）約束條件數也與樣本數相同；而序列最小最優化化(sequential minimal optimization,SMO)演算法是求解SVM對偶問題的一種啟發式演算法，它的思路是：**每次只選擇一個變數優化，而固定住其他變數**，比如選擇$\alpha_1$進行優化，而固定住$\alpha_i,i=2,3,...,N$，但由於我們的問題中有一個約束：$\sum_i^N\alpha_iy_i=0$，需要另外選擇一個$\alpha_2$來配合$\alpha_1$做改變，當兩者中任何一個變數確定後，另外一個也就隨之確定了，比如確定$\alpha_2$後：
$$
\alpha_1=-y_i\sum_{i=2}^N\alpha_iy_i(關係3)
$$
**選擇兩個變數後，如果優化？**
我們在選擇好兩個變數後，如何進行優化呢？比如選擇的$\alpha_1,\alpha_2$，由於剩餘的$\alpha_3,\alpha_4,...,\alpha_N$都視作常量，在$Q^*$中可以忽略，重新整理一下此時的$Q^*$：
$$
\min_{\alpha_1,\alpha_2}\frac{1}{2}\alpha_1^2 x_1^Tx_1+\frac{1}{2}\alpha_2^2x_2^Tx_2+\alpha_1\alpha_2y_1y_2x_1^Tx_2+\frac{1}{2}\alpha_1y_1x_1^T\sum_{i=3}^N\alpha_iy_ix_i+\frac{1}{2}\alpha_2y_2x_2^T\sum_{i=3}^N\alpha_iy_ix_i-\alpha_1-\alpha_2\\
s.t.\alpha_1y_1+\alpha_2y_2=-\sum_{i=3}^N\alpha_iy_i=\eta\\
\alpha_1\geq0,\alpha_2\geq0
$$
這裡求解其實就很簡單了，將關係3帶入，消掉$\alpha_1$後可以發現，優化的目標函式其實是關於$\alpha_2$的二次函式（且開口朝上）：
$$
\min_{\alpha_2}\frac{1}{2}(x_1-x_2)^T(x_1-x_2)\alpha_2^2+(-y_2\eta x_1^Tx_1+y_1\eta x_1^Tx_2+\frac{1}{2}y_2x_2^T\gamma-\frac{1}{2}y_2x_1^T\gamma-1+y_1y_2)\alpha_2\\
s.t.\alpha_2\geq0,y_1(\eta-\alpha_2y_2)\geq0
$$
這裡，$\eta=-\sum_{i=3}^N\alpha_iy_i,\gamma=\sum_{i=3}^N\alpha_iy_ix_i$
所以該問題無約束的最優解為：
$$
\alpha_2^{unc}=-\frac{-y_2\eta x_1^Tx_1+y_1\eta x_1^Tx_2+\frac{1}{2}y_2x_2^T\gamma-\frac{1}{2}y_2x_1^T\gamma-1+y_1y_2}{(x_1-x_2)^T(x_1-x_2)}(公式1)
$$
接下來，我們對上面的表示式做一些優化，大家注意每次迭代時，$\gamma,\eta$都有大量的重複計算（每次僅修改了$\alpha$的兩個變數，剩餘部分其實無需重複計算），而且對於$\alpha_1,\alpha_2$的更新也沒有有效利用它上一階段的取值（記作$\alpha_1^{old},\alpha_2^{old}$）：
我們記：
$$
g(x)=\sum_{i=1}^N\alpha_iy_ix_i^Tx+b
$$
記：
$$
E_i=g(x_i)-y_i
$$
這裡$g(x)$表示模型對$x$的預測值，$E_i$表示預測值與真實值之差，於是我們有：
$$
x_1^T\gamma=g(x_1)-\alpha_1^{old}y_1x_1^Tx_1-\alpha_2^{old}y_2x_2^Tx_1-b^{old}\\
x_2^T\gamma=g(x_2)-\alpha_1^{old}y_1x_1^Tx_2-\alpha_2^{old}y_2x_2^Tx_2-b^{old}
$$
另外：
$$
\eta=\alpha_1^{old}y_1+\alpha_2^{old}y_2
$$
帶入公式1，可得：
$$
\alpha_2^{unc}=\alpha_2^{old}+\frac{y_2(E_1^{old}-E_2^{old})}{\beta}
$$
這裡$\beta=(x_1-x_2)^T(x_1-x_2)$，到這一步，可以發現計算量大大降低，因為$E_1^{old},E_2^{old}$可先快取到記憶體中，但別忘了$\alpha_2$還有約束條件$\alpha_2\geq0,y_1(\eta-\alpha_2y_2)\geq0$，所以需要進一步對它的最優解分情況討論：
當$y_1y_2=1$時，
$$
\alpha_2^{new}=\left\{\begin{matrix}
0 & \alpha_2^{unc}<0\\
\alpha_2^{unc} & 0\leq\alpha_2^{unc}\leq \alpha_1^{old}+\alpha_2^{old}\\
\alpha_1^{old}+\alpha_2^{old} & \alpha_2^{unc}>\alpha_1^{old}+\alpha_2^{old}
\end{matrix}\right.
$$
當$y_1y_2=-1$時，
$$
\alpha_2^{new}=\left\{\begin{matrix}
\alpha_2^{unc} & \alpha_2^{unc}\geq max\{0,\alpha_2^{old}-\alpha_1^{old}\}\\
max\{0,\alpha_2^{old}-\alpha_1^{old}\} & \alpha_2^{unc}< max\{0, \alpha_2^{old}-\alpha_1^{old}\}
\end{matrix}\right.
$$
到這兒，我們可以發現，SMO演算法可以極大的方便$Q^*$的求解，而且是以解析解方式，得到$\alpha_2^{new}$後，由於$\alpha_1^{new}y_1+\alpha_2^{new}y_2=\alpha_1^{old}y_1+\alpha_2^{old}y_2$，可得到$\alpha_1^{new}$的更新公式：
$$
\alpha_1^{new}=\alpha_1^{old}+y_1y_2(\alpha_2^{old}-\alpha_2^{new})
$$
最後，得到$w$的更新公式：
$$
w^{new}=w^{old}+(\alpha_1^{new}-\alpha_1^{old})y_1x_1+(\alpha_2^{new}-\alpha_2^{old})y_2x_2
$$
**對$b$和$E$的更新**
而對於$b$的更新同樣藉助於$\alpha_1,\alpha_2$更新，在更新後，傾向於$\alpha_1^{new}>0,\alpha_2^{new}>0$，還記得前面的互補鬆弛條件吧（條件4），即對於$\alpha_i>0$的情況，必然要有$1-y_i(w^Tx_i+b)=0$成立，即$w^Tx_i+b=y_i$，所以對$(x_1,y_1),(x_2,y_2)$有如下關係：
$$
{w^{new}}^Tx_1+b=y_1(關係4)\\
{w^{new}}^Tx_2+b=y_2(關係5)\\
$$
對關係4和關係5可以分別計算出$b_1^{new}=y_1-{w^{new}}^Tx_1,b_2^{new}=y_2-{w^{new}}^Tx_2$，對$b$的更新，可以取兩者的均值：
$$
b^{new}=\frac{b_1^{new}+b_2^{new}}{2}
$$
接下來，對於$E_1,E_2$的更新就很自然了：
$$
E_1^{new}={w^{new}}^Tx_1+b^{new}-y_1\\
E_2^{new}={w^{new}}^Tx_2+b^{new}-y_2
$$
那接下來還有一個問題，那就是$\alpha_1,\alpha_2$如何選擇的問題
**如何選擇兩個優化變數？**
這可以啟發式選擇，分為兩步：第一步是如何選擇$\alpha_1$，第二步是在選定$\alpha_1$時，如何選擇一個不錯的$\alpha_2$：
**$\alpha_1$的選擇**
選擇$\alpha_1$同感知機模型類似，選擇一個不滿足KKT條件的點$(x_i,y_i)$，即不滿足如下兩種情況之一的點：
$$
\left\{\begin{matrix}
\alpha_i=0\Leftrightarrow y_i(w^Tx_i+b)\geq1\\
\alpha_i>0\Leftrightarrow y_i(w^Tx_i+b)=1
\end{matrix}\right.
$$
**$\alpha_2$的選擇**
對$\alpha_2$的選擇傾向於選擇使其變化儘可能大的點，由前面的更新公式可知是使$\mid E_1^{old}-E_2^{old}\mid$最大的點，所以選擇的兩個點$(x_1,y_1)$和$(x_2,y_2)$會更傾向於選擇異類的點
### 四.程式碼實現
```python
import numpy as np
import matplotlib.pyplot as plt
import copy
import random
import os
os.chdir('../')
from ml_models import utils
%matplotlib inline
```
```python
#定義一個繪製決策邊界以及支援向量的函式（並放到utils中）
def plot_decision_function(X, y, clf, support_vectors=None):
plot_step = 0.02
x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, plot_step),
np.arange(y_min, y_max, plot_step))
Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
plt.contourf(xx, yy, Z, alpha=0.4)
plt.scatter(X[:, 0], X[:, 1], alpha=0.8, c=y, edgecolor='k')
# 繪製支援向量
if support_vectors is not None:
plt.scatter(X[support_vectors, 0], X[support_vectors, 1], s=80, c='none', alpha=0.7, edgecolor='red')
```
```python
"""
硬間隔支援向量機的smo實現,放到ml_models.svm模組
"""
class HardMarginSVM(object):
def __init__(self, epochs=100):
self.w = None
self.b = None
self.alpha = None
self.E = None
self.epochs = epochs
# 記錄支援向量
self.support_vectors = None
def init_params(self, X, y):
"""
:param X: (n_samples,n_features)
:param y: (n_samples,) y_i\in\{0,1\}
:return:
"""
n_samples, n_features = X.shape
self.w = np.zeros(n_features)
self.b = .0
self.alpha = np.zeros(n_samples)
self.E = np.zeros(n_samples)
# 初始化E
for i in range(0, n_samples):
self.E[i] = np.dot(self.w, X[i, :]) + self.b - y[i]
def _select_j(self, best_i):
"""
選擇j
:param best_i:
:return:
"""
valid_j_list = [i for i in range(0, len(self.alpha)) if self.alpha[i] > 0 and i != best_i]
best_j = -1
# 優先選擇使得|E_i-E_j|最大的j
if len(valid_j_list) > 0:
max_e = 0
for j in valid_j_list:
current_e = np.abs(self.E[best_i] - self.E[j])
if current_e > max_e:
best_j = j
max_e = current_e
else:
# 隨機選擇
l = list(range(len(self.alpha)))
seq = l[: best_i] + l[best_i + 1:]
best_j = random.choice(seq)
return best_j
def _meet_kkt(self, w, b, x_i, y_i, alpha_i):
"""
判斷是否滿足KKT條件
:param w:
:param b:
:param x_i:
:param y_i:
:return:
"""
if alpha_i < 1e-7:
return y_i * (np.dot(w, x_i) + b) >= 1
else:
return abs(y_i * (np.dot(w, x_i) + b) - 1) < 1e-7
def fit(self, X, y2, show_train_process=False):
"""
:param X:
:param y2:
:param show_train_process: 顯示訓練過程
:return:
"""
y = copy.deepcopy(y2)
y[y == 0] = -1
# 初始化引數
self.init_params(X, y)
for _ in range(0, self.epochs):
if_all_match_kkt = True
for i in range(0, len(self.alpha)):
x_i = X[i, :]
y_i = y[i]
alpha_i_old = self.alpha[i]
E_i_old = self.E[i]
# 外層迴圈：選擇違反KKT條件的點i
if not self._meet_kkt(self.w, self.b, x_i, y_i, alpha_i_old):
if_all_match_kkt = False
# 內層迴圈，選擇使|Ei-Ej|最大的點j
best_j = self._select_j(i)
alpha_j_old = self.alpha[best_j]
x_j = X[best_j, :]
y_j = y[best_j]
E_j_old = self.E[best_j]
# 進行更新
# 1.首先獲取無裁剪的最優alpha_2
eta = np.dot(x_i - x_j, x_i - x_j)
# 如果x_i和x_j很接近，則跳過
if eta < 1e-3:
continue
alpha_j_unc = alpha_j_old + y_j * (E_i_old - E_j_old) / eta
# 2.裁剪並得到new alpha_2
if y_i == y_j:
if alpha_j_unc < 0:
alpha_j_new = 0
elif 0 <= alpha_j_unc <= alpha_i_old + alpha_j_old:
alpha_j_new = alpha_j_unc
else:
alpha_j_new = alpha_i_old + alpha_j_old
else:
if alpha_j_unc < max(0, alpha_j_old - alpha_i_old):
alpha_j_new = max(0, alpha_j_old - alpha_i_old)
else:
alpha_j_new = alpha_j_unc
# 如果變化不夠大則跳過
if np.abs(alpha_j_new - alpha_j_old) < 1e-5:
continue
# 3.得到alpha_1_new
alpha_i_new = alpha_i_old + y_i * y_j * (alpha_j_old - alpha_j_new)
# 4.更新w
self.w = self.w + (alpha_i_new - alpha_i_old) * y_i * x_i + (alpha_j_new - alpha_j_old) * y_j * x_j
# 5.更新alpha_1,alpha_2
self.alpha[i] = alpha_i_new
self.alpha[best_j] = alpha_j_new
# 6.更新b
b_i_new = y_i - np.dot(self.w, x_i)
b_j_new = y_j - np.dot(self.w, x_j)
if alpha_i_new > 0:
self.b = b_i_new
elif alpha_j_new > 0:
self.b = b_j_new
else:
self.b = (b_i_new + b_j_new) / 2.0
# 7.更新E
for k in range(0, len(self.E)):
self.E[k] = np.dot(self.w, X[k, :]) + self.b - y[k]
# 顯示訓練過程
if show_train_process is True:
utils.plot_decision_function(X, y2, self, [i, best_j])
utils.plt.pause(0.1)
utils.plt.clf()
# 如果所有的點都滿足KKT條件，則中止
if if_all_match_kkt is True:
break
# 計算支援向量
self.support_vectors = np.where(self.alpha > 1e-3)[0]
# 利用所有的支援向量，更新b
self.b = np.mean([y[s] - np.dot(self.w, X[s, :]) for s in self.support_vectors.tolist()])
# 顯示最終結果
if show_train_process is True:
utils.plot_decision_function(X, y2, self, self.support_vectors)
utils.plt.show()
def get_params(self):
"""
輸出原始的係數
:return: w
"""
return self.w, self.b
def predict_proba(self, x):
"""
:param x:ndarray格式資料: m x n
:return: m x 1
"""
return utils.sigmoid(x.dot(self.w) + self.b)
def predict(self, x):
"""
:param x:ndarray格式資料: m x n
:return: m x 1
"""
proba = self.predict_proba(x)
return (proba >= 0.5).astype(int)
```
### 五.效果檢驗
```python
from sklearn.datasets import make_classification
# 生成分類資料
data, target = make_classification(n_samples=100, n_features=2, n_classes=2, n_informative=1, n_redundant=0,
n_repeated=0, n_clusters_per_class=1, class_sep=2.0)
```
```python
plt.scatter(data[:,0],dat,1],c=target)
```