接下來，我們根據6.2.1節中的最後兩個式子進行優化，其中一個是最小化的目標函數，一個是在優化過程中必須遵循的約束條件。不久之前，人們還在使用二次規劃求解工具（quadratic solver）來求解上述最優化問題，這種工具是一種用於在線性約束下優化具有多個變量的二次目標函數的軟件。而這些二次規劃求解工具則需要強大的計算能力支撐，另外在實現上也十分複雜。所有需要做的圍繞優化的事情就是訓練分類器，一旦得到alpha的最優值，我們就得到了分隔超平面（2維平面中就是直線）並能夠將之用於數據分類。
下面我們就開始討論SMO算法，然後給出一個簡化的版本，以便讀者能夠正確理解它的工作流程。後一節將會給出SMO算法的完整版，它比簡化版的運行速度要快很多。
6.3.1 Platt的SMO算法
1996年，John Platt發佈了一個稱為SMO1的強大算法，用於訓練SVM。SMO表示序列最小優化（Sequential Minimal Optimization）。Platt的SMO算法是將大優化問題分解為多個小優化問題來求解的。這些小優化問題往往很容易求解，並且對它們進行順序求解的結果與將它們作為整體來求解的結果是完全一致的。在結果完全相同的同時，SMO算法的求解時間短很多。
SMO算法的目標是求出一系列alpha和b，一旦求出了這些alpha，就很容易計算出權重向量w並得到分隔超平面。
SMO算法的工作原理是：每次循環中選擇兩個alpha進行優化處理。一旦找到一對合適的alpha，那麼就增大其中一個同時減小另一個。這裡所謂的「合適」就是指兩個alpha必須要符合一定的條件，條件之一就是這兩個alpha必須要在間隔邊界之外，而其第二個條件則是這兩個alpha還沒有進行過區間化處理或者不在邊界上。
1. John C. Platt, 「Using Analytic QP and Sparseness to Speed Training of Support Vector Machines」 in Advances in Neural Information Processing Systems 11, M. S. Kearns, S. A. Solla, D. A. Cohn, eds(MIT Press, 1999), 557–63.
6.3.2 應用簡化版SMO算法處理小規模數據集
Platt SMO算法的完整實現需要大量代碼。在接下來的第一個例子中，我們將會對算法進行簡化處理，以便瞭解算法的基本工作思路，之後再基於簡化版給出完整版。簡化版代碼雖然量少但執行速度慢。Platt SMO算法中的外循環確定要優化的最佳alpha對。而簡化版卻會跳過這一部分，首先在數據集上遍歷每一個alpha，然後在剩下的alpha集合中隨機選擇另一個alpha，從而構建alpha對。這裡有一點相當重要，就是我們要同時改變兩個alpha。之所以這樣做是因為我們有一個約束條件：
由於改變一個alpha可能會導致該約束條件失效，因此我們總是同時改變兩個alpha。
為此，我們將構建一個輔助函數，用於在某個區間範圍內隨機選擇一個整數。同時，我們也需要另一個輔助函數，用於在數值太大時對其進行調整。下面的程序清單給出了這兩個函數的實現。讀者可以打開一個文本編輯器將這些代碼加入到svmMLiA.py文件中。
程序清單6-1 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
def selectJrand(i,m):
j=i
while (j==i):
j = int(random.uniform(0,m))
return j
def clipAlpha(aj,H,L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
在testSet.txt文件中保存了圖6-3所給出的數據是。接下來，我們就將在這個文件上應用SMO算法。程序清單6-1中的第一個函數就是我們所熟知的loadDatSet函數，該函數打開文件並對其進行逐行解析，從而得到每行的類標籤和整個數據矩陣。
下一個函數selectJrand有兩個參數值，其中i是第一個alpha的下標，m是所有alpha的數目。只要函數值不等於輸入值i，函數就會進行隨機選擇。
最後一個輔助函數就是clipAlpha，它是用於調整大於H或小於L的alpha值。儘管上述3個輔助函數本身做的事情不多，但在分類器中卻很有用處。
在輸入並保存程序清單6-1中的代碼之後，運行如下命令：
>>> import svmMLiA
>>> dataArr,labelArr = svmMLiA.loadDataSet(\'testSet.txt\')
>>> labelArr
[-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0...
可以看得出來，這裡採用的類別標籤是-1和1，而不是0和1。
上述工作完成之後，就可以使用SMO算法的第一個版本了。
該SMO函數的偽代碼看上去像下面這個樣子：
創建一個alpha向量並將其初始化為0向量
當迭代次數小於最大迭代次數時（外循環）
對數據集中的每個數據向量（內循環）：
如果該數據向量可以被優化：
隨機選擇另外一個數據向量
同時優化這兩個向量
如果兩個向量都不能被優化，退出內循環
如果所有向量都沒被優化，增加迭代數目，繼續下一次循環
程序清單6-2中的代碼是SMO算法的一個有效版本。在Python中，如果某行以符號結束，那麼就意味著該行語句沒有結束並會在下一行延續。下面的代碼當中有很多很長的語句必須要分成多行來寫。因此，下面的程序中使用了多個符號。打開文件svmMLiA.py之後輸入如下程序清單中的代碼。
程序清單6-2 簡化版SMO算法
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose
b = 0; m,n = shape(dataMatrix)
alphas = mat(zeros((m,1)))
iter = 0
while (iter < maxIter):
alphaPairsChanged = 0
for i in range(m):
fXi = float(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
j = selectJrand(i,m)
fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
Ej = fXj - float(labelMat[j])
alphaIold = alphas[i].copy;
alphaJold = alphas[j].copy;
# ❸（以下六行）保證alpha在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
alphas[j] -= labelMat[j]*(Ei - Ej)/eta
alphas[j] = clipAlpha(alphas[j],H,L)
if (abs(alphas[j] - alphaJold) < 0.00001): print \"j not moving enough\"; continue
#❹ 對i進行修改，修改量與j相同，但方向相反
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
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)
if (alphaPairsChanged == 0): iter += 1
else: iter = 0
print \"iteration number: %d\" % iter
return b,alphas
這個函數比較大，或許是我所知道的本書中最大的一個函數。該函數有5個輸入參數，分別是：數據集、類別標籤、常數C、容錯率和退出前最大的循環次數。在本書，我們構建函數時採用了通用的接口，這樣就可以對算法和數據源進行組合或配對處理。上述函數將多個列表和輸入參數轉換成NumPy矩陣，這樣就可以簡化很多數學處理操作。由於轉置了類別標籤，因此我們得到的就是一個列向量而不是列表。於是類別標籤向量的每行元素都和數據矩陣中的行一一對應。我們也可以通過矩陣dataMatIn的shape屬性得到常數m和n。最後，我們就可以構建一個alpha列矩陣，矩陣中元素都初始化為0，並建立一個iter變量。該變量存儲的則是在沒有任何alpha改變的情況下遍歷數據集的次數。當該變量達到輸入值 maxIter時，函數結束運行並退出。
每次循環當中，將alphaPairsChanged先設為0，然後再對整個集合順序遍歷。變量alphaPairsChanged用於記錄alpha是否已經進行優化。當然，在循環結束時就會得知這一點。首先，fXi能夠計算出來，這就是我們預測的類別。然後，基於這個實例的預測結果和真實結果的比對，就可以計算誤差Ei。如果誤差很大，那麼可以對該數據實例所對應的alpha值進行優化。對該條件的測試處於上述程序清單的❶處。在if語句中，不管是正間隔還是負間隔都會被測試。並且在該if語句中，也要同時檢查alpha值，以保證其不能等於0或C。由於後面alpha小於0或大於C時將被調整為0或C，所以一旦在該if語句中它們等於這兩個值的話，那麼它們就已經在「邊界」上了，因而不再能夠減小或增大，因此也就不值得再對它們進行優化了。
接下來，可以利用程序清單6-1中的輔助函數來隨機選擇第二個alpha值，即alpha[j]❷。同樣，可以採用第一個alpha值(alpha[i])的誤差計算方法，來計算這個alpha值的誤差。這個過程可以通過copy的方法來實現，因此稍後可以將新的alpha值與老的alpha值進行比較。Python則會通過引用的方式傳遞所有列表，所以必須明確地告知Python要為alphaIold和alphaJold分配新的內存；否則的話，在對新值和舊值進行比較時，我們就看不到新舊值的變化。之後我們開始計算L和H❸，它們用於將alpha[j]調整到0到C之間。如果L和H相等，就不做任何改變，直接執行continue語句。這在Python中，則意味著本次循環結束直接運行下一次for的循環。
Eta是alpha[j]的最優修改量，在那個很長的計算代碼行中得到。如果eta為0，那就是說需要退出for循環的當前迭代過程。該過程對真實SMO算法進行了簡化處理。如果eta為0，那麼計算新的alpha[j]就比較麻煩了，這裡我們就不對此進行詳細的介紹了。有需要的讀者可以閱讀Platt的原文來瞭解更多的細節。現實中，這種情況並不常發生，因此忽略這一部分通常也無傷大雅。於是，可以計算出一個新的alpha[j]，然後利用程序清單6-1中的輔助函數以及L與H值對其進行調整。
然後，就是需要檢查alpha[j]是否有輕微改變。如果是的話，就退出for循環。然後，alpha[i]和alpha[j]同樣進行改變，雖然改變的大小一樣，但是改變的方向正好相反（即如果一個增加，那麼另外一個減少）❹。在對alpha[i]和alpha[j]進行優化之後，給這兩個alpha值設置一個常數項b❺。
最後，在優化過程結束的同時，必須確保在合適的時機結束循環。如果程序執行到for循環的最後一行都不執行continue語句，那麼就已經成功地改變了一對alpha，同時可以增加alphaPairsChanged的值。在for循環之外，需要檢查alpha值是否做了更新，如果有更新則將iter設為0後繼續運行程序。只有在所有數據集上遍歷maxIter次，且不再發生任何alpha修改之後，程序才會停止並退出while循環。
為瞭解實際效果，可以運行如下命令：
>>> b,alphas = svmMLiA.smoSimple(dataArr, labelArr, 0.6, 0.001, 40)
運行後輸出類似如下結果：
iteration number: 29
j not moving enough
iteration number: 30
iter: 30 i:17, pairs changed 1
j not moving enough
iteration number: 0
j not moving enough
iteration number: 1
上述運行過程需要幾分鐘才會收斂。一旦運行結束，我們可以對結果進行觀察：
>>> b
matrix([[-3.84064413]])
我們可以直接觀察alpha矩陣本身，但是其中的零元素太多。為了觀察大於0的元素的數量，可以輸入如下命令：
>>> alphas[alphas>0]
matrix([[ 0.12735413, 0.24154794, 0.36890208]])
由於SMO算法的隨機性，讀者運行後所得到的結果可能會與上述結果不同。alphas[alphas>0]命令是數組過濾（array filtering）的一個實例，而且它只對NumPy類型有用，卻並不適用於Python中的正則表（regular list）。如果輸入alpha>0，那麼就會得到一個布爾數組，並且在不等式成立的情況下，其對應值為正確的。於是，在將該布爾數組應用到原始的矩陣當中時，就會得到一個NumPy矩陣，並且其中矩陣僅僅包含大於0的值。
為了得到支持向量的個數，輸入：
>>> shape(alphas[alphas>0])
為瞭解哪些數據點是支持向量，輸入：
>>> for i in range(100):
... if alphas[i]>0.0: print dataArr[i],labelArr[i]
得到的結果類似如下：
...
[4.6581910000000004, 3.507396] -1.0
[3.4570959999999999, -0.082215999999999997] -1.0
[6.0805730000000002, 0.41888599999999998] 1.0
在原始數據集上對這些支持向量畫圈之後的結果如圖6-4所示。
圖6-4 在示例數據集上運行簡化版SMO算法後得到的結果，包括畫圈的支持向量與分隔超平面
利用前面的設置，我運行了10次程序並取其平均時間。結果是，這個過程在一台性能較差的筆記本上需要14.5秒。雖然結果看起來並不是太差，但是別忘了這只是一個僅有100個點的小規模數據集而已。在更大的數據集上，收斂時間會變得更長。在下一節中，我們將通過構建完整SMO算法來加快其運行速度。