機器學習：支持向量機2

本文來自同步博客。

P.S. 不知道怎麼顯示數學公式以及更好的排版內容。因此若是以爲文章下面格式亂的話請自行跳轉到上述連接。後續我將再也不對數學公式進行截圖，畢竟行內公式截圖的話排版會很亂。看原博客地址會有更好的體驗。

上一篇文章介紹了機器學習中支持向量機的基本原理，而且在文章末尾介紹了一種利用Python求解二項規劃問題極值的方法。這篇文章我將利用這種方法一步步求解上文中說起的$-\vec{\alpha}-$、$-\vec{w}-$、$-b-$，藉此複習和驗證支持向量機的知識點。

數據

下面看一組測試數據：

data = {
    '+': [
        [1, 7],
        [2, 8],
        [3, 8],
        [2, 6.5]
    ],
    '-': [
        [5, 1],
        [6, -1],
        [7, 3]
    ]
}

數據data是擁有兩種已經分好類的數據，每種類型的數據的元素都是二維向量，能夠在笛卡爾座標系中表示。

依照上一篇文章講述的原理，咱們須要利用這些數據求解一個$-\vec{\alpha}-$向量。也就是咱們須要求解使得二項規劃方程值最小時的$-\vec{\alpha}-$向量：
$$
F(\vec{\alpha}) = \frac{1}{2}\vec{\alpha}_{T}H\vec{\alpha} + \vec{c}\vec{\alpha} + c_0， \vec{y}^{T}\vec{\alpha} = 0， \vec{\alpha} \ge 0
$$

很明顯，在支持向量機中，$-c_0 = 0-$。

參數求解

首先利用輸入的測試data準備上述方程中出現的變量$-H，c，c_0-$。參考下面代碼：

def parseXYC(d):
    X = []
    y = []
    c = []
    for _, v in enumerate(d['+']):
        X.append(np.array(v))
        y.append(1)
        c.append(-1)
    for _, v in enumerate(d['-']):
        X.append(np.array(v))
        y.append(-1)
        c.append(-1)
    return X, y, c, 0

X, y, c, c0 = parseXYC(data)

parseXYC函數把data格式化成$-X, y, c, c_0-$。

而後計算$-H-$矩陣的值。比較簡單，一行代碼就能夠獲得：

H = np.array([y[i] * y[j] * np.dot(X[i], X[j]) for i in range(len(X)) for j in range(len(X))]).reshape(len(X), len(X))

求解$-\vec{\alpha}-$

全部數據都準備好了，接下來就是帶入optimize.minimize函數中計算結果。

這裏有幾個超出本文描述範圍的難點須要簡單說起一下：

1. optimize.minimize函數求解二項規劃使用的SLSQP方法既須要用到二項方程的雅各比導函數，也須要用到約束條件函數的雅各比導函數。不清楚這點致使我在測試過程當中一直沒法求解到正確的值。
2. 不等式約束條件$-\vec{\alpha} \ge 0-$沒法在做爲約束條件參數constraints傳遞給optimize.minimize函數。我猜想是由於我構造的不等式參數是錯誤的，所以沒法讓不等式約束條件生效。我尚沒法解決這個問題，但願瞭解該問題的同窗能留言賜教。做爲一種補救方法，我利用邊界約束參數bounds描述$-\vec{\alpha} \ge 0-$這個不等式。
3. 求解出來的$-\vec{\alpha}-$向量中，部分應該爲0的元素沒法完整精確到0。我觀察測試結果總結出來的精確度應該在1e-16，所以在負16次方這個精確度下的值我都假設它就是0。通過繪圖驗證，我發現這個假設是合理的。

下面開代碼實現：

# 定義二項規劃方程fun及其雅各比方程jac
def fun(x, sign=1.):
    return sign * (0.5 * np.dot(x.T, np.dot(H, x))+ np.dot(c, x) + c0)
def jac(x, sign=1.):
    return sign * (np.dot(x.T, H) + c)

# 定義等式約束條件方程feq及其雅各比方程jeq
def feq(x):
    return np.dot(y, x)
def jeq(x):
    return np.array(y)

# 生成相關參數
diff = 1e-16
bounds = [(0, None) for _ in range(len(y))] # x >= 0
constraints = [{ 'type': 'eq', 'fun': feq, 'jac': jeq }]# y*x = 0
options = { 'ftol': diff, 'disp': True }
guess = np.array([0 for _ in range(len(X))])

# 計算結果
res_cons = optimize.minimize(fun, guess, method='SLSQP', jac=jac, bounds=bounds, constraints=constraints, options=options)
alpha = [ 0 if abs(x - 0) <= diff else x for x in res_cons.x ]

# 輸出結果與校驗y*alpha的值是否爲0
print('raw alpha: ', res_cons.x)
print('fmt alpha: ', alpha)
print('check y*alpha: ', 'is 0'if (abs(np.dot(y, res_cons.x) - 0) < diff ) else 'is not 0')

求解$-\vec{w}-$和$-b-$

# 計算w = sum(xi*yi*Xi)
w = np.sum([ np.array([0, 0]) if alpha[i] == 0 else (alpha[i] * y[i] * X[i]) for i in range(len(alpha))], axis=0)
print('w: ', w)

# 計算b，對support vector有：yi(w*xi + b) = 1，既有：b = 1/yi - w*xi
B = [( 0 if alpha[i] == 0 else ( 1 / y[i] - np.dot(w, X[i]) ) ) for i in range(len(alpha))]
B = list(filter(lambda x: x != 0, B))
b = 0 if len(B) <= 0 else B[0]
print('b: ', b)

至此支持向量機的參數求解過程完畢。

運行結果以下圖所示：

 

繪圖

最後把數據繪製成圖像。

limit = 11
plt.xlim(-2, limit)
plt.ylim(-2, limit)
# 繪製數據點
[plt.scatter(X[i][0],X[i][1], s=100, color=('r' if y[i] > 0 else 'y')) for i in range(len(X))]
# 繪製分割超平面L： wx + b = 0
plt.plot([i for i in range(limit)], [(-b - w[0]*i)/w[1] for i in range(limit)])
# 繪製上下邊： wx + b = 1/-1
plt.plot([i for i in range(limit)], [(1-b - w[0]*i)/w[1] for i in range(limit)])
plt.plot([i for i in range(limit)], [(-1-b - w[0]*i)/w[1] for i in range(limit)])
plt.show()

效果以下圖。其中紅點爲'+'樣本，綠點爲'-'樣本。中間的藍色線爲分類的標準線。邊界線，即紅色線和綠色線分別穿過各自類別中最靠近分類標準線的點。這些點就是支持向量，只有這些向量所對應的$-vec{\alpha}-$份量才爲非零值。

 

 

本文源碼