
機(jī)器學(xué)習(xí)算法與Python實(shí)踐之(四)支持向量機(jī)(SVM)實(shí)現(xiàn)
八、SVM的實(shí)現(xiàn)之SMO算法
終于到SVM的實(shí)現(xiàn)部分了。那么神奇和有效的東西還得回歸到實(shí)現(xiàn)才可以展示其強(qiáng)大的功力。SVM有效而且存在很高效的訓(xùn)練算法,這也是工業(yè)界非常青睞SVM的原因。
前面講到,SVM的學(xué)習(xí)問(wèn)題可以轉(zhuǎn)化為下面的對(duì)偶問(wèn)題:
需要滿(mǎn)足的KKT條件:
也就是說(shuō)找到一組αi可以滿(mǎn)足上面的這些條件的就是該目標(biāo)的一個(gè)最優(yōu)解。所以我們的優(yōu)化目標(biāo)是找到一組最優(yōu)的αi*。一旦求出這些αi*,就很容易計(jì)算出權(quán)重向量w*和b,并得到分隔超平面了。
這是個(gè)凸二次規(guī)劃問(wèn)題,它具有全局最優(yōu)解,一般可以通過(guò)現(xiàn)有的工具來(lái)優(yōu)化。但當(dāng)訓(xùn)練樣本非常多的時(shí)候,這些優(yōu)化算法往往非常耗時(shí)低效,以致無(wú)法使用。從SVM提出到現(xiàn)在,也出現(xiàn)了很多優(yōu)化訓(xùn)練的方法。其中,非常出名的一個(gè)是1982年由Microsoft Research的John C. Platt在論文《Sequential Minimal Optimization: A Fast Algorithm for TrainingSupport Vector Machines》中提出的Sequential Minimal Optimization序列最小化優(yōu)化算法,簡(jiǎn)稱(chēng)SMO算法。SMO算法的思想很簡(jiǎn)單,它將大優(yōu)化的問(wèn)題分解成多個(gè)小優(yōu)化的問(wèn)題。這些小問(wèn)題往往比較容易求解,并且對(duì)他們進(jìn)行順序求解的結(jié)果與將他們作為整體來(lái)求解的結(jié)果完全一致。在結(jié)果完全一致的同時(shí),SMO的求解時(shí)間短很多。在深入SMO算法之前,我們先來(lái)了解下坐標(biāo)下降這個(gè)算法,SMO其實(shí)基于這種簡(jiǎn)單的思想的。
8.1、坐標(biāo)下降(上升)法
假設(shè)要求解下面的優(yōu)化問(wèn)題:
在這里,我們需要求解m個(gè)變量αi,一般來(lái)說(shuō)是通過(guò)梯度下降(這里是求最大值,所以應(yīng)該叫上升)等算法每一次迭代對(duì)所有m個(gè)變量αi也就是α向量進(jìn)行一次性?xún)?yōu)化。通過(guò)誤差每次迭代調(diào)整α向量中每個(gè)元素的值。而坐標(biāo)上升法(坐標(biāo)上升與坐標(biāo)下降可以看做是一對(duì),坐標(biāo)上升是用來(lái)求解max最優(yōu)化問(wèn)題,坐標(biāo)下降用于求min最優(yōu)化問(wèn)題)的思想是每次迭代只調(diào)整一個(gè)變量αi的值,其他變量的值在這次迭代中固定不變。
最里面語(yǔ)句的意思是固定除αi之外的所有αj(i不等于j),這時(shí)W可看作只是關(guān)于αi的函數(shù),那么直接對(duì)αi求導(dǎo)優(yōu)化即可。這里我們進(jìn)行最大化求導(dǎo)的順序i是從1到m,可以通過(guò)更改優(yōu)化順序來(lái)使W能夠更快地增加并收斂。如果W在內(nèi)循環(huán)中能夠很快地達(dá)到最優(yōu),那么坐標(biāo)上升法會(huì)是一個(gè)很高效的求極值方法。
用個(gè)二維的例子來(lái)說(shuō)明下坐標(biāo)下降法:我們需要尋找f(x,y)=x2+xy+y2的最小值處的(x*, y*),也就是下圖的F*點(diǎn)的地方。
假設(shè)我們初始的點(diǎn)是A(圖是函數(shù)投影到xoy平面的等高線(xiàn)圖,顏色越深值越?。?,我們需要達(dá)到F*的地方。那最快的方法就是圖中黃色線(xiàn)的路徑,一次性就到達(dá)了,其實(shí)這個(gè)是牛頓優(yōu)化法,但如果是高維的話(huà),這個(gè)方法就不太高效了(因?yàn)樾枰蠼饩仃嚨哪?,這個(gè)不在這里討論)。我們也可以按照紅色所指示的路徑來(lái)走。從A開(kāi)始,先固定x,沿著y軸往讓f(x, y)值減小的方向走到B點(diǎn),然后固定y,沿著x軸往讓f(x, y)值減小的方向走到C點(diǎn),不斷循環(huán),直到到達(dá)F*。反正每次只要我們都往讓f(x, y)值小的地方走就行了,這樣腳踏實(shí)地,一步步走,每一步都使f(x, y)慢慢變小,總有一天,皇天不負(fù)有心人的。到達(dá)F*也是時(shí)間問(wèn)題。到這里你可能會(huì)說(shuō),這紅色線(xiàn)比黃色線(xiàn)貧富差距也太嚴(yán)重了吧。因?yàn)檫@里是二維的簡(jiǎn)單的情況嘛。如果是高維的情況,而且目標(biāo)函數(shù)很復(fù)雜的話(huà),再加上樣本集很多,那么在梯度下降中,目標(biāo)函數(shù)對(duì)所有αi求梯度或者在牛頓法中對(duì)矩陣求逆,都是很耗時(shí)的。這時(shí)候,如果W只對(duì)單個(gè)αi優(yōu)化很快的時(shí)候,坐標(biāo)下降法可能會(huì)更加高效。
8.2、SMO算法
SMO算法的思想和坐標(biāo)下降法的思想差不多。唯一不同的是,SMO是一次迭代優(yōu)化兩個(gè)α而不是一個(gè)。為什么要優(yōu)化兩個(gè)呢?
我們回到這個(gè)優(yōu)化問(wèn)題。我們可以看到這個(gè)優(yōu)化問(wèn)題存在著一個(gè)約束,也就是
假設(shè)我們首先固定除α1以外的所有參數(shù),然后在α1上求極值。但需要注意的是,因?yàn)槿绻潭é?以外的所有參數(shù),由上面這個(gè)約束條件可以知道,α1將不再是變量(可以由其他值推出),因?yàn)閱?wèn)題中規(guī)定了:
因此,我們需要一次選取兩個(gè)參數(shù)做優(yōu)化,比如αi和αj,此時(shí)αi可以由αj和其他參數(shù)表示出來(lái)。這樣回代入W中,W就只是關(guān)于αj的函數(shù)了,這時(shí)候就可以只對(duì)αj進(jìn)行優(yōu)化了。在這里就是對(duì)αj進(jìn)行求導(dǎo),令導(dǎo)數(shù)為0就可以解出這個(gè)時(shí)候最優(yōu)的αj了。然后也可以得到αi。這就是一次的迭代過(guò)程,一次迭代只調(diào)整兩個(gè)拉格朗日乘子αi和αj。SMO之所以高效就是因?yàn)樵诠潭ㄆ渌麉?shù)后,對(duì)一個(gè)參數(shù)優(yōu)化過(guò)程很高效(對(duì)一個(gè)參數(shù)的優(yōu)化可以通過(guò)解析求解,而不是迭代。雖然對(duì)一個(gè)參數(shù)的一次最小優(yōu)化不可能保證其結(jié)果就是所優(yōu)化的拉格朗日乘子的最終結(jié)果,但會(huì)使目標(biāo)函數(shù)向極小值邁進(jìn)一步,這樣對(duì)所有的乘子做最小優(yōu)化,直到所有滿(mǎn)足KKT條件時(shí),目標(biāo)函數(shù)達(dá)到最?。?
總結(jié)下來(lái)是:
重復(fù)下面過(guò)程直到收斂{
(1)選擇兩個(gè)拉格朗日乘子αi和αj;
(2)固定其他拉格朗日乘子αk(k不等于i和j),只對(duì)αi和αj優(yōu)化w(α);
(3)根據(jù)優(yōu)化后的αi和αj,更新截距b的值;
}
那訓(xùn)練里面這兩三步驟到底是怎么實(shí)現(xiàn)的,需要考慮什么呢?下面我們來(lái)具體分析下:
(1)選擇αi和αj:
我們現(xiàn)在是每次迭代都優(yōu)化目標(biāo)函數(shù)的兩個(gè)拉格朗日乘子αi和αj,然后其他的拉格朗日乘子保持固定。如果有N個(gè)訓(xùn)練樣本,我們就有N個(gè)拉格朗日乘子需要優(yōu)化,但每次我們只挑兩個(gè)進(jìn)行優(yōu)化,我們就有N(N-1)種選擇。那到底我們要選擇哪對(duì)αi和αj呢?選擇哪對(duì)才好呢?想想我們的目標(biāo)是什么?我們希望把所有違法KKT條件的樣本都糾正回來(lái),因?yàn)槿绻袠颖径紳M(mǎn)足KKT條件的話(huà),我們的優(yōu)化就完成了。那就很直觀了,哪個(gè)害群之馬最嚴(yán)重,我們得先對(duì)他進(jìn)行思想教育,讓他盡早回歸正途。OK,那我們選擇的第一個(gè)變量αi就選違法KKT條件最嚴(yán)重的那一個(gè)。那第二個(gè)變量αj怎么選呢?
我們是希望快點(diǎn)找到最優(yōu)的N個(gè)拉格朗日乘子,使得代價(jià)函數(shù)最大,換句話(huà)說(shuō),要最快的找到代價(jià)函數(shù)最大值的地方對(duì)應(yīng)的N個(gè)拉格朗日乘子。這樣我們的訓(xùn)練時(shí)間才會(huì)短。就像你從廣州去北京,有飛機(jī)和綠皮車(chē)給你選,你選啥?(就算你不考慮速度,也得考慮下空姐的感受嘛,別辜負(fù)了她們渴望看到你的期盼,哈哈)。有點(diǎn)離題了,anyway,每次迭代中,哪對(duì)αi和αj可以讓我更快的達(dá)到代價(jià)函數(shù)值最大的地方,我們就選他們?;蛘哒f(shuō),走完這一步,選這對(duì)αi和αj代價(jià)函數(shù)值增加的值最多,比選擇其他所有αi和αj的結(jié)合中都多。這樣我們才可以更快的接近代價(jià)函數(shù)的最大值,也就是達(dá)到優(yōu)化的目標(biāo)了。再例如,下圖,我們要從A點(diǎn)走到B點(diǎn),按藍(lán)色的路線(xiàn)走c2方向的時(shí)候,一跨一大步,按紅色的路線(xiàn)走c1方向的時(shí)候,只能是人類(lèi)的一小步。所以,藍(lán)色路線(xiàn)走兩步就邁進(jìn)了成功之門(mén),而紅色的路線(xiàn),人生曲折,好像成功遙遙無(wú)期一樣,故曰,選擇比努力更重要!
真啰嗦!說(shuō)了半天,其實(shí)就一句話(huà):為什么每次迭代都要選擇最好的αi和αj,就是為了更快的收斂!那實(shí)踐中每次迭代到底要怎樣選αi和αj呢?這有個(gè)很好聽(tīng)的名字叫啟發(fā)式選擇,主要思想是先選擇最有可能需要優(yōu)化(也就是違反KKT條件最嚴(yán)重)的αi,再針對(duì)這樣的αi選擇最有可能取得較大修正步長(zhǎng)的αj。具體是以下兩個(gè)過(guò)程:
1)第一個(gè)變量αi的選擇:
SMO稱(chēng)選擇第一個(gè)變量的過(guò)程為外層循環(huán)。外層訓(xùn)練在訓(xùn)練樣本中選取違法KKT條件最嚴(yán)重的樣本點(diǎn)。并將其對(duì)應(yīng)的變量作為第一個(gè)變量。具體的,檢驗(yàn)訓(xùn)練樣本(xi, yi)是否滿(mǎn)足KKT條件,也就是:
該檢驗(yàn)是在ε范圍內(nèi)進(jìn)行的。在檢驗(yàn)過(guò)程中,外層循環(huán)首先遍歷所有滿(mǎn)足條件0<αj<C的樣本點(diǎn),即在間隔邊界上的支持向量點(diǎn),檢驗(yàn)他們是否滿(mǎn)足KKT條件,然后選擇違反KKT條件最嚴(yán)重的αi。如果這些樣本點(diǎn)都滿(mǎn)足KKT條件,那么遍歷整個(gè)訓(xùn)練集,檢驗(yàn)他們是否滿(mǎn)足KKT條件,然后選擇違反KKT條件最嚴(yán)重的αi。
優(yōu)先選擇遍歷非邊界數(shù)據(jù)樣本,因?yàn)榉沁吔鐢?shù)據(jù)樣本更有可能需要調(diào)整,邊界數(shù)據(jù)樣本常常不能得到進(jìn)一步調(diào)整而留在邊界上。由于大部分?jǐn)?shù)據(jù)樣本都很明顯不可能是支持向量,因此對(duì)應(yīng)的α乘子一旦取得零值就無(wú)需再調(diào)整。遍歷非邊界數(shù)據(jù)樣本并選出他們當(dāng)中違反KKT 條件為止。當(dāng)某一次遍歷發(fā)現(xiàn)沒(méi)有非邊界數(shù)據(jù)樣本得到調(diào)整時(shí),遍歷所有數(shù)據(jù)樣本,以檢驗(yàn)是否整個(gè)集合都滿(mǎn)足KKT條件。如果整個(gè)集合的檢驗(yàn)中又有數(shù)據(jù)樣本被進(jìn)一步進(jìn)化,則有必要再遍歷非邊界數(shù)據(jù)樣本。這樣,不停地在遍歷所有數(shù)據(jù)樣本和遍歷非邊界數(shù)據(jù)樣本之間切換,直到整個(gè)樣本集合都滿(mǎn)足KKT條件為止。以上用KKT條件對(duì)數(shù)據(jù)樣本所做的檢驗(yàn)都以達(dá)到一定精度ε就可以停止為條件。如果要求十分精確的輸出算法,則往往不能很快收斂。
對(duì)整個(gè)數(shù)據(jù)集的遍歷掃描相當(dāng)容易,而實(shí)現(xiàn)對(duì)非邊界αi的掃描時(shí),首先需要將所有非邊界樣本的αi值(也就是滿(mǎn)足0<αi<C)保存到新的一個(gè)列表中,然后再對(duì)其進(jìn)行遍歷。同時(shí),該步驟跳過(guò)那些已知的不會(huì)改變的αi值。
2)第二個(gè)變量αj的選擇:
在選擇第一個(gè)αi后,算法會(huì)通過(guò)一個(gè)內(nèi)循環(huán)來(lái)選擇第二個(gè)αj值。因?yàn)榈诙€(gè)乘子的迭代步長(zhǎng)大致正比于|Ei-Ej|,所以我們需要選擇能夠最大化|Ei-Ej|的第二個(gè)乘子(選擇最大化迭代步長(zhǎng)的第二個(gè)乘子)。在這里,為了節(jié)省計(jì)算時(shí)間,我們建立一個(gè)全局的緩存用于保存所有樣本的誤差值,而不用每次選擇的時(shí)候就重新計(jì)算。我們從中選擇使得步長(zhǎng)最大或者|Ei-Ej|最大的αj。
(2)優(yōu)化αi和αj:
選擇這兩個(gè)拉格朗日乘子后,我們需要先計(jì)算這些參數(shù)的約束值。然后再求解這個(gè)約束最大化問(wèn)題。
首先,我們需要給αj找到邊界L<=αj<=H,以保證αj滿(mǎn)足0<=αj<=C的約束。這意味著αj必須落入這個(gè)盒子中。由于只有兩個(gè)變量(αi, αj),約束可以用二維空間中的圖形來(lái)表示,如下圖:
不等式約束使得(αi,αj)在盒子[0, C]x[0, C]內(nèi),等式約束使得(αi, αj)在平行于盒子[0, C]x[0, C]的對(duì)角線(xiàn)的直線(xiàn)上。因此要求的是目標(biāo)函數(shù)在一條平行于對(duì)角線(xiàn)的線(xiàn)段上的最優(yōu)值。這使得兩個(gè)變量的最優(yōu)化問(wèn)題成為實(shí)質(zhì)的單變量的最優(yōu)化問(wèn)題。由圖可以得到,αj的上下界可以通過(guò)下面的方法得到:
我們優(yōu)化的時(shí)候,αj必須要滿(mǎn)足上面這個(gè)約束。也就是說(shuō)上面是αj的可行域。然后我們開(kāi)始尋找αj,使得目標(biāo)函數(shù)最大化。通過(guò)推導(dǎo)得到αj的更新公式如下:
這里Ek可以看做對(duì)第k個(gè)樣本,SVM的輸出與期待輸出,也就是樣本標(biāo)簽的誤差。
而η實(shí)際上是度量?jī)蓚€(gè)樣本i和j的相似性的。在計(jì)算η的時(shí)候,我們需要使用核函數(shù),那么就可以用核函數(shù)來(lái)取代上面的內(nèi)積。
得到新的αj后,我們需要保證它處于邊界內(nèi)。換句話(huà)說(shuō),如果這個(gè)優(yōu)化后的值跑出了邊界L和H,我們就需要簡(jiǎn)單的裁剪,將αj收回這個(gè)范圍:
最后,得到優(yōu)化的αj后,我們需要用它來(lái)計(jì)算αi:
到這里,αi和αj的優(yōu)化就完成了。
(3)計(jì)算閾值b:
優(yōu)化αi和αj后,我們就可以更新閾值b,使得對(duì)兩個(gè)樣本i和j都滿(mǎn)足KKT條件。如果優(yōu)化后αi不在邊界上(也就是滿(mǎn)足0<αi<C,這時(shí)候根據(jù)KKT條件,可以得到y(tǒng)igi(xi)=1,這樣我們才可以計(jì)算b),那下面的閾值b1是有效的,因?yàn)楫?dāng)輸入xi時(shí)它迫使SVM輸出yi。
同樣,如果0<αj<C,那么下面的b2也是有效的:
如果0<αi<C和0<αj<C都滿(mǎn)足,那么b1和b2都有效,而且他們是相等的。如果他們兩個(gè)都處于邊界上(也就是αi=0或者αi=C,同時(shí)αj=0或者αj=C),那么在b1和b2之間的閾值都滿(mǎn)足KKT條件,一般我們?nèi)∷麄兊钠骄礲=(b1+b2)/2。所以,總的來(lái)說(shuō)對(duì)b的更新如下:
每做完一次最小優(yōu)化,必須更新每個(gè)數(shù)據(jù)樣本的誤差,以便用修正過(guò)的分類(lèi)面對(duì)其他數(shù)據(jù)樣本再做檢驗(yàn),在選擇第二個(gè)配對(duì)優(yōu)化數(shù)據(jù)樣本時(shí)用來(lái)估計(jì)步長(zhǎng)。
(4)凸優(yōu)化問(wèn)題終止條件:
SMO算法的基本思路是:如果說(shuō)有變量的解都滿(mǎn)足此最優(yōu)化問(wèn)題的KKT條件,那么這個(gè)最優(yōu)化問(wèn)題的解就得到了。因?yàn)镵KT條件是該最優(yōu)化問(wèn)題的充分必要條件(證明請(qǐng)參考文獻(xiàn))。所以我們可以監(jiān)視原問(wèn)題的KKT條件,所以所有的樣本都滿(mǎn)足KKT條件,那么就表示迭代結(jié)束了。但是由于KKT條件本身是比較苛刻的,所以也需要設(shè)定一個(gè)容忍值,即所有樣本在容忍值范圍內(nèi)滿(mǎn)足KKT條件則認(rèn)為訓(xùn)練可以結(jié)束;當(dāng)然了,對(duì)于對(duì)偶問(wèn)題的凸優(yōu)化還有其他終止條件,可以參考文獻(xiàn)。
8.3、SMO算法的Python實(shí)現(xiàn)
8.3.1、Python的準(zhǔn)備工作
我使用的Python是2.7.5版本的。附加的庫(kù)有Numpy和Matplotlib。而Matplotlib又依賴(lài)dateutil和pyparsing兩個(gè)庫(kù),所以我們需要安裝以上三個(gè)庫(kù)。前面兩個(gè)庫(kù)還好安裝,直接在官網(wǎng)下對(duì)應(yīng)版本就行。但我找后兩個(gè)庫(kù)的時(shí)候,就沒(méi)那么容易了。后來(lái)發(fā)現(xiàn),其實(shí)對(duì)Python的庫(kù)的下載和安裝可以借助pip工具的。這個(gè)是安裝和管理Python包的工具。
8.3.2、SMO算法的Python實(shí)現(xiàn)
在代碼中已經(jīng)有了比較詳細(xì)的注釋了。不知道有沒(méi)有錯(cuò)誤的地方,如果有,還望大家指正(每次的運(yùn)行結(jié)果都有可能不同,另外,感覺(jué)有些結(jié)果似乎不太正確,但我還沒(méi)發(fā)現(xiàn)哪里出錯(cuò)了,如果大家找到有錯(cuò)誤的地方,還望大家指點(diǎn)下,衷心感謝)。里面我寫(xiě)了個(gè)可視化結(jié)果的函數(shù),但只能在二維的數(shù)據(jù)上面使用。直接貼代碼:
SVM.py
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
#################################################
# SVM: support vector machine
# Author : zouxy
# Date : 2013-12-12
# HomePage : http://blog.csdn.net/zouxy09
# Email : zouxy09@qq.com
#################################################
from numpy import *
import time
import matplotlib.pyplot as plt
# calulate kernel value
def calcKernelValue(matrix_x, sample_x, kernelOption):
kernelType = kernelOption[0]
numSamples = matrix_x.shape[0]
kernelValue = mat(zeros((numSamples, 1)))
if kernelType == 'linear':
kernelValue = matrix_x * sample_x.T
elif kernelType == 'rbf':
sigma = kernelOption[1]
if sigma == 0:
sigma = 1.0
for i in xrange(numSamples):
diff = matrix_x[i, :] - sample_x
kernelValue[i] = exp(diff * diff.T / (-2.0 * sigma**2))
else:
raise NameError('Not support kernel type! You can use linear or rbf!')
return kernelValue
# calculate kernel matrix given train set and kernel type
def calcKernelMatrix(train_x, kernelOption):
numSamples = train_x.shape[0]
kernelMatrix = mat(zeros((numSamples, numSamples)))
for i in xrange(numSamples):
kernelMatrix[:, i] = calcKernelValue(train_x, train_x[i, :], kernelOption)
return kernelMatrix
# define a struct just for storing variables and data
class SVMStruct:
def __init__(self, dataSet, labels, C, toler, kernelOption):
self.train_x = dataSet # each row stands for a sample
self.train_y = labels # corresponding label
self.C = C # slack variable
self.toler = toler # termination condition for iteration
self.numSamples = dataSet.shape[0] # number of samples
self.alphas = mat(zeros((self.numSamples, 1))) # Lagrange factors for all samples
self.b = 0
self.errorCache = mat(zeros((self.numSamples, 2)))
self.kernelOpt = kernelOption
self.kernelMat = calcKernelMatrix(self.train_x, self.kernelOpt)
# calculate the error for alpha k
def calcError(svm, alpha_k):
output_k = float(multiply(svm.alphas, svm.train_y).T * svm.kernelMat[:, alpha_k] + svm.b)
error_k = output_k - float(svm.train_y[alpha_k])
return error_k
# update the error cache for alpha k after optimize alpha k
def updateError(svm, alpha_k):
error = calcError(svm, alpha_k)
svm.errorCache[alpha_k] = [1, error]
# select alpha j which has the biggest step
def selectAlpha_j(svm, alpha_i, error_i):
svm.errorCache[alpha_i] = [1, error_i] # mark as valid(has been optimized)
candidateAlphaList = nonzero(svm.errorCache[:, 0].A)[0] # mat.A return array
maxStep = 0; alpha_j = 0; error_j = 0
# find the alpha with max iterative step
if len(candidateAlphaList) > 1:
for alpha_k in candidateAlphaList:
if alpha_k == alpha_i:
continue
error_k = calcError(svm, alpha_k)
if abs(error_k - error_i) > maxStep:
maxStep = abs(error_k - error_i)
alpha_j = alpha_k
error_j = error_k
# if came in this loop first time, we select alpha j randomly
else:
alpha_j = alpha_i
while alpha_j == alpha_i:
alpha_j = int(random.uniform(0, svm.numSamples))
error_j = calcError(svm, alpha_j)
return alpha_j, error_j
# the inner loop for optimizing alpha i and alpha j
def innerLoop(svm, alpha_i):
error_i = calcError(svm, alpha_i)
### check and pick up the alpha who violates the KKT condition
## satisfy KKT condition
# 1) yi*f(i) >= 1 and alpha == 0 (outside the boundary)
# 2) yi*f(i) == 1 and 0<alpha< C (on the boundary)
# 3) yi*f(i) <= 1 and alpha == C (between the boundary)
## violate KKT condition
# because y[i]*E_i = y[i]*f(i) - y[i]^2 = y[i]*f(i) - 1, so
# 1) if y[i]*E_i < 0, so yi*f(i) < 1, if alpha < C, violate!(alpha = C will be correct)
# 2) if y[i]*E_i > 0, so yi*f(i) > 1, if alpha > 0, violate!(alpha = 0 will be correct)
# 3) if y[i]*E_i = 0, so yi*f(i) = 1, it is on the boundary, needless optimized
if (svm.train_y[alpha_i] * error_i < -svm.toler) and (svm.alphas[alpha_i] < svm.C) or\
(svm.train_y[alpha_i] * error_i > svm.toler) and (svm.alphas[alpha_i] > 0):
# step 1: select alpha j
alpha_j, error_j = selectAlpha_j(svm, alpha_i, error_i)
alpha_i_old = svm.alphas[alpha_i].copy()
alpha_j_old = svm.alphas[alpha_j].copy()
# step 2: calculate the boundary L and H for alpha j
if svm.train_y[alpha_i] != svm.train_y[alpha_j]:
L = max(0, svm.alphas[alpha_j] - svm.alphas[alpha_i])
H = min(svm.C, svm.C + svm.alphas[alpha_j] - svm.alphas[alpha_i])
else:
L = max(0, svm.alphas[alpha_j] + svm.alphas[alpha_i] - svm.C)
H = min(svm.C, svm.alphas[alpha_j] + svm.alphas[alpha_i])
if L == H:
return 0
# step 3: calculate eta (the similarity of sample i and j)
eta = 2.0 * svm.kernelMat[alpha_i, alpha_j] - svm.kernelMat[alpha_i, alpha_i] \
- svm.kernelMat[alpha_j, alpha_j]
if eta >= 0:
return 0
# step 4: update alpha j
svm.alphas[alpha_j] -= svm.train_y[alpha_j] * (error_i - error_j) / eta
# step 5: clip alpha j
if svm.alphas[alpha_j] > H:
svm.alphas[alpha_j] = H
if svm.alphas[alpha_j] < L:
svm.alphas[alpha_j] = L
# step 6: if alpha j not moving enough, just return
if abs(alpha_j_old - svm.alphas[alpha_j]) < 0.00001:
updateError(svm, alpha_j)
return 0
# step 7: update alpha i after optimizing aipha j
svm.alphas[alpha_i] += svm.train_y[alpha_i] * svm.train_y[alpha_j] \
* (alpha_j_old - svm.alphas[alpha_j])
# step 8: update threshold b
b1 = svm.b - error_i - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \
* svm.kernelMat[alpha_i, alpha_i] \
- svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \
* svm.kernelMat[alpha_i, alpha_j]
b2 = svm.b - error_j - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \
* svm.kernelMat[alpha_i, alpha_j] \
- svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \
* svm.kernelMat[alpha_j, alpha_j]
if (0 < svm.alphas[alpha_i]) and (svm.alphas[alpha_i] < svm.C):
svm.b = b1
elif (0 < svm.alphas[alpha_j]) and (svm.alphas[alpha_j] < svm.C):
svm.b = b2
else:
svm.b = (b1 + b2) / 2.0
# step 9: update error cache for alpha i, j after optimize alpha i, j and b
updateError(svm, alpha_j)
updateError(svm, alpha_i)
return 1
else:
return 0
# the main training procedure
def trainSVM(train_x, train_y, C, toler, maxIter, kernelOption = ('rbf', 1.0)):
# calculate training time
startTime = time.time()
# init data struct for svm
svm = SVMStruct(mat(train_x), mat(train_y), C, toler, kernelOption)
# start training
entireSet = True
alphaPairsChanged = 0
iterCount = 0
# Iteration termination condition:
# Condition 1: reach max iteration
# Condition 2: no alpha changed after going through all samples,
# in other words, all alpha (samples) fit KKT condition
while (iterCount < maxIter) and ((alphaPairsChanged > 0) or entireSet):
alphaPairsChanged = 0
# update alphas over all training examples
if entireSet:
for i in xrange(svm.numSamples):
alphaPairsChanged += innerLoop(svm, i)
print '---iter:%d entire set, alpha pairs changed:%d' % (iterCount, alphaPairsChanged)
iterCount += 1
# update alphas over examples where alpha is not 0 & not C (not on boundary)
else:
nonBoundAlphasList = nonzero((svm.alphas.A > 0) * (svm.alphas.A < svm.C))[0]
for i in nonBoundAlphasList:
alphaPairsChanged += innerLoop(svm, i)
print '---iter:%d non boundary, alpha pairs changed:%d' % (iterCount, alphaPairsChanged)
iterCount += 1
# alternate loop over all examples and non-boundary examples
if entireSet:
entireSet = False
elif alphaPairsChanged == 0:
entireSet = True
print 'Congratulations, training complete! Took %fs!' % (time.time() - startTime)
return svm
# testing your trained svm model given test set
def testSVM(svm, test_x, test_y):
test_x = mat(test_x)
test_y = mat(test_y)
numTestSamples = test_x.shape[0]
supportVectorsIndex = nonzero(svm.alphas.A > 0)[0]
supportVectors = svm.train_x[supportVectorsIndex]
supportVectorLabels = svm.train_y[supportVectorsIndex]
supportVectorAlphas = svm.alphas[supportVectorsIndex]
matchCount = 0
for i in xrange(numTestSamples):
kernelValue = calcKernelValue(supportVectors, test_x[i, :], svm.kernelOpt)
predict = kernelValue.T * multiply(supportVectorLabels, supportVectorAlphas) + svm.b
if sign(predict) == sign(test_y[i]):
matchCount += 1
accuracy = float(matchCount) / numTestSamples
return accuracy
# show your trained svm model only available with 2-D data
def showSVM(svm):
if svm.train_x.shape[1] != 2:
print "Sorry! I can not draw because the dimension of your data is not 2!"
return 1
# draw all samples
for i in xrange(svm.numSamples):
if svm.train_y[i] == -1:
plt.plot(svm.train_x[i, 0], svm.train_x[i, 1], 'or')
elif svm.train_y[i] == 1:
plt.plot(svm.train_x[i, 0], svm.train_x[i, 1], 'ob')
# mark support vectors
supportVectorsIndex = nonzero(svm.alphas.A > 0)[0]
for i in supportVectorsIndex:
plt.plot(svm.train_x[i, 0], svm.train_x[i, 1], 'oy')
# draw the classify line
w = zeros((2, 1))
for i in supportVectorsIndex:
w += multiply(svm.alphas[i] * svm.train_y[i], svm.train_x[i, :].T)
min_x = min(svm.train_x[:, 0])[0, 0]
max_x = max(svm.train_x[:, 0])[0, 0]
y_min_x = float(-svm.b - w[0] * min_x) / w[1]
y_max_x = float(-svm.b - w[0] * max_x) / w[1]
plt.plot([min_x, max_x], [y_min_x, y_max_x], '-g')
plt.show()
測(cè)試的數(shù)據(jù)來(lái)自這里。有100個(gè)樣本,每個(gè)樣本兩維,最后是對(duì)應(yīng)的標(biāo)簽,例如:
3.542485 1.977398 -1
3.018896 2.556416 -1
7.551510 -1.580030 1
2.114999 -0.004466 -1
……
測(cè)試代碼中首先加載這個(gè)數(shù)據(jù)庫(kù),然后用前面80個(gè)樣本來(lái)訓(xùn)練,再用剩下的20個(gè)樣本的測(cè)試,并顯示訓(xùn)練后的模型和分類(lèi)結(jié)果。測(cè)試代碼如下:
test_SVM.py
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
#################################################
# SVM: support vector machine
# Author : zouxy
# Date : 2013-12-12
# HomePage : http://blog.csdn.net/zouxy09
# Email : zouxy09@qq.com
#################################################
from numpy import *
import SVM
################## test svm #####################
## step 1: load data
print "step 1: load data..."
dataSet = []
labels = []
fileIn = open('E:/Python/Machine Learning in Action/testSet.txt')
for line in fileIn.readlines():
lineArr = line.strip().split('\t')
dataSet.append([float(lineArr[0]), float(lineArr[1])])
labels.append(float(lineArr[2]))
dataSet = mat(dataSet)
labels = mat(labels).T
train_x = dataSet[0:81, :]
train_y = labels[0:81, :]
test_x = dataSet[80:101, :]
test_y = labels[80:101, :]
## step 2: training...
print "step 2: training..."
C = 0.6
toler = 0.001
maxIter = 50
svmClassifier = SVM.trainSVM(train_x, train_y, C, toler, maxIter, kernelOption = ('linear', 0))
## step 3: testing
print "step 3: testing..."
accuracy = SVM.testSVM(svmClassifier, test_x, test_y)
## step 4: show the result
print "step 4: show the result..."
print 'The classify accuracy is: %.3f%%' % (accuracy * 100)
SVM.showSVM(svmClassifier)
運(yùn)行結(jié)果如下:
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
step 1: load data... 數(shù)據(jù)分析師培訓(xùn)
step 2: training...
---iter:0 entire set, alpha pairs changed:8
---iter:1 non boundary, alpha pairs changed:7
---iter:2 non boundary, alpha pairs changed:1
---iter:3 non boundary, alpha pairs changed:0
---iter:4 entire set, alpha pairs changed:0
Congratulations, training complete! Took 0.058000s!
step 3: testing...
step 4: show the result...
The classify accuracy is: 100.000%
訓(xùn)練好的模型圖:
數(shù)據(jù)分析咨詢(xún)請(qǐng)掃描二維碼
若不方便掃碼,搜微信號(hào):CDAshujufenxi
LSTM 模型輸入長(zhǎng)度選擇技巧:提升序列建模效能的關(guān)鍵? 在循環(huán)神經(jīng)網(wǎng)絡(luò)(RNN)家族中,長(zhǎng)短期記憶網(wǎng)絡(luò)(LSTM)憑借其解決長(zhǎng)序列 ...
2025-07-11CDA 數(shù)據(jù)分析師報(bào)考條件詳解與準(zhǔn)備指南? ? 在數(shù)據(jù)驅(qū)動(dòng)決策的時(shí)代浪潮下,CDA 數(shù)據(jù)分析師認(rèn)證愈發(fā)受到矚目,成為眾多有志投身數(shù) ...
2025-07-11數(shù)據(jù)透視表中兩列相乘合計(jì)的實(shí)用指南? 在數(shù)據(jù)分析的日常工作中,數(shù)據(jù)透視表憑借其強(qiáng)大的數(shù)據(jù)匯總和分析功能,成為了 Excel 用戶(hù) ...
2025-07-11尊敬的考生: 您好! 我們誠(chéng)摯通知您,CDA Level I和 Level II考試大綱將于 2025年7月25日 實(shí)施重大更新。 此次更新旨在確保認(rèn) ...
2025-07-10BI 大數(shù)據(jù)分析師:連接數(shù)據(jù)與業(yè)務(wù)的價(jià)值轉(zhuǎn)化者? ? 在大數(shù)據(jù)與商業(yè)智能(Business Intelligence,簡(jiǎn)稱(chēng) BI)深度融合的時(shí)代,BI ...
2025-07-10SQL 在預(yù)測(cè)分析中的應(yīng)用:從數(shù)據(jù)查詢(xún)到趨勢(shì)預(yù)判? ? 在數(shù)據(jù)驅(qū)動(dòng)決策的時(shí)代,預(yù)測(cè)分析作為挖掘數(shù)據(jù)潛在價(jià)值的核心手段,正被廣泛 ...
2025-07-10數(shù)據(jù)查詢(xún)結(jié)束后:分析師的收尾工作與價(jià)值深化? ? 在數(shù)據(jù)分析的全流程中,“query end”(查詢(xún)結(jié)束)并非工作的終點(diǎn),而是將數(shù) ...
2025-07-10CDA 數(shù)據(jù)分析師考試:從報(bào)考到取證的全攻略? 在數(shù)字經(jīng)濟(jì)蓬勃發(fā)展的今天,數(shù)據(jù)分析師已成為各行業(yè)爭(zhēng)搶的核心人才,而 CDA(Certi ...
2025-07-09【CDA干貨】單樣本趨勢(shì)性檢驗(yàn):捕捉數(shù)據(jù)背后的時(shí)間軌跡? 在數(shù)據(jù)分析的版圖中,單樣本趨勢(shì)性檢驗(yàn)如同一位耐心的偵探,專(zhuān)注于從單 ...
2025-07-09year_month數(shù)據(jù)類(lèi)型:時(shí)間維度的精準(zhǔn)切片? ? 在數(shù)據(jù)的世界里,時(shí)間是最不可或缺的維度之一,而year_month數(shù)據(jù)類(lèi)型就像一把精準(zhǔn) ...
2025-07-09CDA 備考干貨:Python 在數(shù)據(jù)分析中的核心應(yīng)用與實(shí)戰(zhàn)技巧? ? 在 CDA 數(shù)據(jù)分析師認(rèn)證考試中,Python 作為數(shù)據(jù)處理與分析的核心 ...
2025-07-08SPSS 中的 Mann-Kendall 檢驗(yàn):數(shù)據(jù)趨勢(shì)與突變分析的有力工具? ? ? 在數(shù)據(jù)分析的廣袤領(lǐng)域中,準(zhǔn)確捕捉數(shù)據(jù)的趨勢(shì)變化以及識(shí)別 ...
2025-07-08備戰(zhàn) CDA 數(shù)據(jù)分析師考試:需要多久?如何規(guī)劃? CDA(Certified Data Analyst)數(shù)據(jù)分析師認(rèn)證作為國(guó)內(nèi)權(quán)威的數(shù)據(jù)分析能力認(rèn)證 ...
2025-07-08LSTM 輸出不確定的成因、影響與應(yīng)對(duì)策略? 長(zhǎng)短期記憶網(wǎng)絡(luò)(LSTM)作為循環(huán)神經(jīng)網(wǎng)絡(luò)(RNN)的一種變體,憑借獨(dú)特的門(mén)控機(jī)制,在 ...
2025-07-07統(tǒng)計(jì)學(xué)方法在市場(chǎng)調(diào)研數(shù)據(jù)中的深度應(yīng)用? 市場(chǎng)調(diào)研是企業(yè)洞察市場(chǎng)動(dòng)態(tài)、了解消費(fèi)者需求的重要途徑,而統(tǒng)計(jì)學(xué)方法則是市場(chǎng)調(diào)研數(shù) ...
2025-07-07CDA數(shù)據(jù)分析師證書(shū)考試全攻略? 在數(shù)字化浪潮席卷全球的當(dāng)下,數(shù)據(jù)已成為企業(yè)決策、行業(yè)發(fā)展的核心驅(qū)動(dòng)力,數(shù)據(jù)分析師也因此成為 ...
2025-07-07剖析 CDA 數(shù)據(jù)分析師考試題型:解鎖高效備考與答題策略? CDA(Certified Data Analyst)數(shù)據(jù)分析師考試作為衡量數(shù)據(jù)專(zhuān)業(yè)能力的 ...
2025-07-04SQL Server 字符串截取轉(zhuǎn)日期:解鎖數(shù)據(jù)處理的關(guān)鍵技能? 在數(shù)據(jù)處理與分析工作中,數(shù)據(jù)格式的規(guī)范性是保證后續(xù)分析準(zhǔn)確性的基礎(chǔ) ...
2025-07-04CDA 數(shù)據(jù)分析師視角:從數(shù)據(jù)迷霧中探尋商業(yè)真相? 在數(shù)字化浪潮席卷全球的今天,數(shù)據(jù)已成為企業(yè)決策的核心驅(qū)動(dòng)力,CDA(Certifie ...
2025-07-04CDA 數(shù)據(jù)分析師:開(kāi)啟數(shù)據(jù)職業(yè)發(fā)展新征程? ? 在數(shù)據(jù)成為核心生產(chǎn)要素的今天,數(shù)據(jù)分析師的職業(yè)價(jià)值愈發(fā)凸顯。CDA(Certified D ...
2025-07-03