支援向量機（Support Vector Machine）

支援向量機

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

GV0SHYC3S{P{Q4QVB66UN6T.png

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

0T0FN3C5MKTQE`DO{N}QX(C.png

①超平面

介紹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：

https://github.com/GreenArrow2017/MachineLearning