標籤：

Python3《機器學習實戰》學習筆記（八）：支持向量機原理篇之手撕線性SVM

01-24

轉載請註明作者和出處： zhuanlan.zhihu.com/ml-j

機器學習知乎專欄：zhuanlan.zhihu.com/ml-j

CSDN博客專欄：blog.csdn.net/column/de

Github代碼獲取：github.com/Jack-Cherish

Python版本： Python3.x

運行平台： Windows

IDE： Sublime text3

一、前言

說來慚愧，斷更快半個月了，本打算是一周一篇的。感覺SVM瞬間難了不少，推導耗費了很多時間，同時身邊的事情也不少，忙了許久。本篇文章參考了諸多大牛的文章寫成的，對於什麼是SVM做出了生動的闡述，同時也進行了線性SVM的理論推導，以及最後的編程實踐，公式較多，還需靜下心來一點一點推導。

本文出現的所有代碼，均可在我的github上下載，歡迎Follow、Star：github.com/Jack-Cherish

如果對於代碼理解不夠的，可以結合本文，觀看由南京航空航天大學碩士：深度眸，為大家免費錄製的視頻：

  • SVM視頻講解：pan.baidu.com/s/1pL7gt3 密碼：80ge
  • 視頻交流群：ML與DL視頻分享群(678455658)，歡迎提出寶貴意見。

二、什麼是SVM？

SVM的英文全稱是Support Vector Machines，我們叫它支持向量機。支持向量機是我們用於分類的一種演算法。讓我們以一個小故事的形式，開啟我們的SVM之旅吧。

在很久以前的情人節，一位大俠要去救他的愛人，但天空中的魔鬼和他玩了一個遊戲。

魔鬼在桌子上似乎有規律放了兩種顏色的球，說："你用一根棍分開它們？要求：盡量在放更多球之後，仍然適用。"

於是大俠這樣放，乾的不錯？

然後魔鬼，又在桌上放了更多的球，似乎有一個球站錯了陣營。顯然，大俠需要對棍做出調整。

SVM就是試圖把棍放在最佳位置，好讓在棍的兩邊有儘可能大的間隙。這個間隙就是球到棍的距離。

現在好了，即使魔鬼放了更多的球，棍仍然是一個好的分界線。

魔鬼看到大俠已經學會了一個trick(方法、招式)，於是魔鬼給了大俠一個新的挑戰。

現在，大俠沒有棍可以很好幫他分開兩種球了，現在怎麼辦呢？當然像所有武俠片中一樣大俠桌子一拍，球飛到空中。然後，憑藉大俠的輕功，大俠抓起一張紙，插到了兩種球的中間。

現在，從空中的魔鬼的角度看這些球，這些球看起來像是被一條曲線分開了。

再之後，無聊的大人們，把這些球叫做data，把棍子叫做classifier, 找到最大間隙的trick叫做optimization，拍桌子叫做kernelling, 那張紙叫做hyperplane。

更為直觀地感受一下吧(需要翻牆)：youtube.com/watch?

概述一下：

當一個分類問題，數據是線性可分的，也就是用一根棍就可以將兩種小球分開的時候，我們只要將棍的位置放在讓小球距離棍的距離最大化的位置即可，尋找這個最大間隔的過程，就叫做最優化。但是，現實往往是很殘酷的，一般的數據是線性不可分的，也就是找不到一個棍將兩種小球很好的分類。這個時候，我們就需要像大俠一樣，將小球拍起，用一張紙代替小棍將小球進行分類。想要讓數據飛起，我們需要的東西就是核函數(kernel)，用於切分小球的紙，就是超平面。

也許這個時候，你還是似懂非懂，沒關係。根據剛才的描述，可以看出，問題是從線性可分延伸到線性不可分的。那麼，我們就按照這個思路，進行原理性的剖析。

三、線性SVM

先看下線性可分的二分類問題。

上圖中的(a)是已有的數據，紅色和藍色分別代表兩個不同的類別。數據顯然是線性可分的，但是將兩類數據點分開的直線顯然不止一條。上圖的(b)和(c)分別給出了B、C兩種不同的分類方案，其中黑色實線為分界線，術語稱為「決策面」。每個決策面對應了一個線性分類器。雖然從分類結果上看，分類器A和分類器B的效果是相同的。但是他們的性能是有差距的，看下圖：

在"決策面"不變的情況下，我又添加了一個紅點。可以看到，分類器B依然能很好的分類結果，而分類器C則出現了分類錯誤。顯然分類器B的"決策面"放置的位置優於分類器C的"決策面"放置的位置，SVM演算法也是這麼認為的，它的依據就是分類器B的分類間隔比分類器C的分類間隔大。這裡涉及到第一個SVM獨有的概念"分類間隔"。在保證決策面方向不變且不會出現錯分樣本的情況下移動決策面，會在原來的決策面兩側找到兩個極限位置（越過該位置就會產生錯分現象），如虛線所示。虛線的位置由決策面的方向和距離原決策面最近的幾個樣本的位置決定。而這兩條平行虛線正中間的分界線就是在保持當前決策面方向不變的前提下的最優決策面。兩條虛線之間的垂直距離就是這個最優決策面對應的分類間隔。顯然每一個可能把數據集正確分開的方向都有一個最優決策面（有些方向無論如何移動決策面的位置也不可能將兩類樣本完全分開），而不同方向的最優決策面的分類間隔通常是不同的，那個具有「最大間隔」的決策面就是SVM要尋找的最優解。而這個真正的最優解對應的兩側虛線所穿過的樣本點，就是SVM中的支持樣本點，稱為"支持向量"。

1、數學建模

求解這個"決策面"的過程，就是最優化。一個最優化問題通常有兩個基本的因素：1）目標函數，也就是你希望什麼東西的什麼指標達到最好；2）優化對象，你期望通過改變哪些因素來使你的目標函數達到最優。在線性SVM演算法中，目標函數顯然就是那個"分類間隔"，而優化對象則是決策面。所以要對SVM問題進行數學建模，首先要對上述兩個對象（"分類間隔"和"決策面"）進行數學描述。按照一般的思維習慣，我們先描述決策面。

數學建模的時候，先在二維空間建模，然後再推廣到多維。

（1）"決策面"方程

我們都知道二維空間下一條直線的方式如下所示：

現在我們做個小小的改變，讓原來的x軸變成x1，y軸變成x2。

移項得：

將公式向量化得：

進一步向量化，用w列向量和x列向量和標量γ進一步向量化：

其中，向量w和x分別為：

這裡w1=a，w2=-1。我們都知道，最初的那個直線方程a和b的幾何意義，a表示直線的斜率，b表示截距，a決定了直線與x軸正方向的夾角，b決定了直線與y軸交點位置。那麼向量化後的直線的w和r的幾何意義是什麼呢？

現在假設：

可得：

在坐標軸上畫出直線和向量w：

藍色的線代表向量w，紅色的先代表直線y。我們可以看到向量w和直線的關係為垂直關係。這說明了向量w也控制這直線的方向，只不過是與這個直線的方向是垂直的。標量γ的作用也沒有變，依然決定了直線的截距。此時，我們稱w為直線的法向量。

二維空間的直線方程已經推導完成，將其推廣到n為空間，就變成了超平面方程。(一個超平面，在二維空間的例子就是一個直線)但是它的公式沒變，依然是：

不同之處在於：

我們已經順利推導出了"決策面"方程，它就是我們的超平面方程，之後，我們統稱其為超平面方程。

（2）"分類間隔"方程

現在，我們依然對於一個二維平面的簡單例子進行推導。

我們已經知道間隔的大小實際上就是支持向量對應的樣本點到決策面的距離的二倍。那麼圖中的距離d我們怎麼求？我們高中都學過，點到直線的距離距離公式如下：

公式中的直線方程為Ax0+By0+C=0，點P的坐標為(x0,y0)。

現在，將直線方程擴展到多維，求得我們現在的超平面方程，對公式進行如下變形：

這個d就是"分類間隔"。其中||w||表示w的二範數，求所有元素的平方和，然後再開方。比如對於二維平面：

那麼，

我們目的是為了找出一個分類效果好的超平面作為分類器。分類器的好壞的評定依據是分類間隔W=2d的大小，即分類間隔W越大，我們認為這個超平面的分類效果越好。此時，求解超平面的問題就變成了求解分類間隔W最大化的為題。W的最大化也就是d最大化的。

（3）約束條件

看起來，我們已經順利獲得了目標函數的數學形式。但是為了求解w的最大值。我們不得不面對如下問題：

  • 我們如何判斷超平面是否將樣本點正確分類？
  • 我們知道要求距離d的最大值，我們首先需要找到支持向量上的點，怎麼在眾多的點中選出支持向量上的點呢？

上述我們需要面對的問題就是約束條件，也就是說我們優化的變數d的取值範圍受到了限制和約束。事實上約束條件一直是最優化問題里最讓人頭疼的東西。但既然我們已經知道了這些約束條件確實存在，就不得不用數學語言對他們進行描述。但SVM演算法通過一些巧妙的小技巧，將這些約束條件融合到一個不等式裡面。

這個二維平面上有兩種點，我們分別對它們進行標記：

  • 紅顏色的圓點標記為1，我們人為規定其為正樣本；
  • 藍顏色的五角星標記為-1，我們人為規定其為負樣本。

對每個樣本點xi加上一個類別標籤yi：

如果我們的超平面方程能夠完全正確地對上圖的樣本點進行分類，就會滿足下面的方程：

如果我們要求再高一點，假設決策面正好處於間隔區域的中軸線上，並且相應的支持向量對應的樣本點到決策面的距離為d，那麼公式進一步寫成：

上述公式的解釋就是，對於所有分類標籤為1和-1樣本點，它們到直線的距離都大於等於d(支持向量上的樣本點到超平面的距離)。公式兩邊都除以d，就可以得到：

其中，

因為||w||和d都是標量。所上述公式的兩個矢量，依然描述一條直線的法向量和截距。

上述兩個公式，都是描述一條直線，數學模型代表的意義是一樣的。現在，讓我們對wd和γd重新起個名字，就叫它們w和γ。因此，我們就可以說："對於存在分類間隔的兩類樣本點，我們一定可以找到一些超平面面，使其對於所有的樣本點均滿足下面的條件："

上述方程即給出了SVM最優化問題的約束條件。這時候，可能有人會問了，為什麼標記為1和-1呢？因為這樣標記方便我們將上述方程變成如下形式：

正是因為標籤為1和-1，才方便我們將約束條件變成一個約束方程，從而方便我們的計算。

（4）線性SVM優化問題基本描述

現在整合一下思路，我們已經得到我們的目標函數：

我們的優化目標是是d最大化。我們已經說過，我們是用支持向量上的樣本點求解d的最大化的問題的。那麼支持向量上的樣本點有什麼特點呢？

你贊同這個觀點嗎？所有支持向量上的樣本點，都滿足如上公式。如果不贊同，請重看"分類間隔"方程推導過程。

現在我們就可以將我們的目標函數進一步化簡：

因為，我們只關心支持向量上的點。隨後我們求解d的最大化問題變成了||w||的最小化問題。進而||w||的最小化問題等效於

為什麼要做這樣的等效呢？這是為了在進行最優化的過程中對目標函數求導時比較方便，但這絕對不影響最優化問題最後的求解。我們將最終的目標函數和約束條件放在一起進行描述：

這裡n是樣本點的總個數，縮寫s.t.表示"Subject to"，是"服從某某條件"的意思。上述公式描述的是一個典型的不等式約束條件下的二次型函數優化問題，同時也是支持向量機的基本數學模型。

（5）求解準備

我們已經得到支持向量機的基本數學模型，接下來的問題就是如何根據數學模型，求得我們想要的最優解。在學習求解方法之前，我們得知道一點，想用我下面講述的求解方法有一個前提，就是我們的目標函數必須是凸函數。理解凸函數，我們還要先明確另一個概念，凸集。在凸幾何中，凸集(convex set)是在)凸組合下閉合的放射空間的子集。看一幅圖可能更容易理解：

左右量圖都是一個集合。如果集合中任意2個元素連線上的點也在集合中，那麼這個集合就是凸集。顯然，上圖中的左圖是一個凸集，上圖中的右圖是一個非凸集。

凸函數的定義也是如此，其幾何意義表示為函數任意兩點連線上的值大於對應自變數處的函數值。若這裡凸集C即某個區間L，那麼，設函數f為定義在區間L上的函數，若對L上的任意兩點x1，x2和任意的實數λ，λ屬於(0,1)，總有：

則函數f稱為L上的凸函數，當且僅當其上鏡圖（在函數圖像上方的點集）為一個凸集。再看一幅圖，也許更容易理解：

像上圖這樣的函數，它整體就是一個非凸函數，我們無法獲得全局最優解的，只能獲得局部最優解。比如紅框內的部分，如果單獨拿出來，它就是一個凸函數。對於我們的目標函數：

很顯然，它是一個凸函數。所以，可以使用我接下來講述的方法求取最優解。

通常我們需要求解的最優化問題有如下幾類：

  • 無約束優化問題，可以寫為：

  • 有等式約束的優化問題，可以寫為：

  • 有不等式約束的優化問題，可以寫為：

對於第(a)類的優化問題，嘗嘗使用的方法就是費馬大定理(Fermat)，即使用求取函數f(x)的導數，然後令其為零，可以求得候選最優值，再在這些候選值中驗證；如果是凸函數，可以保證是最優解。這也就是我們高中經常使用的求函數的極值的方法。

對於第(b)類的優化問題，常常使用的方法就是拉格朗日乘子法（Lagrange Multiplier) ，即把等式約束h_i(x)用一個係數與f(x)寫為一個式子，稱為拉格朗日函數，而係數稱為拉格朗日乘子。通過拉格朗日函數對各個變數求導，令其為零，可以求得候選值集合，然後驗證求得最優值。

對於第(c)類的優化問題，常常使用的方法就是KKT條件。同樣地，我們把所有的等式、不等式約束與f(x)寫為一個式子，也叫拉格朗日函數，係數也稱拉格朗日乘子，通過一些條件，可以求出最優值的必要條件，這個條件稱為KKT條件。

必要條件和充要條件如果不理解，可以看下面這句話：

  • A的必要條件就是A可以推出的結論
  • A的充分條件就是可以推出A的前提

了解到這些，現在讓我們再看一下我們的最優化問題：

現在，我們的這個對優化問題屬於哪一類？很顯然，它屬於第(c)類問題。在學習求解最優化問題之前，我們還要學習兩個東西：拉格朗日函數和KKT條件。

（6）拉格朗日函數

首先，我們先要從宏觀的視野上了解一下拉格朗日對偶問題出現的原因和背景。

我們知道我們要求解的是最小化問題，所以一個直觀的想法是如果我能夠構造一個函數，使得該函數在可行解區域內與原目標函數完全一致，而在可行解區域外的數值非常大，甚至是無窮大，那麼這個沒有約束條件的新目標函數的優化問題就與原來有約束條件的原始目標函數的優化問題是等價的問題。這就是使用拉格朗日方程的目的，它將約束條件放到目標函數中，從而將有約束優化問題轉換為無約束優化問題。

隨後，人們又發現，使用拉格朗日獲得的函數，使用求導的方法求解依然困難。進而，需要對問題再進行一次轉換，即使用一個數學技巧：拉格朗日對偶。

所以，顯而易見的是，我們在拉格朗日優化我們的問題這個道路上，需要進行下面二個步驟：

  • 將有約束的原始目標函數轉換為無約束的新構造的拉格朗日目標函數
  • 使用拉格朗日對偶性，將不易求解的優化問題轉化為易求解的優化

下面，進行第一步：將有約束的原始目標函數轉換為無約束的新構造的拉格朗日目標函數

公式變形如下：

其中αi是拉格朗日乘子，αi大於等於0，是我們構造新目標函數時引入的係數變數(我們自己設置)。現在我們令：

當樣本點不滿足約束條件時，即在可行解區域外：

此時，我們將αi設置為正無窮，此時θ(w)顯然也是正無窮。

當樣本點滿足約束條件時，即在可行解區域內：

此時，顯然θ(w)為原目標函數本身。我們將上述兩種情況結合一下，就得到了新的目標函數：

此時，再看我們的初衷，就是為了建立一個在可行解區域內與原目標函數相同，在可行解區域外函數值趨近於無窮大的新函數，現在我們做到了。

現在，我們的問題變成了求新目標函數的最小值，即：

這裡用p*表示這個問題的最優值，且和最初的問題是等價的。

接下來，我們進行第二步：將不易求解的優化問題轉化為易求解的優化

我們看一下我們的新目標函數，先求最大值，再求最小值。這樣的話，我們首先就要面對帶有需要求解的參數w和b的方程，而αi又是不等式約束，這個求解過程不好做。所以，我們需要使用拉格朗日函數對偶性，將最小和最大的位置交換一下，這樣就變成了：

交換以後的新問題是原始問題的對偶問題，這個新問題的最優值用d*來表示。而且d*<=p*。我們關心的是d=p的時候，這才是我們要的解。需要什麼條件才能讓d=p呢？

  • 首先必須滿足這個優化問題是凸優化問題。
  • 其次，需要滿足KKT條件。

凸優化問題的定義是：求取最小值的目標函數為凸函數的一類優化問題。目標函數是凸函數我們已經知道，這個優化問題又是求最小值。所以我們的最優化問題就是凸優化問題。

接下里，就是探討是否滿足KKT條件了。

（7）KKT條件

我們已經使用拉格朗日函數對我們的目標函數進行了處理，生成了一個新的目標函數。通過一些條件，可以求出最優值的必要條件，這個條件就是接下來要說的KKT條件。一個最優化模型能夠表示成下列標準形式：

KKT條件的全稱是Karush-Kuhn-Tucker條件，KKT條件是說最優值條件必須滿足以下條件：

  • 條件一：經過拉格朗日函數處理之後的新目標函數L(w,b,α)對x求導為零：
  • 條件二： h_j(x)=0；
  • 條件三： alpha*g_k(k)=0 ；

對於我們的優化問題：

顯然，條件二已經滿足了。另外兩個條件為啥也滿足呢？

這裡原諒我省略一系列證明步驟，感興趣的可以移步這裡：blog.csdn.net/xianlingm

這裡已經給出了很好的解釋。現在，凸優化問題和KKT都滿足了，問題轉換成了對偶問題。而求解這個對偶學習問題，可以分為三個步驟：首先要讓L(w,b,α)關於w和b最小化，然後求對α的極大，最後利用SMO演算法求解對偶問題中的拉格朗日乘子。現在，我們繼續推導。

（8）對偶問題求解

第一步：

根據上述推導已知：

首先固定α，要讓L(w,b,α)關於w和b最小化，我們分別對w和b偏導數，令其等於0，即：

將上述結果帶回L(w,b,α)得到：

從上面的最後一個式子，我們可以看出，此時的L(w,b,α)函數只含有一個變數，即αi。

第二步：

現在內側的最小值求解完成，我們求解外側的最大值，從上面的式子得到

現在我們的優化問題變成了如上的形式。對於這個問題，我們有更高效的優化演算法，即序列最小優化（SMO）演算法。我們通過這個優化演算法能得到α，再根據α，我們就可以求解出w和b，進而求得我們最初的目的：找到超平面，即"決策平面"。

總結一句話：我們為啥使出吃奶的勁兒進行推導？因為我們要將最初的原始問題，轉換到可以使用SMO演算法求解的問題，這是一種最流行的求解方法。為啥用這種求解方法？因為它牛逼啊！

對於上述問題，我們就可以使用SMO演算法進行求解了，但是，SMO演算法又是什麼呢？

2、SMO演算法

現在，我們已經得到了可以用SMO演算法求解的目標函數，但是對於怎麼編程實現SMO演算法還是感覺無從下手。那麼現在就聊聊如何使用SMO演算法進行求解。

（1）Platt的SMO演算法

1996年，John Platt發布了一個稱為SMO的強大演算法，用於訓練SVM。SM表示序列最小化(Sequential Minimal Optimizaion)。Platt的SMO演算法是將大優化問題分解為多個小優化問題來求解的。這些小優化問題往往很容易求解，並且對它們進行順序求解的結果與將它們作為整體來求解的結果完全一致的。在結果完全相同的同時，SMO演算法的求解時間短很多。

SMO演算法的目標是求出一系列alpha和b，一旦求出了這些alpha，就很容易計算出權重向量w並得到分隔超平面。

SMO演算法的工作原理是：每次循環中選擇兩個alpha進行優化處理。一旦找到了一對合適的alpha，那麼就增大其中一個同時減小另一個。這裡所謂的"合適"就是指兩個alpha必須符合以下兩個條件，條件之一就是兩個alpha必須要在間隔邊界之外，而且第二個條件則是這兩個alpha還沒有進進行過區間化處理或者不在邊界上。

（2）SMO演算法的解法

先來定義特徵到結果的輸出函數為：

接著，我們回憶一下原始優化問題，如下：

求導得：

將上述公式帶入輸出函數中：

與此同時，拉格朗日對偶後得到最終的目標化函數：

將目標函數變形，在前面增加一個符號，將最大值問題轉換成最小值問題：

實際上，對於上述目標函數，是存在一個假設的，即數據100%線性可分。但是，目前為止，我們知道幾乎所有數據都不那麼"乾淨"。這時我們就可以通過引入所謂的鬆弛變數(slack variable)，來允許有些數據點可以處於超平面的錯誤的一側。這樣我們的優化目標就能保持仍然不變，但是此時我們的約束條件有所改變：

根據KKT條件可以得出其中αi取值的意義為：

  • 對於第1種情況，表明αi是正常分類，在邊界內部；
  • 對於第2種情況，表明αi是支持向量，在邊界上；
  • 對於第3種情況，表明αi是在兩條邊界之間。

而最優解需要滿足KKT條件，即上述3個條件都得滿足，以下幾種情況出現將會不滿足：

也就是說，如果存在不能滿足KKT條件的αi，那麼需要更新這些αi，這是第一個約束條件。此外，更新的同時還要受到第二個約束條件的限制，即：

因為這個條件，我們同時更新兩個α值，因為只有成對更新，才能保證更新之後的值仍然滿足和為0的約束，假設我們選擇的兩個乘子為α1和α2：

其中， ksi為常數。因為兩個因子不好同時求解，所以可以先求第二個乘子α2的解（α2 new），得到α2的解（α2 new）之後，再用α2的解（α2 new）表示α1的解（α1 new ）。為了求解α2 new ，得先確定α2 new的取值範圍。假設它的上下邊界分別為H和L，那麼有：

接下來，綜合下面兩個條件：

當y1不等於y2時，即一個為正1，一個為負1的時候，可以得到：

所以有：

此時，取值範圍如下圖所示：

當y1等於y2時，即兩個都為正1或者都為負1，可以得到：

所以有：

此時，取值範圍如下圖所示：

如此，根據y1和y2異號或同號，可以得出α2 new的上下界分別為：

這個界限就是編程的時候需要用到的。已經確定了邊界，接下來，就是推導迭代式，用於更新 α值。

我們已經知道，更新α的邊界，接下來就是討論如何更新α值。我們依然假設選擇的兩個乘子為α1和α2。固定這兩個乘子，進行推導。於是目標函數變成了：

為了描述方便，我們定義如下符號：

最終目標函數變為：

我們不關心constant的部分，因為對於α1和α2來說，它們都是常數項，在求導的時候，直接變為0。對於這個目標函數，如果對其求導，還有個未知數α1，所以要推導出α1和α2的關係，然後用α2代替α1，這樣目標函數就剩一個未知數了，我們就可以求導了，推導出迭代公式。所以現在繼續推導α1和α2的關係。注意第一個約束條件：

我們在求α1和α2的時候，可以將α3,α4,...,αn和y3,y4,...,yn看作常數項。因此有：

我們不必關心常數B的大小，現在將上述等式兩邊同時乘以y1，得到(y1y1=1)：

其中γ為常數By1，我們不關心這個值，s=y1y2。接下來，我們將得到的α1帶入W(α2)公式得：

這樣目標函數中就只剩下α2了，我們對其求偏導（注意：s=y1y2，所以s的平方為1，y1的平方和y2的平方均為1）：

繼續化簡，將s=y1y2帶入方程。

我們令：

Ei為誤差項，η為學習速率。

再根據我們已知的公式：

將α2 new繼續化簡得：

這樣，我們就得到了最終需要的迭代公式。這個是沒有經過剪輯是的解，需要考慮約束：

根據之前推導的α取值範圍，我們得到最終的解析解為：

又因為：

消去γ得：

這樣，我們就知道了怎樣計算α1和α2了，也就是如何對選擇的α進行更新。

當我們更新了α1和α2之後，需要重新計算閾值b，因為b關係到了我們f(x)的計算，也就關係到了誤差Ei的計算。

我們要根據α的取值範圍，去更正b的值，使間隔最大化。當α1 new在0和C之間的時候，根據KKT條件可知，這個點是支持向量上的點。因此，滿足下列公式：

公式兩邊同時乘以y1得(y1y1=1)：

因為我們是根據α1和α2的值去更新b，所以單獨提出i=1和i=2的時候，整理可得：

其中前兩項為：

將上述兩個公式，整理得：

同理可得b2 new為：

當b1和b2都有效的時候，它們是相等的，即：

當兩個乘子都在邊界上，則b閾值和KKT條件一致。當不滿足的時候，SMO演算法選擇他們的中點作為新的閾值：

最後，更新所有的α和b，這樣模型就出來了，從而即可求出我們的分類函數。

現在，讓我們梳理下SMO演算法的步驟：

  • 步驟1：計算誤差：

  • 步驟2：計算上下界L和H：

  • 步驟3：計算η：

  • 步驟4：更新αj：

  • 步驟5：根據取值範圍修剪αj：

  • 步驟6：更新αi：

  • 步驟7：更新b1和b2：

  • 步驟8：根據b1和b2更新b：

四、編程求解線性SVM

已經梳理完了SMO演算法實現步驟，接下來按照這個思路編寫代碼，進行實戰練習。

（1）可視化數據集

我們先使用簡單的數據集進行測試，數據集下載地址：github.com/Jack-Cherish

編寫程序可視化數據集，看下它是長什麼樣的：

# -*- coding:UTF-8 -*-import matplotlib.pyplot as pltimport numpy as np"""函數說明:讀取數據Parameters: fileName - 文件名Returns: dataMat - 數據矩陣 labelMat - 數據標籤Author: Jack CuiBlog: http://blog.csdn.net/c406495762Zhihu: https://www.zhihu.com/people/Jack--Cui/Modify: 2017-09-21"""def loadDataSet(fileName): dataMat = []; labelMat = [] fr = open(fileName) for line in fr.readlines(): #逐行讀取，濾除空格等 lineArr = line.strip().split(" ") dataMat.append([float(lineArr[0]), float(lineArr[1])]) #添加數據 labelMat.append(float(lineArr[2])) #添加標籤 return dataMat,labelMat"""函數說明:數據可視化Parameters: dataMat - 數據矩陣 labelMat - 數據標籤Returns: 無Author: Jack CuiBlog: http://blog.csdn.net/c406495762Zhihu: https://www.zhihu.com/people/Jack--Cui/Modify: 2017-09-21"""def showDataSet(dataMat, labelMat): data_plus = [] #正樣本 data_minus = [] #負樣本 for i in range(len(dataMat)): if labelMat[i] > 0: data_plus.append(dataMat[i]) else: data_minus.append(dataMat[i]) data_plus_np = np.array(data_plus) #轉換為numpy矩陣 data_minus_np = np.array(data_minus) #轉換為numpy矩陣 plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1]) #正樣本散點圖 plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1]) #負樣本散點圖 plt.show()if __name__ == "__main__": dataMat, labelMat = loadDataSet("testSet.txt") showDataSet(dataMat, labelMat)

運行程序，查看結果：

這就是我們使用的二維數據集，顯然線性可分。現在我們使用簡化版的SMO演算法進行求解。

（2）簡化版SMO演算法

按照上述已經推導的步驟編寫代碼：

# -*- coding:UTF-8 -*-from time import sleepimport matplotlib.pyplot as pltimport numpy as npimport randomimport types"""函數說明:讀取數據Parameters: fileName - 文件名Returns: dataMat - 數據矩陣 labelMat - 數據標籤Author: Jack CuiBlog: http://blog.csdn.net/c406495762Zhihu: https://www.zhihu.com/people/Jack--Cui/Modify: 2017-09-21"""def loadDataSet(fileName): dataMat = []; labelMat = [] fr = open(fileName) for line in fr.readlines(): #逐行讀取，濾除空格等 lineArr = line.strip().split(" ") dataMat.append([float(lineArr[0]), float(lineArr[1])]) #添加數據 labelMat.append(float(lineArr[2])) #添加標籤 return dataMat,labelMat"""函數說明:隨機選擇alphaParameters: i - alpha m - alpha參數個數Returns: j -Author: Jack CuiBlog: http://blog.csdn.net/c406495762Zhihu: https://www.zhihu.com/people/Jack--Cui/Modify: 2017-09-21"""def selectJrand(i, m): j = i #選擇一個不等於i的j while (j == i): j = int(random.uniform(0, m)) return j"""函數說明:修剪alphaParameters: aj - alpha值 H - alpha上限 L - alpha下限Returns: aj - alpah值Author: Jack CuiBlog: http://blog.csdn.net/c406495762Zhihu: https://www.zhihu.com/people/Jack--Cui/Modify: 2017-09-21"""def clipAlpha(aj,H,L): if aj > H: aj = H if L > aj: aj = L return aj"""函數說明:簡化版SMO演算法Parameters: dataMatIn - 數據矩陣 classLabels - 數據標籤 C - 鬆弛變數 toler - 容錯率 maxIter - 最大迭代次數Returns: 無Author: Jack CuiBlog: http://blog.csdn.net/c406495762Zhihu: https://www.zhihu.com/people/Jack--Cui/Modify: 2017-09-23"""def smoSimple(dataMatIn, classLabels, C, toler, maxIter): #轉換為numpy的mat存儲 dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose() #初始化b參數，統計dataMatrix的維度 b = 0; m,n = np.shape(dataMatrix) #初始化alpha參數，設為0 alphas = np.mat(np.zeros((m,1))) #初始化迭代次數 iter_num = 0 #最多迭代matIter次 while (iter_num < maxIter): alphaPairsChanged = 0 for i in range(m): #步驟1：計算誤差Ei fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b Ei = fXi - float(labelMat[i]) #優化alpha，更設定一定的容錯率。 if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)): #隨機選擇另一個與alpha_i成對優化的alpha_j j = selectJrand(i,m) #步驟1：計算誤差Ej fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b Ej = fXj - float(labelMat[j]) #保存更新前的aplpha值，使用深拷貝 alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy(); #步驟2：計算上下界L和H if (labelMat[i] != labelMat[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 L==H: print("L==H"); continue #步驟3：計算eta eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T if eta >= 0: print("eta>=0"); continue #步驟4：更新alpha_j alphas[j] -= labelMat[j]*(Ei - Ej)/eta #步驟5：修剪alpha_j alphas[j] = clipAlpha(alphas[j],H,L) if (abs(alphas[j] - alphaJold) < 0.00001): print("alpha_j變化太小"); continue #步驟6：更新alpha_i alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j]) #步驟7：更新b_1和b_2 b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T #步驟8：根據b_1和b_2更新b 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 #統計優化次數 alphaPairsChanged += 1 #列印統計信息 print("第%d次迭代 樣本:%d, alpha優化次數:%d" % (iter_num,i,alphaPairsChanged)) #更新迭代次數 if (alphaPairsChanged == 0): iter_num += 1 else: iter_num = 0 print("迭代次數: %d" % iter_num) return b,alphas"""函數說明:分類結果可視化Parameters: dataMat - 數據矩陣 w - 直線法向量 b - 直線解決Returns: 無Author: Jack CuiBlog: http://blog.csdn.net/c406495762Zhihu: https://www.zhihu.com/people/Jack--Cui/Modify: 2017-09-23"""def showClassifer(dataMat, w, b): #繪製樣本點 data_plus = [] #正樣本 data_minus = [] #負樣本 for i in range(len(dataMat)): if labelMat[i] > 0: data_plus.append(dataMat[i]) else: data_minus.append(dataMat[i]) data_plus_np = np.array(data_plus) #轉換為numpy矩陣 data_minus_np = np.array(data_minus) #轉換為numpy矩陣 plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7) #正樣本散點圖 plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7) #負樣本散點圖 #繪製直線 x1 = max(dataMat)[0] x2 = min(dataMat)[0] a1, a2 = w b = float(b) a1 = float(a1[0]) a2 = float(a2[0]) y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2 plt.plot([x1, x2], [y1, y2]) #找出支持向量點 for i, alpha in enumerate(alphas): if abs(alpha) > 0: x, y = dataMat[i] plt.scatter([x], [y], s=150, c="none", alpha=0.7, linewidth_=1.5, edgecolor="red") plt.show()"""函數說明:計算wParameters: dataMat - 數據矩陣 labelMat - 數據標籤 alphas - alphas值Returns: 無Author: Jack CuiBlog: http://blog.csdn.net/c406495762Zhihu: https://www.zhihu.com/people/Jack--Cui/Modify: 2017-09-23"""def get_w(dataMat, labelMat, alphas): alphas, dataMat, labelMat = np.array(alphas), np.array(dataMat), np.array(labelMat) w = np.dot((np.tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas) return w.tolist()if __name__ == "__main__": dataMat, labelMat = loadDataSet("testSet.txt") b,alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40) w = get_w(dataMat, labelMat, alphas) showClassifer(dataMat, w, b)

程序運行結果：

其中，中間的藍線為求出來的分類器，用紅圈圈出的點為支持向量點。

五、總結

  • 本文主要進行了線性SVM的推導，並通過編程實現一個簡化版SMO演算法；
  • 本文的簡化版SMO演算法在選取α的時候，沒有選擇啟發式的選擇方法，並且兩個乘子的計算沒有進行優化，所以演算法比較耗時，下一篇文章會講解相應的優化方法；
  • 本文討論的是線性SVM，沒有使用核函數，下一篇文章將會講解如何應用核函數，將SVM應用於非線性數據集；
  • 如有問題，請留言。如有錯誤，還望指正，謝謝！

PS： 如果覺得本篇本章對您有所幫助，歡迎關注、評論、贊！

參考資料：

  • [1] 五歲小孩也能看懂的SVM：zhihu.com/question/2109
  • [2] 五歲小孩也能看懂的SVM ：reddit.com/r/MachineLea
  • [3] pluskid大牛博客：blog.pluskid.org/?
  • [4] 陳東嶽老師文章：zhuanlan.zhihu.com/p/24
  • [5] 深入理解拉格朗日乘子法和KKT條件：blog.csdn.net/xianlingm
  • [6] 充分條件和必要條件：zhihu.com/question/3046
  • [7] 凸函數：zh.wikipedia.org/wiki/%
  • [8]《機器學習實戰》第6章內容。
  • [9] SVM之SMO演算法：cnblogs.com/zangrunqian

推薦閱讀：

※機器學習演算法實踐-SVM核函數和軟間隔
※機器學習演算法實踐-Platt SMO和遺傳演算法優化SVM
※淺談最優化問題的KKT條件

TAG:机器学习 | SVM | Python |