接下来，我们根据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,labelMatdef selectJrand(i,m): j=i while (j==i): j = int(random.uniform(0,m)) return jdef clipAlpha(aj,H,L): if aj > H: aj = H if L > aj: aj = Lreturn 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: 29j not moving enoughiteration number: 30iter: 30 i:17, pairs changed 1j not moving enoughiteration number: 0j not moving enoughiteration number: 1 上述运行过程需要几分钟才会收敛。一旦运行结束，我们可以对结果进行观察：
>>> bmatrix([[-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算法来加快其运行速度。