機器學習演算法實踐-SVM中的SMO演算法

01-30

前言

前兩篇關於SVM的文章分別總結了SVM基本原理和核函數以及軟間隔原理，本文我們就針對前面推導出的SVM對偶問題的一種高效的優化方法-序列最小優化演算法(Sequential Minimal Optimization, SMO)的原理進行總結並進行相應的Python實現。

坐標上升演算法(Coordinate Ascent)

在SMO演算法之前，還是需要總結下坐標上升演算法，因為SMO演算法的思想與坐標上升演算法的思想類似。

坐標上升演算法每次通過更新多元函數中的一維，經過多次迭代直到收斂來達到優化函數的目的。簡單的講就是不斷地選中一個變數做一維最優化直到函數達到局部最優點。

假設我們需要求解的問題形式為(類似我們SVM的對偶形式):

$max limits_{alpha} W(alpha_{1}, alpha_{2}, …, alpha_{N})$

演算法過程偽碼:

例子

若我們的優化目標為一個二元函數:

$arg min limits_{x_{1}, x_{2}} f(x_{1}, x_{2}) = -x_{1}^{2} - 3x_{2}^{2} + 2x_{1}x_{2} + 6$

我們先給一個 $(x_{1}, x_{2})$ 的初值然後開始迭代。

先固定 $x_1$ ，把 $f$ 看做 $x_2$ 的一元函數求最優值，可以簡單的進行求導獲取解析解:
$frac{partial{f}}{partial{x_{1}}} = -2x_{1} + 2x_{2} = 0 rightarrow x_{1} = x_{2}$
在固定x2x2, 把ff看成x1x1的一元函數求最優值，得到x1x1的解析解:
$frac{partial{f}}{partial{x_{2}}} = -6x_{2} + 2x_{1} rightarrow x_{2} = frac{1}{3}x_{1}$

按照上面兩個過程不斷交替的優化 $x_1$ 和 $x_2$ ，直到函數收斂。

通過下面的圖就可以看出，優化的過程，因為每次只優化一個變數，每次迭代的方向都是沿著坐標軸方向的。

因為每次只是做一維優化，所以每個循環中的優化過程的效率是很高的, 但是迭代的次數會比較多。

序列最小優化演算法(SMO)

SMO演算法介紹

SMO的思想類似坐標上升演算法，我們需要優化一系列的αα的值，我們每次選擇盡量少的 $alpha$ 來優化，不斷迭代直到函數收斂到最優值。

來到SVM的對偶問題上，對偶形式:

$arg max limits_{alpha} sum_{i=1}^{N} alpha_{i} - frac{1}{2}sum_{i=1}^{N}sum_{j=1}^{N}y_{i}y_{j}alpha_{i}alpha_{j}langle x_{i}, x_{j} rangle$

subject to $alpha_{i} ge 0$ , $sum_{i=1}^{N}alpha_{i}y_{i}=0$

其中我們需要對 $(alpha_{1}, alpha_{2}, …, alpha_{N})$ 進行優化，但是這個凸二次優化問題的其他求解演算法的複雜度很高，但是Platt提出的SMO演算法可以高效的求解上述對偶問題，他把原始問題的求解 $N$ 個參數二次規劃問題分解成多個二次規劃問題求解，每個字問題只需要求解2各參數，節省了時間成本和內存需求。

與坐標上升演算法不同的是，我們在SMO演算法中我們每次需要選擇一對變數 $(alpha_i, alpha_j)$ , 因為在SVM中，我們的 $alpha$ 並不是完全獨立的，而是具有約束的:

$sum_{i=1}^{N}alpha_{i}y_{i} = 0$

因此一個 $alpha$ 改變，另一個也要隨之變化以滿足條件。

SMO演算法原理

獲得沒有修剪的原始解

假設我們選取的兩個需要優化的參數為 $alpha_1, alpha_2$ , 剩下的 $alpha_{3}, alpha_{4}, …, alpha_{N}$ 則固定，作為常數處理。將SVM優化問題進行展開就可以得到(把與 $alpha_{1}, alpha_{2}$ 無關的項合併成常數項 $C$ ):

$W(alpha_{1}, alpha_{2}) = alpha_{1} + alpha_{2} - frac{1}{2}K_{1,1}y_{1}^{2}alpha_{1}^{2} - frac{1}{2}K_{2,2}y_{2}^{2}alpha_{2}^{2} - K_{1,2}y_{1}y_{2}alpha_{1}alpha_{2} - y_{1}alpha_{1}sum_{i=3}^{N}alpha_{i}y_{i}K_{i,1} - y_{2}alpha_{2}sum_{i=3}^{N}alpha_{i}y_{i}K_{i, 2} + C$

於是就是一個二元函數的優化:

$arg max limits_{alpha_{1}, alpha_{2}} W(alpha_{1}, alpha_{2})$

根據約束條件 $sum_{i=1}^{N}alpha_{i}y_{i} = 0$ 可以得到 $alpha_1$ 與 $alpha_2$ 的關係:

$alpha_{1}y_{1} + alpha_{2}y_{2} = -sum_{i=3}^{N}alpha_{i}y_{i} = zeta$

兩邊同時乘上 $y_1$ , 由於 $y_{i}y_{i} = 1$ 得到:

$alpha_{1} = zeta y_{1} - alpha_{2}y_{1}y_{2}$

令 $v_{1} = sum_{i=3}^{N}alpha_{i}y_{i}K_{i, 1}$ , $v_{2} = sum_{i=3}^{N}alpha_{i}y_{i}K_{i, 2}$ ， $alpha_1$ 的表達式代入得到:

$W(alpha_{2}) = -frac{1}{2}K_{1, 1}(zeta - alpha_{2}y_{2})^{2} - frac{1}{2}K_{2, 2}alpha_{2}^{2} - y_{2}(zeta - alpha_{2}y_{2})alpha_{2}K_{1, 2} - v_{1}(zeta - alpha_{2}y_{2}) - v_{2}y_{2}alpha_{2} + alpha_{1} + alpha_{2} + C$

後面我們需要對這個一元函數進行求極值， $W$ 對 $alpha$ 的一階導數為0得到:

$frac{partial{W(alpha_{2})}}{partial{alpha_{2}}} = -(K_{1, 1} + K_{2, 2} - 2K_{1, 2})alpha_{2} + K_{1, 1}zeta y_{2} - K_{1, 2}zeta y_{2} + v_{1}y_{2} - v_{2}y_{2} - y_{1}y_{2} + y_{2}^{2} = 0$

下面我們稍微對上式進行下變形，使得 $alpha_{2}^{new}$ 能夠用更新前的 $alpha_{2}^{old}$ 表示，而不是使用不方便計算的 $zeta$ 。

因為SVM對數據點的預測值為: $f(x) = sum_{i=1}^{N}alpha_{i}y_{i} K(x_{i}, x) + b$

則 $v_1$ 以及 $v_2$ 的值可以表示成:

$v_{1} = sum_{i=3}^{N}alpha_{i}y_{i}K_{1, i} = f(x_{1}) - alpha_{1}y_{1}K_{1, 1} - alpha_{2}y_{2}K_{1, 2} - b$

$v_{2} = sum_{i=3}^{N}alpha_{i}y_{i}K_{2, i} = f(x_{2}) - alpha_{1}y_{1}K_{1, 2} - alpha_{2}y_{2}K_{2, 2} - b$

已知 $alpha_{1} = (zeta - alpha_{2}y_{2})y_{2}$ , 可得到:

$v_{1} - v_{2} = f(x_{1}) - f(x_{2}) - K_{1, 1}zeta + K_{1, 2}zeta + (K_{1, 1} + K_{2, 2} - 2K_{1, 2})alpha_{2}y_{2}$

將 $v_{1} - v_{2}$ 的表達式代入到 $frac{partial{W(alpha_{2})}}{partial{alpha_{2}}}$ 中可以得到: $frac{partial{W(alpha_{2})}}{partial{alpha_{2}}} = -(K_{1, 1}) + K_{2, 2} - 2K_{1, 2})alpha_{2}^{new} +(K_{1, 1}) + K_{2, 2} - 2K_{1, 2})alpha_{2}^{old} + y_{2}left[ y_{2} - y_{1} + f(x_{1}) - f(x_{2}) right]$

我們記 $E_i$ 為SVM預測值與真實值的誤差: $E_{i} = f(x_{i}) - y_{i}$

令 $eta = K_{1, 1} + K_{2, 2} - 2K_{1, 2}$ 得到最終的一階導數表達式:

$frac{partial{W(alpha_{2})}}{partial{alpha_{2}}} = -eta alpha_{2}^{new} + eta alpha_{2}^{old} + y_{2}(E_{1} - E_{2}) = 0$

得到:

$alpha_{2}^{new} = alpha_{2}^{old} + frac{y_{2}(E_{1} - E_{2})}{eta}$

這樣我們就得到了通過舊的 $alpha_2$ 獲取新的 $alpha_2$ 的表達式, $alpha_{1}^{new}$ 便可以通過 $alpha_{2}^{new}$ 得到。

對原始解進行修剪

上面我們通過對一元函數求極值的方式得到的最優 $alpha_{i}, alpha_{j}$ 是未考慮約束條件下的最優解，我們便更正我們上部分得到的 $alpha_{2}^{new}$ 為 $alpha_{2}^{new, unclipped}$ , 即:

$alpha_{2}^{new, unclipped} = alpha_{2}^{old} + frac{y_{2}(E_{1} - E_{2})}{eta}$

但是在SVM中我們的 $alpha_i$ 是有約束的，即:

$alpha_{1}y_{1} + alpha_{2}y_{2} = -sum_{i=3}^{N}alpha_{i}y_{i} = zeta$

$0 le alpha_{i} le C$

此約束為方形約束(Bosk constraint), 在二維平面中我們可以看到這是個限制在方形區域中的直線（見下圖）。

(如左圖) 當 $y_{1} ne y_{2}$ 時，線性限制條件可以寫成: $alpha_{1} - alpha_{2} = k$ ，根據 $k$ 的正負可以得到不同的上下界，因此統一表示成:

下界: $L = max(0, alpha_{2}^{old} - alpha_{1}^{old})$
上界: $H = min(C, C + alpha_{2}^{old} - alpha_{1}^{old})$

(如右圖) 當 $y_{1} = y_{2}$ 時，限制條件可寫成: $alpha_{1} + alpha_{2} = k$ , 上下界表示成:

下界: $L = max(0, alpha_{1}^{old} + alpha_{2}^{old} - C)$
上界: $H = min(C, alpha_{2}^{old} + alpha_{1}^{old})$

根據得到的上下界，我們可以得到修剪後的 $alpha_{2}^{new}$ :

$alpha_{2}^{new} = begin{cases} H & alpha_{2}^{new, unclipped} > H alpha_{2}^{new, unclipped} & L le alpha_{2}^{new, unclipped} le H L & alpha_{2}^{new, unclipped} < L end{cases}$

得到了 $alpha_{2}^{new}$ 我們便可以根據 $alpha_{1}^{old}y_{1} + alpha_{2}^{old}y_{2} = alpha_{1}^{new}y_{1} + alpha_{2}^{new}y_{2}$ 得到 $alpha_{1}^{new}$ :

$alpha_{1}^{new} = alpha_{1}^{old} + y_{1}y_{2}(alpha_{2}^{old} - alpha_{2}^{new})$

OK，這樣我們就知道如何將選取的一對 $alpha_{i}, alpha_{j}$ 進行優化更新了。

更新閾值b

當我們更新了一對 $alpha_{i}, alpha_{j}$ 之後都需要重新計算閾值 $b$ ，因為 $b$ 關係到我們 $f(x)$ 的計算，關係到下次優化的時候誤差 $E_i$ 的計算。

為了使得被優化的樣本都滿足KKT條件，

當 $alpha_{1}^{new}$ 不在邊界，即 $0 < alpha_{1}^{new} < C$ , 根據KKT條件可知相應的數據點為支持向量，滿足 $y_{1}(w^{T} + b) = 1$ , 兩邊同時乘上 $y_1$ 得到 $sum_{i=1}^{N}alpha_{i}y_{i}K_{i, 1} + b = y_{1}$ , 進而得到 $b_{1}^{new}$ 的值:

$b_{1}^{new} = y_{1} - sum_{i=3}^{N}alpha_{i}y_{i}K_{i, 1} - alpha_{1}^{new}y_{1}K_{1, 1} - alpha_{2}^{new}y_{2}K_{2, 1}$

其中上式的前兩項可以寫成:

$y_{1} - sum_{i=3}^{N}alpha_{i}y_{i}K_{i, 1} = -E_{1} + alpha_{1}^{old}y_{1}K_{1, 1} + alpha_{2}^{old}y_{2}K_{2, 1} + b^{old}$

當 $0 < alpha_{2}^{new} < C$ , 可以得到bnew2b2new的表達式(推導同上):

$b_{2}^{new} = -E_{2} - y_{1}K_{1, 2}(alpha_{1}^{new} - alpha_{1}^{old}) - y_{2}K_{2, 2}(alpha_{2}^{new} - alpha_{2}^{old}) + b^{old}$

當 $b_1$ 和 $b_2$ 都有效的時候他們是相等的, 即 $b^{new} = b_{1}^{new} = b_{2}^{new}$ 。

當兩個乘子 $alpha_{1}, alpha_{2}$ 都在邊界上，且 $L ne H$ 時， $b1, b2$ 之間的值就是和KKT條件一直的閾值。SMO選擇他們的中點作為新的閾值:

$b^{new} = frac{b_{1}^{new} + b_{2}^{new}}{2}$

簡化版SMO演算法實現

這裡我主要針對SMO中已選取的一對 $alpha_{i}, alpha_{j}$ 值的優化過程進行下Python實現，其中 $alpha_{i}, alpha_{j}$ 的選取直接使用傻瓜的遍歷方式，並使用100數據點進行訓練。

首先是一些輔助函數，用來幫助載入數據，修剪 $alpha$ 的值以及隨機選取 $alpha$

def load_data(filename):n dataset, labels = [], []n with open(filename, r) as f:n for line in f:n x, y, label = [float(i) for i in line.strip().split()]n dataset.append([x, y])n labels.append(label)n return dataset, labelsndef clip(alpha, L, H):n 修建alpha的值到L和H之間.n n if alpha < L:n return Ln elif alpha > H:n return Hn else:n return alphandef select_j(i, m):n 在m中隨機選擇除了i之外剩餘的數n n l = list(range(m))n seq = l[: i] + l[i+1:]n return random.choice(seq)n

為了能在最後繪製SVM分割線，我們需要根據獲取的 $alpha$ ，數據點以及標籤來獲取 $w$ 的值:

def get_w(alphas, dataset, labels):n 通過已知數據點和拉格朗日乘子獲得分割超平面參數wn n alphas, dataset, labels = np.array(alphas), np.array(dataset), np.array(labels)n yx = labels.reshape(1, -1).T*np.array([1, 1])*datasetn w = np.dot(yx.T, alphas)n return w.tolist()n

簡化版SMO演算法的實現,即便沒有添加啟發式的 $alpha$ 選取，SMO演算法仍然有比較多的公式需要實現，我本人按照上文的推導進行實現的時候就因為寫錯了一個下標演算法一直跑不出想要的結果。

此實現主要包含兩重循環，外層循環是控制最大迭代步數，此迭代步數是在每次有優化一對αα之後進行判斷所選取的 $alpha$ 是否已被優化，如果沒有則進行加一，如果連續max_iter步數之後仍然沒有 $alpha$ 被優化，則我們就認為所有的 $alpha$ 基本已經被優化，優化便可以終止了.

def simple_smo(dataset, labels, C, max_iter):n 簡化版SMO演算法實現，未使用啟發式方法對alpha對進行選擇.n :param dataset: 所有特徵數據向量n :param labels: 所有的數據標籤n :param C: 軟間隔常數, 0 <= alpha_i <= Cn :param max_iter: 外層循環最大迭代次數n n dataset = np.array(dataset)n m, n = dataset.shapen labels = np.array(labels)n # 初始化參數n alphas = np.zeros(m)n b = 0n it = 0n def f(x):n "SVM分類器函數 y = w^Tx + b"n # Kernel function vector.n x = np.matrix(x).Tn data = np.matrix(dataset)n ks = data*xn # Predictive value.n wx = np.matrix(alphas*labels)*ksn fx = wx + bn return fx[0, 0]nn while it < max_iter:n pair_changed = 0n for i in range(m):n a_i, x_i, y_i = alphas[i], dataset[i], labels[i]n fx_i = f(x_i)n E_i = fx_i - y_in j = select_j(i, m)n a_j, x_j, y_j = alphas[j], dataset[j], labels[j]n fx_j = f(x_j)n E_j = fx_j - y_jn K_ii, K_jj, K_ij = np.dot(x_i, x_i), np.dot(x_j, x_j), np.dot(x_i, x_j)n eta = K_ii + K_jj - 2*K_ijn if eta <= 0:n print(WARNING eta <= 0)n continuen # 獲取更新的alpha對n a_i_old, a_j_old = a_i, a_jn a_j_new = a_j_old + y_j*(E_i - E_j)/etan # 對alpha進行修剪n if y_i != y_j:n L = max(0, a_j_old - a_i_old)n H = min(C, C + a_j_old - a_i_old)n else:n L = max(0, a_i_old + a_j_old - C)n H = min(C, a_j_old + a_i_old)n a_j_new = clip(a_j_new, L, H)n a_i_new = a_i_old + y_i*y_j*(a_j_old - a_j_new)n if abs(a_j_new - a_j_old) < 0.00001:n #print(WARNING alpha_j not moving enough)n continuen alphas[i], alphas[j] = a_i_new, a_j_newn # 更新閾值bn b_i = -E_i - y_i*K_ii*(a_i_new - a_i_old) - y_j*K_ij*(a_j_new - a_j_old) + bn b_j = -E_j - y_i*K_ij*(a_i_new - a_i_old) - y_j*K_jj*(a_j_new - a_j_old) + bn if 0 < a_i_new < C:n b = b_in elif 0 < a_j_new < C:n b = b_jn else:n b = (b_i + b_j)/2n pair_changed += 1n print(INFO iteration:{} i:{} pair_changed:{}.format(it, i, pair_changed))n if pair_changed == 0:n it += 1n else:n it = 0n print(iteration number: {}.format(it))n return alphas, bn

Ok, 下面我們就用訓練數據對SVM進行優化, 並對最後優化的分割線以及數據點進行可視化

if __main__ == __name__:n # 載入訓練數據n dataset, labels = load_data(testSet.txt)n # 使用簡化版SMO演算法優化SVMn alphas, b = simple_smo(dataset, labels, 0.6, 40)n # 分類數據點n classified_pts = {+1: [], -1: []}n for point, label in zip(dataset, labels):n if label == 1.0:n classified_pts[+1].append(point)n else:n classified_pts[-1].append(point)n fig = plt.figure()n ax = fig.add_subplot(111)n # 繪製數據點n for label, pts in classified_pts.items():n pts = np.array(pts)n ax.scatter(pts[:, 0], pts[:, 1], label=label)n # 繪製分割線n w = get_w(alphas, dataset, labels)n x1, _ = max(dataset, key=lambda x: x[0])n x2, _ = min(dataset, key=lambda x: x[0])n a1, a2 = wn y1, y2 = (-b - a1*x1)/a2, (-b - a1*x2)/a2n ax.plot([x1, x2], [y1, y2])n # 繪製支持向量n for i, alpha in enumerate(alphas):n if abs(alpha) > 1e-3:n x, y = dataset[i]n ax.scatter([x], [y], s=150, c=none, alpha=0.7,n linewidth=1.5, edgecolor=#AB3319)n plt.show()n

優化最後我們可以看到針對100個數據的 $alpha$ 只有少部分是大於零的，即對應的數據點就是支持向量:

為了能直觀的顯示支持向量，我將其標註了出來，最終可視化的效果如下圖:

總結

本文從坐標上升演算法開始介紹，並對SMO演算法的原理進行了簡單的推導，針對SMO演算法中對αα對的優化並使用了Python進行了簡化版的SMO實現，並針對小數據集進行了訓練得到了對應優化後的SVM。

實現代碼以及訓練數據鏈接: https://github.com/PytLab/MLBox/tree/master/svm

參考

優化演算法——坐標上升法
支持向量機系列（5）——SMO演算法解對偶問題
《Machine Learning in Action》