SVM大解密（附程式碼和公式）
來源 | 我i智慧（公眾號ID：AInewworld）
寫在之前
支援向量機（SVM），一個神祕而眾知的名字，在其出來就受到了莫大的追捧，號稱最優秀的分類演算法之一，以其簡單的理論構造了複雜的演算法，又以其簡單的用法實現了複雜的問題，不得不說確實完美。
本系列旨在以基礎化的過程，例項化的形式一探SVM的究竟。曾經也只用過整合化的SVM軟體包，效果確實好。因為眾人皆說原理複雜就對其原理卻沒怎麼研究，最近經過一段時間的研究感覺其原理還是可以理解，這裡希望以一個從懵懂到略微熟知的角度記錄一下學習的過程。其實網路上講SVM演算法的多不勝數，部落格中也有許多大師級博主的文章，寫的也很簡單明瞭，可是在看過之和總是感覺像差點什麼，當然對於那些基礎好的可能一看就懂了，然而對於像我們這些薄基礎的一遍下來也能馬馬虎虎懂，過一兩天後又忘了公式怎麼來的了。比如說在研究SVM之前，你是否聽說過拉格朗日乘子法？你是否知道什麼是對偶問題？你是否瞭解它們是怎麼解決問題的？Ok這些不知道的話，更別說什麼是KKT條件了，哈哈，有沒有說到你的心聲，不用怕，學學就會了。話說像拉格朗日乘子法，在大學裡面學數學的話，不應該沒學過，然你學會了嗎？你知道是幹什麼的嗎？如果那個時候就會了，那你潛質相當高了。作為一個過來的人，將以簡單例項化形式記錄自己的學習過程，力圖幫助新手級學習者少走彎路。
一、關於拉格朗日乘子法和KKT條件
1）關於拉格朗日乘子法
首先來了解拉格朗日乘子法，那麼為什麼需要拉格朗日乘子法？記住，有拉格朗日乘子法的地方，必然是一個組合優化問題。那麼帶約束的優化問題很好說，就比如說下面這個：
這是一個帶等式約束的優化問題，有目標值，有約束條件。那麼想想假設沒有約束條件這個問題是怎麼求解的呢？是不是直接f對各個x求導等於0,，解x就可以了，可以看到沒有約束的話，求導為0，那麼各個x均為0吧，這樣f=0了，最小。但是x都為0不滿足約束條件呀，那麼問題就來了。這裡在說一點的是，為什麼上面說求導為0就可以呢？理論上多數問題是可以的，但是有的問題不可以。如果求導為0一定可以的話，那麼f一定是個凸優化問題，什麼是凸的呢？像下面這個左圖：
凸的就是開口朝一個方向（向上或向下）。更準確的數學關係就是：
注意的是這個條件是對函式的任意x取值。如果滿足第一個就是開口向上的凸，第二個是開口向下的凸。可以看到對於凸問題，你去求導的話，是不是隻有一個極點，那麼他就是最優點，很合理。類似的看看上圖右邊這個圖，很明顯這個條件對任意的x取值不滿足，有時滿足第一個關係，有時滿足第二個關係，對應上面的兩處取法就是，所以這種問題就不行，再看看你去對它求導，會得到好幾個極點。然而從圖上可以看到，只有其中一個極點是最優解，其他的是區域性最優解，那麼當真實問題的時候你選擇那個？說了半天要說啥呢，就是拉格朗日法是一定適合於凸問題的，不一定適合於其他問題，還好我們最終的問題是凸問題。
回頭再來看看有約束的問題，既然有了約束不能直接求導，那麼如果把約束去掉不就可以了嗎？怎麼去掉呢？這才需要拉格朗日方法。既然是等式約束，那麼我們把這個約束乘一個係數加到目標函式中去，這樣就相當於既考慮了原目標函式，也考慮了約束條件，比如上面那個函式，加進去就變為：
這裡可以看到與相乘的部分都為0，所以的取值為全體實數。現在這個優化目標函式就沒有約束條件了吧，既然如此，求法就簡單了，分別對x求導等於0，如下：
把它在帶到約束條件中去，可以看到，2個變數兩個等式，可以求解，最終可以得到,這樣再帶回去求x就可以了。那麼一個帶等式約束的優化問題就通過拉格朗日乘子法完美的解決了。那麼更高一層的，帶有不等式的約束問題怎麼辦？那麼就需要用更一般化的拉格朗日乘子法即KKT條件來解決這種問題了。
2）關於KKT條件
繼續討論關於帶等式以及不等式的約束條件的凸函式優化。任何原始問題約束條件無非最多3種，等式約束，大於號約束，小於號約束，而這三種最終通過將約束方程化簡化為兩類：約束方程等於0和約束方程小於0。再舉個簡單的方程為例，假設原始約束條件為下列所示：
那麼把約束條件變個樣子：
為什麼都變成等號與小於號，方便後面的，反正式子的關係沒有發生任何變化就行了。
現在將約束拿到目標函式中去就變成：
那麼KKT條件的定理是什麼呢？就是如果一個優化問題在轉變完後變成
其中g是不等式約束，h是等式約束（像上面那個只有不等式約束，也可能有等式約束）。那麼KKT條件就是函式的最優值必定滿足下面條件：
這三個式子前兩個好理解，重點是第三個式子不好理解，因為我們知道在約束條件變完後，所有的g(x)<=0，且，然後求和還要為0，無非就是告訴你，要麼某個不等式,要麼其對應的。那麼為什麼KKT的條件是這樣的呢？
假設有一個目標函式，以及它的約束條件，形象的畫出來就如下：
假設就這麼幾個吧，最終約束是把自變數約束在一定範圍，而函式是在這個範圍內尋找最優解。函式開始也不知道該取哪一個值是吧，那就隨便取一個，假設某一次取得自變數集合為x1，發現一看，不滿足約束，然後再換呀換，換到了x2,發現可以了，但是這個時候函式值不是最優的，並且x2使得g1(x)與g2(x)等於0了，而g3(x)還是小於0。
這個時候，我們發現在x2的基礎上再尋找一組更優解要靠誰呢？當然是要靠約束條件g1(x)與g2(x)，因為他們等於0了，很極限呀，一不小心，走錯了就不滿足它們兩了，這個時候我們會選擇g1(x)與g2(x)的梯度方向往下走，這樣才能最大程度的拜託g1(x)與g2(x)=0的命運，使得他們滿足小於0的約束條件對不對。至於這個時候需不需要管g3(x)呢？正常來說管不管都可以，如果管了，也取g3在x2處的梯度的話，因為g3已經滿足了小於0的條件，這個時候在取在x2處的梯度，你能保證它是往好的變了還是往差的變了？答案是都有可能。運氣好，往好的變了，可以更快得到結果，運氣不好，往差的變了，反而適得其反。
那麼如果不管呢？因為g1(x)與g2(x)已經在邊緣了，所以取它的梯度是一定會讓目標函式變好的。綜合來看，這個時候我們就不選g3。那麼再往下走，假設到了自變數優化到了x3，這個時候發現g2(x)與g3(x)等於0，也就是走到邊了，而g1(x)小於0，可變化的空間綽綽有餘，那麼這個時候舉要取g2(x)與g3(x)的梯度方向作為變化的方向，而不用管g1(x)。那麼一直這樣走呀走，最終找到最優解。可以看到的是，上述如果g1(x)、g2(x)=0的話，我們是需要優化它的，又因為他們本身的條件是小於0的，所以最終的公式推導上表明，是要乘以一個正係數作為他們梯度增長的倍數，而那些不需要管的g(x)為了統一表示，這個時候可以將這個係數設定為0，那麼這一項在這一次的優化中就沒有了。那麼把這兩種綜合起來就可以表示為
也即是某次的g(x)在為最優解起作用，那麼它的係數值(可以)不為0。如果某次g(x)沒有為下一次的最優解x的獲得起到作用，那麼它的係數就必須為0，這就是這個公式的含義。
比如上面例子的目標值與約束：
將約束提到函式中有：
此時分別對x1、x2求導數：
而我們還有一個條件就是,那麼也就是：
這樣我們就去討論下，要麼g=0，要麼，這裡兩個g兩個，這樣我們就需要討論四種情況，可能你會說，這是約束條件少的情況，那麼如果有10個約束條件，這樣就有10個g和10個，你去給我討論？多少種組合，不知道，但是換個思路，我們非得去10個一起去討論？機智的學者想到一種方法，考慮到這個條件，那麼我兩個兩個討論不就可以了，比如現在我就討論7，8，讓其他的不變，為什麼選或者至少選兩個討論呢，因為這個式子求和為0，改變一個顯然是不行的，那就改變兩個，你增我就減，這樣和可以為0。再問為什麼不討論3個呢？也可以，這不是麻煩嘛，一個俗語怎麼說來著，三個和尚沒水喝，假設你改變了一個，另外兩個你說誰去減或者加使得和為0，還是兩個都變化一點呢？不好說吧，自然界都是成雙成對的才和諧，沒有成三成四的（有的話也少）。
這裡順便提一下後面會介紹到的內容，就是實現SVM演算法的SMO方法，在哪裡，會有很多，那麼人們怎麼解決的呢，就是隨便選擇兩個去變化，看看結果好的話，就接受，不好的話就捨棄在選擇兩個，如此反覆，後面介紹。
可以看到像這種簡單的討論完以後就可以得到解了。
x1=110/101=1.08;x2=90/101=0.89,那麼它得到結果對不對呢？這裡因為函式簡單，可以在matlab下畫出來，同時約束條件也可以畫出來，那麼原問題以及它的約束面畫出來就如下所示：
這是擷取下來的符合約束要求的目標面
可以看到最優解確實就是上面我們求的那個解。既然簡單的問題可以這樣解，那麼複雜一點的只需要簡單化，照樣可以解，至此KKT條件解這類約束性問題就是這樣，它對後續的SVM求解最優解至關重要。
二、SVM的理論基礎
上節我們探討了關於拉格朗日乘子和KKT條件，這為後面SVM求解奠定基礎，本節希望通俗的細說一下原理部分。
一個簡單的二分類問題如下圖：
我們希望找到一個決策面使得兩類分開，這個決策面一般表示就是W'X+b=0,現在的問題是找到對應的W和b使得分割最好，知道logistic分類 機器學習之logistic迴歸與分類的可能知道，這裡的問題和那裡的一樣，也是找權值。在那裡，我們是根據每一個樣本的輸出值與目標值得誤差不斷的調整權值W和b來求得最終的解的。當然這種求解最優的方式只是其中的一種方式。那麼SVM的求優方式是怎樣的呢？
這裡我們把問題反過來看，假設我們知道了結果，就是上面這樣的分類線對應的權值W和b。那麼我們會看到，在這兩個類裡面，是不是總能找到離這個線最近的點，向下面這樣：
然後定義一下離這個線最近的點到這個分介面（線）的距離分別為d1,d2。那麼SVM找最優權值的策略就是，先找到最邊上的點，再找到這兩個距離之和D，然後求解D的**最大值**，想想如果按照這個策略是不是可以實現最優分類，是的。好了，還是假設找到了這樣一個分介面W'X+b=0,那麼做離它最近的兩類點且平行於分類面，如上面的虛線所示。
好了再假設我們有這兩個虛線，那麼真實的分介面我們認為正好是這兩個分介面的中間線，這樣d1就等於d2了。因為真實的分介面為W'X+b=0，那麼就把兩個虛線分別設定為W'X+b=1和W'X+b=-1,可以看到虛線相對於真實面只是上下移動了1個單位距離，可能會說你怎麼知道正好是一個距離？確實不知道，就假設上下是k個距離吧，那麼假設上虛線現在為W'X+b=k，兩邊同時除k可以吧，這樣上虛線還是可以變成W'X+b=1,同理下虛線也可以這樣，然後他們的中線就是W1'X+b1=0吧，可以看到從k到1，權值無非從w變化到w1,b變到b1,我在讓w=w1,b=b1，不是又回到了起點嗎，也就是說，這個中間無非是一個倍數關係。
所以我們只需要先確定使得上下等於1的距離，再去找這一組權值，這一組權值會自動變化到一定倍數使得距離為1的。
好了再看看D=d1+d2怎麼求吧，假設分介面W'X+b=0，再假設X是兩維的，那麼分介面再細寫出來就是：W1'X1+W2'X2+b=0。上分界線：W1'X1+W2'X2+b=1，這是什麼，兩條一次函式（y=kx+b）的曲線是不是，那麼初中就學過兩直線的距離吧，
這裡W=(w1,w2)，是個向量，||W||為向量的距離，那麼||W||^2=W'W。下介面同理。這樣
,
要使D最大，就要使分母最小，這樣優化問題就變為,乘一個係數0.5沒影響，但是在後面卻有用。
注意的是這可不是一個約束條件，而是對所有的每個樣本xi都有一個這樣的約束條件。轉換到這種形式以後是不是很像上節說到的KKT條件下的優化問題了，就是這個。但是有一個問題，我們說上節的KKT是在凸函式下使用的，那麼這裡的目標函式是不是呢？答案是的，想想W'*W，函式乘出來應該很單一，不能有很多極點，當然也也可以數學證明是的。
好了那樣的話就可以引入拉格朗日乘子法了，優化的目標變為：
然後要求這個目標函式最優解，求導吧，
這兩個公式非常重要，簡直是核心公式。
求導得到這個應該很簡單吧，那我問你為什麼W'W 對w求導是w呢？如果你知道，那麼你很厲害了，反正開始我是一直沒轉過來。其實說起來也很簡單，如果光去看看為什麼求導以後，轉置就沒了，不太好想明白，設想一下假設現在是二維樣本點，也就是最終的W=(w1,w2)，那麼W'W=w1*w1+w2*w2那麼對w1求導就是2w1,對w2就是2w2,這樣寫在一起就是對w求導得到(2w1,2w2)=2w了，然後乘前面一個1/2（這也就是為什麼要加一個1/2），就變成w了。
好了得到上面的兩個公式，再帶回L中把去w和b消掉，你又可能發現，w確實可以消，因為有等式關係，那b怎麼辦？上述對b求導的結果竟然不含有b，上天在開玩笑嗎？其實沒有，雖然沒有b，但是有那個求和為0呀，帶進去你會驚人的發現，b還真的可以消掉，就是因為了那個等式。簡單帶下：
那麼求解最最開始的函式的最小值等價到這一步以後就是求解W的最大值了，因為使用了拉格朗日乘子法後，原問題就變為其對偶問題了，最小變成了最大，至於為什麼，等到詳細研究過對偶問題再來解釋吧。不瞭解的，只需要知道求W的極值即可。整理一下，經過這麼一圈的轉化，最終的問題為：
為什麼有ai >0$,這是上節說到的KKT條件的必須。至此問題來源部分到這。
細心的你肯可能會發現，上述所有的構造等等都是在資料完全線性可分，且分介面完全將兩類分開，那麼如果出現了下面這種情況：
正負兩類的最遠點沒有明顯的分解面，搞不好正類的最遠點反而會跑到負類裡面去了，負類最遠點跑到正類裡面去了，要是這樣的話，你的分介面都找不到，因為你不可能找到將它們完全分開的分介面，那麼這些點在實際情況是有的，就是一些離群點或者噪聲點，因為這一些點導致整個系統用不了。當然如果不做任何處理確實用不了，但是我們處理一下就可以用了。SVM考慮到這種情況，所以在上下分介面上加入鬆弛變數e,認為如果正類中有點到上介面的距離小於e，那麼認為他是正常的點，哪怕它在上介面稍微偏下一點的位置，同理下介面。還是以上面的情況，我們目測下的是理想的分解面應該是下面這種情況：
如果按照這種分會發現4個離群點，他們到自己對應分介面的距離表示如上，理論上講，我們給每一個點都給一個自己的鬆弛變數ei，如果一個分介面求出來了，那麼比較這個點到自己對應的介面（上、下介面）的距離是不是小於這個值，要是小於這個值，就認為這個介面分的可以，比如上面的e3這個點，雖然看到明顯偏離了正軌，但是計算發現它的距離d小於等於我們給的e3，那麼我們說這個分介面可以接受。你可能會說那像上面的e10，距離那麼遠了，他肯定是大於預設給這個點的ei了對吧，確實是這樣的，但是我們還發現什麼？這個點是分對了的點呀，所以你管他大不大於預設值，反正不用調整分介面。需要調整分介面的情況是隻有當類似e3這樣的點的距離大於了e3的時候。
你發現目標函式裡面多了一點東西，而加上這個是合理的，我們在優化的同時，也使得總的鬆弛變數之和最小。常數C決定了鬆弛變數之和的影響程度，如果越大，影響越嚴重，那麼在優化的時候會更多的注重所有點到分介面的距離，優先保證這個和小。好了將問題寫在一起吧：
三、SMO演算法原理與實戰求解
上節我們討論到解SVM問題最終演化為求下列帶約束條件的問題：
問題的解就是找到一組使得W最小。
現在我們來看看最初的約束條件是什麼樣子的：
這是最初的一堆約束條件吧，現在有多少個約束條件就會有多少個
我們知道
某一次迭代後，分類面為粗藍線所示，上下距離為1的分介面如細藍線所示，而理想的分介面如紫虛線所示。那麼我們想想，要想把粗藍線變化到紫虛線，在這一次是哪些個點在起作用？很顯然是界於細藍線邊上以及它們之間的所有樣本點在起作用吧，而對於那些在細藍線之外的點，比如正類的四個圈和反類的三個叉，它們在這一次的分類中就已經分對了，那還考慮它們幹什麼？所以這一次就不用考慮這些分對了的點。那麼我們用數學公式可以看到，對於在這一次就分對了的點，它們滿足什麼關係，顯然，那麼顯然它們的
那麼滿足KKT條件的，我們說如果一個點滿足KKT條件，那麼它就不需要調整，一旦不滿足，就需要調整。由上可知，不滿足KKT條件的也有三種情況：
至此我們可以說，簡單的，線性的，帶有鬆弛條件（可以容錯的）的整個SMO演算法就完了，剩下的就是迴圈，選擇兩個
當然了，這裡面還有些問題就是如何去優化這些步驟，最明顯的就是怎麼去選擇這兩個
先不管那麼多，就先讓他盲目盲目的找吧，設定一個代數，盲目到一定代數停止就ok了，下面就來一個盲目找
我的樣本是這樣的：
程式如下：
%%
% * svm 簡單演算法設計
%
%% 載入資料
% * 最終data格式：m*n，m樣本數，n維度
% * label:m*1 標籤必須為-1與1這兩類
clc
clear
close all
data = load('data_test2.mat');
data = data.data;
train_data = data(1:end-1,:)';
label = data(end,:)';
[num_data,d] = size(train_data);
data = train_data;
%% 定義向量機引數
alphas = zeros(num_data,1);
% 係數
b = 0;
% 鬆弛變數影響因子
C = 0.6;
iter = 0;
max_iter = 40;
%%
while iter < max_iter
alpha_change = 0;
for i = 1:num_data
%輸出目標值
pre_Li = (alphas.*label)'*(data*data(i,:)') + b;
%樣本i誤差
Ei = pre_Li - label(i);
% 滿足KKT條件
if (label(i)*Ei<-0.001 && alphas(i)<C)||(label(i)*Ei>0.001 && alphas(i)>0)
% 選擇一個和 i 不相同的待改變的alpha(2)--alpha(j)
j = randi(num_data,1);
if j == i
temp = 1;
while temp
j = randi(num_data,1);
if j ~= i
temp = 0;
end
end
end
% 樣本j的輸出值
pre_Lj = (alphas.*label)'*(data*data(j,:)') + b;
%樣本j誤差
Ej = pre_Lj - label(j);
%更新上下限
if label(i) ~= label(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));
end
if L==H %上下限一樣結束本次迴圈
continue;end
%計算eta
eta = 2*data(i,:)*data(j,:)'- data(i,:)*data(i,:)' - ...
data(j,:)*data(j,:)';
%儲存舊值
alphasI_old = alphas(i);
alphasJ_old = alphas(j);
%更新alpha(2)，也就是alpha(j)
alphas(j) = alphas(j) - label(j)*(Ei-Ej)/eta;
%限制範圍
if alphas(j) > H
alphas(j) = H;
elseif alphas(j) < L
alphas(j) = L;
end
%如果alpha(j)沒怎麼改變，結束本次迴圈
if abs(alphas(j) - alphasJ_old)<1e-4
continue;end
%更新alpha(1)，也就是alpha(i)
alphas(i) = alphas(i) + label(i)*label(j)*(alphasJ_old-alphas(j));
%更新系數b
b1 = b - Ei - label(i)*(alphas(i)-alphasI_old)*data(i,:)*data(i,:)'-...
label(j)*(alphas(j)-alphasJ_old)*data(i,:)*data(j,:)';
b2 = b - Ej - label(i)*(alphas(i)-alphasI_old)*data(i,:)*data(j,:)'-...
label(j)*(alphas(j)-alphasJ_old)*data(j,:)*data(j,:)';
%b的幾種選擇機制
if alphas(i)>0 && alphas(i)<C
b = b1;
elseif alphas(j)>0 && alphas(j)<C
b = b2;
else
b = (b1+b2)/2;
end
%確定更新了，記錄一次
alpha_change = alpha_change + 1;
end
end
% 沒有實行alpha交換，迭代加1
if alpha_change == 0
iter = iter + 1;
%實行了交換，迭代清0
else
iter = 0;
end
disp(['iter ================== ',num2str(iter)]);
end
%% 計算權值W
W = (alphas.*label)'*data;
%記錄支援向量位置
index_sup = find(alphas ~= 0);
%計算預測結果
predict = (alphas.*label)'*(data*data') + b;
predict = sign(predict);
%% 顯示結果
figure;
index1 = find(predict==-1);
data1 = (data(index1,:))';
plot(data1(1,:),data1(2,:),'+r');
hold on
index2 = find(predict==1);
data2 = (data(index2,:))';
plot(data2(1,:),data2(2,:),'*');
hold on
dataw = (data(index_sup,:))';
plot(dataw(1,:),dataw(2,:),'og','LineWidth',2);
% 畫出分介面，以及b上下正負1的分介面
hold on
k = -W(1)/W(2);
x = 1:0.1:5;
y = k*x + b;
plot(x,y,x,y-1,'r--',x,y+1,'r--');
title(['鬆弛變數範圍C = ',num2str(C)]);
程式中設定了鬆弛變數前的係數C是事先規定的，表明鬆弛變數項的影響程度大小。下面是幾個不同C下的結果：
這是80個樣本點，matlab下還是挺快2秒左右就好了。上圖中，把真實分介面，上下範圍為1的介面，以及那些α不為0的點（綠色標出）都畫了出來，可以看到，C越大，距離越起作用，那麼落在分界線之間的點就越少。同時可以看到，三種情況下，真實的分介面（藍色）都可以將兩種樣本完全分開（我的樣本並沒有重疊，也就是完全是可分的）。
好了，這就是隨機選取α的實驗，第一個α是按順序遍歷所有的α，第二個α是在剩下的α中在隨機選一個。當第二個α選了iter次還沒有發現不滿足KKT條件的，就退出迴圈。同時程式中的KKT條件略有不同，不過是一樣的。下面介紹如何進行啟發式的選取α呢？我們分析分析，比如上一次的一些點的α在0-C之間，也就是這些點不滿足條件需要調整，那麼一次迴圈後，他調整了一點，在下一次中這些點是不是還是更有可能不滿足條件，因為每一次的調整比較不可能完全。而那些在上一次本身滿足條件的點，那麼在下一次後其實還是更有可能滿足條件的。所以在啟發式的尋找α過程中，我們並不是遍歷所有的點的α，而是遍歷那些在0-C之間的α，而0-C反應到點上就是那些屬於邊界之間的點是不是。當某次遍歷在0-C之間找不到α了，那麼我們再去整體遍歷一次，這樣就又會出現屬於邊界之間α了，然後再去遍歷這些α，如此迴圈。那麼在遍歷屬於邊界之間α的時候，因為是需要選兩個α的，第一個可以隨便選，那第二個呢？這裡在用一個啟發式的思想，第1個α選擇後，其對應的點與實際標籤是不是有一個誤差，屬於邊界之間α的所以點每個點都會有一個自己的誤差，這個時候選擇剩下的點與第一個α點產生誤差之差最大的那個點。
程式如下：
%%
% * svm 簡單演算法設計 --啟發式選擇
%
%% 載入資料
% * 最終data格式：m*n，m樣本數，n維度
% * label:m*1 標籤必須為-1與1這兩類
clc
clear
close all
data = load('data_test2.mat');
data = data.data;
train_data = data(1:end-1,:)';
label = data(end,:)';
[num_data,d] = size(train_data);
data = train_data;
%% 定義向量機引數
alphas = zeros(num_data,1);
b = 0;
error = zeros(num_data,2);
tol = 0.001;
C = 0.6;
iter = 0;
max_iter = 40;
%%
alpha_change = 0;
entireSet = 1;%作為一個標記看是選擇全遍歷還是部分遍歷
while (iter < max_iter) && ((alpha_change > 0) || entireSet)
alpha_change = 0;
%% -----------全遍歷樣本-------------------------
if entireSet
for i = 1:num_data
Ei = calEk(data,alphas,label,b,i);%計算誤差
if (label(i)*Ei<-0.001 && alphas(i)<C)||...
(label(i)*Ei>0.001 && alphas(i)>0)
%選擇下一個alphas
[j,Ej] = select(i,data,num_data,alphas,label,b,C,Ei,entireSet);
alpha_I_old = alphas(i);
alpha_J_old = alphas(j);
if label(i) ~= label(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));
end
if L==H
continue;end
eta = 2*data(i,:)*data(j,:)'- data(i,:)*...
data(i,:)' - data(j,:)*data(j,:)';
if eta >= 0
continue;end
alphas(j) = alphas(j) - label(j)*(Ei-Ej)/eta;
%限制範圍
if alphas(j) > H
alphas(j) = H;
elseif alphas(j) < L
alphas(j) = L;
end
if abs(alphas(j) - alpha_J_old) < 1e-4
continue;end
alphas(i) = alphas(i) + label(i)*...
label(j)*(alpha_J_old-alphas(j));
b1 = b - Ei - label(i)*(alphas(i)-alpha_I_old)*...
data(i,:)*data(i,:)'- label(j)*...
(alphas(j)-alpha_J_old)*data(i,:)*data(j,:)';
b2 = b - Ej - label(i)*(alphas(i)-alpha_I_old)*...
data(i,:)*data(j,:)'- label(j)*...
(alphas(j)-alpha_J_old)*data(j,:)*data(j,:)';
if (alphas(i) > 0) && (alphas(i) < C)
b = b1;
elseif (alphas(j) > 0) && (alphas(j) < C)
b = b2;
else
b = (b1+b2)/2;
end
alpha_change = alpha_change + 1;
end
end
iter = iter + 1;
%% --------------部分遍歷(alphas=0~C)的樣本--------------------------
else
index = find(alphas>0 & alphas < C);
for ii = 1:length(index)
i = index(ii);
Ei = calEk(data,alphas,label,b,i);%計算誤差
if (label(i)*Ei<-0.001 && alphas(i)<C)||...
(label(i)*Ei>0.001 && alphas(i)>0)
%選擇下一個樣本
[j,Ej] = select(i,data,num_data,alphas,label,b,C,Ei,entireSet);
alpha_I_old = alphas(i);
alpha_J_old = alphas(j);
if label(i) ~= label(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));
end
if L==H
continue;end
eta = 2*data(i,:)*data(j,:)'- data(i,:)*...
data(i,:)' - data(j,:)*data(j,:)';
if eta >= 0
continue;end
alphas(j) = alphas(j) - label(j)*(Ei-Ej)/eta;
%限制範圍
if alphas(j) > H
alphas(j) = H;
elseif alphas(j) < L
alphas(j) = L;
end
if abs(alphas(j) - alpha_J_old) < 1e-4
continue;end
alphas(i) = alphas(i) + label(i)*...
label(j)*(alpha_J_old-alphas(j));
b1 = b - Ei - label(i)*(alphas(i)-alpha_I_old)*...
data(i,:)*data(i,:)'- label(j)*...
(alphas(j)-alpha_J_old)*data(i,:)*data(j,:)';
b2 = b - Ej - label(i)*(alphas(i)-alpha_I_old)*...
data(i,:)*data(j,:)'- label(j)*...
(alphas(j)-alpha_J_old)*data(j,:)*data(j,:)';
if (alphas(i) > 0) && (alphas(i) < C)
b = b1;
elseif (alphas(j) > 0) && (alphas(j) < C)
b = b2;
else
b = (b1+b2)/2;
end
alpha_change = alpha_change + 1;
end
end
iter = iter + 1;
end
%% --------------------------------
if entireSet %第一次全遍歷了，下一次就變成部分遍歷
entireSet = 0;
elseif alpha_change == 0
%如果部分遍歷所有都沒有找到需要交換的alpha，再改為全遍歷
entireSet = 1;
end
disp(['iter ================== ',num2str(iter)]);
end
%% 計算權值W
W = (alphas.*label)'*data;
%記錄支援向量位置
index_sup = find(alphas ~= 0);
%計算預測結果
predict = (alphas.*label)'*(data*data') + b;
predict = sign(predict);
%% 顯示結果
figure;
index1 = find(predict==-1);
data1 = (data(index1,:))';
plot(data1(1,:),data1(2,:),'+r');
hold on
index2 = find(predict==1);
data2 = (data(index2,:))';
plot(data2(1,:),data2(2,:),'*');
hold on
dataw = (data(index_sup,:))';
plot(dataw(1,:),dataw(2,:),'og','LineWidth',2);
% 畫出分介面，以及b上下正負1的分介面
hold on
k = -W(1)/W(2);
x = 1:0.1:5;
y = k*x + b;
plot(x,y,x,y-1,'r--',x,y+1,'r--');
title(['鬆弛變數範圍C = ',num2str(C)]);
其中的子函式，一個是計算誤差函式，一個是選擇函式如下：
function Ek = calEk(data,alphas,label,b,k)
pre_Li = (alphas.*label)'*(data*data(k,:)') + b;
Ek = pre_Li - label(k);
function [J,Ej] = select(i,data,num_data,alphas,label,b,C,Ei,choose)
maxDeltaE = 0;maxJ = -1;
if choose == 1 %全遍歷---隨機選擇alphas
j = randi(num_data ,1);
if j == i
temp = 1;
while temp
j = randi(num_data,1);
if j ~= i
temp = 0;
end
end
end
J = j;
Ej = calEk(data,alphas,label,b,J);
else %部分遍歷--啟發式的選擇alphas
index = find(alphas>0 & alphas < C);
for k = 1:length(index)
if i == index(k)
continue;
end
temp_e = calEk(data,alphas,label,b,k);
deltaE = abs(Ei - temp_e); %選擇與Ei誤差最大的alphas
if deltaE > maxDeltaE
maxJ = k;
maxDeltaE = deltaE;
Ej = temp_e;
end
end
J = maxJ;
end
至此算是完了，試驗了一下，兩者的效果其實差不多（反而隨機選取的效果更好一點，感覺是因為隨機保證了更多的可能，畢竟隨機選擇包括了你的特殊選擇，但是特殊選擇到後期是特殊不起來的，反而隨機會把那些差一點的選擇出來），但是速度上當樣本小的時候，基本上差不多，但是當樣本大的時候，啟發式的特殊選擇明顯佔優勢了。我試驗了400個樣本點的情況，隨機選擇10多秒把，而啟發式2,3秒就好了。可見效果差不多的情況下，啟發式選擇是首要選擇。
至此兩種方式下的方法都實驗完了。那麼我們看到，前面(三節)所講的一切以及實驗，分類的樣本都是線性樣本，那麼如果來了非線性樣本該如何呢？而SVM的強大之處更在於對非線性樣本的準確劃分。那麼前面的理論對於非線性樣本是否適用？我們又該如何處理非線性樣本呢？請看下節SVM非線性樣本的分類。
四、SVM非線性分類原理實驗
前面幾節我們討論了SVM原理、求解線性分類下SVM的SMO方法。本節將分析SVM處理非線性分類的相關問題。
一般的非線性分類如下左所示（後面我們將實戰下面這種情況）：
可以看到在原始空間中你想用一個直線分類面劃分開來是不可能了，除非圓。而當你把資料點對映一下成右圖所示的情況後，現在資料點明顯看上去是線性可分的，那麼在這個空間上的資料點我們再用前面的SVM演算法去處理，就可以得到每個資料點的分類情況了，而這個分類情況也是我們在低維空間的情況。也就是說，單純的SVM是不能處理非線性問題的，說白了只能處理線性問題，但是來了非線性樣本怎麼辦呢？我們是在樣本上做的文章，我把非線性樣本變成線性樣本，再去把變化後的線性樣本拿去分類，經過這麼一圈，就達到了把非線性樣本分開的目的，所以只看開頭和結尾的話發現，SVM竟然可以分非線性問題，其實呢還是分的線性問題。
現在的問題是如何找到這個對映關係對吧，就比如上面那個情況，我們可以人為計算出這種對映，比如一個樣本點是用座標表示的(x1,x2),它有個類標籤，假設為1，那麼把這個點對映到三維中變成,對每個點我都這麼去對映，假設一個原始點樣本集是這樣的：
然後按照上面那個公式去把每個點對映成3維座標點後，畫出來是這樣的：
可以看到是線性可分的吧，如果還看不清把視角換個角度（右檢視）：
現在能看清楚了吧。那這是二維的點到三維，對映的關係就是上面的那個關係，那如果是三維到四維，四維到N維呢？這個關係你還想去找嗎？理論上是找的到的，但是實際上人工去找你怎麼去找？你怎麼知道資料的對映關係是這樣的是那樣的？不可能知道。然而我們真的需要找到這種關係嗎？答案是不需要的，返回去看看前三節的關於SVM的理論部分可以看到，無論是計算a呀，還是b呀等等，只要涉及到原始資料點的，都是以內積的形式出來的，也就是說是一個點的向量與另一個點的向量相乘的，向量內積出來是一個值。
就拿a來更新來說，如下：
最後也是得到一個值比如C2。既然SVM裡面所有涉及到原始資料的地方都是以向量的形式出現的，那麼我們還需要管它的對映關係嗎？因為它也不需要你去計算說具體到比如說三維以後，三維裡面的三個座標值究竟是多少，他需要的是內積以後的一個結果值。那麼好辦了，我就假設有一個黑匣子，輸入原始資料維度下的兩個座標向量，然後經過黑匣子這麼一圈，出來一個值，這個值我們就認為是高維度下的值。而黑匣子的潛在意義就相當於一個高維對映器一樣。更重要的是我們並不需要知道黑匣子究竟是怎麼對映的，只需要知道它的低緯度下的形式就可以了。常用的黑匣子就是徑向基函式，而這個黑匣子在數學上就叫做核函式，例如徑向基函式的外在形式如下所示：
o是需要預先設定的引數。至於這個黑匣子把初始資料對映到多少維了，誰知道呢，既然是黑匣子，那就是看不到的，上帝給了人類這麼一個黑匣子就已經很夠意思了。可以看到的是原始資料結果黑匣子算了以後，出來就是一個值了，而這個值就認為是高維度下的資料通過內積計算而來的值。當然上帝還留了一個窗戶，就是
好了既然黑匣子是藏著的，那也就只能說這麼多了。有趣的是上帝給的這個黑匣子不止一個，有好幾個，只是上面的那個普遍效果更好而已。基於此，那麼對於上節的SMO演算法，如果拿來求解非線性資料的話，我們只需要將其中對應的內積部分改成核函式的形式即可。一個資料核函式程式如下：
function result = Kernel(data1,data2,sigma)
% data裡面每一行資料是一個樣本(的行向量)
[m1,~] = size(data1);
[m2,~] = size(data2);
result = zeros(m1,m2);
for i = 1:m1
for j = 1:m2
result(i,j) = exp(-norm(data1(i,:)-data2(j,:))/(2*sigma^2));
end
end
有了此核函式，我們用上節的隨機遍歷α
然後把主程式對應的部分用上述核函式代替：
%%
% * svm 簡單演算法設計
%
%% 載入資料
% * 最終data格式：m*n，m樣本數，n維度
% * label:m*1 標籤必須為-1與1這兩類
clc
clear
% close all
data = load('data_test1.mat');
data = data.data;
train_data = data(1:end-1,:)';
label = data(end,:)';
[num_data,d] = size(train_data);
data = train_data;
%% 定義向量機引數
alphas = zeros(num_data,1);
% 係數
b = 0;
% 鬆弛變數影響因子
C = 0.6;
iter = 0;
max_iter = 80;
% 核函式的引數
sigma = 4;
%%
while iter < max_iter
alpha_change = 0;
for i = 1:num_data
%輸出目標值
pre_Li = (alphas.*label)'*Kernel(data,data(i,:),sigma) + b;
%樣本i誤差
Ei = pre_Li - label(i);
% 滿足KKT條件
if (label(i)*Ei<-0.001 && alphas(i)<C)||(label(i)*Ei>0.001 && alphas(i)>0)
% 選擇一個和 i 不相同的待改變的alpha(2)--alpha(j)
j = randi(num_data,1);
if j == i
temp = 1;
while temp
j = randi(num_data,1);
if j ~= i
temp = 0;
end
end
end
% 樣本j的輸出值
pre_Lj = (alphas.*label)'*Kernel(data,data(j,:),sigma) + b;
%樣本j誤差
Ej = pre_Lj - label(j);
%更新上下限
if label(i) ~= label(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));
end
if L==H %上下限一樣結束本次迴圈
continue;end
%計算eta
eta = 2*Kernel(data(i,:),data(j,:),sigma)- ...
Kernel(data(i,:),data(i,:),sigma)...
- Kernel(data(j,:),data(j,:),sigma);
%儲存舊值
alphasI_old = alphas(i);
alphasJ_old = alphas(j);
%更新alpha(2)，也就是alpha(j)
alphas(j) = alphas(j) - label(j)*(Ei-Ej)/eta;
%限制範圍
if alphas(j) > H
alphas(j) = H;
elseif alphas(j) < L
alphas(j) = L;
end
%如果alpha(j)沒怎麼改變，結束本次迴圈
if abs(alphas(j) - alphasJ_old)<1e-4
continue;end
%更新alpha(1)，也就是alpha(i)
alphas(i) = alphas(i) + label(i)*label(j)*(alphasJ_old-alphas(j));
%更新系數b
b1 = b - Ei - label(i)*(alphas(i)-alphasI_old)*...
Kernel(data(i,:),data(i,:),sigma) - label(j)*...
(alphas(j)-alphasJ_old)*Kernel(data(i,:),data(j,:),sigma);
b2 = b - Ej - label(i)*(alphas(i)-alphasI_old)*...
Kernel(data(i,:),data(j,:),sigma)- label(j)*...
(alphas(j)-alphasJ_old)*Kernel(data(j,:),data(j,:),sigma);
%b的幾種選擇機制
if alphas(i)>0 && alphas(i)<C
b = b1;
elseif alphas(j)>0 && alphas(j)<C
b = b2;
else
b = (b1+b2)/2;
end
%確定更新了，記錄一次
alpha_change = alpha_change + 1;
end
end
% 沒有實行alpha交換，迭代加1
if alpha_change == 0
iter = iter + 1;
%實行了交換，迭代清0
else
iter = 0;
end
disp(['iter ================== ',num2str(iter)]);
end
%% 計算權值W
% W = (alphas.*label)'*data;
%記錄支援向量位置
index_sup = find(alphas ~= 0);
%計算預測結果
predict = (alphas.*label)'*Kernel(data,data,sigma) + b;
predict = sign(predict);
%% 顯示結果
figure;
index1 = find(predict==-1);
data1 = (data(index1,:))';
plot(data1(1,:),data1(2,:),'+r');
hold on
index2 = find(predict==1);
data2 = (data(index2,:))';
plot(data2(1,:),data2(2,:),'*');
hold on
dataw = (data(index_sup,:))';
plot(dataw(1,:),dataw(2,:),'og','LineWidth',2);
title(['核函式引數sigma = ',num2str(sigma)]);
下面是幾個不同引數下的結果：
可以看到σ到4以後就分不出來了。綠色的為支援向量，可以看到在σ在0.6到1之間是最少的，結果應該也是最好的。至此SMO實驗非線性樣本完畢。
當今學者已經有非常多的人研究SVM演算法，同時開發了許多開源的程式，這些程式都是經過不斷優化的，效能比起我們這裡自己編的來說要好得多，所以在實際應用中通常都是用他們無私貢獻的軟體包。一個典型的軟體包就是臺灣一個教授團隊的LIBSVM軟體包，那麼你是否想一窺其用法，看看它的效能如何呢？請看下節matlab下LIBSVM的簡單使用。
五、MATLAB下libsvm的簡單使用：分類與迴歸
本節簡單介紹一下libsvm的使用方法。關於libsvm似乎曾經使用過，那個時候主要用libsvm進行簡單的人臉識別實驗。
1）介紹與分類實驗
那麼現在最新版本的libsvm為3.2.0，下載地址如下： http://www.csie.ntu.edu.tw/~cjlin/libsvm/
下載下來的libsvm其實包含好多個平臺的工具箱軟體，c++，matlab，java，python都有。他們的函式使用方法是一樣的。
那麼在下載完以後，點選裡面的matlab下平臺，直接在點選裡面的make.m函式就可以了。正常情況下如果你的matlab含有編譯平臺的話直接就可以執行了，如果沒有，還需要選擇一個平臺 mex -setup 。小提醒一下，這個編譯過程不要在c盤下使用，也就是libsvm先不要放在c盤，涉及到許可權，機器不讓編譯。編譯完後在matlab的設定路徑中新增進去編譯的資料夾及其內容，那麼就可以使用了。正常編譯的過程是這樣的： 在上面的人臉識別實驗中曾經介紹過裡面的主要函式，這裡為了放在一塊，把那裡的拿過來吧：
目前版LIBSVM（3.2.0）在matlab下編譯完後只有四個函式，libsvmread，Libsvmwrite，svmtrain（matlab自帶的工具箱中有一個同名的函式），svmpredict。
libsvmread主要用於讀取資料
這裡的資料是非matlab下的.mat資料，比如說是.txt，.data等等，這個時候需要使用libsvmread函式進行轉化為matlab可識別資料，比如自帶的資料是heart_scale資料，那麼匯入到matlab有兩種方式，一種使用libsvmread函式，在matlab下直接libsvmread(heart_scale);第二種方式為點選matlab的‘匯入資料’按鈕，然後導向heart_scale所在位置，直接選擇就可以了。個人感覺第二種方式超級棒，無論對於什麼資料，比如你在哪個資料庫下下載的資料，如何把它變成matlab下資料呢？因為有的資料libsvmread讀取不管用，但是‘匯入資料’後就可以變成matlab下資料。
libsvmwrite寫函式，就是把已知資料存起來
使用方式為：libsvmwrite(‘filename’,label_vector, instance_matrix);
label_vector是標籤，instance_matrix為資料矩陣（注意這個資料必須是稀疏矩陣，就是裡面的資料不包含沒用的資料（比如很多0），有這樣的資料應該去掉再存）。
svmtrain訓練函式，訓練資料產生模型的
一般直接使用為：model=svmtrain(label,data,cmd); label為標籤，data為訓練資料（資料有講究，每一行為一個樣本的所有資料，列數代表的是樣本的個數），每一個樣本都要對應一個標籤（分類問題的話一般為二分類問題，也就是每一個樣本對應一個標籤）。cmd為相應的命令集合，都有哪些命令呢？很多，-v，-t,-g,-c,等等，不同的引數代表的含義不同，比如對於分類問題，這裡-t就表示選擇的核函式型別，-t=0時線性核。-t=1多項式核，-t=2，徑向基函式（高斯），-t=3，sigmod核函式，新版出了個-t=4，預計算核（還不會用）；-g為核函式的引數係數，-c為懲罰因子係數，-v為交叉驗證的數，預設為5，這個引數在svmtrain寫出來使用與不寫出來不使用的時候，model出來的東西不一樣，不寫的時候，model為一個結構體，是一個模型，可以帶到svmpredict中直接使用，寫出來的時候，出來的是一個訓練模型的準確率，為一個數值。一般情況下就這幾個引數重要些，還有好多其他引數，可以自己參考網上比較全的，因為下面的這種方法的人臉識別就用到這麼幾個引數，其他的就不寫了。
svmpredict訓練函式，使用訓練的模型去預測來的資料型別。
使用方式為：
[predicted_label,accuracy,decision_values/prob_estimates]= svmpredict(testing_label_vector,testing_instance_matrix,model,’libsvm_options’)
或者：
[predicted_label]=svmpredict(testing_label_vector,testing_instance_matrix, model, ‘libsvm_options’)
第一種方式中，輸出為三個引數，預測的型別，準確率，評估值（非分類問題用著），輸入為測試型別（這個可與可無，如果沒有，那麼預測的準確率accuracy就沒有意義了，如果有，那麼就可以通過這個值與預測出來的那個型別值相比較得出準確率accuracy，但是要說明一點的是，無論這個值有沒有，在使用的時候都得加上，即使沒有，也要隨便加上一個型別值，反正你也不管它對不對，這是函式使用所規定的的），再就是輸入資料值，最後是引數值（這裡的引數值只有兩種選擇，-p和-b引數），曾經遇到一個這樣的問題，比如說我在訓練函式中規定了-g引數為0.1，那麼在預測的時候是不是也要規定這個引數呢？當你規定了以後，程式反而錯誤，提醒沒有svmpredict的-g引數，原因是在svmtrain後會出現一個model，而在svmpredict中你已經用了這個model，而這個model中就已經包含了你所有的訓練引數了，所以svmpredict中沒有這個引數，那麼對於的libsvm_options就是-p和-b引數了。對於函式的輸出，兩種方式呼叫的方法不一樣，第一種呼叫把所有需要的資料都呼叫出來了，二第二種呼叫，只呼叫了predicted_label預測的型別，這裡我們可以看到，在單純的分類預測模型中，其實第二種方式更好一些吧，既簡單有實用。
致此，四個函式在分類問題中的介紹大概如此，當然還有很多可以優化的細節就不詳細說了，比如可以再使用那些引數的時候，你如果不規定引數的話，所有的-引數都是使用預設的，預設的就可能不是最好的吧，這樣就涉及到如何去優化這個引數了。
使用就介紹到這裡吧，下面實戰一下，樣本集選擇前面使用的200個非線性樣本集，函式如下：
%%
% * libsvm 工具箱簡單使用
%
%% 載入資料
% * 最終data格式：m*n，m樣本數，n維度
% * label:m*1 標籤為-1與1這兩類
clc
clear
close all
data = load('data_test1.mat');
data = data.data';
%選擇訓練樣本個數
num_train = 80;
%構造隨機選擇序列
choose = randperm(length(data));
train_data = data(choose(1:num_train),:);
gscatter(train_data(:,1),train_data(:,2),train_data(:,3));
label_train = train_data(:,end);
test_data = data(choose(num_train+1:end),:);
label_test = test_data(:,end);
predict = zeros(length(test_data),1);
%% ----訓練模型並預測分類
model = svmtrain(label_train,train_data(:,1:end-1),'-t 2');
% -t = 2 選擇徑向基函式核
true_num = 0;
for i = 1:length(test_data)
% 作為預測，svmpredict第一個引數隨便給個就可以
predict(i) = svmpredict(1,test_data(i,1:end-1),model);
end
%% 顯示結果
figure;
index1 = find(predict==1);
data1 = (test_data(index1,:))';
plot(data1(1,:),data1(2,:),'or');
hold on
index2 = find(predict==-1);
data2 = (test_data(index2,:))';
plot(data2(1,:),data2(2,:),'*');
hold on
indexw = find(predict~=(label_test));
dataw = (test_data(indexw,:))';
plot(dataw(1,:),dataw(2,:),'+g','LineWidth',3);
accuracy = length(find(predict==label_test))/length(test_data);
title(['predict the testing data and the accuracy is :',num2str(accuracy)]);
可以看到，關於svm的部分就那麼一點，其他的都是輔助吧，那麼一個結果如下：
資料人為設定了一些重疊，這個結果算是非常好了。當然對於libsvm函式，裡面還有許多細節，像引數選擇等等，不同的引數結果是不一樣的，這就待你去探究了。
2）迴歸實驗
迴歸問題不像分類問題，迴歸問題相當於根據訓練樣本訓練出一個擬合函式一樣，可以根據這個擬合函式可以來預測給定一個樣本的輸出值。可以看到分類問題輸出的是樣本所屬於的類，而回歸問題輸出的是樣本的預測值。
常用的地方典型的比如股票預測，人口預測等等此類預測問題。
libsvm同樣可以進行迴歸預測，所需要改變的只是裡面的引數設定。檢視libsvm的官網介紹引數詳情如下：
options:
-s svm_type : set type of SVM (default 0)
0 -- C-SVC
1 -- nu-SVC
2 -- one-class SVM
3 -- epsilon-SVR
4 -- nu-SVR
-t kernel_type : set type of kernel function (default 2)
0 -- linear: u'*v
1 -- polynomial: (gamma*u'*v + coef0)^degree
2 -- radial basis function: exp(-gamma*|u-v|^2)
3 -- sigmoid: tanh(gamma*u'*v + coef0)
-d degree : set degree in kernel function (default 3)
-g gamma : set gamma in kernel function (default 1/num_features)
-r coef0 : set coef0 in kernel function (default 0)
-c cost : set the parameter C of C-SVC, epsilon-SVR, and nu-SVR (default 1)
-n nu : set the parameter nu of nu-SVC, one-class SVM, and nu-SVR (default 0.5)
-p epsilon : set the epsilon in loss function of epsilon-SVR (default 0.1)
-m cachesize : set cache memory size in MB (default 100)
-e epsilon : set tolerance of termination criterion (default 0.001)
-h shrinking: whether to use the shrinking heuristics, 0 or 1 (default 1)
-b probability_estimates: whether to train a SVC or SVR model for probability estimates, 0 or 1 (default 0)
-wi weight: set the parameter C of class i to weight*C, for C-SVC (default 1)
可以看到-s svm_type 控制的就是訓練型別，而當-s等於3或4的時候，就是迴歸模型SVR。
-s 3 就是常用的帶懲罰項的 SVR模型，我們用這個實驗。我使用的是libsvm3.2.0工具箱，版本不同可能會帶來呼叫方式的不同。測試實驗的程式碼如下，可能會有一些細節需要自己去探索：
close all;
clear;
clc;
%%
% 生成待迴歸的資料
x = (-1:0.1:1)';
y = -100*x.^3 + x.^2 - x + 1;
% 加點噪聲
y = y+ 20*rand(length(y),1);
%% 採用交叉驗證選擇引數
mse = 10^7;
for log2c = -10:0.5:3,
for log2g = -10:0.5:3,
% -v 交叉驗證引數：在訓練的時候需要，測試的時候不需要，否則出錯
cmd = ['-v 3 -c ', num2str(2^log2c), ' -g ', num2str(2^log2g) , ' -s 3 -p 0.4 -t 3'];
cv = svmtrain(y,x,cmd);
if (cv < mse),
mse = cv; bestc = 2^log2c; bestg = 2^log2g;
end
end
end
%% 訓練--
cmd = ['-c ', num2str(2^bestc), ' -g ', num2str(2^bestg) , ' -s 3 -p 0.4 -n 0.1'];
model = svmtrain(y,x,cmd);
% model
% 利用建立的模型看其在訓練集合上的迴歸效果
% 注意libsvm3.2.0的svmpredict函式必須有三個引數輸出
[py,~,~] = svmpredict(y,x,model);
figure;
plot(x,y,'o');
hold on;
plot(x,py,'g+');
%%
% 進行預測新的x值
%-- 產生[-1 1]的隨機數
testx = -2+(2-(-2))*rand(10,1);
testy = zeros(10,1);% 理論y值無所謂
[ptesty,~,~] = svmpredict(testy,testx,model);
hold on;
plot(testx,ptesty,'r*');
legend('原始資料','迴歸資料','新資料');
grid on;
% title('t=0:線性核')
% title('t=1:多項式核')
% title('t=2:徑向基函式（高斯）')
title('t=3:sigmod核函式')
這裡我隨機生成一個3次函式的隨機資料，測試了幾種不同svm裡面的核函式：
因為我們的資料是由三次函式模擬生成的，所以可以看到，在這種情況下使用線性核t=0時候效果更好，然而實際情況下一般我們也不知道資料的分佈函式，所以在選擇核函式的時候還是需要多實驗，找到最適合自己資料的核函式。
這裡採用了交叉驗證的方式自適應選擇模型中重要的兩個引數，需要注意的是引數的範圍，不要太大，步長可能也需要控制，否則在資料量很大的時候需要執行很久。
至此SVM系列文章就到這裡吧，感謝能看到這裡的朋友~_~。
作者部落格地址：https://blog.csdn.net/on2way
AI科技大本營現招聘AI記者和資深編譯，有意者請將簡歷投至：gulei@csdn.net，期待你的加入！
如果你暫時不能加入營長的隊伍，也歡迎與營長分享你的精彩文章，投稿郵箱：suiling@csdn.net
AI科技大本營讀者群（計算機視覺、機器學習、深度學習、NLP、Python、AI硬體、AI+金融、AI+PM方向）正在招募中，關注AI科技大本營微信公眾號，後臺回覆：讀者群，聯絡營長，新增營長請備註姓名，研究方向。
☟☟☟點選 | 閱讀原文 | 檢視更多精彩內容
相關文章
- 一文讀懂支援向量機SVM（附實現程式碼、公式）公式
- opencv + SVM 程式碼OpenCV
- 學習SVM，這篇文章就夠了！（附詳細程式碼）
- LOAM原始碼分析附公式推導原始碼公式
- 支援向量機（SVM）從原理到python程式碼實現Python
- Diffusion系列 - DDIM 公式推導 + 程式碼 -（三）公式
- 好工具，SWOT分析大解密！解密
- SGDClassifier和LR,SVM的區別
- 通達信漲跌動因指標公式原始碼附圖指標公式原始碼
- SVM
- 機器學習(四)：通俗理解支援向量機SVM及程式碼實踐機器學習
- 如何開發和評估Vanilla LSTM模型？（附程式碼）模型
- 通達信騎牛找馬選股指標公式原始碼附圖指標公式原始碼
- 支援向量機（SVM）和python實現（二）Python
- 基於PCA和SVM的人臉識別PCA
- SVM原理
- 阿里面試題BIO和NIO數量問題附答案和程式碼阿里面試題
- 表白程式碼製作（附原始碼）原始碼
- 32、條件格式和公式公式
- 教你excel表格怎麼提取手機號碼和姓名，不用公式和函式Excel公式函式
- SVM 的核函式選擇和調參函式
- PySpark和SparkSQL基礎：如何利用Python程式設計執行Spark（附程式碼）SparkSQLPython程式設計
- 什麼是Python的迭代器和生成器?（附程式碼）Python
- 小明上學（201812-1/CCF)———附帶思路和完整程式碼
- 《TensorFlow 機器學習方案手冊》（附 pdf 和完整程式碼）機器學習
- 初識MapReduce的應用場景（附JAVA和Python程式碼）JavaPython
- 支援向量機（Support Vector Machine，SVM）—— 線性SVMMac
- RSI+BOLL副圖公式原始碼 通訊達公式編寫公式原始碼
- 教你excel篩選手機號碼公式和excel篩選正常電話號碼Excel公式
- 部落格園 ----LaTex使用教程，內附各種公式，建議收藏和備忘！！公式
- SVM家族（一）
- opencv中SVMOpenCV
- 200 行 Python 程式碼做個換臉程式（附原始碼）Python原始碼
- Python機器學習筆記：SVM（1）——SVM概述Python機器學習筆記
- 另類KDJ指標公式原始碼 2019通達信指標公式指標公式原始碼
- 使用 SVM 和決策樹進行整合學習
- Node和http：一本通【附tcp實現http小程式碼】HTTPTCP
- 低程式碼和無程式碼的區別