机器学习课程笔记【五】- 支持向量机（2）

支持向量机算法是最有效的有监督学习算法之一，为全面掌握，首先学习间隔以及数据划分的思想，然后了解最优间隔分类器（涉及拉格朗日对偶，高维特征空间的核函数），最后介绍SVMs的一种高效实现——SMO算法。
本节为支持向量机部分笔记的第二小节，主要内容包括：结合拉格朗日对偶看理想间隔分类器，核函数推导，正则化与先线性不可分，SMO算法及其python实现代码，调用sklearn库中的svm函数分别完成线性可分和线性不可分的分类任务。

第一小节：机器学习课程笔记【五】- 支持向量机（1）

6. 再看理想间隔分类器

第四小节我们提出了一个找理想间隔分类的原始优化问题：

其中约束条件：

从KKT补充条件中我们知道当且仅当训练样本的函数间隔等于1则αi>0\alpha_i>0αi>0，看下图，最大化的间隔用实线表示，

有最小间隔的点刚好里决策边界最近，这里有三个虚线上的点代表的训练样本的αi\alpha_iαi再优化问题的解中非0。这三个点就是该问题中的支持向量，支持向量的数量会比训练集样本数量小很多。
继续，一个关键的点是我们需要把我们的算法视作内积的形式即

然后为我们的优化问题构建拉格朗日
乘子式：

这里只有αi\alpha_iαi没有βi\beta_iβi，因为这个问题只有一个不等式约束
之后找到这个问题的对偶形式，为此我们首先需最小化

令导数为0有：

也即：

把该上式代回拉格朗日乘子式中，有：

根据原始的拉格朗日乘子式求对b的偏导，令其为0，有：

所以将ω\omegaω代回之后的式子可以写为：

得到下述对偶优化问题：

这里对偶问题与原始问题优化的值相等，而且KKT条件也满足。于是我们可以通过解对偶问题来解原始问题。
对偶问题中，最大化问题需要拟合的参数是αi′s\alpha_i'sαi′s，这样我们就可以反向找到最优的ω\omegaω，(α\alphaα的函数)。找到ω∗\omega^*ω∗之后，可以直接根据原始优化问题得出b的值：

而且根据下式：

我们可以这样计算ωTx+b\omega^Tx+bωTx+b:

因此确定了αi\alpha_iαi之后，需要计算样本之间的内积，既可以得到预测结果。之前的分析中我们已经得出，只有维数不多的支持向量中αi\alpha_iαi非0，其他都为0，这样计算量也就非常小了。
通过以上对优化问题的对偶形式的讨论，我们对优化的问题更加理解，并可以以内积的形式写出全部的算法。下一节，就利用这些属性来应用核函数到分类问题中，对应的算法-支持向量机-在高维空间的性能表现，也会非常优秀。

7. 核函数

回到我们之前的线性回归问题，输入的xxx表示房子的居住面积，我们现在考虑使用x，x2，x3x，x^2，x^3x，x2，x3这三个特征来构造一个立方函数。称原始的输入xxx为输入参数，称经过映射后的量x，x2，x3x，x^2，x^3x，x2，x3为输入特征，定义φ\varphiφ为特征映射函数，这个例子中，我们有：

与使用原始输入参数xxx不同，SVMs使用输入特征φ(x)\varphi(x)φ(x)，这样我们之前的算法中的xxx都应替换为φ(x)\varphi(x)φ(x)。因为原来的算法可以完全写为向量内积的形式，则我们可以把内积都替换为<φ(x),φ(z)><\varphi(x),\varphi(z)><φ(x),φ(z)>。这样，给定特征映射函数φ\varphiφ，定义核函数(Kernel)为：

这样之前使用<x,z><x,z><x,z>的算法都可以替换为K(x,z)K(x,z)K(x,z)。
通常K(x,z)K(x,z)K(x,z)计算代价很低，即使φ(x)\varphi(x)φ(x)计算代价高。在算法中如果使用高效的方法计算K(x,z)K(x,z)K(x,z)，可以在给定φ\varphiφ的情况下，使用SVMs在高维特征空间进行学习，并且不用明确地计算向量φ(x)\varphi(x)φ(x)。
看这样一个例子，假设：

也可以展开计算：

与原始的写法对比：

这样映射函数可以写为

可以看到，直接计算映射函数需要O(n2)O(n^2)O(n2)的时间，而计算和函数只需要O(n)O(n)O(n)的时间复杂度。
继续考虑这样一个核函数

它对应的特征映射为：

参数ccc控制了xix_ixi和xixjx_ix_jxixj之间的权重关系。
更一般地，核函数：

对应于映射到(n+dd)\left(\begin{array}{c}n+d \\d\end{array}\right)(n+dd)维特征空间，但是其本身的计算复杂度，仍为线性。
从另外一个视角来看，如果φ(x)\varphi(x)φ(x)和φ(z)\varphi(z)φ(z)非常接近，那么K(x,z)=φ(x)Tφ(z)K(x,z)=\varphi(x)^T \varphi(z)K(x,z)=φ(x)Tφ(z)会非常大，若φ(x)\varphi(x)φ(x)和φ(z)\varphi(z)φ(z)接近正交，那么K(x,z)=φ(x)Tφ(z)K(x,z)=\varphi(x)^T \varphi(z)K(x,z)=φ(x)Tφ(z)会比较。所以这个K(x,z)K(x,z)K(x,z)也可以视作衡量φ(x)\varphi(x)φ(x)和φ(z)\varphi(z)φ(z)或者说xxx和zzz相似性的手段。
了解了上面这个结论，假设在我们自己的算法中选择了如下所示的核函数：

从相似性衡量的角度来看，当xxx和zzz接近，那么式子接近111，如果xxx和zzz差的多，那么式子接近000，这个核函数应用于SVM的话，我们称为高斯核，并且对应与一个无穷高维的特征映射。更一般地，给定某个函数KKK，我们应该如何断定它是否是一个有效的核函数，也就是说怎么样断定是否存在φ\varphiφ使得K(x,z)=φ(x)φ(z)K(x,z)=\varphi(x)\varphi(z)K(x,z)=φ(x)φ(z)对于所有xxx和zzz成立？
假定K是对应于映射φ\varphiφ的一个有效的核函数，考虑一个有mmm个点的有限即{x(1),...,x(m)}\{x^{(1)},...,x^{(m)}\}{x(1),...,x(m)}，并定义一个m×mm \times mm×m的方阵KKK，且有：

那么这个矩阵也被称为核矩阵，这样如果KKK是一个有效的核函数，那么有：

因此对应的和矩阵是对称矩阵，如果用φk(x)\varphi_k(x)φk(x)表示向量φ(x)\varphi(x)φ(x)的第kkk处坐标，那么对于任意向量zzz都有：

也就是说，对于任意的向量zzz，核矩阵都是正的半正定矩阵。这样，如果核函数KKK有效，那么对应的核矩阵必然为正对称半正定矩阵。这也得出一个核函数KKK有效的一个充要条件即Mercer定理，即核函数有效的充要条件是核矩阵为对称半正定矩阵，反过来，任何半正定函数都可以作为核函数。
这样，如果一个算法可以写为内积的形式那么我们就可以使用一个核函数来替代，这样我们的算法就可以高效地在高维特征空间进行计算。

8. 正则化和线性不可分的情况

上述关于SVM的推导是基于数据线性可分的情况，虽然通过φ\varphiφ将数据映射到高维特征空间会增加数据可分的可能性，但我们无法确保数据可分。而且在一些情况中，找到的分隔超平面因为离群点的存在并不是我们恰好需要的。如下图所示，增加离群点会使得分隔超平面的剧烈变化：

为了使算法对于线性不可分的数据和存在离群点的数据集同样有效，我们将做如下优化，称为L1正则化：

这样一来，就允许了函数间隔小于1。如果某个样本的函数间隔为1−ξi1-\xi_i1−ξi，可以在目标函数中应用一个代价因子使其增加到CξiC\xi_iCξi，参数CCC控制了∣∣ω∣∣2||\omega||^2∣∣ω∣∣2尽可能大和所有样本的函数间隔至少为1这两个优化目标的权重。
和之前一样，构造拉格朗日乘子式：

其中αi′s\alpha_i'sαi′s和ri′sr_i'sri′s是拉格朗日乘子并且严格不小于0，我们通过令拉格朗日函数关于ω\omegaω和bbb的导数分别为0，然后代回化简，得到问题的对偶形式如下：

同样地ω\omegaω也可以写作αi′s\alpha_i'sαi′s的函数，然后使用之前提到的式子进行预测：

加入L1正则化之后，对偶问题的变化也仅仅是0≤αi0 \leq \alpha_i0≤αi变为0≤αi≤C0 \leq \alpha_i \leq C0≤αi≤C，b∗b^*b∗的计算公式对应调整而已。
KKT对偶补充条件为：

这样一来，问题就转化为解对偶问题。

9. SMO算法

SMO(sequential minimal optimization)算法就给出了解对偶问题的一个方案，再学习SMO算法之前，先来看协同上升算法。

9.1 协同上升

考虑这样一个无约束的优化问题：

这里将WWW视为αi′s\alpha_i'sαi′s的函数并且暂且忽略这个问题和SVMs联系。之前我们也已经学习过梯度下降算法和牛顿方法，这里我们学习一种新的方法即协同上升算法：

在内循环中，我们保持除了αi\alpha_iαi之外的所有参数不变，然后根据这个αi\alpha_iαi重新优化WWW，这里给出的协同上升算法的版本根据这样的参数顺序不断重新优化WWW：

算法执行的过程如下：

图中椭圆是我们要优化的二次函数的等值线，每一次优化，协同上升算法迈出的步子总与坐标轴平行。

9.2 SMO

接下来讨论SMO，我们要解决的对偶问题是：

加入有一组αi′s\alpha_i'sαi′s满足约束条件，那么我们保持α2,...,αm\alpha_2,...,\alpha_mα2,...,αm不变，我们针对α1\alpha_1α1来应用协同上升算法，但这样实际上并不能得到一丁点的优化。因为：

或者写为：

也就是说α1\alpha_1α1由其他的α\alphaα确定，我们不能作像协同算法的那种迭代。
但是我们可以同时优化两个α\alphaα，这样就可以满足约束条件，这也就是SMO算法，过程如下：

SMO算法非常高效，因为αi\alpha_iαi和αj\alpha_jαj的更新十分高效。假定我们有一组αi′s\alpha_i'sαi′s满足约束条件，并且我们现在决定保持α3,...,αm\alpha_3,...,\alpha_mα3,...,αm不变，根据α1\alpha_1α1和α2\alpha_2α2来优化WWW，首先有：

因为等式右边固定所以可以写为：

这样可以画出α1\alpha_1α1和α2\alpha_2α2的约束即：

将α1\alpha_1α1写为α2\alpha_2α2的函数：

这样优化目标可以写为：

因为除了α2\alpha_2α2和α1\alpha_1α1之外的其他α\alphaα均为常数，那么优化问题就是一个简单的关于α2\alpha_2α2的二次函数，结合上面的图像有：

得到新的α2\alpha_2α2后代回可以得到α1\alpha_1α1。之后是bbb的更新，因为：

则：

损失可以写为：

因为新的b1b_1b1的前两项可以写为：

最终新的b1b_1b1为：

同理：

更新bbb的规则：

这样我们可以更新这个损失：

还有一个问题是SMO怎样选择要优化的一对参数。
第一个变量的选择称为外层循环，这个变量需要选择训练集中范围KKT条件最严重的样本点，KKT条件即：

一般来说，首先选择违反0<αi<C⇒yig(xi)=10<\alpha_i<C \Rightarrow y_ig(x_i)=10<αi<C⇒yig(xi)=1这个条件最严重的点，比如0<αi<C0<\alpha_i<C0<αi<C对应的点(xi,yi)(x_i,y_i)(xi,yi)，yi(xi)y_i(x_i)yi(xi)远远大于1或者小于1，则这个点就是违反该条件最严重的点。如果0<αi<C0<\alpha_i<C0<αi<C所有得到点都满足KKT条件，在选择违反alphai=0⇒yig(xi)≥1alpha_i=0 \Rightarrow y_ig(x_i) \geq1alphai=0⇒yig(xi)≥1和alphai=C⇒yig(xi)≤1alpha_i=C \Rightarrow y_ig(x_i) \leq 1alphai=C⇒yig(xi)≤1的点。
第二个变量的选择称为内层循环，假设我们在外层循环找到了α1\alpha_1α1，那么第二个变量的选择标准是让∣E1−E2∣|E_1-E_2|∣E1−E2∣有足够大的变化，由于α1\alpha_1α1确定之后，E1E_1E1也就确定了，所以要想∣E1−E2∣|E_1-E_2|∣E1−E2∣最大，只需要E1E_1E1为负时选择最小的EiE_iEi，E1E_1E1为负时，选择最大的E2E_2E2，我们可以将所有的EiE_iEi进行保存避免重复运算。
如果内层循环找到的点不能让目标函数有足够大的下降，可以采用遍历支持向量点来作为α2\alpha_2α2直到目标函数有足够的下降，若仍没有，则需要跳出内层循环重新选择α1\alpha_1α1。
简化版SMO算法的python实现，代码来自《机器学习实战》这本书。

'''
随机数，不等于J
'''
def selectJrand(i,m):j=i #we want to select any J not equal to iwhile (j==i):j = int(random.uniform(0,m))  # 一直在挑选随机数j，直到不等于i，随机数的范围在0~mreturn j  # 返回挑选好的随机数'''
门限函数
'''
def clipAlpha(aj,H,L):  # 最大不能超过H，最小不能低于Lif aj > H:aj = Hif L > aj:aj = Lreturn aj'''
简化版的SMO函数
'''
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):  # 输入数据，标记，常数C，容错率，最大迭代次数dataMatrix = mat(dataMatIn);   # 转换成矩阵labelMat = mat(classLabels).transpose()  # 转换成矩阵，并转置，标记成为一个列向量，每一行和数据矩阵对应m,n = shape(dataMatrix)  # 行，列   b = 0;  # 参数b的初始化alphas = mat(zeros((m,1)))  # 参数alphas是个list，初始化也是全0，大小等于样本数iter = 0  # 当前迭代次数，maxIter是最大迭代次数while (iter < maxIter):  # 当超过最大迭代次数，推出alphaPairsChanged = 0  # 标记位，记录alpha在该次循环中，有没有优化for i in range(m):  # 第i个样本fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b  # 第i样本的预测类别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)  # 随机选择第j个样本fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b  # 样本j的预测类别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"; continueeta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].Tif eta >= 0: print "eta>=0"; continuealphas[j] -= labelMat[j]*(Ei - Ej)/etaalphas[j] = clipAlpha(alphas[j],H,L)  # 门限函数阻止alpha_j的修改量过大#如果修改量很微小if (abs(alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; continue# alpha_i的修改方向相反alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#update i by the same amount as j#the update is in the oppostie direction# 为两个alpha设置常数项bb1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].Tb2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].Tif (0 < alphas[i]) and (C > alphas[i]): b = b1elif (0 < alphas[j]) and (C > alphas[j]): b = b2else: b = (b1 + b2)/2.0# 说明alpha已经发生改变alphaPairsChanged += 1print "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)#如果没有更新，那么继续迭代；如果有更新，那么迭代次数归0，继续优化if (alphaPairsChanged == 0): iter += 1else: iter = 0print "iteration number: %d" % iter# 只有当某次优化更新达到了最大迭代次数，这个时候才返回优化之后的alpha和breturn b,alphas

当然我们也可以用现成的，对于线性可分的情况

import numpy as np
import pylab as pl   #画图用
from sklearn import svm#随机生成两组二位数据
np.random.seed(0)#使每次产生随机数不变
X = np.r_[np.random.randn(20,2)-[2,2],np.random.randn(20,2)+[2,2]]#注意这里np.r_[],而不是np.r_（）我都打错了，会报错TypeError: 'RClass' object is not callable
#np.r_是按列连接两个矩阵，就是把两矩阵上下相加，要求列数相等，np.c_是按行连接两个矩阵，就是把两矩阵左右相加，要求行数相等Y = [0] * 20+[1] * 20#Python原来可以这么简单的创建重复的列表呀clf=svm.SVC(kernel='linear')
clf.fit(X,Y)w=clf.coef_[0]
a=-w[0]/w[1]
xx=np.linspace(-5,5)#产生-5到5的线性连续值，间隔为1
yy=a*xx-(clf.intercept_[0])/w[1]  #clf.intercept_[0]是w3.即为公式a1*x1+a2*x2+w3中的w3。(clf.intercept_[0])/w[1]即为直线的截距#得出支持向量的方程
b=clf.support_vectors_[0]
yy_down=a*xx+(b[1]-a*b[0])#(b[1]-a*b[0])就是简单的算截距
b=clf.support_vectors_[-1]
yy_up=a*xx+(b[1]-a*b[0])print("w:",w) #打印出权重系数
print("a:",a) #打印出斜率
print("suport_vectors_:",clf.support_vectors_)#打印出支持向量
print("clf.coef_:",clf.coef_)                  #打印出权重系数，还是w#这个就是画出来而已。很简单，也不太常用，都用matplotlib去了。不多说了
pl.plot(xx,yy,'k-')
pl.plot(xx,yy_down,'k--')
pl.plot(xx,yy_up,'k--')pl.scatter(clf.support_vectors_[:,0],clf.support_vectors_[:,0],s=80,facecolors='none')
pl.scatter(X[:,0],X[:,1],c=Y,cmap=pl.cm.Paired)pl.axis('tight')
pl.show()

线性不可分：

from time import time   #对程序运行时间计时用的
import logging           #打印程序进展日志用的
import matplotlib.pyplot as plt  #绘图用的from sklearn.model_selection import train_test_split
from sklearn.datasets import fetch_lfw_people
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix
from sklearn.decomposition import PCA
from sklearn.svm import SVClogging.basicConfig(level=logging.INFO,format='%(asctime)s %(message)s')
lfw_people=fetch_lfw_people(min_faces_per_person=70,resize=0.4)  #名人的人脸数据集、n_samples,h,w=lfw_people.images.shape  #多少个实例，h,w高度，宽度值X=lfw_people.data   #特征向量矩阵
n_feature=X.shape[1]#每个人有多少个特征值Y=lfw_people.target
target_names=lfw_people.target_names
n_classes=target_names.shape[0]     #多少类
print("Total dataset size")
print("n_samples:",n_samples)
print("n_feature:",n_feature)
print("n_classes:",n_classes)X_train,X_test,Y_train,Y_test=train_test_split(X,Y,test_size=0.25)  #选取0.25的测试集#降维
n_components=150  #PCA算法中所要保留的主成分个数n，也即保留下来的特征个数n
print("Extracting the top %d eigenfaces from %d faces" % (n_components,X_train.shape[0]))
t0=time()
pca=PCA(svd_solver='randomized',n_components=n_components,whiten=True).fit(X_train)#训练一个pca模型print("Train PCA in %0.3fs" % (time()-t0))eigenfaces  = pca.components_.reshape((n_components,h,w))  #提取出来特征值之后的矩阵print("Prijecting the input data on the eigenfaces orthonarmal basis")
t0=time()
X_train_pca = pca.transform(X_train)     #将训练集与测试集降维
X_test_pca = pca.transform(X_test)
print("Done PCA in %0.3fs" % (time()-t0))#终于到SVM训练了
print("Fiting the classifier to the training set")
t0=time()
param_grid ={'C':[1e3,5e3,1e4,5e4,1e5],#C是对错误的惩罚'gamma':[0.0001,0.0005,0.001,0.005,0.01,0.1],}#gamma核函数里多少个特征点会被使用}#对参数尝试不同的值
clf = GridSearchCV(SVC(kernel='rbf'),param_grid)
clf=clf.fit(X_train_pca,Y_train)
print("Done Fiting in %0.3fs" % (time()-t0))print("Best estimotor found by grid search:")
print(clf.best_estimator_)print("Predicting people's names on the test set")
t0=time()
Y_pred = clf.predict(X_test_pca)
print("done Predicting in %0.3fs" % (time()-t0))print(classification_report(Y_test,Y_pred,target_names=target_names)) #生成一个小报告呀
print(confusion_matrix(Y_test,Y_pred,labels=range(n_classes)))#这个也是，生成的矩阵的意思是有多少个被分为此类。#把分类完的图画出来12个。#这个函数就是画图的
def plot_gallery(images,titles,h,w,n_row=3,n_col=4):plt.figure(figsize=(1.8*n_col,2.4*n_row))plt.subplots_adjust(bottom=0,left=.01,right=.99,top=.90,hspace=.35)for i in range(n_row*n_col):plt.subplot(n_row,n_col,i+1)plt.imshow(images[i].reshape((h,w)),cmap=plt.cm.gray)plt.title(titles[i],size=12)plt.xticks(())plt.yticks(())#这个函数是生成一个固定格式的字符串的
def title(y_pred,y_test,target_names,i):pred_name=target_names[y_pred[i]].rsplit(' ',1)[-1]true_name = target_names[y_test[i]].rsplit(' ', 1)[-1]return "predicted: %s\n true: %s" %(pred_name,true_name)predicted_titles=[title(Y_pred,Y_test,target_names,i) for i in range(Y_pred.shape[0])] #这个for循环的用法很简介plot_gallery(X_test,predicted_titles,h,w)eigenfaces_titles=["eigenface %d " % i for i in range(eigenfaces.shape[0])]
plot_gallery(eigenfaces,eigenfaces_titles,h,w)plt.show()

欢迎扫描二维码关注微信公众号深度学习与数学 [每天获取免费的大数据、AI等相关的学习资源、经典和最新的深度学习相关的论文研读，算法和其他互联网技能的学习，概率论、线性代数等高等数学知识的回顾]