支持向量機(svm)的想法與前面介紹的感知機模型相似,找一個超平面將正負樣本分開,但svm的想法要更深刻了一步,它要求正負樣本中離超平面最近的點的距離要儘量的大,因此svm模型建模能夠分爲兩個子問題:python
(1)分的對:怎麼能讓超平面將正負樣本分的開;
(2)分的好:怎麼能讓距離超平面最近的點的距離儘量的大。算法
對於第一個子問題:將樣本分開,與感知機模型同樣,咱們也能夠定義模型目標函數爲:緩存
因此對每對樣本\((x,y)\),只要知足\(y\cdot (w^Tx+b)>0\),即表示模型將樣本正確分開了app
對於第二個子問題:怎麼能讓離超平面最近的點的距離儘量的大,對於這個問題,又能夠拆解爲兩個小問題:dom
(1)怎麼度量距離?
(2)距離超平面最近的點如何定義?函數
距離的度量很簡單,可使用高中時代就知道的點到面的距離公式:優化
距離超平面最近的點,咱們能夠強制定義它爲知足\(|w^Tx+b|=1\)的點(注意,正負樣本都要知足),爲何能夠這樣定義呢?咱們能夠反過來看,一個訓練好的模型能夠知足:(1)要使得正負樣本距離超平面最近的點的距離都儘量大,那麼這個距離必然要相等,(2)參數\(w,b\)能夠等比例的變化,而不會影響到模型自身,因此\(|w^Tx+b|=1\)天然也能夠知足,因此這時最近的點的距離能夠表示爲:spa
同時第一個子問題的條件要調整爲\(y\cdot(w^Tx+b)\geq1\),而\(\max d^*\)能夠等價的表示爲\(\min \frac{1}{2}w^Tw\),因此svm模型的求解能夠表述爲以下優化問題:3d
對於上面優化問題的求解每每轉化爲對其對偶問題的求解,首先,構造其拉格朗日函數:code
這時,原優化問題(設爲\(P\))就等價於:
這裏簡單說明一下爲何等價,首先看裏面\(\max\)那一層
對每一個樣本都有約束條件\(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)\),因此有:
而最優解天然也知足原始的約束條件,即:
由條件1,2,3,咱們能夠得出更強地約束條件:
證實也很簡單,由條件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\)。
進一步的,能夠推導出這樣的關係:
因此條件4有個很形象的稱呼:互補鬆弛條件,而對於知足關係1的樣本,也有個稱呼,叫支持向量
好的,咱們繼續看svm的對偶問題(設爲\(Q\))的定義:
很幸運,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,b\),可得對偶問題關於\(\alpha\)的表達式:
顯然,等價於以下優化問題(設爲\(Q^*\)):
該問題是關於\(\alpha\)的凸二次規劃(QP)問題,能夠經過一些優化計算包(好比cvxopt)直接求解最優的\(\alpha^*\),再由條件5,可知:
而關於\(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條件的解就是最優解。
關於對偶問題(\(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,\alpha_2\),因爲剩餘的\(\alpha_3,\alpha_4,...,\alpha_N\)都視做常量,在\(Q^*\)中能夠忽略,從新整理一下此時的\(Q^*\):
這裏求解其實就很簡單了,將關係3帶入,消掉\(\alpha_1\)後能夠發現,優化的目標函數實際上是關於\(\alpha_2\)的二次函數(且開口朝上):
這裏,\(\eta=-\sum_{i=3}^N\alpha_iy_i,\gamma=\sum_{i=3}^N\alpha_iy_ix_i\)
因此該問題無約束的最優解爲:
接下來,咱們對上面的表達式作一些優化,你們注意每次迭代時,\(\gamma,\eta\)都有大量的重複計算(每次僅修改了\(\alpha\)的兩個變量,剩餘部分其實無需重複計算),並且對於\(\alpha_1,\alpha_2\)的更新也沒有有效利用它上一階段的取值(記做\(\alpha_1^{old},\alpha_2^{old}\)):
咱們記:
記:
這裏\(g(x)\)表示模型對\(x\)的預測值,\(E_i\)表示預測值與真實值之差,因而咱們有:
另外:
帶入公式1,可得:
這裏\(\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\)時,
當\(y_1y_2=-1\)時,
到這兒,咱們能夠發現,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}\)的更新公式:
最後,獲得\(w\)的更新公式:
對\(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)\)有以下關係:
對關係4和關係5能夠分別計算出\(b_1^{new}=y_1-{w^{new}}^Tx_1,b_2^{new}=y_2-{w^{new}}^Tx_2\),對\(b\)的更新,能夠取二者的均值:
接下來,對於\(E_1,E_2\)的更新就很天然了:
那接下來還有一個問題,那就是\(\alpha_1,\alpha_2\)如何選擇的問題
如何選擇兩個優化變量?
這能夠啓發式選擇,分爲兩步:第一步是如何選擇\(\alpha_1\),第二步是在選定\(\alpha_1\)時,如何選擇一個不錯的\(\alpha_2\):
\(\alpha_1\)的選擇
選擇\(\alpha_1\)同感知機模型相似,選擇一個不知足KKT條件的點\((x_i,y_i)\),即不知足以下兩種狀況之一的點:
\(\alpha_2\)的選擇
對\(\alpha_2\)的選擇傾向於選擇使其變化儘量大的點,由前面的更新公式可知是使\(\mid E_1^{old}-E_2^{old}\mid\)最大的點,因此選擇的兩個點\((x_1,y_1)\)和\((x_2,y_2)\)會更傾向於選擇異類的點
import numpy as np import matplotlib.pyplot as plt import copy import random import os os.chdir('../') from ml_models import utils %matplotlib inline
#定義一個繪製決策邊界以及支持向量的函數(並放到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')
""" 硬間隔支持向量機的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)
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)
plt.scatter(data[:,0],dat,1],c=target)
<matplotlib.collections.PathCollection at 0x19e700bb390>
#訓練 svm = HardMarginSVM() svm.fit(data, target)
plot_decision_function(data, target, svm, svm.support_vectors)
#可視化訓練過程,建議在pycharm中運行,notebook會生成不少張圖片 # svm = HardMarginSVM() # svm.fit(data, target,show_train_process=True)
你們能夠將上面的代碼多運行幾回,能夠發現若是有異常點等狀況出現時(即線性不可分時),模型訓練的結果會很難看,後面小節將會對這種狀況作處理,教模型如何去「容忍」這些很差看的點,或者巧妙地經過座標映射的方式將低維數據映射到高維空間進而能夠線性可分
我的以爲更可能是爲了引入核技巧,由於對偶問題進行計算時,有關於兩個點內積的計算:\(x_i^Tx_j\),這能夠方便的用核函數替代\(\kappa(x_i,x_j)\),便於處理非線性可分的狀況