Python機器學習(十二)支持向量機算法


1. 解決什么問題?

最基本的應用是數據分類,特別是對於非線性不可分數據集。支持向量機不僅能對非線性可分數據集進行分類,對於非線性不可分數據集的也可以分類

(我認為這才是支持向量機的真正魅力所在,因為現實場景中,樣本數據往往是非線性不可分的)。

現實場景一 :樣本數據大部分是線性可分的,但是只是在樣本中含有少量噪聲或特異點,去掉這些噪聲或特異點后線性可分 => 用支持向量機的軟間隔方法進行分類;

現實場景二:樣本數據完全線性不可分 => 引入核函數,將低維不可分的非線性數據集轉化為高維可分的數據集,用支持向量機的軟間隔方法進行分類;

 

 

一定要明白一點:對分類器准確性有影響只是樣本中的支持向量,因此,其他樣本在計算分類器的權重矩陣時可以直接過濾掉,大大節省運行時間。

目標:獲取離邊界最近樣本點到超平面最遠距離

 

 

2. SVM的思想演化?

博客 https://www.cnblogs.com/zhizhan/p/4430253.html 關於SVM的演化說得很透徹,也很形象,這里借用一下。

2.1 硬間隔支持向量機

SVM中最關鍵的思想之一就是引入和定義了“間隔”這個概念。這個概念本身很簡單,以二維空間為例,就是點到分類直線之間的距離。假設直線為y=wx+b,那么只要使所有正分類點到該直線的距離與所有負分類點到該直線的距離的總和達到最大,這條直線就是最優分類直線。這樣,原問題就轉化為一個約束優化問題,可以直接求解。這叫做硬間隔最大化,得到的SVM模型稱作硬間隔支持向量機

2.2 軟間隔支持向量機

但是新問題出現了,在實際應用中,我們得到的數據並不總是完美的線性可分的,其中可能會有個別噪聲點,他們錯誤的被分類到了其他類中。如果將這些特異的噪點去除后,可以很容易的線性可分。但是,我們對於數據集中哪些是噪聲點卻是不知道的,如果以之前的方法進行求解,會無法進行線性分開。是不是就沒辦法了呢?假設在y=x+1直線上下分為兩類,若兩類中各有對方的幾個噪點,在人的眼中,仍然是可以將兩類分開的。這是因為在人腦中是可以容忍一定的誤差的,仍然使用y=x+1直線分類,可以在最小誤差的情況下進行最優的分類。同樣的道理,我們在SVM中引入誤差的概念,將其稱作“松弛變量”。通過加入松弛變量,在原距離函數中需要加入新的松弛變量帶來的誤差,這樣,最終的優化目標函數變成了兩個部分組成:距離函數和松弛變量誤差。這兩個部分的重要程度並不是相等的,而是需要依據具體問題而定的,因此,我們加入權重參數C,將其與目標函數中的松弛變量誤差相乘,這樣,就可以通過調整C來對二者的系數進行調和。如果我們能夠容忍噪聲,那就把C調小,讓他的權重降下來,從而變得不重要;反之,我們需要很嚴格的噪聲小的模型,則將C調大一點,權重提升上去,變得更加重要。通過對參數C的調整,可以對模型進行控制。這叫做軟間隔最大化,得到的SVM稱作軟間隔支持向量機

2.3 非線性支持向量機

之前的硬間隔支持向量機和軟間隔支持向量機都是解決線性可分數據集或近似線性可分數據集的問題的。但是如果噪點很多,甚至會造成數據變成了線性不可分的,那該怎么辦?最常見的例子是在二維平面笛卡爾坐標系下,以原點(0,0)為圓心,以1為半徑畫圓,則圓內的點和圓外的點在二維空間中是肯定無法線性分開的。但是,學過初中幾何就知道,對於圓圈內(含圓圈)的點:x^2+y^2≤1,圓圈外的則x^2+y^2>1。我們假設第三個維度:z=x^2+y^2,那么在第三維空間中,可以通過z是否大於1來判斷該點是否在圓內還是圓外。這樣,在二維空間中線性不可分的數據在第三維空間很容易的線性可分了。這就是非線性支持向量機

實際中,對某個實際問題函數來尋找一個合適的空間進行映射是非常困難的,幸運的是,在計算中發現,我們需要的只是兩個向量在新的映射空間中的內積結果,而映射函數到底是怎么樣的其實並不需要知道。這一點不太好理解,有人會問,既然不知道映射函數,那怎么能知道映射后在新空間中的內積結果呢?答案其實是可以的。這就需要引入了核函數的概念。核函數是這樣的一種函數:仍然以二維空間為例,假設對於變量x和y,將其映射到新空間的映射函數為φ,則在新空間中,二者分別對應φ(x)和φ(y),他們的內積則為<φ(x),φ(y)>。我們令函數Kernel(x,y)=<φ(x),φ(y)>=k(x,y),可以看出,函數Kernel(x,y)是一個關於x和y的函數!而與φ無關!這是一個多么好的性質!我們再也不用管φ具體是什么映射關系了,只需要最后計算Kernel(x,y)就可以得到他們在高維空間中的內積,這樣就可以直接帶入之前的支持向量機中計算!

核函數不是很好找到,一般是由數學家反向推導出來或拼湊出來的。現在知道的線性核函數有多項式核函數、高斯核函數等。其中,高斯核函數對應的支持向量機是高斯徑向基函數(RBF),是最常用的核函數。

 

RBF核函數可以將維度擴展到無窮維的空間,因此,理論上講可以滿足一切映射的需求。為什么會是無窮維呢?我以前都不太明白這一點。后來老師講到,RBF對應的是泰勒級數展開,在泰勒級數中,一個函數可以分解為無窮多個項的加和,其中,每一個項可以看做是對應的一個維度,這樣,原函數就可以看做是映射到了無窮維的空間中。這樣,在實際應用中,RBF是相對最好的一個選擇。當然,如果有研究的話,還可以選用其他核函數,可能會在某些問題上表現更好。但是,RBF是在對問題不了解的情況下,對最廣泛問題效果都很不錯的核函數。因此,使用范圍也最廣。

3. SVM原理

3.1 目標函數

求使幾何間隔最大的分離超平面

                

以及相應的分類決策函數

                 

也就是求出w,b

3.2 幾何間隔定義

3.3 學習的對偶算法

為了求解線性可分支持向量機的最優化問題,將它作為原始最優化問題,應用拉格朗日對偶性,通過求解對偶問題得到原始問題最優解。

3.4 最大軟間隔分類器(線性可分支持向量學習的最優化問題)

這就是線性可分支持向量機的對偶算法,這樣做的優點,一是對偶問題往往更容易求解;二是自然引入核函數,進而推廣到非線性分類問題。         

 

 的求解就轉化為  的求解。上述公式就是對不同變量求偏導,具體推導過程見《統計學習方法》第7章 或者相關博客。

3.5 KKT條件

 

 

上圖的總結有點小問題,分類不能合並,正確的分類描述應該是:

在超平面上或誤分的樣本點是不能正確分類的。

4. SMO算法(sequential minimal optimization 序列最小最優化)

高效實現支持向量機的算法是SMO算法,其基本思路是:如果所有變量的解都滿足此最優化的KKT條件,那么這個最優化問題的解就得到了。因為KKT條件是該優化問題的充分必要條件。否則,

選擇兩個變量,固定其他變量,針對這兩個變量構建一個二次規划問題,這個二次規划問題關於這兩個變量的解應該更接近原始二次規划問題的解,因為這會使原始二次規划問題的目標函數值變得更小。

重要的是,這時子問題可以通過解析方法求解,這樣可以大大提高整個算法的計算速度。

子問題有兩個變量,一個是違反KKT條件最嚴重的那一個,另一個由約束條件自動確定。如此,SMO算法將原問題不斷分解為子問題並對子問題求解,進而達到求解原問題的目的。

並且   推出

4.1  范圍

 假定不固定,其他拉格朗日因子固定(也就是常量),得到

  常量 , 為變量, 為因變量,L,H的范圍就是  的 范圍

 

4.2 計算拉格朗日乘子 

第一個  公式如下,推導過程見《統計學習方法》第7章或者相關博客。

     

     

選擇一個違背KKT的數據項,根據下面結論得到優化后 

 

s表示clip_image101,從而可以得到另一個優化后的 ,同樣也需要進行約束。

4.3 計算b1,b2

在每次完成兩個變量的優化后,都要重新計算閾值b。具體推導見《統計學習方法》

 

5. 代碼實現

5.1 輸入數據

這里使用兩個原始數據文件 trainingData.txt,testData.txt。

trainingData.txt

復制代碼
  1 -0.214824    0.662756    -1.000000
  2 -0.061569    -0.091875    1.000000
  3 0.406933    0.648055    -1.000000
  4 0.223650    0.130142    1.000000
  5 0.231317    0.766906    -1.000000
  6 -0.748800    -0.531637    -1.000000
  7 -0.557789    0.375797    -1.000000
  8 0.207123    -0.019463    1.000000
  9 0.286462    0.719470    -1.000000
 10 0.195300    -0.179039    1.000000
 11 -0.152696    -0.153030    1.000000
 12 0.384471    0.653336    -1.000000
 13 -0.117280    -0.153217    1.000000
 14 -0.238076    0.000583    1.000000
 15 -0.413576    0.145681    1.000000
 16 0.490767    -0.680029    -1.000000
 17 0.199894    -0.199381    1.000000
 18 -0.356048    0.537960    -1.000000
 19 -0.392868    -0.125261    1.000000
 20 0.353588    -0.070617    1.000000
 21 0.020984    0.925720    -1.000000
 22 -0.475167    -0.346247    -1.000000
 23 0.074952    0.042783    1.000000
 24 0.394164    -0.058217    1.000000
 25 0.663418    0.436525    -1.000000
 26 0.402158    0.577744    -1.000000
 27 -0.449349    -0.038074    1.000000
 28 0.619080    -0.088188    -1.000000
 29 0.268066    -0.071621    1.000000
 30 -0.015165    0.359326    1.000000
 31 0.539368    -0.374972    -1.000000
 32 -0.319153    0.629673    -1.000000
 33 0.694424    0.641180    -1.000000
 34 0.079522    0.193198    1.000000
 35 0.253289    -0.285861    1.000000
 36 -0.035558    -0.010086    1.000000
 37 -0.403483    0.474466    -1.000000
 38 -0.034312    0.995685    -1.000000
 39 -0.590657    0.438051    -1.000000
 40 -0.098871    -0.023953    1.000000
 41 -0.250001    0.141621    1.000000
 42 -0.012998    0.525985    -1.000000
 43 0.153738    0.491531    -1.000000
 44 0.388215    -0.656567    -1.000000
 45 0.049008    0.013499    1.000000
 46 0.068286    0.392741    1.000000
 47 0.747800    -0.066630    -1.000000
 48 0.004621    -0.042932    1.000000
 49 -0.701600    0.190983    -1.000000
 50 0.055413    -0.024380    1.000000
 51 0.035398    -0.333682    1.000000
 52 0.211795    0.024689    1.000000
 53 -0.045677    0.172907    1.000000
 54 0.595222    0.209570    -1.000000
 55 0.229465    0.250409    1.000000
 56 -0.089293    0.068198    1.000000
 57 0.384300    -0.176570    1.000000
 58 0.834912    -0.110321    -1.000000
 59 -0.307768    0.503038    -1.000000
 60 -0.777063    -0.348066    -1.000000
 61 0.017390    0.152441    1.000000
 62 -0.293382    -0.139778    1.000000
 63 -0.203272    0.286855    1.000000
 64 0.957812    -0.152444    -1.000000
 65 0.004609    -0.070617    1.000000
 66 -0.755431    0.096711    -1.000000
 67 -0.526487    0.547282    -1.000000
 68 -0.246873    0.833713    -1.000000
 69 0.185639    -0.066162    1.000000
 70 0.851934    0.456603    -1.000000
 71 -0.827912    0.117122    -1.000000
 72 0.233512    -0.106274    1.000000
 73 0.583671    -0.709033    -1.000000
 74 -0.487023    0.625140    -1.000000
 75 -0.448939    0.176725    1.000000
 76 0.155907    -0.166371    1.000000
 77 0.334204    0.381237    -1.000000
 78 0.081536    -0.106212    1.000000
 79 0.227222    0.527437    -1.000000
 80 0.759290    0.330720    -1.000000
 81 0.204177    -0.023516    1.000000
 82 0.577939    0.403784    -1.000000
 83 -0.568534    0.442948    -1.000000
 84 -0.011520    0.021165    1.000000
 85 0.875720    0.422476    -1.000000
 86 0.297885    -0.632874    -1.000000
 87 -0.015821    0.031226    1.000000
 88 0.541359    -0.205969    -1.000000
 89 -0.689946    -0.508674    -1.000000
 90 -0.343049    0.841653    -1.000000
 91 0.523902    -0.436156    -1.000000
 92 0.249281    -0.711840    -1.000000
 93 0.193449    0.574598    -1.000000
 94 -0.257542    -0.753885    -1.000000
 95 -0.021605    0.158080    1.000000
 96 0.601559    -0.727041    -1.000000
 97 -0.791603    0.095651    -1.000000
 98 -0.908298    -0.053376    -1.000000
 99 0.122020    0.850966    -1.000000
100 -0.725568    -0.292022    -1.000000
復制代碼

 

testData.txt

復制代碼
  1 3.542485    1.977398    -1
  2 3.018896    2.556416    -1
  3 7.551510    -1.580030    1
  4 2.114999    -0.004466    -1
  5 8.127113    1.274372    1
  6 7.108772    -0.986906    1
  7 8.610639    2.046708    1
  8 2.326297    0.265213    -1
  9 3.634009    1.730537    -1
 10 0.341367    -0.894998    -1
 11 3.125951    0.293251    -1
 12 2.123252    -0.783563    -1
 13 0.887835    -2.797792    -1
 14 7.139979    -2.329896    1
 15 1.696414    -1.212496    -1
 16 8.117032    0.623493    1
 17 8.497162    -0.266649    1
 18 4.658191    3.507396    -1
 19 8.197181    1.545132    1
 20 1.208047    0.213100    -1
 21 1.928486    -0.321870    -1
 22 2.175808    -0.014527    -1
 23 7.886608    0.461755    1
 24 3.223038    -0.552392    -1
 25 3.628502    2.190585    -1
 26 7.407860    -0.121961    1
 27 7.286357    0.251077    1
 28 2.301095    -0.533988    -1
 29 -0.232542    -0.547690    -1
 30 3.457096    -0.082216    -1
 31 3.023938    -0.057392    -1
 32 8.015003    0.885325    1
 33 8.991748    0.923154    1
 34 7.916831    -1.781735    1
 35 7.616862    -0.217958    1
 36 2.450939    0.744967    -1
 37 7.270337    -2.507834    1
 38 1.749721    -0.961902    -1
 39 1.803111    -0.176349    -1
 40 8.804461    3.044301    1
 41 1.231257    -0.568573    -1
 42 2.074915    1.410550    -1
 43 -0.743036    -1.736103    -1
 44 3.536555    3.964960    -1
 45 8.410143    0.025606    1
 46 7.382988    -0.478764    1
 47 6.960661    -0.245353    1
 48 8.234460    0.701868    1
 49 8.168618    -0.903835    1
 50 1.534187    -0.622492    -1
 51 9.229518    2.066088    1
 52 7.886242    0.191813    1
 53 2.893743    -1.643468    -1
 54 1.870457    -1.040420    -1
 55 5.286862    -2.358286    1
 56 6.080573    0.418886    1
 57 2.544314    1.714165    -1
 58 6.016004    -3.753712    1
 59 0.926310    -0.564359    -1
 60 0.870296    -0.109952    -1
 61 2.369345    1.375695    -1
 62 1.363782    -0.254082    -1
 63 7.279460    -0.189572    1
 64 1.896005    0.515080    -1
 65 8.102154    -0.603875    1
 66 2.529893    0.662657    -1
 67 1.963874    -0.365233    -1
 68 8.132048    0.785914    1
 69 8.245938    0.372366    1
 70 6.543888    0.433164    1
 71 -0.236713    -5.766721    -1
 72 8.112593    0.295839    1
 73 9.803425    1.495167    1
 74 1.497407    -0.552916    -1
 75 1.336267    -1.632889    -1
 76 9.205805    -0.586480    1
 77 1.966279    -1.840439    -1
 78 8.398012    1.584918    1
 79 7.239953    -1.764292    1
 80 7.556201    0.241185    1
 81 9.015509    0.345019    1
 82 8.266085    -0.230977    1
 83 8.545620    2.788799    1
 84 9.295969    1.346332    1
 85 2.404234    0.570278    -1
 86 2.037772    0.021919    -1
 87 1.727631    -0.453143    -1
 88 1.979395    -0.050773    -1
 89 8.092288    -1.372433    1
 90 1.667645    0.239204    -1
 91 9.854303    1.365116    1
 92 7.921057    -1.327587    1
 93 8.500757    1.492372    1
 94 1.339746    -0.291183    -1
 95 3.107511    0.758367    -1
 96 2.609525    0.902979    -1
 97 3.263585    1.367898    -1
 98 2.912122    -0.202359    -1
 99 1.731786    0.589096    -1
100 2.387003    1.573131    -1
復制代碼

 

5.2 SMO算法實現

定義數據結構體optStruct,用於緩存,提高運行速度。SMO算法具體實現如下(mySVMMLiA.py)

每個方法的作用,以及每行代碼的作用,我都做了詳細的注解,希望對大家的理解有幫助。

 # -*- coding: utf-8 -*-
 """
 Created on Tue Sep  4 16:58:16 2018
 支持向量機代碼實現
 SMO(Sequential Minimal Optimization)最小序列優化
 @author: weixw
 """
 import numpy as np
 #核轉換函數(一個特征空間映射到另一個特征空間,低維空間映射到高維空間)
 #高維空間解決線性問題,低維空間解決非線性問題
 #線性內核 = 原始數據矩陣(100*2)與原始數據第一行矩陣轉秩乘積(2*1) =>(100*1)
 #非線性內核公式:k(x,y) = exp(-||x - y||**2/2*(e**2))
 #1.原始數據每一行與原始數據第一行作差, 
 #2.平方   
 def kernelTrans(dataMat, rowDataMat, kTup):
     m,n=np.shape(dataMat)
     #初始化核矩陣 m*1
     K = np.mat(np.zeros((m,1)))
     if kTup[0] == 'lin': #線性核
         K = dataMat*rowDataMat.T
     elif kTup[0] == 'rbf':#非線性核
         for j in range(m):
             #xi - xj
             deltaRow = dataMat[j,:] - rowDataMat
             K[j] = deltaRow*deltaRow.T
         #1*m m*1 => 1*1
         K = np.exp(K/(-2*kTup[1]**2))
     else: raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
     return K
         
 #定義數據結構體,用於緩存,提高運行速度
 class optStruct:
     def __init__(self, dataSet, labelSet, C, toler, kTup):
         self.dataMat = np.mat(dataSet) #原始數據,轉換成m*n矩陣
         self.labelMat = np.mat(labelSet).T #標簽數據 m*1矩陣
         self.C = C #懲罰參數,C越大,容忍噪聲度小,需要優化;反之,容忍噪聲度高,不需要優化;
                    #所有的拉格朗日乘子都被限制在了以C為邊長的矩形里
         self.toler = toler #容忍度
         self.m = np.shape(self.dataMat)[0] #原始數據行長度
         self.alphas = np.mat(np.zeros((self.m,1))) # alpha系數,m*1矩陣
         self.b = 0 #偏置
         self.eCache = np.mat(np.zeros((self.m,2))) # 保存原始數據每行的預測值
         self.K = np.mat(np.zeros((self.m,self.m))) # 核轉換矩陣 m*m
         for i in range(self.m):
             self.K[:,i] = kernelTrans(self.dataMat, self.dataMat[i,:], kTup)
             
 #計算原始數據第k項對應的預測誤差  1*m m*1 =>1*1  
 #oS:結構數據
 #k: 原始數據行索引           
 def calEk(oS, k):
     #f(x) = w*x + b 
     fXk = float(np.multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
     Ek = fXk - float(oS.labelMat[k])
     return Ek
 
 #在alpha有改變都要更新緩存
 def updateEk(oS, k):
     Ek = calEk(oS, k)
     oS.eCache[k] = [1, Ek]
     
 
 #第一次通過selectJrand()隨機選取j,之后選取與i對應預測誤差最大的j(步長最大)
 def selectJ(i, oS, Ei):
     #初始化
     maxK = -1  #誤差最大時對應索引
     maxDeltaE = 0 #最大誤差
     Ej = 0 # j索引對應預測誤差
     #保存每一行的預測誤差值 1相對於初始化為0的更改
     oS.eCache[i] = [1,Ei]
     #獲取數據緩存結構中非0的索引列表(先將矩陣第0列轉化為數組)
     validEcacheList = np.nonzero(oS.eCache[:,0].A)[0]
     #遍歷索引列表,尋找最大誤差對應索引
     if len(validEcacheList) > 1:
         for k in validEcacheList:
             if k == i:
                 continue
             Ek = calEk(oS, k)
             deltaE = abs(Ei - Ek)
             if(deltaE > maxDeltaE):
                 maxK = k
                 maxDeltaE = deltaE
                 Ej = Ek
         return maxK, Ej
     else:
         #隨機選取一個不等於i的j
         j = selectJrand(i, oS.m)
         Ej = calEk(oS, j)
     return j,Ej
 
 #隨機選取一個不等於i的索引          
 def selectJrand(i, m):
     j = i
     while (j == i):
        j = int(np.random.uniform(0, m))
     return j
 
 #alpha范圍剪輯
 def clipAlpha(aj, L, H):
     if aj > H:
         aj = H
     if aj < L:
         aj = L
     return aj
 
 #從文件獲取特征數據,標簽數據
 def loadDataSet(fileName):
     dataSet = []; labelSet = []
     fr = open(fileName)
     for line in fr.readlines():
         #分割
         lineArr = line.strip().split('\t')
         dataSet.append([float(lineArr[0]), float(lineArr[1])])
         labelSet.append(float(lineArr[2]))
     return dataSet, labelSet
 
 #計算 w 權重系數
 def calWs(alphas, dataSet, labelSet):
     dataMat = np.mat(dataSet)
     #1*100 => 100*1
     labelMat = np.mat(labelSet).T
     m, n = np.shape(dataMat)    
     w = np.zeros((n, 1))    
     for i in range(m):
         w += np.multiply(alphas[i]*labelMat[i], dataMat[i,:].T)        
     return w
 #計算原始數據每一行alpha,b,保存到數據結構中,有變化及時更新       
 def innerL(i, oS):
     #計算預測誤差
     Ei = calEk(oS, i)
     #選擇第一個alpha,違背KKT條件2
     #正間隔,負間隔
     if ((oS.labelMat[i] * Ei < -oS.toler) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i] * Ei > oS.toler) and (oS.alphas[i] > 0)):
         #第一次隨機選取不等於i的數據項,其后根據誤差最大選取數據項
         j, Ej = selectJ(i, oS, Ei)
         #初始化,開辟新的內存
         alphaIold = oS.alphas[i].copy()
         alphaJold = oS.alphas[j].copy()
         #通過 a1y1 + a2y2 = 常量
         #    0 <= a1,a2 <= C 求出L,H
         if oS.labelMat[i] != oS.labelMat[j]:
             L = max(0, oS.alphas[j] - oS.alphas[i])
             H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
         else:
             L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
             H = min(oS.C, oS.alphas[j] + oS.alphas[i])
         if L == H : 
             print ("L == H")
             return 0
         #內核分母 K11 + K22 - 2K12
         eta = oS.K[i, i] + oS.K[j, j] - 2.0*oS.K[i, j]
         if eta <= 0:
             print ("eta <= 0")
             return 0
         #計算第一個alpha j
         oS.alphas[j] += oS.labelMat[j]*(Ei - Ej)/eta
         #修正alpha j的范圍
         oS.alphas[j] = clipAlpha(oS.alphas[j], L, H)
         #alpha有改變,就需要更新緩存數據
         updateEk(oS, j)
         #如果優化后的alpha 與之前的alpha變化很小,則舍棄,並重新選擇數據項的alpha
         if (abs(oS.alphas[j] - alphaJold) < 0.00001):
             print ("j not moving enough, abandon it.")
             return 0
         #計算alpha對的另一個alpha i
         # ai_new*yi + aj_new*yj = 常量
         # ai_old*yi + ai_old*yj = 常量 
         # 作差=> ai = ai_old + yi*yj*(aj_old - aj_new)
         oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
         #alpha有改變,就需要更新緩存數據
         updateEk(oS, i)
         #計算b1,b2
         # y(x) = w*x + b => b = y(x) - w*x
         # w = aiyixi(i= 1->N求和)
         #b1_new = y1_new - (a1_new*y1*k11 + a2_new*y2*k21 + ai*yi*ki1(i = 3 ->N求和 常量))
         #b1_old = y1_old - (a1_old*y1*k11 + a2_old*y2*k21 + ai*yi*ki1(i = 3 ->N求和 常量))
         #作差=> b1_new = b1_old + (y1_new - y1_old) - y1*k11*(a1_new - a1_old) - y2*k21*(a2_new - a2_old)
         # => b1_new = b1_old + Ei - yi*(ai_new - ai_old)*kii - yj*(aj_new - aj_old)*kij      
         #同樣可推得 b2_new = b2_old + Ej - yi*(ai_new - ai_old)*kij - yj*(aj_new - aj_old)*kjj
         bi = oS.b - Ei - oS.labelMat[i]*(oS.alphas[i] - alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j] - alphaJold)*oS.K[i,j]
         bj = oS.b - Ej - oS.labelMat[i]*(oS.alphas[i] - alphaIold)*oS.K[i,j] - oS.labelMat[j]*(oS.alphas[j] - alphaJold)*oS.K[j,j]
         #首選alpha i,相對alpha j 更准確
         if (0 < oS.alphas[i]) and (oS.alphas[i] < oS.C):
             oS.b = bi
         elif (0 < oS.alphas[j]) and (oS.alphas[j] < oS.C):
             oS.b = bj
         else:
             oS.b = (bi + bj)/2.0
         return 1
     else:
         return 0
     
 #完整SMO核心算法,包含線性核核非線性核,返回alpha,b
 #dataSet 原始特征數據
 #labelSet 標簽數據
 #C 凸二次規划參數
 #toler 容忍度
 #maxInter 循環次數
 #kTup 指定核方式
 #程序邏輯:
 #第一次全部遍歷,遍歷后根據alpha對是否有修改判斷,
 #如果alpha對沒有修改,外循環終止;如果alpha對有修改,則繼續遍歷屬於支持向量的數據。
 #直至外循環次數達到maxIter
 #相比簡單SMO算法,運行速度更快,原因是:
 #1.不是每一次都全量遍歷原始數據,第一次遍歷原始數據,
 #如果alpha有優化,就遍歷支持向量數據,直至alpha沒有優化,然后再轉全量遍歷,這是如果alpha沒有優化,循環結束;
 #2.外循環不需要達到maxInter次數就終止;
 def smoP(dataSet, labelSet, C, toler, maxInter, kTup = ('lin', 0)):
     #初始化結構體類,獲取實例
     oS = optStruct(dataSet, labelSet, C, toler, kTup)
     iter = 0
     #全量遍歷標志
     entireSet = True
     #alpha對是否優化標志
     alphaPairsChanged = 0
     #外循環 終止條件:1.達到最大次數 或者 2.alpha對沒有優化
     while (iter < maxInter) and ((alphaPairsChanged > 0) or (entireSet)):
         alphaPairsChanged = 0
         #全量遍歷 ,遍歷每一行數據 alpha對有修改,alphaPairsChanged累加
         if entireSet:
             for i in range(oS.m):
                 alphaPairsChanged += innerL(i, oS)
                 print ("fullSet, iter: %d i:%d, pairs changed %d" %(iter, i, alphaPairsChanged))
             iter += 1
         else:
             #獲取(0,C)范圍內數據索引列表,也就是只遍歷屬於支持向量的數據
             nonBounds = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
             for i in nonBounds:
                 alphaPairsChanged += innerL(i, oS)
                 print ("non-bound, iter: %d i:%d, pairs changed %d" %(iter, i, alphaPairsChanged))
             iter += 1
         #全量遍歷->支持向量遍歷
         if entireSet:
             entireSet = False
         #支持向量遍歷->全量遍歷
         elif alphaPairsChanged == 0:
             entireSet = True
         print ("iteation number: %d"% iter)
         print ("entireSet :%s"% entireSet)
         print ("alphaPairsChanged :%d"% alphaPairsChanged)
     return oS.b,oS.alphas
 
 #繪制支持向量
 def drawDataMap(dataArr,labelArr,b,alphas):
     import matplotlib.pyplot as plt
     #alphas.A>0 獲取大於0的索引列表,只有>0的alpha才對分類起作用
     svInd=np.nonzero(alphas.A>0)[0]           
      #分類數據點
     classified_pts = {'+1':[],'-1':[]}
     for point,label in zip(dataArr,labelArr):
         if label == 1.0:
             classified_pts['+1'].append(point)
         else:
             classified_pts['-1'].append(point)
     fig = plt.figure()
     ax = fig.add_subplot(111)
     #繪制數據點
     for label,pts in classified_pts.items():
         pts = np.array(pts)
         ax.scatter(pts[:, 0], pts[:, 1], label = label)
     #繪制分割線
     w = calWs(alphas, dataArr, labelArr)
     #函數形式:max( x ,key=lambda a : b )        #    x可以是任何數值,可以有多個x值
     #先把x值帶入lambda函數轉換成b值,然后再將b值進行比較
     x1, _=max(dataArr, key=lambda x:x[0])
     x2, _=min(dataArr, key=lambda x:x[0])    
     a1, a2 = w
     y1, y2 = (-b - a1*x1)/a2, (-b - a1*x2)/a2
     #矩陣轉化為數組.A
     ax.plot([x1, x2],[y1.A[0][0], y2.A[0][0]])
     
     #繪制支持向量
     for i in svInd:
         x, y= dataArr[i]        
         ax.scatter([x], [y], s=150, c ='none', alpha=0.7, linewidth=1.5, edgecolor = '#AB3319')
     plt.show()
     
      #alpha>0對應的數據才是支持向量,過濾不是支持向量的數據
     sVs= np.mat(dataArr)[svInd] #get matrix of only support vectors
     print ("there are %d Support Vectors.\n" % np.shape(sVs)[0])
     
 #訓練結果    
 def getTrainingDataResult(dataSet, labelSet, b, alphas, k1=1.3):
     datMat = np.mat(dataSet)
     #100*1
     labelMat = np.mat(labelSet).T
     #alphas.A>0 獲取大於0的索引列表,只有>0的alpha才對分類起作用
     svInd=np.nonzero(alphas.A>0)[0]
     sVs=datMat[svInd]
     labelSV = labelMat[svInd];
     m,n = np.shape(datMat)
     errorCount = 0
     for i in range(m):
         kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
         # y(x) = w*x + b => b = y(x) - w*x
         # w = aiyixi(i= 1->N求和)
         predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b
         if np.sign(predict)!=np.sign(labelSet[i]): errorCount += 1
     print ("the training error rate is: %f" % (float(errorCount)/m))
     
 def getTestDataResult(dataSet, labelSet, b, alphas, k1=1.3):
     datMat = np.mat(dataSet)
     #100*1
     labelMat = np.mat(labelSet).T
     #alphas.A>0 獲取大於0的索引列表,只有>0的alpha才對分類起作用
     svInd=np.nonzero(alphas.A>0)[0]
     sVs=datMat[svInd]
     labelSV = labelMat[svInd];
     m,n = np.shape(datMat)
     errorCount = 0
     for i in range(m):
         kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
         # y(x) = w*x + b => b = y(x) - w*x
         # w = aiyixi(i= 1->N求和)
         predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b
         if np.sign(predict)!=np.sign(labelSet[i]): errorCount += 1    
     print ("the test error rate is: %f" % (float(errorCount)/m))  
     

 

5.3 測試代碼(testMySVMMLiA.py)

 

 # -*- coding: utf-8 -*-
 """
 Created on Wed Sep  5 15:22:26 2018
 
 @author: weixw
 """
 
 
 import mySVMMLiA as sm
 
 #通過訓練數據計算 b, alphas
 dataArr,labelArr = sm.loadDataSet('trainingData.txt')
 b, alphas = sm.smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', 0.10))
 sm.drawDataMap(dataArr,labelArr,b,alphas)
 sm.getTrainingDataResult(dataArr, labelArr, b, alphas, 0.10)
 dataArr1,labelArr1 = sm.loadDataSet('testData.txt')
 #測試結果
 sm.getTestDataResult(dataArr1, labelArr1, b, alphas, 0.10)

通過訓練數據計算出 b, 權重矩陣,從而分類超平面和決策分類函數就明確了,然后測試數據以決策分類函數進行預測。

 這里采用高斯核RBF。

5.4 運行結果

“j not moving enough, abandon it”  表示數據項對應的  和   非常接近,不需要優化;

“fullSet” 表示全量數據遍歷;

 “non-bound” 表示非邊界遍歷,也就是只遍歷屬於支持向量的數據項。

 

另外我將支持向量的數據項繪制出來了,這樣更直觀。

 

 

可以看出,有77個支持向量,訓練差錯率是0,測試差錯率6%。

不要讓懶惰占據你的大腦,不要讓妥協拖垮了你的人生。青春就是一張票,能不能趕上時代的快車,你的步伐就掌握在你的腳下。

     

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM