Python實現SVM（Support Vector Machine）

1.SVM概念

支持向量機即 Support Vector Machine，簡稱 SVM 。SVM模型的主要思想是在樣本特征空間上找到最佳的分離超平面（二維是線）使得訓練集上正負樣本間隔最大，這個約束使得在感知機的基礎上保證可以找到一個最好的分割分離超平面（也就是說感知機會有多個解）。SVM是用來解決二分類問題的有監督學習算法，在引入了核方法之后SVM也可以用來解決非線性問題。
一般SVM有下面三種：
一、線性可分支持向量機
硬間隔支持向量機：硬間隔最大化（hard margin maximization），當訓練數據線性可分時，可通過硬間隔最大化學得一個線性可分支持向量機。說的白話一點就是兩團數據分的很開，沒有你中有我我中有你的情況，畫一條線完全可以將數據分開，在數據點不多的情況下我們可以直接通過解析解求得w和b。

[圖1]
二、線性支持向量機
軟間隔支持向量機：軟間隔最大化（soft margin maximization ），當訓練數據近似線性可分時，可通過軟間隔最大化學得一個線性支持向量機。兩團數據分的很開，存在少量樣本點輕微越界的情況，我們可以容忍少部分樣本分錯，保證得到的模型魯棒性最高，泛化能力更強。

[圖2]
三、非線性支持向量機
當訓練數據線性不可分時，可通過核方法以及軟間隔最大化學得一個非線性支持向量機。

[圖3]
SVM 一直被認為是效果最好的現成可用的分類算法之一，即使是現在也常應用在圖像識別、文本識別、生物科學和其他科學領域，在金融量化分類中一直是較為理想的，穩定的分類模型。“現成可用”其實是很重要的，因為一直以來學術界和工業界甚至只是學術界里做理論的和做應用的之間，都有一種“鴻溝”，有些很精妙或者很復雜的算法，在抽象出來的模型里很完美，然而在實際問題上卻顯得很脆弱，效果很差甚至全軍覆沒。而 SVM 則非常穩定，理論嚴謹，可以線性分類，也可以高維度非線性分類，工業應用穩定可靠。

2.線性可分支持向量機

給定訓練樣本集 $D=\{(\vec{x_1},y_1),(\vec{x_2},y_2),\dots,(\vec{x_n},y_n)\},y_i\in\{+1,-1\}$， $i$表示第 $i$個樣本， $n$表示樣本容量。分類學習最基本的想法就是基於訓練集 $D$在特征空間中找到一個最佳划分超平面將正負樣本分開，而SVM算法解決的就是如何找到最佳超平面的問題。超平面可通過如下的線性方程來描述：

$\vec{w}^T\vec{x}+b=0$

其中 $\vec{w}^T$ 表示法向量，決定了超平面的方向； $b$表示偏移量，決定了超平面與原點之間的距離。 對於訓練數據集 $D$假設找到了最佳超平面 $\vec{w^*}\vec{x}+b^*=0$，定義決策分類函數:

$f(\vec{x})=sign(\vec{w^*}\vec{x}+b^*)$

該分類決策函數也稱為線性可分支持向量機，其中 $\vec{w^*}$及 $b^*$為最佳分割超平面所對應的參數， $sign$是取符號的意思，所以SVM模型只輸出正負符號，我們在sklearn等包中的概率等是根據其他算法估算出來的，不是來自SVM算法本身。
在測試時對於線性可分支持向量機可以用一個樣本離划分超平面的距離來表示分類預測的可靠程度，如果樣本離划分超平面越遠則對該樣本的分類越可靠，反之就不那么可靠。

[圖2]，沒錯還是圖2…
空間中超平面可記為 $(\vec{w},b)$，根據點到平面的距離公式，空間中任意點 $\vec{x}$ 到超平面 $(\vec{w},b)$的距離可寫為：

$r=\frac{\vec{w}\vec{x}+b}{||\vec{w}||}$

假設超平面 $(\vec{w},b)$能將訓練樣本正確分類，那么對於正樣本一側的任意一個樣本 $(\vec{x_i},y_i)\in{D}$，應該需要滿足該樣本點 $\vec{x}$到超平面的法向量 $\vec{w^*}$的投影到原點的距離大於一定值 $c$的時候使得該樣本點被預測為正樣本一類，即存在數值 $c$使得當 $\vec{w}^T\vec{x_i}>c$， $y_i=+1$。 $\vec{w}^T\vec{x_i}>c$又可寫為 $\vec{w}^T\vec{x_i}+b>0$。在訓練的時候我們要求限制條件更嚴格點以使最終得到的分類器魯棒性更強，所以我們要求 $\vec{w}^T\vec{x_i}+b>1$。也可以寫為大於其它距離，但都可以通過同比例縮放 $\vec{w}$和 $b$來使得使其變為1，因此為計算方便這里直接選擇1。同樣對於負樣本應該有 $\vec{w}^T\vec{x_i}+b<-1$時 $y_i=-1$。即：

$\begin{cases} \vec{w}^T\vec{x_i}+b\geq+1, & y_i=+1 \\ \vec{w}^T\vec{x_i}+b\leq-1, & y_i=-1 \end{cases}$

亦即：

$y_i(\vec{w}^T\vec{x_i}+b)\geq+1$

如圖2所示，距離最佳超平面 $\vec{w}\vec{x}+b=0$最近的幾個訓練樣本點使上式中的等號成立，它們被稱為“支持向量”（support vector）。記超平面 $\vec{w}\vec{x}+b=+1$和 $\vec{w}\vec{x}+b=-1$之間的距離為 $γ$，該距離又被稱為“間隔”（margin），SVM的核心之一就是想辦法將“間隔” $γ$最大化。下面我們推導一下 $γ$與哪些因素有關：

我們初中就已經學過了兩條平行線間距離的求法，例如 $ax+by=c1$和 $ax+by=c2$，那他們的距離是 $\frac{|c2-c1|}{\sqrt{(a^2+b^2)}}$。注意的是，這里的x和y都表示二維坐標（即代表特征里的x1及x2）。而用w來表示就是 $H1:w_1x_1+w_2x_2=+1$和 $H2:w_1x_1+w_2x_2=-1$，那H1和H2的距離就是 $\frac{|1+1|}{\sqrt{(w_1^2+w_2^2)}}=\frac{2}{||w||}$。也就是w的模的倒數的兩倍。也就是說，我們需要最大化 $\gamma$：

$\gamma=\frac{2}{||\vec{w}||}$

也就是說使兩類樣本距離最大的因素僅僅和最佳超平面的法向量（系數w）有關。 要找到具有“最大間隔”（maximum margin）的最佳超平面，就是找到能滿約束的參數 $\vec{w}$、 $b$使得 $γ$最大，即：

$\begin{cases} \max_{\vec{w},b}\frac{2}{||\vec{w}||} \\ s.t.\quad y_i(\vec{w}^T\vec{x_i}+b)\geq+1, i=1,2,\dots,n \end{cases}$

最大化這個距離等價於最小化 $||\vec{w}||$，我們給最小化公式來個平方，乘以1/2都是等價的:

$\begin{cases} \min_{\vec{w},b}\frac{1}{2}||\vec{w}||^2 \\ s.t.\quad y_i(\vec{w}^T\vec{x_i}+b)\geq+1, i=1,2,\dots,n \end{cases}$

2.1使用最優化包直接求解w和b

那么問題來了，怎么求這個 $\min\frac{1}{2}||\vec{w}||^2$公式的最小化呢？先來個簡單粗暴的方法，針對數據點較少的情況，我們用最優化包scipy.optimize.fmin_l_bfgs_b直接求解（代碼參考來源：余凱_西二旗民工）：

# 定義目標函數，w是未知數據，args是已知數，求w的最優解
def func(w,*args):
    X,Y,c=args
    yp=np.dot(X,w)  # w*x，注意已經將b作為x的第一列進行了計算，w*x就是我們現在預測的y
    idx=np.where(yp*Y<1)[0] # 找到分錯的數據索引位置
    e=yp[idx]-Y[idx]        # y預測值-y真實值，誤差
    cost=np.dot(e,e)+c*np.dot(w,w)  # 平方和損失，c：學習率，加w的二范式懲罰
    grand=2*(np.dot(X[idx].T,e)+c*w)# 梯度下降？？
    return cost,grand
if __name__ == '__main__':
    ##1、SVM直接求參數值#############################################################
    print('\n1、SVM直接求參數值，開始')
    # 生成數據，《統計學習方法》李航，P103，例7.1
    dataSet=np.array([[3,3,1],[4,3,1],[1,1,-1]]) # ,[0,0,-1],[0,1,-1]
    m, n = dataSet.shape
    x=dataSet[:,:-1]
    y=dataSet[:,-1] #.reshape((-1,1))
    # 數據定義
    X=np.append(np.ones([x.shape[0],1]),x,1) # x新增一列全1值，作為截距b
    Y=y
    c=0.001 # 學習率
    w=np.zeros(X.shape[1]) # 初始化一組w系數，全0，也可隨機產生：np.random.rand(X.shape[1])
    # bfgs_b方法求最優化問題
    REF=fmin_l_bfgs_b(func,x0=w,args=(X,Y,c),approx_grad=False) #x0=np.random.rand(X.shape[1]) [0,0,0]
    # 采用scipy.optimize其他包夜可以求得
    REF2=fmin_tnc(func,x0=w,args=(X,Y,c),approx_grad=False)
    # 求得最優化計算后的w
    w=REF[0].round(2)           # 取得w值
    print('w:',w[1:],'b:',w[0]) # 與《統計學習方法》李航，P103，例7.1計算結果一致
    # 畫圖
    plotResult(w)
    print('\n1、SVM直接求參數值，結束')

來看看求解結果，貌似還不錯基本接近李航老師的結果（結果進行了四舍五入）：

[圖4]

2.2使用Dual優化

2.1節中直接使用最優化包求解可能並不高效，迭代過程也不可見，我們可以將該凸二次規划問題通過拉格朗日對偶性來求解。
對於上式 $\min\frac{1}{2}||\vec{w}||^2$，我們引入帶約束的拉格朗日乘子 $α_i≥0$，得到拉格朗日函數：

$L(\vec{w},b,\vec{\alpha})=\frac{1}{2}||\vec{w}||^2-\sum_{i=1}^{n}{\alpha_i(y_i(\vec{w}^Tx_i+b)-1)}$

需要注意公式后半部分求和，它指的是我們有m行數據，就會有m個 $\alpha_i$，個人理解分錯的數據 $(\vec{x_i},y_i)$， ${\alpha_i(y_i(\vec{w}^Tx_i+b)-1)}$值不接近0，分對的數據此項值無限接近0，相當於變相求誤差。

對上式 $L(\vec{w},b,\vec{\alpha})$分別對 $\vec{w}$和 $b$求極值，也就是 $L(\vec{w},b,\vec{\alpha})$對 $\vec{w}$和 $b$的梯度為0： $∂L/∂w=0$和 $∂L/∂b=0$，還需要滿足 $α>=0$。求解這里導數為0的式子可以得到：

$\begin{cases} \vec{w}=\sum_{i=1}^{n}{\alpha_iy_i\vec{x_i}}\\ \sum_{i=1}^{n}{\alpha_iy_i}=0 \end{cases}$

有沒有發現上式非常驚艷，求得 $α$后，帶一下公式 $\vec{w}$就求出來了！將上式帶入 $L(\vec{w},b,\vec{\alpha})$函數替換掉 $\vec{w}$和 $b$，得到以下結果：

$\begin{cases} \max._W(\alpha)={\sum_{i=1}^{n}\alpha_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_jy_iy_j\vec{x_i}^T\vec{x_j}}\\ s.t.\quad \alpha_i\geq0,i=1,2,\dots,n \\ \quad\quad\quad \sum_{i=1}^{n}\alpha_iy_i=0 \end{cases}$

這個就是dual problem（如果我們知道 $α$，我們就知道了 $w$。反過來，如果我們知道w，也可以知道 $α$）。這時候我們就變成了求對 $α$的極大，即是關於對偶變量 $α$的優化問題（沒有了變量 $w，b$，只有 $α$）。當求解得到最優的 $α^*$后，就可以同樣代入到上面的公式，導出 $w^*$和 $b^*$了，最終得出分離超平面和分類決策函數。也就是訓練好了SVM。那來一個新的樣本 $x$后，就可以這樣分類了：

$f(\vec{x})=sign(\vec{w}^T\vec{x}+b)$
$=sign((\sum_{i=1}^{n}{\alpha_iy_i\vec{x_i}})*x+b)$
$=sign(\sum_{i=1}^{n}{\alpha_iy_i}(\vec{x_i}*x)+b)$

在這里，其實很多的 $α_i$都是0，也就是說w只是一些少量樣本的線性加權值。這種“稀疏”的表示實際上看成是KNN的數據壓縮的版本。也就是說，以后新來的要分類的樣本首先根據 $w$和 $b$做一次線性運算，然后看求的結果是大於0還是小於0來判斷正例還是負例。現在有了 $α_i$，我們不需要求出 $w$，只需將新來的樣本和訓練數據中的所有樣本做內積和即可。那有人會說，與前面所有的樣本都做運算是不是太耗時了？其實不然，我們從KKT條件中得到，只有支持向量的 $α_i$不為0，其他情況 $α_i$都是0。因此，我們只需求新來的樣本和支持向量的內積，然后運算即可。這種寫法為下面要提到的核函數（kernel）(就是將 $\vec{x_i}*x$引入核函數)做了很好的鋪墊。

2.3Python實現smoSimple核心代碼

# 輸入變量：x、y、c：常數c、toler：容錯率、maxIter：最大循環次數
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    #dataMatIn, classLabels, C, toler, maxIter=dataArr,lableArr,0.6,0.001,40
    dataMatrix = np.mat(dataMatIn)             # 數據x轉換為matrix類型
    labelMat = np.mat(classLabels).transpose() # 標簽y轉換為matrix類型，轉換為一列
    b = 0                                      # 截距b
    m,n = np.shape(dataMatrix)                 # 數據x行數、列數
    alphas = np.mat(np.zeros((m,1)))           # 初始化alpha，有多少行數據就產生多少個alpha
    iter = 0                                   # 遍歷計數器
    while (iter < maxIter):
        #print( "iteration number: %d" % iter)
        alphaPairsChanged = 0                  # 記錄alpha是否已被優化，每次循環都重置
        for i in range(m):                     # 按行遍歷數據，類似隨機梯度下降
            # i=0
            fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b # 預測值y，g(x)函數，《統計學習方法》李航P127，7.104
            Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions # 誤差，Ei函數，P127，7.105
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                # 找第一個alphas[i]，找到第一個滿足判斷條件的，判斷負間隔or正間隔，並且保證0<alphas<C
                j = selectJrand(i,m)            # 隨機找到第二個alphas[j]
                fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b # 計算預測值
                Ej = fXj - float(labelMat[j])   # 計算alphas[j]誤差
                alphaIold = alphas[i].copy()    # 記錄上一次alphas[i]值
                alphaJold = alphas[j].copy()    # 記錄上一次alphas[j]值
                if (labelMat[i] != labelMat[j]):# 計算H及L值，《統計學習方法》李航，P126
                    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
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                # 《統計學習方法》李航P127，7.107，這里的eta與李航的一致，這里乘了負號
                if eta >= 0: 
                    #print("eta>=0")
                    continue
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta     # 《統計學習方法》李航P127，7.107，更新alphas[j]
                alphas[j] = clipAlpha(alphas[j],H,L)       # alphas[j]調整大於H或小於L的alpha值
                if (abs(alphas[j] - alphaJold) < 0.00001): # 調整后過小，則不更新alphas[i]
                    #print( "j not moving enough")
                    continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j]) #更新alphas[i]，《統計學習方法》李航P127，7.109
                # 更新b值，《統計學習方法》李航P130，7.115，7.116
                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
                if (0 < alphas[i]) and (C > alphas[i]): # 判斷符合條件的b
                    b = b1
                elif (0 < alphas[j]) and (C > alphas[j]):
                    b = b2
                else: 
                    b = (b1 + b2)/2.0
                alphaPairsChanged += 1
                #print( "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
        if (alphaPairsChanged == 0):
            iter += 1
        else:
            iter = 0
    return b,alphas

根據alpha求解w， $\vec{w}=\sum_{i=1}^{n}{\alpha_iy_i\vec{x_i}}$：

def calcWs(alphas,dataArr,labelArr):
    w=sum( np.array(alphas) * np.array(labelArr.reshape((-1,1))) * np.array(np.array(dataArr))  )
    return w

線性可分支持向量機運行結果：

[圖5]

3.線性支持向量機

3.1軟間隔支持向量機

[圖6]
如上圖，直線是我們學習獲得的超平面，但有一個圓點離群了，我們是否應該不考慮這個點，獲得新的虛線所示的超平面，這個虛線所示超平面泛化能力比原有超平面要好很多。 依此想法我們引出了軟間隔支持向量機，不一定分類完全正確的，超平面就是最好，樣本數據身線性不可分的情況也可以使用。軟間隔支持向量機容許在一些樣本上出錯，為此引入了“軟間隔”（soft margin）的概念

[圖3]，沒錯還是圖3…
具體來說，硬間隔支持向量機要求所有的樣本均被最佳超平面正確划分，而軟間隔支持向量機允許某些樣本點不滿足間隔大於等於1的條件 $y_i(\vec{w}\vec{x_i}+b)\geq1$，當然在最大化間隔的時候也要限制不滿足間隔大於等於1的樣本的個數使之盡可能的少。於是我們引入一個懲罰系數 $C>0$，並對每個樣本點 $(\vec{x_i},y_i)$引入一個松弛變量（slack variables） $ξ≥0$:

$\begin{cases} \min_{\vec{w},b}(\frac{1}{2}||\vec{w}||^2+C\sum_{i=1}^{n}\xi_i) \\ s.t.\quad y_i(\vec{w}^T\vec{x_i}+b)\geq1-\xi_i \quad ,i=1,2,\dots,n\\ \quad\quad \xi_i\geq0 \end{cases}$

利用拉格朗日乘子法可得到上式的拉格朗日函數：

$L(\vec{w},b,\vec{\alpha},\vec{\xi},\vec{\mu})=\frac{1}{2}||\vec{w}||^2+C\sum_{i=1}^{n}\xi_i-\sum_{i=1}^{n}\alpha_i(y_i(\vec{w}^T\vec{x_i}+b)-1+\xi_i)-\sum_{i=1}^{n}\mu_i\xi_i$

將上式帶入 $L(\vec{w},b,\vec{\alpha})$函數替換掉 $\vec{w}$和 $b$，得到以下結果：

$\begin{cases} \max_{\vec\alpha}\sum_{i=1}^{n}\alpha_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_jy_iy_j\vec{x_i}^T\vec{x_j}\\ s.t. \quad\quad\sum_{i=1}^{n}\alpha_iy_i=0\quad\quad,i=1,2,\dots,n\\ \quad\quad\quad 0\leq\alpha_i\leq C \end{cases}$

對比軟間隔支持向量機的對偶問題和硬間隔支持向量機的對偶問題可發現二者的唯一差別就在於對偶變量的約束不同，軟間隔支持向量機對對偶變量的約束是 $0≤α_i≤C$，硬間隔支持向量機對對偶變量的約束是 $0≤α_i$。

軟間隔支持向量機，KKT條件要求：
$\begin{cases} \alpha_i\geq0,\mu_i\geq0\\ y_i(\vec{w}\vec{x}+b)-1+\xi_i\geq0\\ \alpha_i(y_i(\vec{w}\vec{x}+b)-1+\xi_i)=0\\ \xi_i\geq0,\mu_i\xi_i=0 \end{cases}$

同硬間隔支持向量機類似，對任意訓練樣本 $(\vec{x_i},y_i)$，總有 $α_i=0$或 $y_i(\vec{w}\vec{x}+b-1+\xi_i)$，若αi=0αi=0，則該樣本不會對最佳決策面有任何影響；若 $α_i>0$則必有 $y_i(\vec{w}\vec{x}+b)=1-\xi_i$，也就是說該樣本是支持向量。若 $α_i<C$則 $μ_i>0$進而有 $ξ_i=0$，即該樣本處在最大間隔邊界上；若 $α_i=C$則 $μ_i=0$此時如果 $xi_i≤1$則該樣本處於最大間隔內部，如果 $ξ_i>1$則該樣本處於最大間隔外部即被分錯了。由此也可看出，軟間隔支持向量機的最終模型僅與支持向量有關。

4.非線性支持向量機

4.1核函數

現實任務中原始的樣本空間D中很可能並不存在一個能正確划分兩類樣本的超平面。工程學中經常會遇到的問題，無法找到一個超平面將兩類樣本進行很好的划分。 這種情況我們引入核函數（ kernel function ）來解決，將 $(\vec{x_i}*x)$引入核函數。

[圖7]
令 $ϕ(\vec{x} )$表示將樣本點 $\vec{x}$映射后的特征向量，類似於線性可分支持向量機中的表示方法，在特征空間中划分超平面所對應的模型可表示為:

$f(\vec{x})=\vec{w}^Tx+b$

其中 $\vec{w}$和 $b$是待求解的模型參數。

$\begin{cases} \min_{\vec{w},b}\frac{1}{2}||\vec{w}||^2\\ s.t. \quad y_i(\vec{w}^T\phi(\vec{x})+b)\geq1\quad,i=1,2,\dots,n \end{cases}$

其拉格朗日對偶問題是:

$\begin{cases} \max_\alpha{\sum_{i=1}^{n}\alpha_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_jy_iy_j\phi(\vec{x_i}^T)\phi(\vec{x_j})}\\ s.t.\quad \alpha_i\geq0\quad\quad\quad,i=1,2,\dots,n \\ \quad\quad\quad \sum_{i=1}^{n}\alpha_iy_i=0 \end{cases}$

需要計算 $ϕ(\vec{x_i}^T)$ $ϕ(\vec{x_j})$，即樣本映射到特征空間之后的內積，由於特征空間可能維度很高，甚至可能是無窮維，因此直接計算 $ϕ(\vec{x_i}^T)$ $ϕ(\vec{x_j})$通常是很困難的，在上文中我們提到其實我們根本不關心單個樣本的表現，只關心特征空間中樣本間兩兩的乘積，因此我們沒有必要把原始空間的樣本一個個地映射到特征空間中，只需要想法辦求解出樣本對應到特征空間中樣本間兩兩的乘積即可。為了解決該問題可設想存在核函數：

$\kappa(\vec{x_i},\vec{x_j})=\phi(\vec{x_i}^T)\phi(\vec{x_j})$

也就是說 $\vec{x_i}$與 $\vec{x_j}$在特征空間的內積等於它們在原始空間中通過函數κ(⋅,⋅)計算的結果，這給求解帶來很大的方便。

$\begin{cases} \max_\alpha{\sum_{i=1}^{n}\alpha_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_jy_iy_j\kappa<\vec{x_i},\vec{x_j}>}\\ s.t.\quad \alpha_i\geq0\quad\quad\quad,i=1,2,\dots,n \\ \quad\quad\quad \sum_{i=1}^{n}\alpha_iy_i=0 \end{cases}$

同樣的我們只關心在高維空間中樣本之間兩兩點乘的結果而不關心樣本是如何變換到高維空間中去的。求解后即可得到:

$f(\vec{x})=\vec{w}^T\phi(\vec{x})+b=\sum_{i=1}^{n}\alpha_iy_i\phi(\vec{x})^T\phi(\vec{x})+b=\sum_{i=1}^{n}\alpha_iy_i\kappa<\vec{x_i},\vec{x_j}>+b$

剩余的問題同樣是求解 $α_i$，然后求解 $\vec{w}$和 $b$即可得到最佳超平面。

常用的核函數：

[圖8]

核函數python實現:

# 核函數
def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space
    m,n = np.shape(X)
    K = np.mat(np.zeros((m,1)))
    if kTup[0]=='lin': K = X * A.T   #linear kernel
    elif kTup[0]=='rbf':
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow*deltaRow.T
        K = np.exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlab
    else: raise NameError('Houston We Have a Problem -- \
    That Kernel is not recognized')
    return K

4.2SMO算法

SMO算法的思想和梯度下降法的思想差不多。唯一不同的是，SMO是一次迭代優化兩個α而不是一個。為什么要優化兩個呢？
SVM的學習問題可以轉化為下面的對偶問題：

$\begin{cases} \max_\alpha{\sum_{i=1}^{n}\alpha_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_jy_iy_j\kappa<\vec{x_i},\vec{x_j}>}\\ s.t.\quad \alpha_i\geq0\quad\quad\quad,i=1,2,\dots,n \\ \quad\quad\quad \sum_{i=1}^{n}\alpha_iy_i=0 \end{cases}$

這個優化問題，我們可以看到這個優化問題存在着一個約束，也就是：

$\sum_{i=1}^{n}\alpha_iy_i=0$

假設我們首先固定除 $α_1$以外的所有參數，然后在 $α_1$上求極值。但需要注意的是，因為如果固定 $α_1$以外的所有參數，由上面這個約束條件可以知道， $α_1$將不再是變量（可以由其他值推出），因為問題中規定了：

$α_1y^{(1)}=\sum_{i=2}^{n}\alpha_iy_i$

因此，我們需要一次選取兩個參數做優化，比如 $α_i$和 $α_j$，此時 $α_i$可以由 $α_j$和其他參數表示出來。這樣回代入W中，W就只是關於 $α_j$的函數了，這時候就可以只對 $α_j$進行優化了。在這里就是對 $α_j$進行求導，令導數為0就可以解出這個時候最優的 $α_j$了。然后也可以得到αi。這就是一次的迭代過程，一次迭代只調整兩個拉格朗日乘子 $α_i$和 $α_j$。SMO之所以高效就是因為在固定其他參數后，對一個參數優化過程很高效（對一個參數的優化可以通過解析求解，而不是迭代。雖然對一個參數的一次最小優化不可能保證其結果就是所優化的拉格朗日乘子的最終結果，但會使目標函數向極小值邁進一步，這樣對所有的乘子做最小優化，直到所有滿足KKT條件時，目標函數達到最小）。

總結下來是：
重復下面過程直到收斂{
（1）選擇兩個拉格朗日乘子 $α_i$和 $α_j$；
（2）固定其他拉格朗日乘子 $α_k$(k不等於i和j)，只對 $α_i$和 $α_j$優化 $w(α)$;
（3）根據優化后的 $α_i$和 $α_j$，更新截距b的值；
}

4.2.1選擇 $α_i$和 $α_j$：

我們現在是每次迭代都優化目標函數的兩個拉格朗日乘子 $α_i$和 $α_j$，然后其他的拉格朗日乘子保持固定。如果有N個訓練樣本，我們就有N個拉格朗日乘子需要優化，但每次我們只挑兩個進行優化，我們就有N(N-1)種選擇。那到底我們要選擇哪對 $α_i$和 $α_j$呢？選擇哪對才好呢？想想我們的目標是什么？我們希望把所有違法KKT條件的樣本都糾正回來，因為如果所有樣本都滿足KKT條件的話，我們的優化就完成了。那就很直觀了，哪個害群之馬最嚴重，我們得先對他進行思想教育，讓他盡早回歸正途。OK，我們選擇的第一個變量αi就選違法KKT條件最嚴重的那一個。第二個變量 $α_j$選擇步長最大的那一個。
1）第一個變量αi的選擇：
SMO稱選擇第一個變量的過程為外層循環。外層訓練在訓練樣本中選取違法KKT條件最嚴重的樣本點。並將其對應的變量作為第一個變量。具體的，檢驗訓練樣本 $(x_i, y_i)$是否滿足KKT條件:

$a_i=0 {\Leftrightarrow}y_ig(x_i)>=1$
$0<a_i<C{\Leftrightarrow}y_ig(x_i)=1$
$a_i=C{\Leftrightarrow}y_ig(x_i)<=1$
$其中，g(x_i)=\sum_{i=1}^Na_jy_j\kappa<\vec{x_i},\vec{x_j}>+b$

該檢驗是在 $ε$范圍內進行的。在檢驗過程中，外層循環首先遍歷所有滿足條件 $0<α_j<C$的樣本點，即在間隔邊界上的支持向量點，檢驗他們是否滿足KKT條件，然后選擇違反KKT條件最嚴重的 $α_i$。如果這些樣本點都滿足KKT條件，那么遍歷整個訓練集，檢驗他們是否滿足KKT條件，然后選擇違反KKT條件最嚴重的 $α_i$。

優先選擇遍歷非邊界數據樣本，因為非邊界數據樣本更有可能需要調整，邊界數據樣本常常不能得到進一步調整而留在邊界上。由於大部分數據樣本都很明顯不可能是支持向量，因此對應的α乘子一旦取得零值就無需再調整。遍歷非邊界數據樣本並選出他們當中違反KKT 條件為止。當某一次遍歷發現沒有非邊界數據樣本得到調整時，遍歷所有數據樣本，以檢驗是否整個集合都滿足KKT條件。如果整個集合的檢驗中又有數據樣本被進一步進化，則有必要再遍歷非邊界數據樣本。這樣，不停地在遍歷所有數據樣本和遍歷非邊界數據樣本之間切換，直到整個樣本集合都滿足KKT條件為止。以上用KKT條件對數據樣本所做的檢驗都以達到一定精度 $ε$就可以停止為條件。如果要求十分精確的輸出算法，則往往不能很快收斂。

對整個數據集的遍歷掃描相當容易，而實現對非邊界 $α_i$的掃描時，首先需要將所有非邊界樣本的 $α_i$值（也就是滿足 $0<α_i<C$）保存到新的一個列表中，然后再對其進行遍歷。同時，該步驟跳過那些已知的不會改變的 $α_i$值。

2）第二個變量 $α_j$的選擇：
在選擇第一個 $α_i$后，算法會通過一個內循環來選擇第二個 $α_j$值。因為第二個乘子的迭代步長大致正比於 $|E_i-E_j|$，所以我們需要選擇能夠最大化 $|E_i-E_j|$的第二個乘子（選擇最大化迭代步長的第二個乘子）。在這里，為了節省計算時間，我們建立一個全局的緩存用於保存所有樣本的誤差值，而不用每次選擇的時候就重新計算。我們從中選擇使得步長最大或者 $|E_i-E_j|$最大的 $α_j$。

4.2.2優化 $α_i$和 $α_j$：

選擇這兩個拉格朗日乘子后，我們需要先計算這些參數的約束值。然后再求解這個約束最大化問題。
首先，我們需要給 $α_j$找到邊界 $L<=α_j<=H$，以保證 $α_j$滿足 $0<=α_j<=C$的約束。這意味着 $α_j$必須落入這個盒子中。由於只有兩個變量 $(α_i, α_j)$，約束可以用二維空間中的圖形來表示，如下圖：

[圖9]
不等式約束使得 $(α_i, α_j)$在盒子[0, C]x[0, C]內，等式約束使得 $(α_i, α_j)$在平行於盒子[0, C]x[0, C]的對角線的直線上。因此要求的是目標函數在一條平行於對角線的線段上的最優值。這使得兩個變量的最優化問題成為實質的單變量的最優化問題。由圖可以得到， $α_j$的上下界可以通過下面的方法得到：

$If\quad y^{(i)}{\neq}y^{(j)} L=max(0,a_j-a_i),H=min(c,c+a_j-a_i)$
$If\quad y^{(i)}=y^{(j)} L=max(0,a_i+a_j-C),H=min(c,a_i+a_j)$

我們優化的時候， $α_j$必須要滿足上面這個約束。也就是說上面是 $α_j$的可行域。然后我們開始尋找 $α_j$，使得目標函數最大化。通過推導得到 $α_j$的更新公式如下：

$a_j:=a_j-\frac{y^{(j)}(E_i-E_j)}{\eta}$

這里 $E_k$可以看做對第k個樣本，SVM的輸出與期待輸出，也就是樣本標簽的誤差。

$E_k=f(x^{(k)})-y^{(k)}$
$\eta=2<x^{(i)},x^{(j)}>-<x^{(i)},x^{(i)}>-<x^{(j)},x^{(j)}>$

而η實際上是度量兩個樣本i和j的相似性的。在計算η的時候，我們需要使用核函數，那么就可以用核函數來取代上面的內積。

得到新的 $α_j$后，我們需要保證它處於邊界內。換句話說，如果這個優化后的值跑出了邊界L和H，我們就需要簡單的裁剪，將 $α_j$收回這個范圍：
$α_j:=\begin{cases} H{\quad} if{\quad}a_j>H \\ a_j{\quad}if{\quad}L<=a_j<=H\\ L {\quad}if{\quad}a_j<L \end{cases}$

最后，得到優化的 $α_j$后，我們需要用它來計算 $α_i$：

$a_i:=a_i+y^{(i)}y^{(j)}(a_j^{old}-a_j)$

到這里， $α_i$和 $α_j$的優化就完成了。

4.2.3計算閾值b：

優化 $α_i$和 $α_j$后，我們就可以更新閾值b，使得對兩個樣本i和j都滿足KKT條件。如果優化后 $α_i$不在邊界上（也就是滿足 $0<α_i<C$，這時候根據KKT條件，可以得到 $y_ig_i(x_i)=1$，這樣我們才可以計算 $b$），那下面的閾值 $b1$是有效的，因為當輸入 $x_i$時它迫使SVM輸出 $y_i$。
$b1:=b-E_i-y^{(i)}(a_i-a_i^{(old)})<x^{(i)},x^{(i)}>-y^{(j)}(a_j-a_j^{(old)}<x^{(i)},x^{(j)}>$
同樣，如果 $0<α_j<C$，那么下面的b2也是有效的：
$b2=b-E_j-y^{(i)}(a_i-a_i^{(old)})<x^{(i)},x^{(j)}>-y^{(j)}(a_j-a_j^{(old)})<x^{(j)},x^{(j)}>$
如果 $0<α_i<C$和 $0<α_j<C$都滿足，那么b1和b2都有效，而且他們是相等的。如果他們兩個都處於邊界上（也就是 $α_i$=0或者 $α_i=C$，同時 $α_j=0$或者 $α_j=C$），那么在b1和b2之間的閾值都滿足KKT條件，一般我們取他們的平均值 $b=(b1+b2)/2$。所以，總的來說對b的更新如下：
$b:=\begin{cases} b1{\quad} if{\quad}0<a_i<C \\ b2{\quad}if{\quad}0<a_j<C\\ (b1+b2)/2 {\quad} otherwise \end{cases}$
每做完一次最小優化，必須更新每個數據樣本的誤差，以便用修正過的分類面對其他數據樣本再做檢驗，在選擇第二個配對優化數據樣本時用來估計步長。

4.2.4凸優化問題終止條件：

SMO算法的基本思路是：如果說有變量的解都滿足此最優化問題的KKT條件，那么這個最優化問題的解就得到了。因為KKT條件是該最優化問題的充分必要條件（證明請參考文獻）。所以我們可以監視原問題的KKT條件，所以所有的樣本都滿足KKT條件，那么就表示迭代結束了。但是由於KKT條件本身是比較苛刻的，所以也需要設定一個容忍值，即所有樣本在容忍值范圍內滿足KKT條件則認為訓練可以結束；當然了，對於對偶問題的凸優化還有其他終止條件，可以參考文獻。注意：SMO算法只是求解SVM的其中一種方法，大家在學習SVM過程中不必過於糾結SMO算法編程實現，只需要將理論思想理解即可。

4.2.5Python實現完整SMO算法核心代碼

尋找第2個步長最大的alphas[j]：

# 尋找第2個步長最大的alphas[j]
def selectJ(i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]  #set valid #choose the alpha that gives the maximum delta E
    validEcacheList = np.nonzero(oS.eCache[:,0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E
            if k == i: continue #don't calc for i, waste of time
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej
    else:   #in this case (first time around) we don't have any valid eCache values
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

SMO主函數:

# SMO主函數
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #full Platt SMO
    oS = optStruct(np.mat(dataMatIn),np.mat(classLabels).transpose(),C,toler, kTup)
    iter = 0
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:   #go over all
            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:#go over non-bound (railed) alphas
            nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i,oS)
                #print( "non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            iter += 1
        if entireSet: entireSet = False #toggle entire set loop
        elif (alphaPairsChanged == 0): entireSet = True  
        #print( "iteration number: %d" % iter)
    return oS.b,oS.alphas

SMO完整代碼輸出結果：

[圖10]

5.完整源代碼

# -*- coding: utf-8 -*-
""" @Time : 2018/11/25 19:42 @Author : hanzi5 @Email : [e-mail] @File : SVM.py @Software: PyCharm """
import numpy as np
#import pandas as pd
from scipy.optimize import fmin_l_bfgs_b
from scipy.optimize import fmin_tnc
#from scipy.optimize import fmin_bfgs 
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.patches import Circle

matplotlib.rcParams['font.family']='SimHei'  # 用來正常顯示中文
plt.rcParams['axes.unicode_minus']=False  # 用來正常顯示負號

# 定義目標函數，w是未知數據，args是已知數，求w的最優解
def func(w,*args):
    X,Y,c=args
    yp=np.dot(X,w)  # w*x，注意已經將b作為x的第一列進行了計算，w*x就是我們現在預測的y
    idx=np.where(yp*Y<1)[0] # 找到分錯的數據索引位置
    e=yp[idx]-Y[idx]        # y預測值-y真實值，誤差
    cost=np.dot(e,e)+c*np.dot(w,w)  # 平方和損失，c：學習率，加w的二范式懲罰
    grand=2*(np.dot(X[idx].T,e)+c*w)# 梯度下降？？
    return cost,grand

def plotResult(w):
    margin=2/np.sqrt(np.dot(w[1:3],w[1:3]))
    plot_x=np.append(np.min(x,0)[0]-0.2,np.max(x,0)[0]+0.2)
    plot_y=-(plot_x*w[1]+w[0])/w[2]
    plt.figure()
    pos=(Y==1) # 正類
    neg=(y==-1) # 負類
    plt.plot(x[pos][:,0],x[pos][:,1],"r+",label="正類")
    plt.plot(x[neg][:,0],x[neg][:,1],"bo",label="負類")
    plt.plot(plot_x,plot_y,"r-",label="分割超平面")
    plt.plot(plot_x,plot_y+margin/2,"g-.",label="")
    plt.plot(plot_x,plot_y-margin/2,"g-.",label="")
    plt.xlabel('x1')
    plt.ylabel('x2')
    plt.title('SVM Demo')
    plt.legend()
    plt.show()

# 簡易smo算法開始####################################################################
# 隨機選擇第2個alpha
def selectJrand(i,m):
    j=i #we want to select any J not equal to i
    while (j==i):
        j = int(np.random.uniform(0,m))
    return j

# 調整大於H或小於L的alpha值
def clipAlpha(aj,H,L):
    if aj > H: 
        aj = H
    if L > aj:
        aj = L
    return aj

# 公共函數，根據公式求w，簡易smo算法及完整smo算法通用
def calcWs(alphas,dataArr,labelArr):
    w=sum( np.array(alphas) * np.array(labelArr.reshape((-1,1))) * np.array(np.array(dataArr))  )
    return w


# 輸入變量：x、y、c：常數c、toler：容錯率、maxIter：最大循環次數
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    #dataMatIn, classLabels, C, toler, maxIter=dataArr,lableArr,0.6,0.001,40
    dataMatrix = np.mat(dataMatIn)             # 數據x轉換為matrix類型
    labelMat = np.mat(classLabels).transpose() # 標簽y轉換為matrix類型，轉換為一列
    b = 0                                      # 截距b
    m,n = np.shape(dataMatrix)                 # 數據x行數、列數
    alphas = np.mat(np.zeros((m,1)))           # 初始化alpha，有多少行數據就產生多少個alpha
    iter = 0                                   # 遍歷計數器
    while (iter < maxIter):
        #print( "iteration number: %d" % iter)
        alphaPairsChanged = 0                  # 記錄alpha是否已被優化，每次循環都重置
        for i in range(m):                     # 按行遍歷數據，類似隨機梯度下降
            # i=0
            fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b # 預測值y，g(x)函數，《統計學習方法》李航P127，7.104
            Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions # 誤差，Ei函數，P127，7.105
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                # 找第一個alphas[i]，找到第一個滿足判斷條件的，判斷負間隔or正間隔，並且保證0<alphas<C
                j = selectJrand(i,m)            # 隨機找到第二個alphas[j]
                fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b # 計算預測值
                Ej = fXj - float(labelMat[j])   # 計算alphas[j]誤差
                alphaIold = alphas[i].copy()    # 記錄上一次alphas[i]值
                alphaJold = alphas[j].copy()    # 記錄上一次alphas[j]值
                if (labelMat[i] != labelMat[j]):# 計算H及L值，《統計學習方法》李航，P126
                    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
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                # 《統計學習方法》李航P127，7.107，這里的eta與李航的一致，這里乘了負號
                if eta >= 0: 
                    #print("eta>=0")
                    continue
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta     # 《統計學習方法》李航P127，7.107，更新alphas[j]
                alphas[j] = clipAlpha(alphas[j],H,L)       # alphas[j]調整大於H或小於L的alpha值
                if (abs(alphas[j] - alphaJold) < 0.00001): # 調整后過小，則不更新alphas[i]
                    #print( "j not moving enough")
                    continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j]) #更新alphas[i]，《統計學習方法》李航P127，7.109
                # 更新b值，《統計學習方法》李航P130，7.115，7.116
                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
                if (0 < alphas[i]) and (C > alphas[i]): # 判斷符合條件的b
                    b = b1
                elif (0 < alphas[j]) and (C > alphas[j]):
                    b = b2
                else: 
                    b = (b1 + b2)/2.0
                alphaPairsChanged += 1
                #print( "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
        if (alphaPairsChanged == 0):
            iter += 1
        else:
            iter = 0
    return b,alphas

# 畫圖
def plot_smoSimple(dataArrWithAlpha,b,w):
    type1_x1 = []
    type1_x2 = []
    type2_x1 = []
    type2_x2 = []
    dataSet=dataArrWithAlpha 
    # 取兩類x1及x2值畫圖
    type1_x1=dataSet[dataSet[:,-2]==-1][:,:-2][:,0].tolist() 
    type1_x2=dataSet[dataSet[:,-2]==-1][:,:-2][:,1].tolist()
    type2_x1=dataSet[dataSet[:,-2]==1][:,:-2][:,0].tolist()
    type2_x2=dataSet[dataSet[:,-2]==1][:,:-2][:,1].tolist()
    
    # 畫點
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(type1_x1,type1_x2, marker='s', s=90)
    ax.scatter(type2_x1,type2_x2, marker='o', s=50, c='red')
    plt.title('Support Vectors Circled')
    
    # 獲取支持向量值，畫橢圓
    dataVectors=dataArrWithAlpha[dataArrWithAlpha[:,-1]>0]
    for d in dataVectors:
        circle = Circle(d[0:2], 0.5, facecolor='none', edgecolor=(0,0.8,0.8), linewidth=3, alpha=0.5)
        ax.add_patch(circle)
        
    # 畫分割超平面
    b=b1.getA()[0][0] # 獲得傳入的b
    w0= w[0]#0.8065
    w1= w[1]#-0.2761
    x = np.arange(-2.0, 12.0, 0.1)
    y = (-w0*x - b)/w1
    ax.plot(x,y)
    ax.axis([-2,12,-8,6])
    plt.show()
# 簡易smo算法結束####################################################################

# 完整smo算法開始####################################################################
class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters 
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = np.shape(dataMatIn)[0]
        self.alphas = np.mat(np.zeros((self.m,1)))
        self.b = 0
        self.eCache = np.mat(np.zeros((self.m,2))) #first column is valid flag
        self.K = np.mat(np.zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

# 核函數
def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space
    m,n = np.shape(X)
    K = np.mat(np.zeros((m,1)))
    if kTup[0]=='lin': K = X * A.T   #linear kernel
    elif kTup[0]=='rbf':
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow*deltaRow.T
        K = np.exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlab
    else: raise NameError('Houston We Have a Problem -- \
    That Kernel is not recognized')
    return K

# 計算誤差
def calcEk(oS, k):
    fXk = float(np.multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek

# 尋找第2個步長最大的alphas[j]
def selectJ(i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]  #set valid #choose the alpha that gives the maximum delta E
    validEcacheList = np.nonzero(oS.eCache[:,0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E
            if k == i: continue #don't calc for i, waste of time
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej
    else:   #in this case (first time around) we don't have any valid eCache values
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

# 計算誤差存入緩存中
def updateEk(oS, k):#after any alpha has changed update the new value in the cache
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]

# 內循環，尋找第2個步長最大的alphas[j]
def innerL(i, oS):
    Ei = calcEk(oS, i)
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        j,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
        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
        eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel
        if eta >= 0: 
            #print( "eta>=0")
            return 0
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS, j) #added this for the Ecache
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): 
            #print( "j not moving enough")
            return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
        updateEk(oS, i) #added this for the Ecache #the update is in the oppostie direction
        b1 = 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]
        b2 = 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]
        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

# SMO主函數
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #full Platt SMO
    oS = optStruct(np.mat(dataMatIn),np.mat(classLabels).transpose(),C,toler, kTup)
    iter = 0
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:   #go over all
            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:#go over non-bound (railed) alphas
            nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i,oS)
                #print( "non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            iter += 1
        if entireSet: entireSet = False #toggle entire set loop
        elif (alphaPairsChanged == 0): entireSet = True  
        #print( "iteration number: %d" % iter)
    return oS.b,oS.alphas

# 測試Rbf數據
def testRbf(dataArrTrain,labelArrTrain,dataArrTest,labelArrTest,k1=1.3):
    b,alphas = smoP(dataArrTrain, labelArrTrain, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important
    datMat=np.mat(dataArrTrain); labelArrTrain = np.mat(labelArrTrain).transpose()
    svInd=np.nonzero(alphas.A>0)[0]
    sVs=datMat[svInd] #get matrix of only support vectors
    labelSV = labelArrTrain[svInd];
    #print( "there are %d Support Vectors" % np.shape(sVs)[0])
    m,n = np.shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b
        if np.sign(predict)!=np.sign(labelArrTrain[i]): errorCount += 1
    print( "the training error rate is: %f" % (float(errorCount)/m))
    errorCount = 0
    datMat=np.mat(dataArrTest)
    #labelMat = np.mat(labelArrTest).transpose()
    m,n = np.shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b
        if np.sign(predict)!=np.sign(labelArrTest[i]): errorCount += 1    
    print( "the test error rate is: %f" % (float(errorCount)/m)   )
    return b,alphas

# 畫圖，rbf核函數數據
def plot_smoCompletion ():
    xcord0 = []; ycord0 = []; xcord1 = []; ycord1 = []
    fw = open('D:/python_data/testSetRBF2.txt', 'w')#generate data
    fig = plt.figure()
    ax = fig.add_subplot(111)
    xcord0 = []; ycord0 = []; xcord1 = []; ycord1 = []
    for i in range(100):
        [x,y] = np.random.uniform(0,1,2)
        xpt=x*np.cos(2.0*np.pi*y); ypt = x*np.sin(2.0*np.pi*y)
        if (x > 0.5):
            xcord0.append(xpt); ycord0.append(ypt)
            label = -1.0
        else:
            xcord1.append(xpt); ycord1.append(ypt)
            label = 1.0
        fw.write('%f\t%f\t%f\n' % (xpt, ypt, label))
    ax.scatter(xcord0,ycord0, marker='s', s=90)
    ax.scatter(xcord1,ycord1, marker='o', s=50, c='red')
    plt.title('Non-linearly Separable Data for Kernel Method')
    plt.show()
    fw.close()
# 完整smo算法結束####################################################################

if __name__ == '__main__':
    ##1、SVM直接求參數值#############################################################
    print('\n1、SVM直接求參數值，開始')
    # 生成數據，《統計學習方法》李航，P103，例7.1
    dataSet=np.array([[3,3,1],[4,3,1],[1,1,-1]]) # ,[0,0,-1],[0,1,-1]
    m, n = dataSet.shape
    x=dataSet[:,:-1]
    y=dataSet[:,-1] #.reshape((-1,1))
    # 數據定義
    X=np.append(np.ones([x.shape[0],1]),x,1) # x新增一列全1值，作為截距b
    Y=y
    c=0.001 # 學習率
    w=np.zeros(X.shape[1]) # 初始化一組w系數，全0，也可隨機產生：np.random.rand(X.shape[1])
    # bfgs_b方法求最優化問題
    REF=fmin_l_bfgs_b(func,x0=w,args=(X,Y,c),approx_grad=False) #x0=np.random.rand(X.shape[1]) [0,0,0]
    # 采用scipy.optimize其他包夜可以求得
    REF2=fmin_tnc(func,x0=w,args=(X,Y,c),approx_grad=False)
    # 求得最優化計算后的w
    w=REF[0].round(2)           # 取得w值
    print('w:',w[1:],'b:',w[0]) # 與《統計學習方法》李航，P103，例7.1計算結果一致
    # 畫圖
    plotResult(w)
    print('\n1、SVM直接求參數值，結束')
    
    ##2、SVM簡易SMO算法#############################################################
    print('\n2、SVM簡易SMO算法，開始')
    fileIn = 'D:/python_data/testSet.txt'
    #dataSet=pd.read_table(fileIn,names=['x1','x2','y']).values
    dataSet=np.loadtxt(fileIn)
    dataArr=dataSet[:,:-1] # x
    labelArr=dataSet[:,-1] # y
    b1,alphas1=smoSimple(dataArr,labelArr,0.6,0.001,50) # 輸入變量：x、y、c：常數c、toler：容錯率、maxIter：最大循環次數
    dataArrWithAlpha1=np.array(np.concatenate((dataSet,alphas1),axis=1)) # 把alphas1與原始數據合並
    w1=calcWs(alphas1,dataArr,labelArr)                                    # 根據alpha求w
    print('b:',b1,'\nw:',w1,'\ndata，alphas，支撐向量:\n',dataArrWithAlpha1[dataArrWithAlpha1[:,-1]>0] )# 注意這里的篩選方式與pd.DataFrame篩選方式一致，array類型的才可以這樣寫，np.ndarray及np.matrix類型不可以使用
    plot_smoSimple(dataArrWithAlpha1,b1,w1)  # 畫圖
    print('2、SVM簡易SMO算法，結束')
    
    ##3、SVM完整SMO算法#############################################################
    print('\n3、SVM完整SMO算法，開始')
    dataSetTrain = np.loadtxt('D:/python_data/testSetRBF.txt')
    dataSetTest  = np.loadtxt('D:/python_data/testSetRBF2.txt')
    # 訓練集
    dataArrTrain=dataSetTrain[:,:-1] # 訓練集x
    labelArrTrain=dataSetTrain[:,-1] # 訓練集y
    # 測試集
    dataArrTest=dataSetTest[:,:-1]   # 測試集x
    labelArrTest=dataSetTest[:,-1]   # 測試集y
    # 調用主函數
    b2,alphas2=testRbf(dataArrTrain,labelArrTrain,dataArrTest,labelArrTest,k1=1.3)
    w2=calcWs(alphas2,dataArrTrain,labelArrTrain)                                    # 根據alpha求w
    dataArrWithAlpha2=np.array(np.concatenate((dataSetTrain,alphas2),axis=1)) # 把alphas1與原始數據合並
    print('b:',b1,'\nw:',w1,'\ndata，alphas，支撐向量:\n',dataArrWithAlpha2[dataArrWithAlpha2[:,-1]>0] )# 注意這里的篩選方式與pd.DataFrame篩選方式一致，array類型的才可以這樣寫，np.ndarray及np.matrix類型不可以使用
    plot_smoCompletion() # 畫圖，訓練集
    print('3、SVM完整SMO算法，結束')

數據，testSet.txt

3.542485	1.977398	-1
3.018896	2.556416	-1
7.551510	-1.580030	1
2.114999	-0.004466	-1
8.127113	1.274372	1
7.108772	-0.986906	1
8.610639	2.046708	1
2.326297	0.265213	-1
3.634009	1.730537	-1
0.341367	-0.894998	-1
3.125951	0.293251	-1
2.123252	-0.783563	-1
0.887835	-2.797792	-1
7.139979	-2.329896	1
1.696414	-1.212496	-1
8.117032	0.623493	1
8.497162	-0.266649	1
4.658191	3.507396	-1
8.197181	1.545132	1
1.208047	0.213100	-1
1.928486	-0.321870	-1
2.175808	-0.014527	-1
7.886608	0.461755	1
3.223038	-0.552392	-1
3.628502	2.190585	-1
7.407860	-0.121961	1
7.286357	0.251077	1
2.301095	-0.533988	-1
-0.232542	-0.547690	-1
3.457096	-0.082216	-1
3.023938	-0.057392	-1
8.015003	0.885325	1
8.991748	0.923154	1
7.916831	-1.781735	1
7.616862	-0.217958	1
2.450939	0.744967	-1
7.270337	-2.507834	1
1.749721	-0.961902	-1
1.803111	-0.176349	-1
8.804461	3.044301	1
1.231257	-0.568573	-1
2.074915	1.410550	-1
-0.743036	-1.736103	-1
3.536555	3.964960	-1
8.410143	0.025606	1
7.382988	-0.478764	1
6.960661	-0.245353	1
8.234460	0.701868	1
8.168618	-0.903835	1
1.534187	-0.622492	-1
9.229518	2.066088	1
7.886242	0.191813	1
2.893743	-1.643468	-1
1.870457	-1.040420	-1
5.286862	-2.358286	1
6.080573	0.418886	1
2.544314	1.714165	-1
6.016004	-3.753712	1
0.926310	-0.564359	-1
0.870296	-0.109952	-1
2.369345	1.375695	-1
1.363782	-0.254082	-1
7.279460	-0.189572	1
1.896005	0.515080	-1
8.102154	-0.603875	1
2.529893	0.662657	-1
1.963874	-0.365233	-1
8.132048	0.785914	1
8.245938	0.372366	1
6.543888	0.433164	1
-0.236713	-5.766721	-1
8.112593	0.295839	1
9.803425	1.495167	1
1.497407	-0.552916	-1
1.336267	-1.632889	-1
9.205805	-0.586480	1
1.966279	-1.840439	-1
8.398012	1.584918	1
7.239953	-1.764292	1
7.556201	0.241185	1
9.015509	0.345019	1
8.266085	-0.230977	1
8.545620	2.788799	1
9.295969	1.346332	1
2.404234	0.570278	-1
2.037772	0.021919	-1
1.727631	-0.453143	-1
1.979395	-0.050773	-1
8.092288	-1.372433	1
1.667645	0.239204	-1
9.854303	1.365116	1
7.921057	-1.327587	1
8.500757	1.492372	1
1.339746	-0.291183	-1
3.107511	0.758367	-1
2.609525	0.902979	-1
3.263585	1.367898	-1
2.912122	-0.202359	-1
1.731786	0.589096	-1
2.387003	1.573131	-1

testSetRBF.txt

-0.214824	0.662756	-1.000000
-0.061569	-0.091875	1.000000
0.406933	0.648055	-1.000000
0.223650	0.130142	1.000000
0.231317	0.766906	-1.000000
-0.748800	-0.531637	-1.000000
-0.557789	0.375797	-1.000000
0.207123	-0.019463	1.000000
0.286462	0.719470	-1.000000
0.195300	-0.179039	1.000000
-0.152696	-0.153030	1.000000
0.384471	0.653336	-1.000000
-0.117280	-0.153217	1.000000
-0.238076	0.000583	1.000000
-0.413576	0.145681	1.000000
0.490767	-0.680029	-1.000000
0.199894	-0.199381	1.000000
-0.356048	0.537960	-1.000000
-0.392868	-0.125261	1.000000
0.353588	-0.070617	1.000000
0.020984	0.925720	-1.000000
-0.475167	-0.346247	-1.000000
0.074952	0.042783	1.000000
0.394164	-0.058217	1.000000
0.663418	0.436525	-1.000000
0.402158	0.577744	-1.000000
-0.449349	-0.038074	1.000000
0.619080	-0.088188	-1.000000
0.268066	-0.071621	1.000000
-0.015165	0.359326	1.000000
0.539368	-0.374972	-1.000000
-0.319153	0.629673	-1.000000
0.694424	0.641180	-1.000000
0.079522	0.193198	1.000000
0.253289	-0.285861	1.000000
-0.035558	-0.010086	1.000000
-0.403483	0.474466	-1.000000
-0.034312	0.995685	-1.000000
-0.590657	0.438051	-1.000000
-0.098871	-0.023953	1.000000
-0.250001	0.141621	1.000000
-0.012998	0.525985	-1.000000
0.153738	0.491531	-1.000000
0.388215	-0.656567	-1.000000
0.049008	0.013499	1.000000
0.068286	0.392741	1.000000
0.747800	-0.066630	-1.000000
0.004621	-0.042932	1.000000
-0.701600	0.190983	-1.000000
0.055413	-0.024380	1.000000
0.035398	-0.333682	1.000000
0.211795	0.024689	1.000000
-0.045677	0.172907	1.000000
0.595222	0.209570	-1.000000
0.229465	0.250409	1.000000
-0.089293	0.068198	1.000000
0.384300	-0.176570	1.000000
0.834912	-0.110321	-1.000000
-0.307768	0.503038	-1.000000
-0.777063	-0.348066	-1.000000
0.017390	0.152441	1.000000
-0.293382	-0.139778	1.000000
-0.203272	0.286855	1.000000
0.957812	-0.152444	-1.000000
0.004609	-0.070617	1.000000
-0.755431	0.096711	-1.000000
-0.526487	0.547282	-1.000000
-0.246873	0.833713	-1.000000
0.185639	-0.066162	1.000000
0.851934	0.456603	-1.000000
-0.827912	0.117122	-1.000000
0.233512	-0.106274	1.000000
0.583671	-0.709033	-1.000000
-0.487023	0.625140	-1.000000
-0.448939	0.176725	1.000000
0.155907	-0.166371	1.000000
0.334204	0.381237	-1.000000
0.081536	-0.106212	1.000000
0.227222	0.527437	-1.000000
0.759290	0.330720	-1.000000
0.204177	-0.023516	1.000000
0.577939	0.403784	-1.000000
-0.568534	0.442948	-1.000000
-0.011520	0.021165	1.000000
0.875720	0.422476	-1.000000
0.297885	-0.632874	-1.000000
-0.015821	0.031226	1.000000
0.541359	-0.205969	-1.000000
-0.689946	-0.508674	-1.000000
-0.343049	0.841653	-1.000000
0.523902	-0.436156	-1.000000
0.249281	-0.711840	-1.000000
0.193449	0.574598	-1.000000
-0.257542	-0.753885	-1.000000
-0.021605	0.158080	1.000000
0.601559	-0.727041	-1.000000
-0.791603	0.095651	-1.000000
-0.908298	-0.053376	-1.000000
0.122020	0.850966	-1.000000
-0.725568	-0.292022	-1.000000

testSetRBF2.txt

-0.117432	-0.772533	-1.000000
0.424996	0.685671	-1.000000
-0.366447	-0.032439	1.000000
-0.456538	-0.326062	-1.000000
-0.753378	-0.060207	-1.000000
-0.737534	0.087643	-1.000000
-0.093421	-0.126775	1.000000
-0.351430	0.926426	-1.000000
0.012625	-0.190257	1.000000
-0.351591	0.693200	-1.000000
0.388269	-0.852720	-1.000000
0.142043	0.318140	1.000000
0.264801	-0.401142	1.000000
0.136152	-0.090225	1.000000
0.974734	0.016945	-1.000000
0.244563	0.488831	-1.000000
-0.280123	-0.056987	1.000000
0.091556	-0.081223	1.000000
-0.090611	-0.016636	1.000000
-0.813602	-0.266067	-1.000000
-0.000300	0.000043	1.000000
0.242741	0.303703	1.000000
-0.817048	-0.008500	-1.000000
0.529506	0.798931	-1.000000
0.524553	-0.567558	-1.000000
0.381955	0.001980	1.000000
-0.602431	0.473769	-1.000000
0.231215	-0.234673	1.000000
-0.749613	-0.345047	-1.000000
0.224492	0.182258	1.000000
0.172463	-0.513747	-1.000000
0.147603	0.910376	-1.000000
0.734239	-0.286891	-1.000000
0.154865	-0.548473	-1.000000
-0.147829	0.009083	1.000000
-0.017590	-0.044190	1.000000
0.440672	0.691269	-1.000000
-0.357727	0.550873	-1.000000
0.457356	0.133797	1.000000
0.232346	-0.239350	1.000000
0.478807	0.852203	-1.000000
-0.026059	0.040957	1.000000
0.775044	-0.460803	-1.000000
-0.267704	-0.306558	1.000000
-0.855806	0.110270	-1.000000
0.480211	-0.521624	-1.000000
0.765252	-0.047894	-1.000000
-0.644040	0.586377	-1.000000
-0.110076	0.719103	-1.000000
0.133619	0.524928	-1.000000
-0.173815	0.732079	-1.000000
0.942163	0.221746	-1.000000
-0.687198	0.047586	-1.000000
0.630384	0.281265	-1.000000
0.038130	0.010993	1.000000
-0.008068	-0.118929	1.000000
0.534217	0.111413	-1.000000
-0.325657	0.273019	1.000000
-0.015575	0.736060	-1.000000
-0.333588	-0.932616	-1.000000
0.348827	0.863883	-1.000000
0.442462	0.389002	-1.000000
-0.101086	0.128123	1.000000
0.097092	0.321638	1.000000
0.125921	0.621098	-1.000000
-0.357218	0.095364	1.000000
0.803109	0.081817	-1.000000
-0.514764	-0.111062	-1.000000
0.047659	-0.176539	1.000000
-0.007704	0.044554	1.000000
-0.472410	-0.592451	-1.000000
-0.197716	-0.791531	-1.000000
-0.515508	-0.339354	-1.000000
0.642045	0.251044	-1.000000
-0.065945	-0.162388	1.000000
-0.744879	-0.045031	-1.000000
0.209621	0.472639	-1.000000
-0.544970	-0.167653	-1.000000
-0.431730	0.679799	-1.000000
0.542874	-0.724989	-1.000000
0.074219	0.805839	-1.000000
0.451146	-0.375435	-1.000000
0.743191	0.606076	-1.000000
0.529628	-0.195778	-1.000000
0.368259	0.348836	-1.000000
-0.933910	-0.214162	-1.000000
-0.119933	-0.518656	-1.000000
-0.627876	0.753061	-1.000000
-0.116855	-0.127113	1.000000
0.061212	-0.019597	1.000000
0.897120	-0.198540	-1.000000
-0.469099	0.232863	-1.000000
-0.212642	0.637247	-1.000000
0.012848	0.035170	1.000000
0.652035	-0.129829	-1.000000
-0.001256	0.719979	-1.000000
0.921166	-0.204044	-1.000000
-0.844109	-0.536146	-1.000000
0.428299	-0.813288	-1.000000
0.731897	-0.242729	-1.000000

參考資料：
1、《機器學習實戰》Peter Harrington著
2、《機器學習》西瓜書，周志華著
3、 斯坦福大學公開課 ：機器學習課程
4、機器學習視頻，鄒博
5、《統計學習方法》李航
6、SVM GarryLau
7、機器學習算法與Python實踐之（二、三、四） zouxy09