支持向量機（Support Vector Machine）

支持向量機

linear regression ， perceptron learning algorithm ， logistics regression都是分類器，咱們可使用這些分類器作線性和非線性的分類，好比下面的一個問題：

這裏的每一條線都是能夠把這個平面分開的，支持向量機要作的就是要在這些能夠選擇的直線中選擇一條最好的直線來做爲分類的直線。再給一個簡單的解釋，好比下面的三個圖片，圓圈區域越大，說明這條直線對這些點放錯的容忍度就越高：

##①超平面 介紹SVM以前，先來看超平面的概念： 其實超平面就是用於分割當前維度的一個空間。好比一維能夠用一個點來進行分割，二維用一條線來進行分割，那麼這些點和線就叫作「超」平面。加雙引號是由於他們其實並非正在的超平面，由於超平面是要求大於三維的。 因此四維空間裏面的超平面就是三維。好比在二維的空間裏的超平面就是一條直線：

三維裏面的超平面：

(其實這裏的應該都不能叫超平面，由於超平面是三維以及三維以上的) 咱們把a ， b ， c看作是W0 , W1 , W2...，把x , y , z看作是x1 , x2 , x3，那麼就有：

而W向量就是這個平面的法向量，咱們要求的就是法向量，證實一下： 以二維爲例：

相減獲得：

而[ (X_0 - X_1) , (Y_0 , Y_1)]這個點就在這個平面上，因此獲得了：

因此W就是平面的法向量，這上面的X表明的是這個平面而不是點。

##②函數間隔的最大化 剛剛說到支持向量機也不是找超平面了，而是找最好的超平面，也就是對於點的犯錯的容忍度越大越好，其實就是函數間隔越大越好：

右邊的明顯要好過左邊的，由於左邊的可犯錯空間大啊。因此咱們要尋找的就是最大最肥的超平面——hperplane。

函數的間隔最大，**其實就是距離直線距離最短的點離直線的距離要最大。**因此先要知道直線的距離怎麼求： 首先咱們假設X0在平面上的投影是X1，則確定有法向量W垂直於X0X1,：

又由於：

右由於X1在平面上的，前半部分就是-b了，能夠寫成：

和上面那條式子相等就獲得：

這就是咱們要求的點到直線的距離了。而若是這個hperplane是正確的話，那麼全部點的分類都是對的，那麼咱們就默認他是對的，因而有：

這裏能夠相乘的條件是，咱們默認label正確的是1錯誤的是-1，若是你的錯誤是0正確是1的話公式是不一樣的。乘上一個Y首先是能夠去掉絕對值，使得函數變得可微，另外乘上以後函數值的絕對值也不會有變化，使得求解更加方便。 因此，最後的咱們的優化目標就是這樣了：

裏面的minimize是指找到距離hperplane最小距離的點，最外面就是挑選一個最好的W,b使得這個距離最小的點距離hperplane是最大的。

##③目標函數的化簡 對於上面的式子，注意看到裏面的那個式子：

舉一個例子：

咱們代入(4,3)這個點，獲得19，彷佛這個數字太大了，咱們不想要他這麼大，咱們兩邊同時除去19，這個時候咱們的超平面就變成了：

數字是邊了，可是這個超平面仍是在這個位置，因此能夠認爲超平面是沒有變化的，這就證實了咱們上面那個式子：

是能夠經過對w,b進行放縮而把左邊的結果放大的無限多倍的。既然這樣，那這個東西留着有什麼意義，直接放縮到1就能夠了，因而咱們把他放縮到1，也就是最小值是1。其實等於1都是差很少的，由於最小值以後都是1，因而就是有了：

那麼target fomula就能夠變成這樣：

放縮以後並非就完了，這個是要加入當作條件來使用的。對於這個：

事實上咱們不太喜歡化簡這樣的，找最大化的不就是找最小化的W嗎？找最小化的W不就是找最小化的W*W^T嗎？不如再加個1/2? 因此問題就變成了：

爲何要加上1/2呢？實際上是爲了後面求導的時候方便化簡的，可是這樣對結果又沒有什麼影響。而W變成平方其實就是用上凸優化，由於平方以後就個凸函數了，這樣變換一樣對於最優結果沒有任何影響。 因此最後要優化的結果：

##④Dual problem and KKT condiction 對於上述有條件的最優化問題，天然就要用上lagrange乘子法了。

右邊的約束條件是要小於等於0的，α ≥ 0，只不過前面是符號因此轉一下而已。 到這裏，其實只要把等式扔到Quadratic Programming裏面直接計算就行了。下面的步驟其實都是圍繞解決這個優化問題展開的。

---

---

先在這停頓一下，咱們考慮一下：

####⑴SVM的機器學習可行性問題： 首先先來觀察一下這個式子，感受似曾相識。他和L2 regularization很像，對於L2 regularization，首先是先要計算Ein的最小值，所謂的Ein就是該模型在當前的訓練數據上犯錯誤的指望值。而後再正則化，因此L2是Minimizing Ein and Regularized L2 Paradigms；而支持向量機正好相反，他是先假設我這個平面是分類正確的，而後minimize W方：

後面咱們會用這個結論的。 回顧一下機器學習要解決的問題： ①Ein ≈ Eout Ein剛剛解釋過了，Eout就是這個model在全局所犯的錯誤，Ein ≈ Eout就是要求這個model是能夠反映全局的，若是不能反映，那就是過擬合了。 ②Ein ≈ 0 這個就是訓練錯誤要接近於0，在這一步就容易發生過擬合的現象了。 而Ein ≈ Eout，也就是泛化能力是被VC dimension限制的，也就是說，越複雜的模型他的VC dimension越複雜。也就是VC bound右邊的Ω會很大，VC bound就會很大，致使Ein 遠遠小於Eout了，由於複雜的模型意味着更加小的Ein。 再提一下，VC dimension就是break point - 1獲得的。 若是是這樣的話，那麼正常來講SVM的VC dimension也會很大啊，由於他的W是和數據維度相關的，數據多少維，那麼W就多少個，而W表明的是自由度，一般也就表明這VC dimension， 可是SVM的效果仍是很好，爲何呢？是什麼東西限制着SVM的VC dimension？ 咱們來看一個例子： 在一個圓上，有三個點，你想找到一條能夠分開的直線，能夠獲得VC dimension是3（以前有同窗看到在一個圓上分類他的VC dimension是無限的的，這是由於有無數多個點給你玩，這裏就三個點，無限你又用不了，因此就只能是三個了啦），可是若是加上限制條件，這條線寬是5，那麼VC dimension就是0了，由於沒有地方塞進去。因此若是是large margin，VC dimension ≤ 3的。如圖：

因此，large margin就是SVM的VCdimension的限制條件，致使它的分類效果很好，VC dimension小了天然泛化能力就行了，這裏就解決了Ein ≈ Eout的問題，Ein ≈ 0這就是咱們後面要用凸優化來解決的問題了。

---

---

回到正題：

如何優化咱們的target function？ 上面的討論咱們已經獲得了VC dimension ≤ d + 1，**咱們悲觀估計一下，就算SVM的VC dimension是d + 1了，加上d就是數據維度，加上的1是偏置值b。那麼若是數據維度很大的話，計算複雜度是很高的，另外，如今咱們所研究的SVM仍是linear separable的，若是是來個nonlinear transform，數據維度就更加大了，再加上通常數據數量都是很的，時間會很長。**因此咱們要想一個方法來把數據維度d轉移到其餘地方去或者之間丟了。 而Daul problem剛好就能夠解決。 回到origin target function： 咱們須要最小化：

對原函數作一些變換：

當這個點是違反了條件的時候，那麼約束條件就會 > 0，α > 0，再maximum α那麼就是無限大了，這個時候就不多是最小值，不可能取他，由於是無限大了。 當這個點是不違反的時候，那麼約束條件就會 < 0，α > 0，再maximum α約束條件就是0，minimize w，b以後就仍是原來的target function。 因此變換以後：

變換以後的問題 == origin target function

---

---

再次停頓一下，考慮一下KKT條件是什麼：

####⑵KKT 條件的引出 對於上述裝換過的target function，有以下性質：

那麼對於任何的條件都會有左邊≥右邊的。左邊再加上一個w,b取最小：

一樣是大於右邊的全部狀況，那麼天然了，我右邊再加上一個取α的最大值也是能夠的：

而在右邊的咱們把條件minimize和maximum調換過的式子就叫作Daul Problem。**因此，原問題是≥對偶問題的。**那麼有沒有什麼辦法能夠把原問題的求解轉換成對偶問題的求解呢？

---

---

####⑶KKT 條件的簡單證實

對偶的意思就是存在一個最優的解使得兩邊的等式成立。因此咱們假設有一個W和B是最優的，那麼有：

而最後能夠看到求出來的解正是咱們要求的f(W)原目標函數，而原式子：

代進去也將是這個結果，由於maximum以後ag(x) = 0，因此本質上這個函數仍是求f(W)的最小值。因此對偶式子和原式在結果上是沒有差異的。根據上面的式子，咱們本質就是要求

的最小值，固然這裏的W,B要替換成原來的變量w，b了。求最小值天然就是求梯度爲0了，因此w,b梯度爲0的條件就有了。還有一個ag(x) = 0的條件，這個實際上是前提條件：

咱們以前說這個式子是等同目標函數的，既然要等同天然是要把後面的g(x)和h(x)消去啊！而h(x) = 0原本就消去了，而g(x) < 0，求最大必然就ag(x) = 0了，由於只有這個條件，才能消去後面的ag(x)把這個minimum maximum式子變成minimumf(w)的式子。因此再加上先前的幾個拉格朗日條件就組成了KKT條件了。 因此KKT condition就是：

最後的幾個條件實際上是lagrange乘子法的條件。

---

---

回到正題，既然咱們知道了能夠利用KKT condition把origin target function轉換到daul problem 來求解，那麼上面這個問題咱們嘗試用KKT條件求解一下： 首先對w，b求偏導：

把結果代回到dual problem：

因此最後咱們的target function就變成了這樣。 最後咱們能夠用QP對這個問題進行求解，求出了α以後，咱們隨便取一個α是非0的，也就是＞0的解，這時候利用α*g(x) = 0的條件獲得

b就求出來了，對於w直接代換上面的公式就行了。 而當α>0，由上面的公式能夠獲得這個點就剛恰好是在邊界上，而這些點就叫作support vector，支持向量的點。咱們的擬合直線也將會由着些點肯定，其餘不是support vector的點α就是0。

---

---

又停頓一下，咱們對這個式子思考一下：

####⑷爲何咱們須要dual problem 其實這個最優問題用普通的QP求解也是能夠的，可是若是咱們的數據維度很大，而通過feature transform以後維度就好更大了，這樣就會致使VC dimension會增大，特別是後面用的多項式核，RBF核等等。通過對偶問題KKT條件的變換以後，咱們的目標式子：

轉換成對偶問題以後，變量個數是N個，約束條件也是N個，於VC dimension就沒有了關係，從某種意義上是簡化了計算複雜度。其實計算複雜度仍是沒有變，只是把維度的計算提高到了變量之間點的內積罷了。將原始SVM轉化爲對偶問題，本意是在非線性變化，進行特徵轉換後，若是d’很大，爲了簡化計算，消除d’的影響。進一步引入Kernel SVM，根本上解決上述問題。注意了，這裏只是從某個角度看確實是消除了d維度的影響，實際上並無消失，只是轉移到了計算內積裏面而已。

---

---

回到正題，咱們的target function：

至於這個α怎麼求，後面會用SMO算法求解。 到這裏linear SVM就算結束了。

這就是分類函數。

---

---

再停頓一下，**什麼是支持向量點，爲何非支持向量的點α = 0？**這裏僅僅思考linear SVM，若是是soft margin又不同了。

####⑸支持向量

若是是支持向量，他的function margin是1；而對於很多支持向量的點，function margin > 1，因此右邊是負數，爲了知足最大，因此α只能爲0了，因此非支持向量的點α就是0。

---

---

##⑤kernel Support Vector Machine 回到正題，剛剛只是講了linear SVM，是對於linear separable有效而已，若是是linear inseparable呢？好比一個圓形，這樣就玩不了。記得以前linear regression和logistics regression講到過一個feature transform，若是是非線性的咱們能夠映射到其餘維度進行解決，好比最多見的polynomial transform，可是這樣問題來了，剛剛不是才把維度d轉移到內積嗎？（用dual problem的KKT condition）在來個feature transform那就是φ(x0)φ(x1)了，維度就更大了。

好比polynomial：

二項式的是這樣的，注意到中間好像多了一個X1Xd，這是爲了後面計算方便而已。兩個作內積：

能夠看到，最後的轉換就只和原始空間有關係而已，對於轉換隻後的z空間的維度沒有關係。好比x空間是2維的，爲了解決nonlinear problem，咱們映射到了z空間，在z空間裏面維度確定會比在x空間的原始維度要大，而最後用z空間作內積咱們就只須要拿x空間的原始維度就行了，由於咱們能夠先內積再升維，而不是先升維再內積。

這種就叫作核函數了。 最後的分類函數用kernel function替代：

剛剛所講的就是核函數的一種—— polynomial kernel function

加上幾個參數，γ就是它的參數了，最後化簡一下：

雖然都是二次轉換，對應到同一個z空間。可是，若是他們的γ係數不一樣，內積就會不同，那麼就表明有不一樣的距離，最終可能會獲得不一樣的SVM margin。因此，係數不一樣，可能會獲得不一樣的hperplane。 看一下γ係數對於hperplane的影響：

使用高階的polynomial kernel的話，獲得的Support Vector數量不會太多，分類面不會太複雜，防止過擬合。一樣也避開了對升維以後維度的依賴。

接下來介紹另一種kernel function——Gaussion kernel function 剛剛介紹的Q階多項式是有限維度的，若是是無限維度的能不能經過kernel來簡化計算？？有一個無限維的kernel function——Gaussion kernel

這和咱們以前見的有些不一樣，只是去掉了下面的方差而已，方差是定值沒有什麼太大的影響。逆推看看它的維度是多少：

推出來後面的維度是無限個（中間用的是Taylor展開，由於e的特殊求導性質能夠簡化）。

分類函數就出來了。 可是核函數的過擬合仍是有一點嚴重的：

γ對於核函數的影響有點大。若是取值很大的話最後就會造成一個一個的小圈圈把那些點圈起來。

---

---

又得停頓一下，思考一下核函數的意義以及他們之間的對比：

####⑹Comparison of Kernels Polynomial Kernel的hyperplanes是由多項式曲線構成。優勢：階數能夠靈活設置，更貼近實際分佈；缺點：當Q很到的時候，若是kernel裏面的值算出來是<1，那就基本接近0了，大於1就會變得很大，增長計算複雜度。並且參數過多，難以選擇合適的值。

Gaussan Kernel的優勢是邊界更加複雜多樣，能最準確地區分數據樣本，數值計算K值波動較小，並且只有一個參數，容易選擇；缺點是因爲特徵轉換到無限維度中，w沒有求解出來，計算速度要低於linear kernel，並且可能會發生過擬合。mysterious——no w；slower；too powerful。

以前說過經過對偶問題，咱們的把數據維度轉移到了內積上，因此從某一方面來看咱們確實是作到了簡化計算複雜度，可是實際上內積仍是屬於一個很大的計算。 因此核函數的功能之一，就是簡化計算，把升維和計算內積合在了一塊兒，減小計算複雜度。把計算步驟結合在了一塊兒，以前是先映射再計算內積，如今是一塊兒作了。核函數的功能之二，就是能夠很好的計算兩個樣本點的類似性，即內積。既然是表明類似性，咱們可不可使用其餘的核函數呢？或者本身建立一個，好比歐氏距離，餘弦距離等等？答案是不行。 先來看一下kernel的矩陣：

這有點像以前的協方差矩陣，只是沒有減去均值，因此對稱半正定是基本性質了。因此天然，咱們本身建立或選擇的時候也要選擇 ①symmetric對稱②positive semi-definite 半正定。這也是核函數有效性的判斷。

---

---

回到正題，剛剛只是講了一下對核函數的理解。 ##⑥Soft-Margin Support Vector Machine 上面應用到的Gaussion Kernel貌似仍是會出現過擬合，並且仍是蠻嚴重的，這說明large margin已經限制不了Gaussion kernel了，咱們須要找其餘方法來處理這個問題。 以前有一個比較簡單的算法——perceptron learning algorithm 這個算法對於nonlinear problem有一個很好的處理方式，咱們不要求必定要分類正確，咱們只要求找到一個錯誤最少的分類就能夠了。因此他的function是這樣：

不正確的就加個1，最小爲止。SVM也能夠用這種方法來限制。

加上一個條件，C參數就是對於這些錯誤懲罰度是多少，條件也變了，正確的≥ 1，錯誤的無論他。無論是小錯仍是大錯。 整合一下：

這個式子其實沒有什麼用， 首先不是線性的，用不了二次規劃，更不要說對偶這些了，其次大錯小錯都是同等對待，connot distinguish small error and large error。 對於上述方案繼續修正： 咱們採用一個ξ做爲一個犯錯程度，程度越大，懲罰越大。懲罰就是這個式子數值會變大，而後SVM要花更多的力氣去處理。

接下來就是對偶問題的推導，和以前的hard其實差很少的，lagrange 乘子法加對偶條件：

一樣，KKT條件：

C - α = β 因此有：0 < α < C 其餘的基本一致：w求導爲0：

b求導：

接下來就是求b了：

求b的公式裏面有一個矛盾的地方，就是咱們要求b首先得要求出來ξ的值，可是ξ的值也只有b的公式能夠求的處理，因此這就有一個雞生蛋蛋生雞的問題。因此咱們口語去掉這個ξ。咱們剛剛用到的是拉格朗日乘子法，後面的β(-ξ)是一個仿射函數，仿射函數有β(-ξ) = 0的性質，因此把β代換一下就獲得了上圖的公式。那麼去掉ξ就是等於0了，那麼就只有C-α不等於0纔有啊，因此當這個α ∈ （0 ， C）的時候就有ξ爲0，然後面咱們會講到當α∈（0，C）的時候這個點實際上是支持向量的點。這樣就能夠求出了b。 接下來看看C取值：

直接從我之前在CSDN裏面寫過的拷貝過來了。

接下來看一下一個比較重要的東西： physical significance of α

爲何βξ = 0？緣由和前一個公式是同樣的，由於要取最大值，因此這裏要等於0，β ≥ 0，而實際公式是negative ξ，因此乘上去要是0纔能有最大；第二，若是不是等於0就不等因而原問題的求解了，不等於0就無故端多了一個inequality，和原問題不對等了。以後才能進行daul problem的轉換。 咱們主要是從上面這兩個公式來看當α取值不一樣的時候對應的物理意義。

當α = 0，得ξ = 0，這個點就是沒有放錯的點，由於ξ = 0，不須要容忍。而α = 0，因此不是支持向量機的點，因此表明的就是在bound外而且分類正確的點。

當α∈（0，C），仍是獲得ξ = 0，這時候就不同了，尚未錯誤的點，可是第一條式子括號裏面等於0了，意味着就是在bound上的點，那麼就是支持向量點了。

當α = C，不能肯定ξ是否是0了，

，表示就是錯了多少，這種有兩種狀況，一種是分類正確了，可是距離太近；或者是分類錯了。 當α > C，不存在的，上面都限制了。

理一下整個思路。 ①找到最好的hperplane，最寬的那個。 ②獲得target function ③發現feature transform以後維度對於計算機複雜度有很大影響，用dual problem轉移到內積處理 ④轉移以後發現仍是複雜度在的，引出了kernel function ⑤發現kernel function仍是有overfitting的狀況，因而又引入了soft margin

---

---

在講SMO算法以前，先講一下對於error function的理解：

####⑺對於SVM error function的理解 咱們把SVM換一種形式。對於ξ，其實他是每個點距離邊界有多遠，一種是violating margin，即不知足y(wTz + b) ≥ 1，那麼ξ就能夠表示成：1 - y(wTz + b) > 0。第二種狀況就是not violating margin，即這個點在邊界以外，就是知足上述公式了，這個時候ξ就是0，咱們整合一下： ξ = max ( 1 - y(wTz + b) , 0 )，代換進原來的支持向量機公式：

這個就是支持向量機的error function，先預判了Ein = 0，也就是全對的狀況，前面有說到。 這個function有點像咱們以前所學的L2 lost function：

這和logistics regression的L2範式的cost function很類似。

其實就差很少是同樣的，沒有什麼差異，可是既然是相同的爲何不用這種方法呢？兩個緣由，一個是這種無條件的最優化問題沒法經過QP解決，即對偶推導和kernel都沒法使用；另外一個是這種形式中包含的max()項可能形成函數並非到處可導，這種狀況難以用微分方法解決。

對比發現，L2 regularization和soft margin SVM形式是同樣的，兩個式子λ和C是互相對應的。soft marginSVM裏面的large margin就對應着L2 regularization裏面的short w，都是讓hypothesis set能夠簡單點。λ和C也是互相對應，λ大，w就小，正則化的程度就越大；C小，Ein就大，響應這個margin也會打，因此增大C和減少λ是一個意思，因此large margin等同於regularization，都是防止過擬合做用的。

若是是按照咱們以前的err0/1，正確爲1，錯誤就是0，那麼有：

能夠看到SVM他是大於err0/1的，根據VC bound理論是能夠用來代替err0/1分類的。 後面再加上logic function的cost function：

而這個幾乎就是和L2-regularized logistic regression同樣的。Logistic Regression和Soft-Margin SVM都是在最佳化err0/1的上界而已。能夠看出，求解regularized logistic regression的問題等同於求解soft-margin SVM的問題。

####⑻損失函數 常見的損失函數： err0/1：

此時soft margin就是這樣了，大於0就是1小於就是0。 不敏感損失函數 —— hinge lost function

還有對數損失函數交叉熵等等。logistics用的是交叉熵，SVM就是用的hinge lost function。支持向量機就是一個結構風險最小化的近似實現，結構風險至關於指望風險(Eout)的一個上界，它是經驗風險（Ein）和置信區間(Ω模型複雜度)的和，經驗風險依賴於決策函數f的選取，可是置信區間是，F的VC維的增函數，二者是矛盾的。矛盾體如今：當VC維數變大的時候能夠選到更好的f使得經驗風險比較小，可是此時的置信區間比較大。這就是對應了VC bound理論。還好去聽了臺灣大學林軒宇老師課程，對這些機器學習理論基礎有了解。

---

---

回到正題，開始SMO算法。 ##⑦SMO算法 target function：

剛剛咱們知道怎麼求w,b，可是那是在知道了α的前提下，如今就來求α。 基本思路： 選擇兩個變量，固定其餘變量，針對兩個變量構建一個二次規劃問題。每次針對兩個變量來求解目標函數的最小值，求解完後，繼續尋找新的變量求目標函數，在每次尋找新α的過程當中，目標函數將進一步獲得優化，直到全部的αi更新完了。而對於α的選取，一個是違反KKT條件最嚴重的那一個，另外一個由約束條件自動肯定。

首先，假設咱們選取了兩個變量α1，α2，固定其餘變量以後：

因此只要求出α2，α1就知道了。

原目標函數化簡以後：

K11指的就是x1和本身自己作核函數。因爲咱們已經固定了除了α1和α2，因此天然其餘的常量咱們能夠去掉了，不如優化w+1，和優化w是同樣的，去掉固定常數項就留下了上圖的公式。

別忘了條件，條件是後面求解的關鍵。 首先咱們要獲得α1，α2的範圍。

由着兩個約束條件限制。

因此有：

因此當咱們更新了α以後，咱們還要根據範圍剪輯α才能夠。

咱們假設：

剪輯範圍：

再假設一個定值，也就是i = 3開始求和的：

目標式子：

用上面的vi代換以後：

求α2的話天然是求導了：

爲0獲得：

代入獲得：

這裏的化簡有點麻煩：

手動證實一下。 用假設替換一下上面的式子：

就能夠了。 SMO算法有兩個要點：①α1的選擇，違反KKT最嚴重的條件②α2的選擇策略

---

---

很重要的問題，變量要怎麼選擇，後面會有例子證實。

####⑼變量的選擇方式 SMO稱選擇第1個變量的過程爲外層循環。外層循環在訓練樣本中選取違反KKT條件最嚴重的樣本點，Violation of the most serious sample of KKT conditions。我第一次看這東西是懵逼的。可是仔細想一下，就是檢測哪個樣本是沒有知足KKT的條件：

首先遍歷全部0 < α < C的樣本點，看看是否是知足的，若是沒有載變量全部的。檢測是否知足KKT。因此在SMO迭代的兩個步驟中，只要α中有一個違背了KKT條件，這一輪迭代完成後，目標函數的值必然會增大。Generally speaking，KKT條件違背的程度越大，迭代後的優化效果越明顯，增幅越大。 α1選完了天然就是選擇第二個α了，第二個變量的選擇叫作內存循環，咱們這裏先用普通隨機選擇，看看效果如何。

---

---

##⑧算法實現——version 1 首先是導入各類各樣的包和一個工具了：

import numpy as np  
import matplotlib.pyplot as plt  
import random  
import seaborn as sea  
import pandas as pd  


def get_positive_and_negative():  
  dataSet = pd.read_csv('Datas/LogiReg_data.txt', names=['V1', 'V2', 'Class'])  
  dataSet.Class[dataSet.Class == 0] = -1  
  dataSet = dataSet[60 : 80]  
  positive = dataSet[dataSet['Class'] == 1]  
  negative = dataSet[dataSet['Class'] == -1]  
  return positive , negative , dataSet  


def show_picture(positive , negative):  
  columns = ['V1', 'V2']  
  fig, ax = plt.subplots(figsize=(10, 5))  
  ax.scatter(positive[columns[0]], positive[columns[1]], s=30, c="b", marker="o", label="class 1")  
  ax.scatter(negative[columns[0]], negative[columns[1]], s=30, c="r", marker="x", label="class -1")  
  ax.legend()  
  ax.set_xlabel('V1')  
  ax.set_ylabel('V3')  
  plt.show()  

def load_data_set():  
  _ , _ , file = get_positive_and_negative()  
  orig_data = file.as_matrix()  
  cols = orig_data.shape[1]  
  data_mat = orig_data[ : , 0 : cols-1]  
  label_mat = orig_data[ : , cols-1 : cols]  
  return  data_mat , label_mat  

positive , negative , data = get_positive_and_negative()  
show_picture(positive , negative)  
print(data)  
複製代碼

第一個是獲得正負樣本，而後顯示，最後一個是加載數據，數據隨便找一個就行了。

positive , negative , data = get_positive_and_negative()  
show_picture(positive , negative)  
複製代碼

最後調用一些看看這些點是什麼：

還有一些是對α的限制和一下工具函數：

''''' Generate a random number '''  
def select_jrand(i , m):  
  j = i  
  while(j == i):  
      j = int(random.uniform(0 , m))  
  return j  
  pass  

''''' restraint the α '''  
def clip_alpha(aj , H , L):  
  if aj > H:  
      aj = H  
  elif aj < L:  
      aj = L  
  return aj  
  pass  
複製代碼

接下來就是實現支持向量機了：

def SVM(data_mat , class_label , C , tolar , max_iter):  

  data_mat = np.mat(data_mat)  
  label_mat = np.mat(class_label)  
  b = 0  
  m , n = np.shape(data_mat)  
  alphas = np.zeros((m , 1))  
  iter = 0  

  while iter < max_iter:  
      #做爲迭代變化量 
      alpha_pairs_changed = 0  
      #做爲第一個a 
      for i in range(m):  
          WT_i = np.dot(np.multiply(alphas , label_mat).T , data_mat)  
          f_xi = float(np.dot(WT_i , data_mat[i , :].T)) + b  
          Ei = f_xi - float(label_mat[i])  
          if ((label_mat[i]*Ei < -tolar) and (alphas[i] < C)) or ((label_mat[i]*Ei > tolar) and (alphas[i] > 0)):  
              j = Tools.select_jrand(i , m)  
              WT_j = np.dot(np.multiply(alphas , label_mat).T , data_mat)  
              f_xj  = float(np.dot(WT_j , data_mat[j , :].T)) + b  
              Ej = f_xj - float(label_mat[j])  
              alpha_iold = alphas[i].copy()  
              alpha_jold = alphas[j].copy()  

              if (label_mat[i] != label_mat[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 H == L:  
                  continue  

              eta = 2.0 * data_mat[i, :] * data_mat[j, :].T - data_mat[i, :] * data_mat[i, :].T - data_mat[j, :] * data_mat[j, :].T  
              if eta >= 0: continue  
              alphas[j] = (alphas[j] - label_mat[j]*(Ei - Ej))/eta  
              alphas[j] = Tools.clip_alpha(alphas[j], H, L)  
              if (abs(alphas[j] - alpha_jold) < 0.00001):  
                  continue  
              alphas[i] = alphas[i] + label_mat[j]*label_mat[i]*(alpha_jold - alphas[j])  


              b1 = b - Ei + label_mat[i]*(alpha_iold - alphas[i])*np.dot(data_mat[i,:], data_mat[i,:].T) +\  
              label_mat[j]*(alpha_jold - alphas[j])*np.dot(data_mat[i,:], data_mat[j,:].T)  
              b2 = b - Ej + label_mat[i]*(alpha_iold - alphas[i])*np.dot(data_mat[i,:], data_mat[j,:].T) +\  
              label_mat[j]*(alpha_jold - alphas[j])*np.dot(data_mat[j,:], data_mat[j,:].T)  
              if (0 < alphas[i]) and (C > alphas[i]):  
                  b = b1  
              elif (0 < alphas[j]) and (C > alphas[j]):  
                  b = b2  
              else:  
                  b = (b1 + b2)/2.0  
              print(b)  
              alpha_pairs_changed += 1  
              pass  
      if alpha_pairs_changed == 0:  
          iter += 1  
      else:  
          iter = 0  

  support_x = []  
  support_y = []  
  class1_x = []  
  class1_y = []  
  class01_x = []  
  class01_y = []  
  for i in range(m):  
      if alphas[i] > 0.0:  
          support_x.append(data_mat[i, 0])  
          support_y.append(data_mat[i, 1])  
  for i in range(m):  
      if label_mat[i] == 1:  
          class1_x.append(data_mat[i, 0])  
          class1_y.append(data_mat[i, 1])  
      else:  
          class01_x.append(data_mat[i, 0])  
          class01_y.append(data_mat[i, 1])  
  w_best = np.dot(np.multiply(alphas, label_mat).T, data_mat)  
  fig, ax = plt.subplots(figsize=(10, 5))  
  ax.scatter(support_x, support_y, s=100, c="y", marker="v", label="support_v")  
  ax.scatter(class1_x, class1_y, s=30, c="b", marker="o", label="class 1")  
  ax.scatter(class01_x, class01_y, s=30, c="r", marker="x", label="class -1")  
  lin_x = np.linspace(0, 100)  
  lin_y = (-float(b) - w_best[0, 0] * lin_x) / w_best[0, 1]  
  plt.plot(lin_x, lin_y, color="black")  
  ax.legend()  
  ax.set_xlabel("factor1")  
  ax.set_ylabel("factor2")  
  plt.show()  
  return b , alphas  
datamat , labelmat = dataSet.load_data_set()  
b, alphas = SVM(datamat , labelmat , 0.6 , 0.001 , 10)  
print(b , alphas)  
複製代碼

首先傳入的後面幾個參數分別是懲罰力度，容忍度。比較重要的應該是這一句：

if ((label_mat[i]*Ei < -tolar) and (alphas[i] < C)) or ((label_mat[i]*Ei > tolar) and (alphas[i] > 0)):  
複製代碼

這句話翻譯過去就是yg(x) < 1 - ξ或者是y(g(x)) > 1+ξ。若是是小於，則這個點是離hperplane比較近，這時候這個點應該是等於C纔對的；若是是大於了，也就是遠大於邊界了，那就是離邊界很遠了，可是α又大於0，離邊界遠意味着不是支持向量，因此α應該是0，因此能夠改變。 後面的那些就是依據公式來的：

每一條都是對應公式寫的。 最後就是打印了。

效果：

能夠看到是極度不穩定。這是幾個月前我實現的，後來如今我又從新實現了一個，用了一些改進方法。 爲何會不穩定，我總結了幾個緣由： ①沒有緩存，更新慢，迭代次數不夠 ②對於α2的選取沒有很好的採起策略 ③對於n，也就是更新公式：

我沒有判斷是否是大於0的。n是什麼東西呢？

他要是小於0意味着這個kernel matrix就不是半正定的了，K11 + K22 < 2K12；另外，這個n實際上是：

的二階導數，小於0就不是凸函數了，哪來的凸優化。因此應該是更新的時候遇到這些狀況致使不穩定的。 基於上面的缺點更換策略。 ##⑨算法實現——version 2 首先要改變的是加上一個緩存，用來保存Ei的值，使得計算更塊。**其次就是α2的選擇策略，在優化過程當中，會經過最大化步長的方式來得到第二個alpha值。第二步優化爲，數據集全程掃描策略與在非邊界alpha對中進行更新策略交替進行。**對於n，會進行判斷是否是大於0，在這裏是用-號的，因此n與咱們表達式上的是想反方向，因此是大於0。 首先仍是工具：

''' load data and define some tool function '''
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import random

def loadDataSet(filename):
  ''' :param filename: :return dataset and label: '''

  dataset = []
  label = []
  fr = open(filename)
  for line in fr.readlines():
      lineArr = line.strip().split('\t')
      dataset.append( [np.float32(lineArr[0]) , np.float32(lineArr[1])] )
      label.append(np.float32(lineArr[2]))
  return dataset , label
  pass

''' select alpha2 randomly '''
def selectAlphaTwo(i , m):
  ''' :param i: :param m: :return: '''
  j = i
  while(j == i):
      j = int(random.uniform(0 , m))
  return j

def rangeSelectionForAlpha(aj , H , L):
  if aj > H:
      aj = H
  if L > aj:
      aj = L
  return aj
  pass

''' calculate Ei '''
def calculateEi(os , k):
  fxk = float(np.multiply(os.alphas, os.labels).T * (os.x * os.x[k, :].T)) + os.b
  Ek = fxk - float(os.labels[k])
  return Ek

''' put the Ei into the cache when calculate Ei '''
def selectj(i , os , Ei):
  maxk = -1
  maxDeltaE = 0
  Ej = 0
  os.eCache[i] = [1 , Ei]
  validEachlist = np.nonzero(os.eCache[: , 0].A)[0]
  if (len(validEachlist) > 1):
      for k in validEachlist:
          if k == i:
              continue
          Ek = calculateEi(os , k)
          deltaE = np.abs(Ei - Ek)
          if deltaE > maxDeltaE:
              maxk = k
              maxDeltaE = deltaE
              Ej = Ek
      return maxk , Ej
      pass
  else:
      j = selectAlphaTwo(i , os.m)
      Ej = calculateEi(os , j)
  return j , Ej
  pass

''' draw picture '''
def drawDataset(data , label , x = None , y = None , line = True , alphas = None , kernel = True):
  index_one = []
  index_negative_one = []
  for i in range(100):
      if label[i] == 1:
          index_one.append(data[i])
      else:
          index_negative_one.append(data[i])
  index_one = np.matrix(index_one)
  index_negative_one = np.matrix(index_negative_one)
  plt.scatter(index_one[ : , 0].tolist() , index_one[: , 1].tolist() , c = 'r' , marker='<' , label = 'class equal one')
  plt.scatter(index_negative_one[: , 0].tolist() , index_negative_one[: , 1].tolist() , c = 'b' , marker='x' , label = 'class equal negative one')
  if line == True:
      plt.plot(x , y)
      pass

  ''' draw the support vector,the point which the α not equal zero '''
  if line == True or kernel == True:
      a1 = []
      for i in range(len(alphas)):
          a = alphas[i]
          if a != 0:
             a1.append(data[i])
      a1 =  np.matrix(a1)
      print('The number of the support vector : ' , len(a1))
      plt.scatter(a1[: , 0].tolist(),a1[: , 1].tolist(), s=150, c='none', alpha=0.7,
                     linewidth=1.5, edgecolor='#AB3319' , label = 'support vector')

  plt.legend()
  plt.xlabel('X axis')
  plt.ylabel('Y axis')
  plt.show()

def updateEk(os,k):
  Ek = calculateEi(os,k)
  os.eCache[k]=[1,Ek]

if __name__ == '__main__':
  data , label = loadDataSet('../Data/testSetRBF.txt')
  drawDataset(data , label , line=False ,kernel=False)


複製代碼

SMO算法惟一的一個類：

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import random
import KernelTransform
class optStruct:
  def __init__(self , dataMat , labels , C , toler):
      self.x = dataMat
      self.labels = labels
      self.C = C
      self.toler = toler
      self.m = np.shape(dataMat)[0]
      self.alphas = np.mat(np.zeros((self.m , 1)))
      self.b = 0
      self.eCache = np.mat(np.zeros((self.m , 2)))
      self.K = np.mat(np.zeros((self.m , self.m)))
      for i in range(self.m):
          self.K[: , i] = KernelTransform.kernelTrans(self.x , self.x[i , :] , kTup=('rbf' , 1.2))
      pass

if __name__ == '__main__':
  os = optStruct([1,2] , [3,4] , 1,1)
  a = os.alphas.tolist()[0][0] -  os.alphas.tolist()[1][0]
  print(max(1.0 , a))


複製代碼

須要解釋的應該只有selectj()了，這個是經過計算最大不長來選擇α2的。 首先咱們假設最大不長是-1，由於相減有絕對值不多是negative；os.eCache是咱們的緩存的Ei，先把Ei存進去，1,表示這個數字不是0，這一步就是獲得這個緩存裏面全部有效（不爲0）的Ei。判斷獲得的列表是否是有東西，沒有就隨機選擇了。仍是再解釋一下爲何要這個創建表格吧！ 咱們在選擇第一個α1的時候，選擇的是在邊界外的點，也就是非邊界的點。 優先選擇遍歷非邊界數據樣本，由於非邊界數據樣本更有可能須要調整，邊界數據樣本經常不能獲得進一步調整而留在邊界上。因爲大部分數據樣本都很明顯不多是支持向量，所以對應的α乘子一旦取得零值就無需再調整。遍歷非邊界數據樣本並選出他們當中違反KKT 條件爲止。當某一次遍歷發現沒有非邊界數據樣本獲得調整時，遍歷全部數據樣本，以檢驗是否整個集合都知足KKT條件。若是整個集合的檢驗中又有數據樣本被進一步進化，則有必要再遍歷非邊界數據樣本。這樣，不停地在遍歷全部數據樣本和遍歷非邊界數據樣本之間切換，直到整個樣本集合都知足KKT條件爲止。以上用KKT條件對數據樣本所作的檢驗都以達到必定精度ε就能夠中止爲條件。若是要求十分精確的輸出算法，則每每不能很快收斂。因此在echa中緩存的第一次選出的α，由於咱們選出來的就是非邊界上的點，α2選擇的時候繼續在上面遍歷，雖然緩存是存了Ei，可是這個Ei不能直接用，由於那個是舊的值。因此α的迭代策略就是非邊界和全局選取兩種交替進行了。

以後就是正式的算法了：

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import random
import Tool
import smo_class
import KernelTransform
def innerL(i ,os):
  Ei = Tool.calculateEi(os , i)
  if ((os.labels[i]*Ei < -os.toler) and
      (os.alphas[i] < os.C)) or ((os.labels[i]*Ei > os.toler) and
                                 (os.alphas[i] > 0)):
      j , Ej = Tool.selectj(i , os , Ei)
      alphaIold = os.alphas[i].copy()
      alphaJold = os.alphas[j].copy()
      if (os.labels[i] != os.labels[j]):
          L = max(0 , os.alphas[j] - os.alphas[i])
          H = min(os.C , os.C + np.array(os.alphas)[j] - np.array(os.alphas)[i])
      else:
          L = max(0 , os.alphas[j] + os.alphas[i] - os.C)
          H = min(os.C , np.array(os.alphas)[j] + np.array(os.alphas)[i])
      if L == H:
          return 0
      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('η> 0，the kernel matrix is not semi-positive definite')
          return 0
      os.alphas[j] -= os.labels[j]*(Ei - Ej)/eta
      os.alphas[j] = Tool.rangeSelectionForAlpha(os.alphas[j] , H , L)
      Tool.updateEk(os , j)

      if (abs(os.alphas[j] - alphaJold) < 0.00001):
          print("j not moving enough")
          return 0
      os.alphas[i] += os.labels[j] * os.labels[i] * (alphaJold - os.alphas[j])
      Tool.updateEk(os , i)
      b1 = os.b - Ei - os.labels[i] * (os.alphas[i] - alphaIold) * \
           os.x[i, :] * os.x[i, :].T - os.labels[j] * \
           (os.alphas[j] - alphaJold) * os.x[i, :] * os.x[j, :].T
      b2 = os.b - Ej - os.labels[i] * (os.alphas[i] - alphaIold) * \
           os.x[i, :] * os.x[j, :].T - os.labels[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
      return 1
  else:
      return 0

def smo(data,labels,C = 0.6,toler = 0.001,maxIter = 40 , kernel = True):
  oS = smo_class.optStruct(np.mat(data),np.mat(labels).transpose(),C,toler)
  iter =0
  entireSet  = True
  alphaPairsChanged = 0
  while(iter < maxIter) and ((alphaPairsChanged >0) or (entireSet)):
      alphaPairsChanged = 0
      if entireSet:
          for i in range(oS.m):
              if kernel == True:
                  alphaPairsChanged += KernelTransform.innerL(i,oS)
              else:
                  alphaPairsChanged += innerL(i, oS)
          print("fullSet,iter: %d i: %d,pairs changed %d" %\
              (iter,i,alphaPairsChanged))
          iter +=1
      else:
          # 兩個元素乘積非零，每兩個元素作乘法[0,1,1,0,0]*[1,1,0,1,0]=[0,1,0,0,0]
          nonBoundIs = np.nonzero((oS.alphas.A > 0)*(oS.alphas.A < C))[0]
          for i in nonBoundIs:
              alphaPairsChanged += innerL(i,oS)
              print("nou-bound,iter: %d i:%d,pairs changed %d" % (iter,i,alphaPairsChanged))
          iter +=1
      # entireSet 控制交替的策略選擇
      if entireSet:
          entireSet = False
      # 必須有alpha對進行更新
      elif(alphaPairsChanged == 0):
          entireSet = True
      print("iteration number：%d" % iter)
  return oS.b,oS.alphas
複製代碼

entireSet就是交換策略的標誌。貌似沒有什麼好說的。 以後就是執行函數這些了：

import Tool
import SMO
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import KernelTransform
''' calculate w and draw the picture, the variable which the α not equal zero , we call support vector '''
def calculateW(alphas , data , labels):
  x = np.mat(data)
  label = np.mat(labels).transpose()
  m , n = np.shape(x)
  w = np.zeros((n , 1))
  for i in range(m):
      w += np.multiply(alphas[i] * label[i] , x[i , :].T)
  return w
  pass

if __name__ == '__main__':
  data, label = Tool.loadDataSet('../Data/testSet.txt')
  b,alphas = SMO.smo(data , label , kernel=False)
  w = calculateW(alphas , data , label)
  x = np.arange(0 , 11)
  print(w)
  y = (-b - w[0]*x)/w[1]
  Tool.drawDataset(data , label , x , y.tolist()[0] , line=True , alphas=alphas)

  data, label = Tool.loadDataSet('../Data/testSetRBF.txt')
  b, alphas = SMO.smo(data, label,kernel=True ,maxIter=100)
  svInd = np.nonzero(alphas.A > 0)[0]
  Tool.drawDataset(data, label,  line=False, alphas=alphas)







複製代碼

有一個是kernel function的，先不用管。 效果：

圈起來的是支持向量點，好不少了。 ##⑩算法實現——version 3 kernel function加上，先看看原來的數據：

須要改的其實就是內積就能夠了，處處看看哪裏有內積就改改他，修改事後的innel和smo：

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import Tool
def kernelTrans(X,A,kTup):
  m,n = np.shape(X)
  K = np.mat(np.zeros((m,1)))
  if kTup[0]=='lin':
      K = X*A.T
  elif kTup[0] =='rbf':
      for j in range(m):
          deltRow = X[j,:]-A
          K[j] = deltRow*deltRow.T
      K = np.exp(K/(-1*kTup[1]**2))
  return K

''' update the innel function '''
def innerL(i ,os):
  Ei = calculateEi(os , i)
  if ((os.labels[i]*Ei < -os.toler) and
      (os.alphas[i] < os.C)) or ((os.labels[i]*Ei > os.toler) and
                                 (os.alphas[i] > 0)):
      j , Ej = Tool.selectj(i , os , Ei)
      alphaIold = os.alphas[i].copy()
      alphaJold = os.alphas[j].copy()
      if (os.labels[i] != os.labels[j]):
          L = max(0 , os.alphas[j] - os.alphas[i])
          H = min(os.C , os.C + np.array(os.alphas)[j] - np.array(os.alphas)[i])
      else:
          L = max(0 , os.alphas[j] + os.alphas[i] - os.C)
          H = min(os.C , np.array(os.alphas)[j] + np.array(os.alphas)[i])
      if L == H:
          return 0
      eta = 2.0 * os.K[i, j] - os.K[i, i] - os.K[j, j]
      if eta >= 0:
          print('η> 0，the kernel matrix is not semi-positive definite')
          return 0
      os.alphas[j] -= os.labels[j]*(Ei - Ej)/eta
      os.alphas[j] = Tool.rangeSelectionForAlpha(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.labels[j] * os.labels[i] * (alphaJold - os.alphas[j])
      updateEk(os , i)
      b1 = os.b - Ei - os.labels[i] * (os.alphas[i] - alphaIold) * \
           os.K[i , i] - os.labels[j] * \
           (os.alphas[j] - alphaJold) *  os.K[i , j]
      b2 = os.b - Ej - os.labels[i] * (os.alphas[i] - alphaIold) * \
           os.K[i , j] - os.labels[j] * \
           (os.alphas[j] - alphaJold) * os.K[j , j]
      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
      return 1
  else:
      return 0

''' updata the Ei '''
def calculateEi(os , k):
  fxk = float(np.multiply(os.alphas, os.labels).T * os.K[:, k] + os.b)
  Ek = fxk - float(os.labels[k])
  return Ek
def updateEk(os,k):
  Ek = calculateEi(os,k)
  os.eCache[k]=[1,Ek]
複製代碼

剛剛那個執行函數其實已經包括了kernel的，因此直接就能夠看到效果了：

用的是Gaussion kernel，不知道怎麼作擬合，就把支持向量點圈出來就行了。 最後附上全部代碼GitHub： github.com/GreenArrow2…