SVM筆記

1.前言

SVM（Support Vector Machine）是一種尋求最大分類間隔的機器學習方法，廣泛應用於各個領域，許多人把SVM當做首選方法，它也被稱之為最優分類器，這是為什么呢？這篇文章將系統介紹SVM的原理、推導過程及代碼實踐。

2.初識SVM

首先我們先來看看SVM做的是什么樣的事，我們先看下面一張圖

圖中有三個分類實例，都將數據正確分類，我們直觀上看，會覺得圖中第三個效果會比較好，這是為什么呢？個人覺得人的直觀感受更偏向於數據均勻對稱的結構。當然，這只是直觀感受，我們從專業的角度出發，第三種情況更優的原因在於它能容忍更多的數據噪聲（tolerate more data noise），從而這個系統更加穩健（Robust），如下圖所示，灰色的圈圈代表容忍噪聲的大小門限，第三個門限最大，它比前面兩個更優，而SVM就是尋求最大門限，構成分離數據正確且抗噪聲強的穩健超平面（圖中分類線是直線，一般數據維度是高維，對應的分類平面就是超平面）。

3.線性SVM

好了，現在我們知道SVM是那個抗噪聲門限最大的超平面（一般說的是間隔margin最大，意思一樣），那我們怎樣找到這個超平面呢？我們先來分析，這個超平面需要滿足那些條件：

• (1) 需要滿足將數據都分類正確，以上圖為例，超平面為\(w^Tx+b=0\),超平面以上的數據Label為+1，超平面以下的數據Label為-1，則分類正確需要滿足\(y_n(w^Tx+b)>0\) ( \(y_n\)是數據Label )
• (2) 這個超平面關於數據集的最大間隔等於超平面到所有數據集中的樣本點的距離的最小值，即\(margin(w,b)=\underset{n=1,\cdots,N}{min}distance(x_n,w,b)\)

根據上面的條件，我們可以對SVM建立一個最優化模型，目標是最大化間隔，但是這個最優化模型是有限制條件的，即上面提的兩個條件.

\[\underset{w,b}{max} \space\space\space margin(w,b) \\ subject\ to\space\space\space every\space\space\ y_n(w^Tx+b)>0 \\ \space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space margin(w,b)=\underset{n=1,\cdots,N}{min}distance(x_n,w,b) \]

這樣看起來不知道怎么求，那我們現在再進一步地簡化這個模型，先簡化樣本點到超平面的距離\(distance(x_n,w,b)\)，先看下面這一張圖

可以取超平面任意點\(x'\)和\(x''\)，它們滿足\(w^Tx'+b=0\)和\(w^Tx''+b=0\),將兩式相減，得到\(w^T(x'-x'')=0\)，這個說明什么呢？由於\(x'-x''\)表示平面上的任意向量，而\(w^T\)和它們點乘為0，說明\(w^T\)是這個平面的法向量,那距離就好辦了，樣本點\(x\)與超平面任一點\(x'\)組成向量，根據我們學過的數學知識知道，平面外一點到平面的距離等於改點與平面內任一點組成的向量點乘上法向量，然后除以法向量的模，即

\[distance(x_n,w,b)= \left|\frac{w^T(x_n-x')}{\left ||w\right||} \right|=\left | \frac{w^Tx_n+b}{\left ||w\right||} \right |= \frac { y_n(w^Tx_n+b)}{\left ||w\right||} \]

現在，我們把得到的結果帶到原來的模型當中去,得到一個簡化模型

\[\underset{w,b}{max} \space\space\space margin(w,b) \\ subject\ to\space\space\space every\space\space\ y_n(w^Tx+b)>0 \\ \space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space margin(w,b)=\underset{n=1,\cdots,N}{min}\frac { y_n(w^Tx_n+b)}{\left ||w\right||} \]

現在看看，好像這個模型還是挺復雜的，不知道怎么解，我們在來簡化簡化，\(\frac { y_n(w^Tx_n+b)}{\left ||w\right||}\)這一項是有些特別的，同倍放大\(w\)和\(b\),這項距離（書上一般寫的是幾何間隔）是不會變的，我們下面做出了一個強硬的化簡：

\[\underset{n=1,\cdots,N}{min}\space\space\space y_n(w^Tx_n+b) = 1 \]

上面的表達式\(y_n(w^Tx_n+b)\)書上叫做函數間隔，為什么可以這樣做呢？實際上我們求距離是用幾何間隔，縮放函數間隔不影響幾何間隔，舉個例子吧，樣本點\(x_1\)是離超平面最近的點，但是\(y_n(w^Tx_1+b)=3.24\),我們把\(w\)和\(b\)同時放縮3.24倍，\(y_n(w^Tx_1+b)=1\),幾何間隔和以前沒變的一樣，因為上下是同時縮放的，我們再來看看超平面\(w^Tx_n+b=0\)，也沒什么變化，\(2x+2=0\)和\(x+1=0\)還是代表同一個超平面.這樣縮放后，\(margin(w,b)=\frac {1}{\left ||w\right||}\),而且參數\(b\)不見了，可以去掉b好了，最小函數間隔為1，所以約束條件也發生了變化。這樣簡化后，我們再來看看模型變成啥樣了

\[\underset{b,w}{max} \space\space\space \frac {1}{\left ||w\right||} \\ subject\ to\space\space\space every\space\space\ y_n(w^Tx_n+b)\geq1 \]

比之前看着簡單多了吧，但是這樣變換有沒有問題呢？其實有一個充要轉必要問題，不夠直觀，原始條件（充要問題）是\(y_n(w^Tx_n+b)>0\)和\(\underset{n=1,\cdots,N}{min} \space y_n(w^Tx_n+b)=1\)，我們把它合起來了，變成\(y_n(w^Tx_n+b)\geq1\)(必要問題),后面這個條件其實沒有前面兩個條件嚴謹，后面條件不能保證最小值就等於\(1\),大於1也同樣滿足，但其實這里必要問題是等價於充要問題的，為什么呢？我們假設最優\((w_1,b_1)\)在\(y_n(w_1^Tx_n+b_1)=1.127\)上，滿足\(y_n(w^Tx_n+b)\geq1\)，可是我們可以找到更優的\((w_2,b_2)=(\frac{w_1}{1.127},\frac{b_1}{1.127})\)，同樣滿足\(y_n(w^Tx_n+b)\geq1\)，所以不管怎樣，最優\((w,b)\)一定會在\(y_n(w^Tx_n+b)=1\)上，所以原問題是等價於簡化后的模型的。
為了習慣表示，我們再把最大化問題轉化成最小化問題（例如我們習慣用梯度下降法最小化代價函數一樣），將\(\underset{w}{max} \space \frac {1}{\left ||w\right||}\)變成\(\underset{w}{min}\frac{1}{2}w^Tw\)（\(\left ||w\right||=\sqrt{w^Tw}\),這里去掉了根號和加了一個\(\frac{1}{2}\)，加\(\frac{1}{2}\)是為了后面求導方便，不影響最優化問題）,所以原模型也可以再進一步化簡：

\[\underset{b,w}{min}\space \frac{1}{2}w^Tw \\ subject\ to\space\space\space \space\ y_n(w^Tx_n+b)\geq1\space for\space all \space n \]

好了，我們已經得到非常精簡的模型，那么我們現在該想想，要用什么方法去得到這個最優解呢？梯度下降法可不可以呢？好像不太適用，梯度下降法沒有約束條件，那怎么辦？幸運的是，這是一個二次規划(quadratic programming)問題,簡稱QP問題，QP問題需要符合兩個條件：

• (1) 一個是目標函數是關於自變量的二次函數，而且函數需要是凸（convex）函數，顯然，\(\frac{1}{2}w^Tw\)是二次函數，而且確實是凸函數。
• (2) 約束條件必須是關於自變量的線性約束條件，可以看出最后的模型約束條件確實是線性約束條件。

下面舉個小例子 $$X=\left[ \begin{matrix}0 & 0\2 & 2 \2 & 0\3 & 0
\end{matrix}\right]， y=\left[ \begin{matrix}-1 \-1 \1 \1\end{matrix}\right]$$
\(X\)是二維平面上的點，分別是\((0,0),(2,2),(2,0),(3,0)\),對應\(y\)中的\(label(-1,-1,1,1)\),現在用上面最優化的模型進行求解超平面，帶入所有數據點，可得到如下限制條件:

\[ -b\geq1 \space\space\space\space(i)\\ -2w_1-2w_2-b\geq1\space\space\space\space(ii)\\ 2w_1+b\geq1 \space\space\space\space(iii)\\ 3w_1+b\geq1\space\space\space\space(iv)\]

結合(i)~(iv)，得到\(w_1\geq1\)和\(w_2\leq-1\),從而得到\(\frac{1}{2}w^Tw=\frac{1}{2}(w_1^2+w_2^2)\geq1\)，所以最優解\((w_1,w_2,b)=(1,-1,-1)\)，最后可畫出如下圖像

四個點我們算起來可能很容易，但40000個點呢，手算就不切實際了，自己編工作量太大，而且規范性差，還好，前輩們建立很多解優化問題的庫，例如libsvm，cvxpy等等。一般求解二次規划函數都可以用以下格式表示：

\[optimal\space u\leftarrow QP(q,p,A,c) \\ \underset{u}{min} \space\space\space\space \frac{1}{2}u^TQu+p^Tu \\ subject\space to\space\space A^Tu\geq c_n \space\space for\space all\space samples \]

將上面的模型帶入這個優化函數，參數分別如下：

\[objective \space function: u = \left [ \begin{matrix}b\\w \end{matrix}\right ];Q=\left[\begin{matrix}0&0_d^T\\0_d&I_d\end{matrix}\right];p=\left[\begin{matrix}0_{d+1}\end{matrix}\right]\\constraints:A=y_n\left[\begin{matrix}1&x_n\end{matrix}\right];c_n=1 \]

下面用matlab(使用quadprog函數)和python(基於cvxpy)分別生成一個線性可分數據集的分類器
MATLAB版

% 實現線性可分SVM
clear all ; close all;
file = importdata('testSet.txt');
[m,n] =size(file);
x_1 = file(:,1);
x_2 = file(:,2);
y_n = file(:,3);
file_1 = sortrows(file,3);%根據第三列進行排序
CountNegative1 = sum(y_n==-1);
x_3 = file_1(1:CountNegative1,1); %當y_n 是-1的時候，存對應的x_1
x_4 = file_1(1:CountNegative1,2); %當y_n 是-1的時候，存對應的x_2
x_5 = file_1(CountNegative1+1:m,1); %當y_n 是1的時候，存對應的x_1
x_6 = file_1(CountNegative1+1:m,2); %當y_n 是1的時候，存對應的x_2
figure;
scatter(x_3,x_4,'red','fill');
hold on;
scatter(x_5,x_6,'blue','fill');
H = [0 0 0;0 1 0;0 0 1];
p = [0;0;0];
A = [ones(CountNegative1,1),file_1(1:CountNegative1,1:2);
    -ones(m-CountNegative1,1),- file_1(CountNegative1+1:m,1:2)];
c = -ones(m,1);
u = quadprog(H,p,A,c); %求出參數[b,w]
x = linspace(2,8,400);
y = -(u(1)/u(3)+u(2)/u(3)*x);
plot(x,y,'green');
sv_index = find(abs(u'*[ones(m,1),file_1(:,1:2)]')-1<1e-5);
[m1,n1]=size(sv_index);
for i = 1:n1
    x_sv =file_1(sv_index(i),1);
    y_sv =file_1(sv_index(i),2);
    plot(x_sv,y_sv,'o','Markersize',20);
end

這里程序不多講，懂上面原理，理解程序還是挺簡單的，函數不懂的請看MATLAB的幫助文檔，這里用的testSet.txt其實是機器學習實戰SVM一節的數據集(下載鏈接),最后運行的結果如下圖所示(那些用圈圈圈起來的點叫支持向量Support vecter,是離超平面最近的點)

python版

import cvxpy as cvx
import numpy as np

f = open('testSet.txt')

X1_Neg1 = []
X1_Pos1 = []
X2_Neg1 = []
X2_Pos1 = []

constraints = []

for line in f.readlines():
    lineSplit = line.strip().split('\t')
    xAxis = float(lineSplit[0])
    yAxis = float(lineSplit[1])
    label = int(lineSplit[2])
    if label == 1:
        X1_Pos1.append(xAxis)
        X2_Pos1.append(yAxis)
    else:
        X1_Neg1.append(xAxis)
        X2_Neg1.append(yAxis)

u = cvx.Variable(3)
H = np.array([[0,0,0],[0,1,0],[0,0,1]]) 

objective = cvx.Minimize(1.0/2*cvx.quad_form(u,H))

for X_Pos1 in zip(X1_Pos1,X2_Pos1):
    temp1 = np.zeros([3,1])
    temp1[0] = 1
    temp1[1] = X_Pos1[0]
    temp1[2] = X_Pos1[1]
    constraints.append(temp1.T*u >= 1)
for X_Neg1 in zip(X1_Neg1,X2_Neg1):
    temp2 = np.zeros([3,1])
    temp2[0] = 1
    temp2[1] = X_Neg1[0]
    temp2[2] = X_Neg1[1]
    constraints.append(-(temp2.T*u) >= 1)
    
prob = cvx.Problem(objective,constraints)
prob.solve() 
print('optimal var:\n\r',u.value)

這里圖就不畫了，和上面一樣。

4.非線性SVM

線性SVM是有很大的局限性的，並不是所有數據都能用線性SVM分類，以下圖為例，你很難用線性SVM進行分類達到高准確，線性SVM最高只能達到75%准確度。

這時候就需要用到一個知識————特征轉換(feature transform)
我們將二維數據映射到三維上，轉換關系如下

\[X_1 = x_1^2 \\ X_2 = x_2^2 \\ X_3 = \sqrt{2}x_1x_2 \]

原始數據集分布就變成下圖分布

我們就可以在當前三維分布利用線性SVM進行分類，如圖所示

看起來分離得很好，我們再將平面映射回初始的二維空間中，看看決策邊界時什么樣子

可以看到決策邊界不是線性的了，還是兩條歪曲的曲線構成的，這就是非線性SVM的一個典型例子，它將數據映射到更高維空間，然后在高維進行線性SVM，最后映射回到原來的狀態，得到復雜的決策邊界————曲面形(原因在於特征轉換是非線性的，所以映射回來，高維線性變成低維非線性)。
我們再來看看非線性SVM模型

\[\underset{b,w}{min}\space \frac{1}{2}w^Tw \\ subject\ to\space\space\space \space\ y_n(w^T z_n+b)\geq1\space for\space all \space n \]

這樣看起來和線性SVM的模型一樣，但其實有區別，約束條件中的\(z_n\)是\(x_n\)經過特征轉換得到的，由於\(z_n\)是高維的，\(w\)也會發生相應的變化，具體要看\(z_n\)的維度(和\(z_n\)維度相同)，我們利用優化庫去解這個二次規划問題的方法和之前提到的線性SVM方法幾乎一樣，假設\(z_n\)的維度是\(\overset{\sim}{d}\)，我們可以列出下面的表達式

• \(Q=\left[\begin{matrix}0&0_\overset{\sim}{d}^T\\0_\overset{\sim}{d}&I_\overset{\sim}{d}\end{matrix}\right];p=\left[\begin{matrix}0_{\overset{\sim}{d}+1}\end{matrix}\right]\\A=y_n\left[\begin{matrix}1&z_n^T\end{matrix}\right];c_n=1\)
• $optimal\space \left [ \begin{matrix}b\w \end{matrix}\right ]\leftarrow QP(Q,p,A,c) $
• \(g_{svm}(x)=sign(w^Tz_n+b)\)

看起來非線性SVM已經介紹完了，其實這里遺留了一個問題，我們現在的變量維度是\(\overset{\sim}{d}+1\)，限制條件的個數為N(N為樣本數)，我們要特別注意這個\(\overset{\sim}{d}\)，它的維度會非常高，以64x64x3的圖片為例，它映射后的\(z_n\)是百萬級的，這樣\(w\)也是百萬級的，這無疑是巨大的計算量，前輩們對這一問題，提出了一個非常好的解決方案————引入拉格朗日對偶，將原始問題轉換對偶問題，從而極大簡化了計算量

Original SVM \(\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\Longrightarrow\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\) Equivalent SVM
\(\overset{\sim}{d}+1\space Variables\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space N+1\space Variables\)
\(N\space\space\space\space constraints\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space N\space\space\space\space constraints\)

我們接下來來看看這個轉換是如何做到的，首先這里要引入一個拉格朗日乘子法，將原始帶有限制條件問題，轉換成沒有限制條件的問題，原始問題經過拉格朗日乘子法后，得到如下表達式(其中\(\alpha_n\)是拉格朗日乘子)：

\[L(w,b,\alpha)=\underbrace{\frac{1}{2}w^Tw}_{objective}+\sum_{n=0}^N \alpha_n\underbrace{(1-y_n(w^Tz_n+b))}_{constraint} \]

這樣，原始問題就可以轉化成求解\(\underset{b,w}{min}(\underset{all \space \alpha_n\geq0}{max}L(w,b,\alpha))\)，那么為什么可以這樣轉換呢？我們可以先分析里面\(\underset{all \space \alpha_n\geq0}{max}L(w,b,\alpha)\)這一項，里面最大化的自變量只有\(\alpha_n\)，我們現在可以分情況討論了

• 如果\((w,b)\)使得\((1-y_n(w^Tz_n+b)>0\),由於是取最大值且\(\alpha_n\geq0\)，那么\(\alpha_n(1-y_n(w^Tz_n+b)\)會趨近於無窮大，這樣我們在最外層的最小化會將它淘汰掉(最小化無窮大沒有意義)
• 如果\((w,b)\)使得\((1-y_n(w^Tz_n+b)<0\)，那么\(\alpha_n(1-y_n(w^Tz_n+b)\)會趨近於0，我們的最大值就是\(\frac{1}{2}w^Tw\)，這樣外層最小化才可以進行.
• 所以這里\((1-y_n(w^Tz_n+b)<0\)是隱藏在最大化中，也就是原始約束條件並沒有消失，只是換了一種方式。

好了，現在我們已經將原始問題轉化成求解\(\underset{b,w}{min}(\underset{all \space \alpha_n\geq0}{max}L(w,b,\alpha))\)，那我們在進行下一步地操作吧，我們知道對於任意\(\alpha'\),有$$\underset{b,w}{min}(\underset{all \space \alpha_n\geq0}{max}L(w,b,\alpha))\geq\underset{b,w}{min}L(w,b,\alpha')$$
這是因為最大值大於等於任意取值，然后我們再往右邊加一個最大化，得

\[\underset{b,w}{min}(\underset{all\space\alpha_n\geq0}{max}L(w,b,\alpha))\geq\underset{all\space\alpha_n'\geq0}{max}(\underset{b,w}{min}L(w,b,\alpha')) \]

加上后依然成立，這是因為最好的\(\alpha'\)也是任意中的一個。
然后，我們是可以把\(\geq\)換成\(=\)，這是由於原問題滿足以下三個條件

• 原問題是一個凸的二次規划問題
• 原問題是有解的
• 限制條件是線性約束條件

滿足上面三個條件，就可以將\(\geq\)的弱對偶轉換成\(=\)的強對偶，具體證明我也不清楚，這是前輩們推導出來了，知道就行了，好了，這就完成了拉格朗日對偶，如下

\[\underbrace{\underset{b,w}{min}(\underset{all\space\alpha_n\geq0}{max}L(w,b,\alpha))}_{original\space SVM}=\underbrace{\underset{all\space\alpha_n\geq0}{max}(\underset{b,w}{min}L(w,b,\alpha))}_{Lagrange\space dual\space SVM} \]

我們接下來用對偶SVM來求解原問題，對偶SVM如下

\[\underset{all\space\alpha_n\geq0}{max}(\underset{b,w}{min}\underbrace{\frac{1}{2}w^Tw+\sum_{n=1}^N \alpha_n(1-y_n(w^Tz_n+b))}_{L(w,b,\alpha)}) \]

里面的最小化是沒有約束條件的，約束條件在最大化上，於是我們可以取梯度為0來求極值,由\(\frac{\partial L(w,b,\alpha)}{\partial b}\)=0，得\(\sum_{n=1}^N\alpha_ny_n=0\)，代入對偶SVM，化簡得

\[\underset{all\space\alpha_n\geq0,\sum_{n=1}^N\alpha_ny_n=0}{max}(\underset{b,w}{min}\frac{1}{2}w^Tw+\sum_{n=1}^N \alpha_n(1-y_n(w^Tz_n))) \]

再由\(\frac{\partial L(w,b,\alpha)}{\partial w}\)=0，得\(w=\sum_{n=1}^N\alpha_ny_nz_n\)，代入對偶SVM，化簡得$$\underset{all\space\alpha_n\geq0,\sum_{n=1}^N\alpha_ny_n=0,w=\sum_{n=0}^N\alpha_ny_nz_n}{max}(-\frac{1}{2}{\left|\sum_{n=1}^N\alpha_ny_nz_n\right|}^2+\sum_{n=1}^N \alpha_n)$$
將上面的表達式再變換一下，得

\[\underset{\alpha}{min}\space\space \frac{1}{2}\sum_{n=1}^N\sum_{m=1}^N\alpha_n\alpha_my_ny_mz_n^Tz_m-\sum_{n=1}^N\alpha_n\\subject\space to \space\space \sum_{n=1}^Ny_n\alpha_n=0;\alpha_n\geq0(for\space n=1,2,\cdots,N) \]

好了，到現在，我們終於把原始的\(\overset{\sim}{d}+1\space Variables\)轉換成\(N+1\space Variables\)。

我們根據上面的推導知道，原始和對偶\(w,b,\alpha\)之間是滿足一些條件的，Karush，Kuhn，Tucker三位科學家總結了它們之間的關系，包含四個條件：

primal-dual optimal(w,b,\(\alpha\))
(1) 原問題需滿足\(y_n(w^Tx_n+b)\geq1\)
(2) 對偶問題\(\alpha_n\geq0\)
(3) 對偶問題最優解\(\sum_{n=1}^N\alpha_ny_n=0,w=\sum_{n=0}^N\alpha_ny_nz_n\)
(4) \(\alpha_n(1-y_n(w^Tx_n+b))\)=0

后人命名為KKT條件(Karush-Kuhn-Tucker，以三位科學家命名)，根據KKT條件，我們就可以通過\(\alpha\)求解\((w,b)\)，求\(w\)顯而易見，這里說一下\(b\)，根據KKT第四個條件，當\(\alpha\)的一個分量\(\alpha_j\)大於0時，對應的點必須滿足\(1-y_n(w^Tx_n+b)=0\)，從而可以求出\(b\)
下面進行代碼實踐
matlab版

%非線性SVM
clear all ; close all;
data = csvread('test_nonlinear_SVM.csv');
[m,n] = size(data);
figure;
scatter(data(1:300,1),data(1:300,2),'red','fill');
hold on;
scatter(data(301:600,1),data(301:600,2),'blue','fill');
%feature transfrom x3 = x1^2+x2^2(x1,x2不變)
data_z = [data(:,1),data(:,2),data(:,1).^2+data(:,2).^2,data(:,3)];  
H = (kron(data_z(:,4),ones(1,3)).*data_z(:,1:3))*(data_z(:,1:3)'.*kron(data_z(:,4),ones(1,3))');
f = -ones(m,1);
A = -eye(m);
b = zeros(m,1);
Aeq = data_z(:,4)';
beq = 0;
alpha = quadprog(H,f,A,b,Aeq,beq);%每個參數的具體意義請參考quadprog函數
w = (alpha.*data_z(:,4))'*data_z(:,1:3);
b = -1-w*data_z(377,1:3)'; %第377個點的alpha大於0
x = linspace(-7.6,7.6,400); 
y1 = sqrt((-(b+w(1)*x+w(3)*x.^2)/w(3))+((w(2)^2)/(4*w(3))))-w(2)/(2*w(3)); 
y2 = -sqrt((-(b+w(1)*x+w(3)*x.^2)/w(3))+((w(2)^2)/(4*w(3))))-w(2)/(2*w(3));
plot(x,y1);
plot(x,y2);

這里cmd markdown有一個高亮bug,代碼單引號與markdown的單引號沖突，不用管，直接復制用就可以了，test_nonlinear_SVM.csv和生成它的matlab代碼請在這個百度雲鏈接下載，這里說一下y1和y2，我們知道決策邊界是\(w^Tz+b=0\)，而\(z=(x1,x2,x3)\)，x1已經用linspace生成，畫圖，我們需要求x2，根據方程\(w^T(x1,x2,x1^2+x2^2)+b=0\)，我們可以求出x2的兩個解(即y1和y2)，最后運行結果如下

import numpy as np
import cvxpy as cvx

tmp = np.loadtxt(open("test_nonlinear_SVM.csv","rb"),delimiter=',',skiprows=0)
X,Y = tmp[:,0:2],tmp[:,2]

X_z = np.array([X[:,0],X[:,1],X[:,0]**2+X[:,1]**2])

alpha = cvx.Variable(X.shape[0],1)
H = np.dot((X_z*np.kron(Y,np.ones([3,1]))).T,(X_z*np.kron(Y,np.ones([3,1]))))
f = -np.ones([X.shape[0],1])

objective = cvx.Minimize(1.0/2*cvx.quad_form(alpha,H)+cvx.sum_entries(cvx.mul_elemwise(f,alpha)))
constraints = [alpha>=0,cvx.sum_entries(cvx.mul_elemwise(Y,alpha))==0]
prob = cvx.Problem(objective,constraints)

prob.solve()
print('alpha_value:',alpha.value)

w = np.dot((np.array(alpha.value)*Y.reshape((X.shape[0],1))).T,X_z.T)
b = Y[np.argmax(alpha.value)]-np.dot(w,X_z[0:3,np.argmax(alpha.value)].reshape(3,1))
print('w:',w,'\r\n','b:',b)

圖就不畫了，和matlab的差不多，大家可以將結果與matlab的最優w和最優b進行比較

5.kernel SVM

核SVM的引入，將大大減少上面非線性SVM的計算量，我們知道，根據上面非線性SVM問題，需要計算\(Q=y_ny_mz_n^Tz_m\)，而我們要注意一下，\(z_n\)的維度通常情況是遠遠大於\(x_n\)的維度的，甚至可以是無限大(上面實例為了簡單，畫圖方便，差的維度看不太出來)，當轉換到非常高的維度時，\(z_nz_m\)的內積計算量會是非常大的，這時發現了一種新的方法，使用它，大大減少了計算量，這個方法叫核方法。

以二次多項式變換為例，原始\(x_n\)的維度是\(d\),變換后\(z_n\)的維度是\(d^2\)

\[z_n=\phi(x_n)=(1,x_1,x_2,\cdots,x_n,x_1^2,x_1x_2,\cdots,x_1x_d,x_2^2,\cdots,x_2x_d,\cdots,x_d^2)$$有什么辦法去簡化$z_nz_m$的內積呢？我們先把二次多項式的內積表達式寫出來,看能不能進行簡化 $$zz^T=\phi(x)\phi(x')=1+\sum_{i=1}^dx_ix_i'+\sum_{i=1}^d\sum_{j=1}^dx_ix_jx_i'x_j'=1+\sum_{i=1}^dx_ix_i'+\sum_{i=1}^dx_ix_i'\sum_{j=1}^dx_jx_j'\\\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space\space=1+x^Tx'+(x^Tx')(x^Tx')$$我們發現這個表達式可以用原始$x_n$的內積表示，我們將這個問題簡化成求原始$x_n$的內積，然后套上一個核函數，就得到轉換后的$z_n$的內積(*計算量從原來的$\omicron(d^2)$轉化成了$\omicron(d)$*)，上面過程可以用下面一個表達式表示 $$zz^T=K_{\phi}(x,x')=1+x^Tx'+(x^Tx')(x^Tx')\]

所以
\(Q=y_ny_mK_{\phi}(x_n,x_m)\)
\(optimal\space b=y_s-w^Tz_s=y_s-(\sum_{n=1}^Na_ny_nz_n)z_s=y_s-\sum_{n=1}^Na_ny_n(K_{\phi}(z_n,z_s))\)
\(g_{svm}(x)=sign(\sum_{n=1}^Na_ny_n(K_{\phi}(x_n,x))+b)\)
一般的二次項核(更常用)可以表達為如下形式

\[K_2(x,x')=(1+\gamma x^Tx')^2\space with\space \gamma>0$$上面例子的$\gamma$等於0.5，不同的$\gamma$會導致不同的決策邊界，也會出現不同的支持向量，下圖表示不同的$\gamma$導致的不同結果 ![1.jpg-621.1kB][11] 通用**多項式核**可以如下表示 $$K_Q(x,x')=(\zeta+\gamma x^Tx')^Q\space with\space\zeta\geq0\space \gamma>0$$多項式核的優點 > (1)比線性核(xx')強大，限制更小 > (2)Q值具有強控制性 多項式核的缺陷 > (1)數值計算困難，$(\zeta+\gamma x^Tx')>1$時，$(\zeta+\gamma x^Tx')^Q$會非常大，當$(\zeta+\gamma x^Tx')<1$時,$(\zeta+\gamma x^Tx')^Q$會非常小，趨近於0 > (2)三個參數($\zeta,\gamma,Q$)很不好選 下面介紹**高斯核**，也叫RBF kernel，高斯核是可以表示無限維度的，下圖是原理推導 ![image_1c12bvfgn1rbj1svg1cjiok6alb22.png-74kB][12] 一般的高斯核可以表示為如下表達式$$K(x,x')=exp(-\gamma\|x-x'\|^2)$$高斯核的**優點** > (1)比多項式核和線性核更強大 > (2)數值計算困難程度比多項式核小 > (3)只有一個參數，比較好選擇 高斯核的**缺陷** > (1)沒有w(無限維度的原因) > (2)運算速度比線性慢 > (3)太強大，容易過擬合，如下圖 ![image_1c12d9lfr1t64fsn1o1o1vvr8vn2f.png-57.4kB][13] 給定一個函數K，我們如何怎樣判斷它是否是一個有效的核函數呢？ 給定m個訓練樣本$\{x^{(1)},x^{(2)},\cdots,x_{(m)}\}$，每一個$x^{(i)}$對應一個特征向量。那么，我們可以將任意兩個$x^{(i)}$和$x^{(j)}$帶入K中，計算得到。i可以從1到m，j可以從1到m，這樣可以計算出m*m的核函數矩陣。如果假設K是有效地核函數，則有$$K_{ij}=K(x^{(i)},x^{(j)})=\phi(x^{(i)})^T\phi(x^{(j)})=\phi(x^{(j)})^T\phi(x^{(i)})=K(x^{(j)},x^{(i)})=K_{ji}$$可見，矩陣K一定是**對稱**的矩陣。 下面來看看另一個條件，對於任意向量z，有 ![image_1c12klkbqc0p1le532f7m1hlv2s.png-26.2kB][14] 由此知道矩陣K需要是**半正定**的，綜合起來**K是有效的核函數 <==> 核函數矩陣K是對稱半正定的。** **python實踐**(僅供參考，cvx.quad_form(alpha,Q)不遵守DCP，我用非線性SVM算出的Q沒問題，這個就有問題，具體原因不清楚) ```python import numpy as np import h5py import cvxpy as cvx def load_dataset(): train_dataset = h5py.File('datasets/train_catvnoncat.h5', "r") train_set_x_orig = np.array(train_dataset["train_set_x"][:]) # your train set features train_set_y_orig = np.array(train_dataset["train_set_y"][:]) # your train set labels test_dataset = h5py.File('datasets/test_catvnoncat.h5', "r") test_set_x_orig = np.array(test_dataset["test_set_x"][:]) # your test set features test_set_y_orig = np.array(test_dataset["test_set_y"][:]) # your test set labels classes = np.array(test_dataset["list_classes"][:]) # the list of classes train_set_y_orig = train_set_y_orig.reshape((1, train_set_y_orig.shape[0])) test_set_y_orig = test_set_y_orig.reshape((1, test_set_y_orig.shape[0])) for i in range(train_set_y_orig.shape[1]): if(train_set_y_orig[0,i]==0): train_set_y_orig[0,i]=-1 for i in range(test_set_y_orig.shape[1]): if(test_set_y_orig[0,i]==0): test_set_y_orig[0,i]=-1 return train_set_x_orig, train_set_y_orig, test_set_x_orig, test_set_y_orig, classes def main(): train_set_x_orig, train_set_y, test_set_x_orig, test_set_y, classes = load_dataset() train_set_x_orig=train_set_x_orig.reshape((train_set_x_orig.shape[0],-1)) alpha = cvx.Variable(train_set_x_orig.shape[0],1) f = -np.ones([train_set_x_orig.shape[0],1]) K2 = (np.ones([train_set_x_orig.shape[0],train_set_x_orig.shape[0]])+1.0/2*np.dot(train_set_x_orig,train_set_x_orig.T))**2 y_m=np.kron(train_set_y,np.ones([train_set_x_orig.shape[0],1])) Q=np.array(y_m.T*K2*y_m) objective = cvx.Minimize(1.0/2*cvx.quad_form(alpha,Q)+cvx.sum_entries(cvx.mul_elemwise(f,alpha))) constraints = [alpha>=0,cvx.sum_entries(cvx.mul_elemwise(train_set_y.T,alpha))==0] prob = cvx.Problem(objective,constraints) prob.solve() if __name__ == '__main__': main() ``` ## 6.軟間隔SVM 有的時候，我們強行要把所有的樣本點都分隔正確，會產生過擬合。如下圖所示 ![image_1c12lu9hdsph70edlgbu164s39.png-12kB][15] 這個時候，我們可以放低標准，容許分類器有一定的錯誤，像下圖一樣 ![image_1c12m1iejsek1s0d6ak7131c0243.png-8.6kB][16] 這樣做的好處是使分類器准確地泛化了樣本的分布，這樣的SVM分類器叫軟間隔SVM，那么怎樣實現它呢？ 首先，我們需要對原始問題進行修改 ![image_1c12o6sqp13gk1ahh1t4hcceoic9.png-34.3kB][17] 這里引入了變量C，它是用來平衡大間隔和容忍噪聲的，C越大，對噪聲容忍度越小，希望分類錯誤的點越少越好，反之，C越小，則希望我們的大間隔越大越好 將上面兩個條件合起來，得到如下問題 ![2.jpg-10.3kB][18] 但這里出現兩個問題了，**一個是【$y_n\neq sign(w^Tz_n+b)$】非線性的**，這樣導致這個問題不是QP問題，對偶和核函數都用不了，第二個是這樣的問題**不能區分小錯誤和大錯誤**，比如，一個分類錯誤樣本離邊界0.001和離邊界1000，這樣分類器認為這兩種是一樣的，這顯然也不合適。 於是，我們引入$\varepsilon_n$參數，表示**錯誤的點離間隔邊界的距離**，用它替換條件中的【$y_n\neq sign(w^Tz_n+b)$】，這樣限制條件變成了線性的，同時，最小化中的錯誤個數計數也被錯誤程度取代，問題也變成了QP問題，皆大歡喜! **soft-margin SVM:** $$\underset{b,w,\varepsilon_n}{min}\space \frac{1}{2}w^Tw+C\sum_{i=1}^N\varepsilon_n \\ subject\ to\space\space\space \space\ y_n(w^Tz_n+b)\geq1-\varepsilon_n\space and\space \varepsilon_n\geq0for\space all \space n\]

現在我們來看看這個QP問題的變量和條件，變量有\(\overset{\sim}{d}+1+N\)(\(w_n\)有\(\overset{\sim}{d}\),\(b\)有1個,\(\varepsilon_n\)有\(N\)個)，限制條件有\(2N\)個
接下來就是用拉格朗日對偶來去掉\(\overset{\sim}{d}\)，過程詳細解說就不寫了，原理和上面對偶SVM類似，這里附上過程







最后得到N個變量和2N+1個限制條件的QP問題。

參考：
1 林軒田——機器學習技法課程
2 http://blog.csdn.net/myarrow/article/details/51261971