2. 程式人生 > >機器學習實戰（五）支援向量機SVM（Support Vector Machine）

機器學習實戰（五）支援向量機SVM（Support Vector Machine）

阿新 • • 發佈：2018-11-02

目錄

0. 前言

1. 尋找最大間隔

2. 拉格朗日乘子法和KKT條件

3. 鬆弛變數

4. 帶鬆弛變數的拉格朗日乘子法和KKT條件

5. 序列最小優化SMO（Sequential Minimal Optimization）

6. 高斯核函式（Gaussian Kernel）

7. 實戰案例

7.1. 簡化版SMO演算法案例

7.2. 完整版SMO演算法案例

7.3. 高斯核函式案例

學習完機器學習實戰的支援向量機，簡單的做個筆記。文中部分描述屬於個人消化後的理解，僅供參考。

本篇綜合了先前的文章，如有不理解，可參考：

吳恩達機器學習（十）支援向量機

所有程式碼和資料可以訪問 我的 github

如果這篇文章對你有一點小小的幫助，請給個關注喔~我會非常開心的~

0. 前言

支援向量機（Support Vector Machine）是一類監督學習的二分類演算法。

通過構建超平面，將資料分割開來，超平面一邊的資料屬於同一類別。

支援向量是指距離分割超平面最近的那些點，支援向量機目的是使得支援向量到超平面的間隔最大。

  • 優點：泛化錯誤率低，結果易於解釋
  • 對引數調節和核函式的選擇敏感
  • 適用資料型別：數值型和標稱型資料

如下圖所示（圖源：百度百科）：

H 為分割超平面，在超平面上下方各建立一個介面 H_1 和 H_2 ，在 H_1 之上的為一類，在 H_2 之下的為一類。

SVM的目的是為了最大化間隔，即 H_1 到 H_2 的距離。 

1. 尋找最大間隔

在SVM中，定義正類 y=1 ，定義反類 y=-1 ，這是為了方便後面計算。

通過 w^Tx+b=0 擬合數據的分界線，即超平面，則可定義上介面為 w^Tx+b=1 ，下介面為 w^Tx+b=-1

 。因可以通過改變係數的倍數表示距離，所以可粗略定義為 1 。

上介面以上的點表示為 w^Tx+b\geqslant 1,\ y=1 ，下介面以下的點表示為 w^Tx+b\leqslant -1,\ y=-1 ，所以綜合可得 y(w^Tx+b)\geqslant 1 。

需最大化上介面到下介面的距離：

d=\frac{\left|c_2-c_1\right|}{\sqrt{w_1^2+...+w_n^2}}=\frac{2}{\left\|w\right\|}

\max\ d\Rightarrow \min\ \frac{\left\|w\right\|}{2}\Rightarrow \min\ \frac{\left\|w\right\|^2}{2}\Rightarrow \min\ \frac{1}{2}w^Tw

所以，可得目標函式和約束條件：

\left\{\begin{align*} &\min\ \frac{1}{2}w^Tw\\ &s.t.\ y^{(i)}(w^Tx^{(i)}+b)\geqslant 1 \end{align*}\right.

2. 拉格朗日乘子法和KKT條件

  • 拉格朗日乘子法：求在約束條件下，函式的極值問題
  • KKT條件：若約束條件為不等式，則需滿足一定條件

在拉格朗日乘子法中，引入了其他的引數，並將原問題轉換為對偶問題（求最小值轉換為求最大值）。

首先修改目標函式和約束條件為：

\left\{\begin{align*} &\min\ \frac{1}{2}w^Tw\\ &s.t.\ y^{(i)}(w^Tx^{(i)}+b)\geqslant 1 \end{align*}\right. \Rightarrow \left\{\begin{align*} &\min\ \frac{1}{2}w^Tw\\ &s.t.\ 1-y^{(i)}(w^Tx^{(i)}+b)\leqslant 0\end{align*}\right.

應用拉格朗日乘子法：

\begin{align*} L(w,b,\alpha) &= \frac{1}{2}w^Tw+\alpha_1g_1(x)+...+\alpha_ng_n(x) \\ &= \frac{1}{2}w^Tw -\alpha_1(y^{(1)}(w^Tx^{(1)}+b)-1)-...-\alpha_n(y^{(n)}(w^Tx^{(n)}+b)-1) \\ &= \frac{1}{2}w^Tw -\sum_{i=1}^{m}\alpha_iy^{(i)}(w^Tx^{(i)}+b)+\sum_{i=1}^{m}\alpha_i \end{align*}

應滿足的KKT條件為：

\begin{align*} &\frac{\partial L}{\partial w} = w-\sum_{i=1}^{m}\alpha_iy^{(i)}x^{(i)}=0 \Rightarrow w=\sum_{i=1}^{m}\alpha_iy^{(i)}x^{(i)} \\ &\frac{\partial L}{\partial b} = -\sum_{i=1}^{m}\alpha_iy^{(i)}=0 \Rightarrow \sum_{i=1}^{m}\alpha_iy^{(i)}=0 \\ &\alpha \geqslant 0 \end{align*}

將KKT條件代回拉格朗日乘子法：

\begin{align*} L(w,b,\alpha) &= \frac{1}{2}(\sum_{i=1}^{m}\alpha_iy^{(i)}x^{(i)})^T(\sum_{j=1}^{m}\alpha_jy^{(j)}x^{(j)}) -\sum_{i=1}^{m}\alpha_iy^{(i)}((\sum_{j=1}^{m}\alpha_jy^{(j)}x^{(j)})^Tx^{(i)}+b)+\sum_{i=1}^{m}\alpha_i \\ &= -\frac{1}{2}\sum_{i,j=1}^{m}\alpha_iy^{(i)}\alpha_jy^{(j)}(x^{(i)})^Tx^{(j)}+b\sum_{i=1}^{m}\alpha_iy^{(i)}+\sum_{i=1}^{m}\alpha_i \\ &= \sum_{i=1}^{m}\alpha_i-\frac{1}{2}\sum_{i,j=1}^{m}\alpha_iy^{(i)}\alpha_jy^{(j)}(x^{(i)})^Tx^{(j)} \end{align*}

因轉換為對偶問題，所以此時目標函式和約束條件：

\left\{\begin{align*} &\max\ w(\alpha)=L(w,b,\alpha)\\ &s.t.\ \alpha \geqslant 0 \\ &s.t.\ \sum_{i=1}^{m}a_iy^{(i)}=0 \end{align*}\right.

3. 鬆弛變數

有時，資料並非 100\% 線性可分，有一部分資料會在上下介面之間，或者處於錯誤的類別。

此時可增加一個鬆弛變數 \epsilon _i\ (\epsilon _i \geqslant 0) ，允許有資料點位於錯誤的一側。

定義上介面以上為 w^Tx+b\geqslant 1-\epsilon _i,\ y=1 ，下介面以下為 w^Tx+b\leqslant -1+\epsilon _i,\ y=-1 ，可綜合為 y(w^Tx+b)\geqslant 1-\epsilon _i

所以，可得目標函式和約束條件，其中 C 是控制權重：

\left\{\begin{align*} &\min\ \frac{1}{2}w^Tw+C\sum_{i=1}^{m}\epsilon_i\\ &s.t.\ y^{(i)}(w^Tx^{(i)}+b)\geqslant 1-\epsilon _i \end{align*}\right.

4. 帶鬆弛變數的拉格朗日乘子法和KKT條件

同理，根據拉格朗日乘子法和KKT條件，可得：

L(w,b,\alpha) = \sum_{i=1}^{m}\alpha_i-\frac{1}{2}\sum_{i,j=1}^{m}\alpha_iy^{(i)}\alpha_jy^{(j)}(x^{(i)})^Tx^{(j)}

所以，目標函式和約束條件為：

\left\{\begin{align*} &\max\ w(\alpha)=L(w,b,\alpha)\\ &s.t.\ C\geqslant \alpha \geqslant 0 \\ &s.t.\ \sum_{i=1}^{m}a_iy^{(i)}=0 \end{align*}\right.

5. 序列最小優化SMO（Sequential Minimal Optimization）

SMO演算法的目標是求出一系列的 \alpha 和 b ，然後通過一開始的KKT條件，得出權重 w ：

w=\sum_{i=1}^{m}\alpha_iy^{(i)}x^{(i)}

SMO演算法的工作原理是，每次迴圈都選擇兩個 \alpha ，同時對兩個 \alpha 進行優化。

因需滿足 \sum_{i=1}^{m}a_iy^{(i)}=0 ，所以增大其中一個 \alpha 的同時，減小另一個 \alpha 。

6. 高斯核函式（Gaussian Kernel）

若資料是非線性可分，就無法使用線性的SVM進行解決。

可通過核函式，將低維空間中的資料對映到高維空間中，可將一個在低維空間中的非線性問題轉換為高維空間下的線性問題求解。

高斯核函式的定義如下，其中 \sigma 表示到達率（函式值跌落到 0 的速度引數），\sigma 越小，支援向量越多：

k(x,l)=\exp(-\frac{\left\|x-l\right\|^2}{2\sigma^2})

其中，\large l 稱為標記點，\large l^{(1)},l^{(2)},...,l^{(m)}，每一個標記點與每一個樣本資料在空間中位於相同位置。所以有：

  • 如果 x 與 l 相近 \Rightarrow f\approx exp(0)\approx1
  • 如果 x 與 l 相隔遠 \Rightarrow f=exp(-\infty )\approx0

超平面的定義如下：

\begin{align*} w^Tx+b &= (\sum_{i=1}^{m}\alpha_iy^{(i)}x^{(i)})^Tx+b \\ &= \sum_{i=1}^{m}\alpha_iy^{(i)}(x^{(i)})^Tx+b \end{align*}

可使用核函式，替換向量的內積：

\begin{align*} w^Tx+b &= (\sum_{i=1}^{m}\alpha_iy^{(i)}x^{(i)})^Tx+b \\ &= \sum_{i=1}^{m}\alpha_iy^{(i)}(x^{(i)})^Tx+b \\ &= \sum_{i=1}^{m}\alpha_iy^{(i)}k(x,x^{(i)})+b \end{align*}

高斯核函式的SVM流程可表示為：

  1. 給定資料集 (x^{(1)},y^{(1)}),...,(x^{(m)},y^{(m)})
  2. SMO演算法計算 \alpha 和 b
  3. 高斯核函式計算 k(x,x^{(1)}),...,k(x,x^{(m)})
  4. 代入 w^Tx+b ，判斷大於 0 則屬於正類，小於 0 則屬於反類

7. 實戰案例

以下將展示書中案例的程式碼段，所有程式碼和資料可以在github中下載：

7.1. 簡化版SMO演算法案例

# coding:utf-8
from numpy import *

"""
簡化版SMO演算法案例
"""


# 載入資料集
def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat, labelMat


# 根據第i個alpha，從m中選擇另一個不同的alpha
def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j


# 設定alpha的上界和下界
def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


# 簡化版本SMO演算法，並沒有確定優化的最佳alpha對，而是隨機選擇
# dataMatIn 資料集
# classLabels 類別標籤
# C 常數C
# toler 容錯率
# maxIter 最大迭代次數
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    b = 0
    m, n = shape(dataMatrix)
    # alphas的維度為m*1，m為資料集大小
    alphas = mat(zeros((m, 1)))
    iter = 0
    while (iter < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            # 預測第i個向量屬於的類別
            fXi = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b
            # 判斷這個預測的誤差
            Ei = fXi - float(labelMat[i])
            # 滿足KKT條件
            if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or ((labelMat[i] * Ei > toler) and (alphas[i] > 0)):
                # 選擇另一個向量j
                j = selectJrand(i, m)
                # 預測第j個向量屬於的類別
                fXj = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b
                # 判斷這個預測的誤差
                Ej = fXj - float(labelMat[j])
                # 建立alpha i 和alpha j的副本
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                # 判斷兩個類別標籤是否相等
                # 設定最大值為0，最小值為C
                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
                # 計算最優修改量
                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
                # 對alpha j作調整
                alphas[j] -= labelMat[j] * (Ei - Ej) / eta
                alphas[j] = clipAlpha(alphas[j], H, L)
                # 如果alpha j只是輕微的調整
                if abs(alphas[j] - alphaJold) < 0.00001:
                    print("j not moving enough")
                    continue
                # alpha a改變大小與alpha j一樣，但是方向相反
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])
                # 調整常數項b
                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 = 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))
        # 如果沒有對任何向量進行修改，則迭代次數+1，否則迭代次數清零
        # 只有在所有迭代中，都沒有對alpha進行修改，才能說明所有alpha已經最優
        if alphaPairsChanged == 0:
            iter += 1
        else:
            iter = 0
        print("iteration number: %d" % iter)
    return b, alphas


if __name__ == '__main__':
    dataArr, labelArr = loadDataSet('testSet.txt')
    b, alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)
    print(b)
    print(alphas)

7.2. 完整版SMO演算法案例

# coding:utf-8
from numpy import *

"""
完整版SMO演算法案例
"""


# 載入資料集
def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat, labelMat


# 建立儲存各引數的類
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        self.eCache = mat(zeros((self.m, 2)))


# 根據第i個alpha，從m中選擇另一個不同的alpha
def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j


# 設定alpha的上界和下界
def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


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


# 內迴圈尋找第二個alpha
# 啟發式方法，最大化步長尋找第二個alpha
def selectJ(i, oS, Ei):
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    # 誤差全域性快取賦值
    # 第0個元素表示是否賦值過
    # 第1個元素表示對應誤差
    oS.eCache[i] = [1, Ei]
    # 獲取所有賦值過的誤差全域性快取
    validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
    # 因不能選擇當前i,所以需判斷數量是否大於1
    # 若是則最大化步長, 若否則隨機選擇alpha j
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:
            if k == i:
                continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            # 選擇最大步長的
            if (deltaE > maxDeltaE):
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej


# 更新誤差全域性快取
def updateEk(oS, k):
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]


# 完整版SMO，內迴圈優化
def innerL(i, oS):
    # alpha_i的誤差
    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)):
        # 選擇alpha_j
        j, Ej = selectJ(i, oS, Ei)
        # 記錄原始的alpha
        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.X[i, :] * oS.X[j, :].T - oS.X[i, :] * oS.X[i, :].T - \
              oS.X[j, :] * oS.X[j, :].T
        if eta >= 0:
            print("eta>=0")
            return 0
        # 更新alpha_j
        oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        # 更新全域性快取
        updateEk(oS, j)
        if abs(oS.alphas[j] - alphaJold) < 0.00001:
            print("j not moving enough")
            return 0
        # 更新alpha_i
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
        # 更新全域性快取
        updateEk(oS, i)
        # 更新常數b
        b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[i, :].T - oS.labelMat[j] * (
                oS.alphas[j] - alphaJold) * oS.X[i, :] * oS.X[j, :].T
        b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[j, :].T - oS.labelMat[j] * (
                oS.alphas[j] - alphaJold) * oS.X[j, :] * oS.X[j, :].T
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
            oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]):
            oS.b = b2
        else:
            oS.b = (b1 + b2) / 2.0
        return 1
    else:
        return 0


# 完整版SMO，外迴圈
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
    oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler)
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        # 完整遍歷
        if entireSet:
            for i in range(oS.m):
                alphaPairsChanged += innerL(i, oS)
                print("fullSet, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged))
            iter += 1
        # 非邊界值遍歷
        else:
            nonBoundIs = 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
        elif alphaPairsChanged == 0:
            entireSet = True
        print("iteration number: %d" % iter)
    return oS.b, oS.alphas


# 根據alpha計算權重係數w
def calcWs(alphas, dataArr, classLabels):
    X = mat(dataArr)
    labelMat = mat(classLabels).transpose()
    m, n = shape(X)
    w = zeros((n, 1))
    for i in range(m):
        w += multiply(alphas[i] * labelMat[i], X[i, :].T)
    return w


if __name__ == '__main__':
    dataArr, labelArr = loadDataSet('testSet.txt')
    b, alphas = smoP(dataArr, labelArr, 0.6, 0.001, 40)
    w = calcWs(alphas, dataArr, labelArr)
    print(w)

7.3. 高斯核函式案例

# coding:utf-8
from numpy import *

"""
高斯核函式案例
"""


# 載入資料集
def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat, labelMat


# 核函式轉換
def kernelTrans(X, A, kTup):
    m, n = shape(X)
    K = mat(zeros((m, 1)))
    # 線性核心
    if kTup[0] == 'lin':
        K = X * A.T
    # 高斯核心
    elif kTup[0] == 'rbf':
        for j in range(m):
            deltaRow = X[j, :] - A
            K[j] = deltaRow * deltaRow.T
        K = exp(K / (-1 * kTup[1] ** 2))
    else:
        raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
    return K


# 儲存引數資料結構
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler, kTup):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0
        self.eCache = mat(zeros((self.m, 2)))
        self.K = mat(zeros((self.m, self.m)))
        for i in range(self.m):
            self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)


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


# 根據第i個alpha，從m中選擇另一個不同的alpha
def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j


# 設定alpha的上界和下界
def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


# 內迴圈尋找第二個alpha
# 啟發式方法，最大化步長尋找第二個alpha
def selectJ(i, oS, Ei):
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    # 誤差全域性快取賦值
    # 第0個元素表示是否賦值過
    # 第1個元素表示對應誤差
    oS.eCache[i] = [1, Ei]
    # 獲取所有賦值過的誤差全域性快取
    validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
    # 因不能選擇當前i,所以需判斷數量是否大於1
    # 若是則最大化步長, 若否則隨機選擇alpha j
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:
            if k == i:
                continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            # 選擇最大步長的
            if (deltaE > maxDeltaE):
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej


# 更新誤差全域性快取
def updateEk(oS, k):
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]


# 完整版SMO，內迴圈優化
def innerL(i, oS):
    # alpha_i的誤差
    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)):
        # 選擇alpha_j
        j, Ej = selectJ(i, oS, Ei)
        # 記錄原始的alpha
        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]
        if eta >= 0:
            print("eta>=0")
            return 0
        # 更新alpha_j
        oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        # 更新全域性快取
        updateEk(oS, j)
        if abs(oS.alphas[j] - alphaJold) < 0.00001:
            print("j not moving enough")
            return 0
        # 更新alpha_i
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
        # 更新全域性快取
        updateEk(oS, i)
        # 更新常數b
        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)):
    oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup)
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        # 完整遍歷
        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:
            nonBoundIs = 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
        elif alphaPairsChanged == 0:
            entireSet = True
        print("iteration number: %d" % iter)
    return oS.b, oS.alphas


# 文字轉換為向量
def img2vector(filename):
    returnVect = zeros((1, 1024))
    fr = open(filename)
    for i in range(32):
        lineStr = fr.readline()
        for j in range(32):
            returnVect[0, 32 * i + j] = int(lineStr[j])
    return returnVect


# 載入資料
def loadImages(dirName):
    from os import listdir
    hwLabels = []
    trainingFileList = listdir(dirName)
    m = len(trainingFileList)
    trainingMat = zeros((m, 1024))
    for i in range(m):
        fileNameStr = trainingFileList[i]
        fileStr = fileNameStr.split('.')[0]
        classNumStr = int(fileStr.split('_')[0])
        if classNumStr == 9:
            hwLabels.append(-1)
        else:
            hwLabels.append(1)
        trainingMat[i, :] = img2vector('%s/%s' % (dirName, fileNameStr))
    return trainingMat, hwLabels


# 測試函式
def testDigits(kTup=('rbf', 10)):
    dataArr, labelArr = loadImages('trainingDigits')
    b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    # 獲取所有支援向量
    svInd = nonzero(alphas.A > 0)[0]
    sVs = datMat[svInd]
    labelSV = labelMat[svInd]
    print("there are %d Support Vectors" % shape(sVs)[0])
    m, n = shape(datMat)
    errorCount = 0
    for i in range(m):
        # 構建核函式
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        # 利用核函式進行預測
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]): errorCount += 1
    print("the training error rate is: %f" % (float(errorCount) / m))
    dataArr, labelArr = loadImages('testDigits')
    errorCount = 0
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    m, n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
        if sign(predict) != sign(labelArr[i]): errorCount += 1
    print("the test error rate is: %f" % (float(errorCount) / m))


if __name__ == '__main__':
    testDigits()

如果這篇文章對你有一點小小的幫助，請給個關注喔~我會非常開心的~

 

 

 

相關推薦

機器學習實戰（五）支援向量機SVM（Support Vector Machine）

目錄 0. 前言 1. 尋找最大間隔 2. 拉格朗日乘子法和KKT條件 3. 鬆弛變數 4. 帶鬆弛變數的拉格朗日乘子法和KKT條件 5. 序列最小優化SMO（Sequential Minimal Optimiz

[完]機器學習實戰 第六章 支援向量機（Support Vector Machine）

[參考] 機器學習實戰（Machine Learning in Action） 本章內容 支援向量機（Support Vector Machine）是最好的現成的分類器，“現成”指的是分類器不加修改即可直接使用。基本形式的SVM分類器就可得到低錯

機器學習實戰-55:支援向量機分類演算法(Support Vector Machine)

支援向量機分類演算法 支援向量機(Support Vector Machine)分類演算法屬於監督學習演算法。常用分類演算法包括：邏輯迴歸(Logistic Regression, LR)、K最近鄰(k-Nearest Neighbor, KNN)、樸素貝葉斯

機器學習實戰第六章支援向量機照葫蘆畫瓢演算法實踐

支援向量機簡要介紹 一些概念： 1.分隔超平面：在二維中直觀來說就是將資料集分隔開來的直線，三維中則是一個平面。觸類旁通。 2.超平面：分類的決策邊界，分佈在超平面一側的所有資料都屬於某個類別，另一側屬於另一個。 3.支援向量：離分隔超平面最近的那些

Spark機器學習系列之13： 支援向量機SVM

C−SVM基本公式推導過程 下面摘抄一小部分內容（不考慮推導細節的話，基本上能理解C-SVM方法推導的整個流程）. 我們用一個超平面劃分圖中對圖中的兩類資料進行分類，超平面寫成f(x)=wTx+b=0,線上性可分的情況下，我們能找到一

機器學習(四)：通俗理解支援向量機SVM及程式碼實踐

[上一篇文章](https://mp.weixin.qq.com/s/cEbGM0_Lrt8elfubxSF9jg)我們介紹了使用邏輯迴歸來處理分類問題，本文我們講一個更強大的分類模型。本文依舊側重程式碼實踐，你會發現我們解決問題的手段越來越豐富，問題處理起來越來越簡單。 支援向量機(Support V

吳恩達機器學習（第十三章）---支援向量機SVM

一、優化目標 邏輯迴歸中的代價函式：  畫出兩種情況下的函式影象可得： y=1: 我們找一條折線來近似表示這個函式影象 y=0：    我們用這兩條折線來近似表示原來的曲線函式可得新的代價函式（假設-log(h(x))為,-log(1

機器學習與深度學習系列連載： 第一部分 機器學習（九）支援向量機2（Support Vector Machine）

另一種視角定義SVM：hinge Loss +kennel trick SVM 可以理解為就是hingle Loss和kernel 的組合 1. hinge Loss 還是讓我們回到二分類的問題，為了方便起見，我們y=1 看做是一類，y=-1 看做是另一類

機器學習（三）——支援向量機SVM

SVM探討 目錄 SVM探討 SVM演算法 根據處理問題的複雜度，SVM 可由簡到繁分為三種： 線性可分支援向量機：硬間隔最大化。 線性支援向量機：資料分佈近似線性可分，可通過軟間隔最大化(懲罰因子，鬆弛變數)來線性分隔樣本點。 非線性支援向量機：通過核函式提升

機器學習——支援向量機SVM（一）

在之前做數學建模的時候就有使用過支援向量機的有關知識，但是對於支援向量機的理解一直不是很深刻，尤其是數學推導部分還是存在一些問題。在最近看周志華西瓜書的過程中同樣發現這一問題，西瓜書中對支援向量機的講解部分不是很詳細，所以我又查找了其他的資料。 支援向量機是一種原創性（非組

機器學習之路（一）---支援向量機SVM

##什麼是支援向量機(SVM)？ 支援向量機 (SVM) 是一個類分類器，正式的定義是一個能夠將不同類樣本在樣本空間分隔的超平面。 換句話說，給定一些標記(label)好的訓練樣本 (監督式學習), SVM演算法輸出一個最優化的分隔超平面。 如何來界定一個超

機器學習——支援向量機SVM（Support Vector Machine）（下）

1、SVM演算法特徵 （1）訓練好的模型的演算法複雜度是由支援向量的個數決定的，而不是由資料的維度決定。所以，SVM不太容易產生overfitting。 （2）SVM訓練出來的模型完全依賴於支援向量（

資料探勘十大演算法——支援向量機SVM（五）：線性支援迴歸

首先感謝“劉建平pinard”的淵博知識以及文中詳細準確的推導！！！ 支援向量機原理SVM系列文章共分為5部分： （一）線性支援向量機 （二）線性支援向量機的軟間隔最大化模型 （三）線性不可分支援向量機與核函式 （四）SMO演算法原理 （五）線性支援迴歸

【深度學習基礎-06】支援向量機SVM（下）-線性不可分

 1 SVM的優點                                      

【深度學習基礎-05】支援向量機SVM（上）-線性可分

Support Vector Machine 目錄 1背景 2 機器學習的一般框架 3 什麼是超平面 4 線性可區分（linear separatable）和線性不可區分（linear inseparatable） 5 如何計算超平面以及舉例 1背景 Vladim

學習筆記——支援向量機svm（2）對偶問題

回顧 上一篇已經講到了將找一個最好的分割超平面轉化為上面的公式並且用二次規劃來求解的問題。但是還是存在一個問題，這個演算法維度（w的維度）不是跟資料量相關的，而是跟資料內在的vc維度有關的，所以當

機器學習 scikit-learn3 模型實踐 - 支援向量機和決策樹

支援向量機和決策樹 - 目錄 1 簡介 1.1 程式碼下載 1.2 程式碼使用方法 3 核心程式碼說明 3.1 模型配置 3.2 模型訓練 3.3 輸出結果 3.3.1 Linea

支援向量機SVM（三）：基於核函式的非線性SVM

前言 線性分類SVM是一種線性分類非常有效的方法，若分類問題是非線性，需要考慮對資料進行 空間變換，將非線性分類問題轉變為線性分類問題，使非線性SVM的學習轉變為線性SVM的學習。 若使用對映ϕ\phiϕ將例項從原空間X\mathcal XX（一般為歐式空

資料探勘十大演算法——支援向量機SVM（一）：線性支援向量機

首先感謝“劉建平pinard”的淵博知識以及文中詳細準確的推導！！！ 本文轉自“劉建平pinard”，原網址為：http://www.cnblogs.com/pinard/p/6097604.html。 支援向量機原理SVM系列文章共分為5部分： （一）線性支援向量機

【機器學習基礎】軟間隔支援向量機

引言 在上一小節中，我們介紹了核支援向量機。於是，不管是簡單的問題還是複雜的問題，我們都可以做得到。 然而，像高斯核的方法太複雜了，可能造成過擬合的問題。導致過擬合的現象的原因有可能是你選擇特徵轉換太強大了，導致無法用最大間隔的方法控制模型的複雜度，還有一個