本文深入探讨了支持向量机(SVM)的基本原理，包括最大间隔分隔数据、寻找最大间隔、SMO优化算法以及核函数的应用。通过实例展示了如何使用SVM进行手写数字识别，强调了核函数参数调整对手写数字识别准确率的影响，并探讨了不同核函数在不同类型数据上的表现。此外，还讨论了SVM在复杂数据上的优势，如非线性可分数据的处理。
目录
6.1 基于最大间隔分隔数据
支持向量机
优点：泛化错误率低，计算开销不大，结果易解释。
缺点：对参数调节和核函数的选择敏感，原始分类器不加修改仅适用于处理二类问题。
适用数据类型：数值型和标称型数据。
线性可分：数据点分布可被一条直线分开。
分隔超平面：将数据集分割开的直线。
超平面：具有高维的数据集，分隔数据的也是高维的对象。分布在超平面两侧的所有数据是不同类别的。
我们希望找到离分隔超平面最近的点，确保它们离分割面的距离尽可能远。这里点到分割面的距离被称为间隔。我们希望间隔尽可能地大。
支持向量就是离分隔超平面最近的那些点。接下来需要最大化支持向量到分割面的距离，需要找到此问题的优化求解方法。
6.2 寻找最大间隔
分隔超平面的形式可以写成，计算点A到分隔超平面的距离，就要计算点到分割面的法线或垂线长度
。
6.2.1 分类器求解的优化问题
输入数据给分类器会输出一个类别标签,这相当于一个类似于Sigmoid的函数在作用。下面将使用类似海维赛德阶跃函数（即单位阶跃函数）的函数对作用得到
，其中当u<o时f(u)输出-1，反之则输出+1。
当计算数据点到分隔面的距离并确定分隔面的放置位置时，间隔通过来计算，这时就能体现出-1和+1类的好处了。如果数据点处于正方向（即+1类）并且离分隔超平面很远的位置时，
会是一个很大的正数，同时
也会是一个很大的正数。而如果数据点处于负方向(-1类）并且离分隔超平面很远的位置时，此时由于类别标签为-1，则
仍然是一个很大的正数。
目前的目标就是找出分类器定义中的w和b。为此，我们必须找到具有最小间隔的数据点，而这些数据点也就是前面提到的支持向量。一旦找到具有最小间隔的数据点，我们就需要对该间隔最大化。这就可以写作：
直接求解十分困难，因此对问题进行优化，给定一些约束条件然后求其最优值。约束条件为。对于这类问题，可以引用拉格朗日乘子进行求解。通过引入拉格朗日乘子，我们可以基于约束条件来表述原来的问题。由于这里的约束条件都是基于数据点的，因此我们就可以将超平面写成数据点的形式。于是，优化目标函数最后可以写成：
其约束条件为：
我们默认了数据必须100%线性可分，而我们知道几乎所有数据都并非如此。这时我们就可以通过引入所谓松弛变量，来允许有些数据点可以处于分割面的错误一侧。这样我们的优化目标就能保持仍然不变，但是此时新的约束条件则变为：
常数c用于控制“最大化间隔”和“保证大部分点的函数间隔小于1.0”这两个目标的权重。在优化算法的实现代码中，常数c是一个参数，因此我们就可以通过调节该参数得到不同的结果。一旦求出了所有的alpha，那么分隔超平面就可以通过这些alpha来表达。
6.2.2 SVM应用的一般框架
6.3 SMO高效优化算法
6.3.1 Platt的SMO算法
SMO表示序列最小优化。Platt的SMO算法是将最大优化问题分解为多个小优化问题来求解的。这些小优化问题往往很容易求解，并且对它们进行顺序求解的结果与将它们作为整体来求解的结果是完全一致的。在结果完全相同的同时，SMO算法的求解时间短很多。
SMO算法的目标是求出一系列alpha和b，一旦求出了这些alpha，就很容易计算出权重向量w并得到分隔超平面。
SMO算法的工作原理是：每次循环中选择两个alpha进行优化处理。一旦找到一对合适的alpha，那么就增大其中一个同时减小另一个。这里所谓的“合适”就是指两个alpha必须要符合一定的条件，条件之一就是这两个alpha必须要在间隔边界之外，而其第二个条件则是这两个alpha还没有进行过区间化处理或者不在边界上。
6.3.2 应用简化版SMO算法处理小规模数据集
简化版算法首先在数据集上遍历每一个alpha，然后在剩下的alpha集合中随机选择另一个alpha，从而构建alpha对。我们要同时改变两个alpha，原因是有一个约束条件：
由于改变一个alpha可能会导致该约束条件失效，因此我们总是同时改变两个alpha。
我们将构建一个辅助函数，用于在某个区间范围内随机选择一个整数。同时，我们也需要另一个辅助函数，用于在数值太大时对其进行调整。下面给出了这两个函数的实现，将其添加到svmMLiA.py 文件中。
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 selectJrand(i,m):
j=i #we want to select any J not equal to 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第一个函数就是我们所熟知的 loadDatSet() 函数，该函数打开文件并对其进行逐行解析，从而得到每行的类标签和整个数据矩阵。
下一个函数 selectJrand() 有两个参数值，其中i是第一个alpha的下标，m是所有alpha的数目。只要函数值不等于输入值i，函数就会进行随机选择。
最后一个辅助函数就是 clipAlpha( )，它是用于调整大于H或小于r的alpha值。
import svmMLiA
dataArr,labelArr = svmMLiA.loadDataSet('testSet.txt')
print(labelArr)
可以看出，采用的类别标签是-1和1，而不是0和1。
SMO函数的伪代码大致如下；
下面的代码是SMO算法的一个有效版本。在Python中，如果某行以\符号结束，那么就意味着该行语句没有结束并会在下一行延续。下面的代码当中有很多很长的语句必须要分成多行来写。因此，下面的程序中使用了多个\符号。打开文件 svmMLiA.py 输入如下代码。
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])#if checks if an example violates KKT conditions
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
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();
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
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#update i by the same amount as j
#the update is in the oppostie direction
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变量。该变量存储的则是在没有任何alphai改变的情况下遍历数据集的次数。当该变量达到输入值maxIter时，函数结束运行并退出。
每次循环当中，将alphaPairsChanged先设为0，然后再对整个集合顺序遍历。变量alphaPairsChanged用于记录alpha是否已经进行优化。当然，在循环结束时就会得知这一点。首先，fXi能够计算出来，这就是我们预测的类别。然后，基于这个实例的预测结果和真实结果的比对，就可以计算误差Ei。如果误差很大，那么可以对该数据实例所对应的alpha值进行优化。在if语句中，不管是正间隔还是负间隔都会被测试。并且在该if语句中，也要同时检查alpha值，以保证其不能等于0或C。由于后面alpha小于0或大于C时将被调整为0或C，所以一旦在该if语句中它们等于这两个值的话，那么它们就已经在“边界”上了，因而不再能够减小或增大，因此也就不值得再对它们进行优化了。`
接下来，可以利用辅助函数来随机选择第二个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为o，那就是说需要退出for循环的当前迭代过程。该过程对真实SMO算法进行了简化处理。如果eta为0，那么计算新的alpha[j]就比较麻烦了，这里我们就不对此进行详细的介绍了。现实中，这种情况并不常发生，因此忽略这一部分通常也无伤大雅。于是，可以计算出一个新的alpha[j]，然后利用程序中的辅助函数以及L与H值对其进行调整。
然后，就是需要检查alpha[j]是否有轻微改变。如果是的话，就退出for循环。然后，alpha[i]和alpha[j]同样进行改变，虽然改变的大小一样，但是改变的方向正好相反（即如果一个增加，那么另外一个减少)。在对alpha[i]和alpha[j]进行优化之后，给这两个a1pha值设置一个常数项b。
最后，在优化过程结束的同时，必须确保在合适的时机结束循环。如果程序执行到for循环的最后一行都不执行continue语句，那么就已经成功地改变了一堆alpha，同时可以增加alphaPairsChanged的值。在for循环之外，需要检查alpha值是否做了更新，如果有更新则将iter设为0后继续运行程序。只有在所有数据集上遍历maxIter次，且不再发生任何alpha修改之后，程序才会停止并退出while循环。
测试函数效果：
import svmMLiA
dataArr,labelArr = svmMLiA.loadDataSet('testSet.txt')
b,alphas = svmMLiA.smoSimple(dataArr,labelArr,0.6,0.001,40)
print(b)
print(alphas[alphas>0])
alpha矩阵本身零元素太多，使用代码对元素进行过滤。由于SMO算法的随机性，运行后得到实际结果可能会不相同。alphas[alphas>0]命令是数组过滤（array filtering) 的一个实例，而且它只对NumPy类型有用,却并不适用于Python中的正则表 (regular list)。如果输入alpha>0，那么就会得到一个布尔数组，并且在不等式成立的情况下，其对应值为正确的。于是，在将该布尔数组应用到原始的矩阵当中时，就会得到一个NumPy矩阵，并且其中矩阵仅仅包含大于0的值。
得到支持向量的个数，以及了解哪些数据点是支持向量，输入：
import svmMLiA
from numpy import *
dataArr,labelArr = svmMLiA.loadDataSet('testSet.txt')
b,alphas = svmMLiA.smoSimple(dataArr,labelArr,0.6,0.001,40)
shape(alphas[alphas>0])
for i in range(100):
if alphas[i]>0.0:
print(dataArr[i],labelArr[i]) 
在原始数据集上对这些支持向量画圈之后的结果如下所示：
6.4 利用完整Platt SMO算法加速优化
Platt SMO算法是通过一个外循环来选择第一个alpha值的，并且其选择过程会在两种方式之间进行交替：一种方式是在所有数据集上进行单遍扫描，另一种则是在非边界alpha中实现单遍扫描。而所谓非边界alpha指的就是那些不等于边界0或C的alpha值。对整个数据集的扫描相当容易，而实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后再对这个表进行遍历。同时，该步骤会跳过那些已知的不会改变的alpha值。
在选择第一个alpha值后，算法会通过一个内循环来确定第二个alpha值。在优化过程中，会通过最大化步长的方式来获得第二个alpha值。在简化版SMO算法中，我们会在选择j之后计算错误率Ej。但在此，我们会建立一个全局的缓存用于保存误差值，并从中选择使得步长或者说Ei-Ej最大的alpha值。
下面的程序包含一个用于清理代码的数据结构和3个用于对E进行缓存的辅助函数。
class optStruct:
def __init__(self, dataMatIn, classLabels, C, toler, kTup): # Initialize the structure with the parameters
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))) # first column is valid flag
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
def selectJ(i, oS, Ei): # this is the second choice -heurstic, and calcs Ej
maxK = -1;
maxDeltaE = 0;
Ej = 0
oS.eCache[i] = [1, Ei] # set valid #choose the alpha that gives the maximum delta E
validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
if (len(validEcacheList)) > 1:
for k in validEcacheList: # loop through valid Ecache values and find the one that maximizes delta E
if k == i: continue # don't calc for i, waste of time
Ek = calcEk(oS, k)
deltaE = abs(Ei - Ek)
if (deltaE > maxDeltaE):
maxK = k;
maxDeltaE = deltaE;
Ej = Ek
return maxK, Ej
else: # in this case (first time around) we don't have any valid eCache values
j = selectJrand(i, oS.m)
Ej = calcEk(oS, j)
return j, Ej
def updateEk(oS, k): # after any alpha has changed update the new value in the cache
Ek = calcEk(oS, k)
oS.eCache[k] = [1, Ek]首要的事情就是建立一个数据结构来保存所有的重要值,而这个过程可以通过一个对象来完成。这里使用对象的目的并不是为了面向对象的编程，而只是作为一个数据结构来使用对象。在将值传给函数时，我们可以通过将所有数据移到一个结构中来实现，这样就可以省掉手工输入的麻烦了。而此时，数据就可以通过一个对象来进行传递。实际上，当完成其实现时，可以很容易通过Python的字典来完成。但是在访问对象成员变量时，这样做会有更多的手工输入操作，对比一下myobject.X和myobject['X']就可以知道这一点。为达到这个目的,需要构建一个仅包含init方法的optStruct类。该方法可以实现其成员变量的填充。除了增加了一个m×2的矩阵成员变量eCache之外，这些做法和简化版SMO一模一样。eCache的第一列给出的是eCache是否有效的标志位，而第二列给出的是实际的E值。
对于给定的alpha值，第一个辅助函数 calcEk() 能够计算E值并返回。以前，该过程是采用内嵌的方式来完成的,但是由于该过程在这个版本的SMO算法中出现频繁,这里必须要将其单独拎出来。
下一个函数 selectJ() 用于选择第二个alpha或者说内循环的alpha值。回想一下，这里的目标是选择合适的第二个alpha值以保证在每次优化中采用最大步长。该函数的误差值与第一个alpha值Ei和下标i有关。首先将输入值Ei在缓存中设置成为有效的。这里的有效(valid)意味着它已经计算好了。在eCaahe中，代码nonzero(oS.eCache[ : ,0].A) [0]构建出了一个非零表。NumPy函数nonzero()返回了一个列表，而这个列表中包含以输入列表为目录的列表值，当然读者可以猜得到，这里的值并非零。nonzera ( )语句返回的是非零z值所对应的alpha值，而不是E值本身。程序会在所有的值上进行循环并选择其中使得改变最大的那个值。如果这是第一次循环的话，那么就随机选择一个alpha值。当然，也存在有许多更复杂的方式来处理第一次循环的情况，而上述做法就能够满足我们的目的。
程序最后一个辅助函数是 updateEk( )，它会计算误差值并存人缓存当中。在对alpha值进行优化之后会用到这个值。
程序代码本身的作用并不大，但是当和优化过程及外循环组合在一起时，就能组成强大的SMO算法。
下面的代码是用于寻找决策边界的优化例程：
def innerL(i, oS):
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)):
j,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand
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] #changed for kernel
if eta >= 0: print("eta>=0"); return 0
oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
updateEk(oS, j) #added this for the Ecache
if (abs(oS.alphas[j] - alphaJold) < 0.00001): print ("j not moving enough"); return 0
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
updateEk(oS, i) #added this for the Ecache #the update is in the oppostie direction
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上面的代码几乎和 smoSimple() 函数一模一样，但是这里的代码已经使用了自己的数据结构。该结构在参数oS中传递。第二个重要的修改就是使用 SelectJ() 而不是 selectJrand() 来选择第二个alpha的值。最后，在alpha值改变时更新Ecache。下面将给出把上述过程打包在一起的代码片段。这就是选择第一个alpha值的外循环。
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): #full Platt SMO
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: #go over all
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:#go over non-bound (railed) alphas
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 #toggle entire set loop
elif (alphaPairsChanged == 0): entireSet = True
print "iteration number: %d" % iter
return oS.b,oS.alphas程序给出的是完整版的Platt SMO算法，其输入和函数 smoSimple() 完全一样。函数一开始构建一个数据结构来容纳所有的数据，然后需要对控制函数退出的一些变量进行初始化。整个代码的主体是while循环,这与 smoSimple() 有些类似,但是这里的循环退出条件更多一些。当迭代次数超过指定的最大值，或者遍历整个集合都未对任意alpha对进行修改时，就退出循环。这里的maxIter变量和函数 smoSimple() 中的作用有一点不同，后者当没有任何alpha发生改变时会将整个集合的一次遍历过程计成一次迭代，而这里的一次迭代定义为一次循环过程，而不管该循环具体做了什么事。此时，如果在优化过程中存在波动就会停止，因此这里的做法优于 smoSimple() 函数中的计数方法。
while循环的内部与 smoSimple() 中有所不同，一开始的for循环在数据集上遍历任意可能的alpha。我们通过调用 innerL() 来选择第二个alpha，并在可能时对其进行优化处理。如果有任意一对alpha值发生改变，那么会返回1。第二个for循环遍历所有的非边界alpha值，也就是不在边界0或c上的值。
接下来，我们对for循环在非边界循环和完整遍历之间进行切换，并打印出迭代次数。最后程序将会返回常数b和alpha值。
import svmMLiA
from numpy import *
dataArr,labelArr = svmMLiA.loadDataSet('testSet.txt')
b,alphas = svmMLiA.smoP(dataArr,labelArr,0.6,0.001,40)
与 smoSimple() 函数相比，同样的数据集上，运行速度明显要更快。
常数c给出的是不同优化问题的权重。常数c一方面要保障所有样例的间隔不小于1.0，另一方面又要使得分类间隔要尽可能大，并且要在这两方面之间平衡。如果c很大，那么分类器将力图通过分隔超平面对所有的样例都正确分类。这种优化的运行结果如下图所示。与之前的图相比，会发现下图的支持向量更多。如果回想一下，就会记得前图实际来自于简化版算法，该算法是通过随机的方式选择alpha对的。这种简单的方式也可以工作，但是效果却不如完整版本好，后者覆盖了整个数据集。读者可能还认为选出的支持向量应该始终最接近分隔超平面。给定c的设置，图中画圈的支持向量就给出了满足算法的一种解。如果数据集非线性可分，就会发现支持向量会在超平面附近聚集成团。
计算alpha值之后，可以基于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上述代码中最重要的部分是for循环，虽然在循环中实现的仅仅是多个数的乘积。已知大部分alpha值为0，而非零alpha所对应的也就是支持向量。虽然上述for循环遍历了数据集中的所有数据，但是最终起作用只有支持向量。由于对w计算毫无作用，所以数据集的其他数据点也就会很容易地被舍弃。
import svmMLiA
from numpy import *
dataArr,labelArr = svmMLiA.loadDataSet('testSet.txt')
b,alphas = svmMLiA.smoP(dataArr,labelArr,0.6,0.001,40)
ws = svmMLiA.calcWs(alphas,dataArr,labelArr)
datMat = mat(dataArr)
print(ws)
print(datMat[0] * mat(ws) + b)
print(labelArr[0])
print(datMat[2] * mat(ws) + b)
print(labelArr[2])
print(datMat[1] * mat(ws) + b)
print(labelArr[1])
6.5 在复杂数据上应用核函数
考虑下图给出的数据，显而易见，在该数据中存在某种可以识别的模式。而我们能否像线性情况一样，利用强大的工具来捕捉数据中的这种模式？接下来，我们就要使用一种称为核函数的工具将数据转换成易于分类器理解的形式。
6.5.1 利用核函数将数据映射到高维空间
上图数据点处于一个圆中，人类的大脑能够意识到这一点。然而，对于分类器而言，它只能识别分类器的结果是大于0还是小于0。如果只在x和y轴构成的坐标系中插入直线进行分类的话，我们并不会得到理想的结果。我们或许可以对圆中的数据进行某种形式的转换，从而得到某些新的变量来表示数据。在这种表示情况下，我们就更容易得到大于0或者小于0的测试结果。在这个例子中，我们将数据从一个特征空间转换到另一个特征空间。在新空间下，我们可以很容易利用已有的工具对数据进行处理。数学家们喜欢将这个过程称之为从一个特征空间到另一个特征空间的映射。在通常情况下，这种映射会将低维特征空间映射到高维空间。
这种从某个特征空间到另一个特征空间的映射是通过核函数来实现的。读者可以把核函数想象成一个包装器（wrapper）或者是接口 (interface)，它能把数据从某个很难处理的形式转换成为另一个较容易处理的形式。如果上述特征空间映射的说法听起来很让人迷糊的话，那么可以将它想象成为另外一种距离计算的方法。前面我们提到过距离计算的方法。距离计算的方法有很多种，不久我们也将看到，核函数一样具有多种类型。经过空间转换之后，我们可以在高维空间中解决线性问题，这也就等价于在低维空间中解决非线性问题。
SVM优化中一个特别好的地方就是，所有的运算都可以写成内积 (inner product，也称点积) 的形式。向量的内积指的是两个向量相乘，之后得到单个标量或者数值。我们可以把内积运算替换成核函数，而不必做简化处理。将内积替换成核函数的方式被称为核技巧 (kernel trick) 或者核“变电”(kernel substation)。
6.5.2 径向基核函数
径向基函数是SVM中常用的一个核函数。径向基函数是一个采用向量作为自变量的函数，能够基于向零距离运算输出一个标量。这个距离可以是从<0,0>向量或者其他向量开始计算的距离。下面是径向基函数的高斯版本：
其中σ是用户定的用于确定到达率 (reach) 或者说函数值跌落到0的速度参数。
该高斯核函数将数据从其特征空间映射到更高维的空间，具体来说这里是映射到一个无穷维的空间。
在svmMLiA.py文件中添加一个函数并稍作修改：
class optStruct:
def __init__(self, dataMatIn, classLabels, C, toler, kTup): # Initialize the structure with the parameters
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))) # first column is valid flag
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 kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space
m,n = shape(X)
K = mat(zeros((m,1)))
if kTup[0]=='lin': K = X * A.T #linear kernel
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)) #divide in NumPy is element-wise not matrix like Matlab
else: raise NameError('Houston We Have a Problem -- \
That Kernel is not recognized')
return KoptStruct引入了一个新变量kTup，kTup是一个包含核函数信息的元组。在初始化方法结束时，矩阵K先被构建，然后再通过调用函数 kernelTrans() 进行填充。全局的k值只需计算一次。然后需要核函数时就可以对它进行调用，省去了很多冗余的计算开销。
kernelTrans函数由3个输入参数：两个数值型变量和一个元组。元组kTup给出的是核函数的信息。元组的第一个参数是描述所用核函数类型的一个字符串，其他2个参数都是核函数可能需要的可选参数。该函数首先构建出一个列向量，然后检查元组以确定核函数的类型。
在线性核函数的情况下，内积计算在“所有数据集”和“数据集中的一行”这两个输入之间展开。在径向基核函数的情况下，在for循环中对于矩阵的每个元素计算高斯函数的值。而在for循环结束之后，我们将计算过程应用到整个向量上去。值得一提的是，在NumPy矩阵中，除法符号意味着对矩阵元素展开计算而不像在MATLAB中一样计算矩阵的逆。
如果遇到一个无法识别的元组，程序就会抛出异常，因为在这种情况下不希望程序再继续运行，这一点相当重要。
为了使用核函数，先期的两个函数 innerL() 和 calcEK() 的代码需要做些修改。修改的结果如下：
def innerL(i, oS):
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)):
j,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand
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] #changed for kernel
if eta >= 0: print "eta>=0"; return 0
oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
updateEk(oS, j) #added this for the Ecache
if (abs(oS.alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; return 0
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
updateEk(oS, i) #added this for the Ecache #the update is in the oppostie direction
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
def calcEkK(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 Ek6.5.3 在测试中使用核函数
接下来将构建一个能进行有效分类的分类器，该分类器使用了径向基核函数。前面提到的径向基函数有一个用户定义的输入σ。首先，我们需要确定它的大小，然后利用该核函数构建出一个分类器。
def testRbf(k1=1.3):
dataArr,labelArr = loadDataSet('testSetRBF.txt')
b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important
datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
svInd=nonzero(alphas.A>0)[0]
sVs=datMat[svInd] #get matrix of only support vectors
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,:],('rbf', k1))
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 = loadDataSet('testSetRBF2.txt')
errorCount = 0
datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
m,n = shape(datMat)
for i in range(m):
kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
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)) 上述代码只有一个可选的输入参数，该输入参数是高斯径向基函数中的一个用户定义变量。整个代码主要是由以前定义的函数集合构成的。首先，程序从文件中读人数据集，然后在该数据集上运行Platt SMO算法，其中核函数的类型为'rbf ' 。
优化过程结束后，在后面的矩阵数学运算中建立了数据的矩阵副本，并且找出那些非零的alpha值，从而得到所需要的支持向量；同时，也就得到了这些支持向量和alpha的类别标签值。这些值仅仅是需要分类的值。
整个代码中最重要的是for循环开始的那两行，它们给出了如何利用核函数进行分类。首先利用结构初始化方法中使用过的kernelTrans( )函数，得到转换后的数据。然后，再用其与前面的alpha及类别标签值求积。其中需要特别注意的另一件事是，在这几行代码中,是如何做到只需要支持向量数据就可以进行分类的。除此之外，其他数据都可以直接舍弃。
与第一个for循环相比，第二个for循环仅仅只有数据集不同，后者采用的是测试数据集。
import svmMLiA
svmMLiA.testRbf()
三次使用了不同的k1参数，测试错误率、训练错误率、支持向量个数均不同。下图为k1非常小（=0.1）时的结果。
上图共有100个数据点，其中的85个为支持向量。优化算法发现，必须使用这些支持向量才能对数据进行正确分类。这可能给人径向基函数到达率太小的直觉。我们可以通过增加σ来观察错误率的变化情况，如下图：
两张图，后者只有27个支持向量，少了许多。同时观察函数testRbf输出值，可以发现测试错误率也在下降。该数据集在这个设置的某处存在着最优值。如果降低σ，那么训练错误率就会降低，但是测试错误率却会上升。
支持向量的数目存在一个最优值。SVM的优点在于它能对数据进行高效分类。如果支持向量太少，就可能会得到一个很差的决策边界；如果支持向量太多，也就相当于每次都利用整个数据集进行分类，这种分类方法称为k近邻。
我们可以对SMO算法中的其他设置进行随意地修改或者建立新的核函数。
6.6 示例：手写识别问题回顾
使用第2章中的一些代码和SMO算法，可以构建一个系统去测试手写数字上的分类器。将knn.py中的 img2vector() 函数复制到svmMLiA.py，再加入如下代码：
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) #load the training set
m = len(trainingFileList)
trainingMat = zeros((m,1024))
for i in range(m):
fileNameStr = trainingFileList[i]
fileStr = fileNameStr.split('.')[0] #take off .txt
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)) 函数 loadImages() 是作为前面kNN.py中的handwritingclassTest()的一部分出现的。它已经被重构为自身的一个函数。其中仅有的一个大区别在于，在kNN.py中代码直接应用类别标签，而同支持向量机一起使用时，类别标签为-1或者+1。因此，一旦碰到数字9，则输出类别标签-1，否则输出+1。本质上，支持向量机是一个二类分类器，其分类结果不是+1就是-1。基于SVM构建多类分类器已有很多研究和对比了。由于这里我们只做二类分类，因此除了1和9之外的数字都被去掉了。
下一个函数 testDigits() 并不是全新的函数，它和 testRbf( ) 的代码几乎一样，唯一的大区别就是它调用了 loadImages() 函数来获得类别标签和数据。另一个细小的不同是现在这里的函数元组kTup是输人参数，而在 testRbf() 中默认的就是使用rbf核函数。如果对于函数 testDigits() 不增加任何输入参数的话，那么kTup的默认值就是(‘rbf’ ,10)。
import svmMLiA
svmMLiA.testDigits(('rbf',20))
取不同σ值，并尝试线性核函数，得到结果如下：
上表给出的结果表明，当径向基核函数中的参数a取10左右时，就可以得到最小的测试错误率。该参数值比前面例子中的取值大得多，而前面的测试错误率在1.3左右。为什么差距如此之大?原因就在于数据的不同。在手写识别的数据中，有1024个特征，而这些特征的值有可能高达1.0。
你可能注意到了一个有趣的现象，即最小的训练错误率并不对应于最小的支持向量数目。另一个值得注意的就是，线性核函数的效果并不是特别的糟糕。可以以牺牲线性核函数的错误率来换取分类速度的提高。
扫一扫
1279
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
