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

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

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

6. 再看理想间隔分类器

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

  其中约束条件：

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

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

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

  这里只有${\alpha }_i$没有${\beta }_i$，因为这个问题只有一个不等式约束
  之后找到这个问题的对偶形式，为此我们首先需最小化

  令导数为0有：

  也即：

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

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

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

  得到下述对偶优化问题：

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

  而且根据下式：

  我们可以这样计算${\omega }^{T}x+b$:

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

7. 核函数

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

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

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

  也可以展开计算：

  与原始的写法对比：

  这样映射函数可以写为

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

  它对应的特征映射为：

  参数$c$控制了$x_i$和$x_i{x}_j$之间的权重关系。
  更一般地，核函数：

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

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

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

  因此对应的和矩阵是对称矩阵，如果用${\phi }_k(x)$表示向量$\phi (x)$的第$k$处坐标，那么对于任意向量$z$都有：

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

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

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

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

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

  其中${\alpha }_i^{\prime }s$和$r_i^{\prime }s$是拉格朗日乘子并且严格不小于0，我们通过令拉格朗日函数关于$\omega$和$b$的导数分别为0，然后代回化简，得到问题的对偶形式如下：

  同样地$\omega$也可以写作${\alpha }_i^{\prime }s$的函数，然后使用之前提到的式子进行预测：

  加入L1正则化之后，对偶问题的变化也仅仅是$0\le {\alpha }_i$变为$0\le {\alpha }_i\le C$，$b^{\ast }$的计算公式对应调整而已。
  KKT对偶补充条件为：

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

9. SMO算法

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

9.1 协同上升

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

  这里将$W$视为${\alpha }_i^{\prime }s$的函数并且暂且忽略这个问题和SVMs联系。之前我们也已经学习过梯度下降算法和牛顿方法，这里我们学习一种新的方法即协同上升算法：

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

  算法执行的过程如下：

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

9.2 SMO

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

  加入有一组${\alpha }_i^{\prime }s$满足约束条件，那么我们保持${\alpha }_2,...,{\alpha }_m$不变，我们针对${\alpha }_1$来应用协同上升算法，但这样实际上并不能得到一丁点的优化。因为：

  或者写为：

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

  SMO算法非常高效，因为${\alpha }_i$和${\alpha }_j$的更新十分高效。假定我们有一组${\alpha }_i^{\prime }s$满足约束条件，并且我们现在决定保持${\alpha }_3,...,{\alpha }_m$不变，根据${\alpha }_1$和${\alpha }_2$来优化$W$，首先有：

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

  这样可以画出${\alpha }_1$和${\alpha }_2$的约束即：

  将${\alpha }_1$写为${\alpha }_2$的函数：

  这样优化目标可以写为：

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

  得到新的${\alpha }_2$后代回可以得到${\alpha }_1$。之后是$b$的更新，因为：

  则：

  损失可以写为：

  因为新的$b_1$的前两项可以写为：

  最终新的$b_1$为：

  同理：

  更新$b$的规则：

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

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

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

'''
随机数，不等于J
'''
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))  # 一直在挑选随机数j，直到不等于i，随机数的范围在0~m
    return j  # 返回挑选好的随机数
 
 
'''
门限函数
'''
def clipAlpha(aj,H,L):  # 最大不能超过H，最小不能低于L
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return 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"; 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)  # 门限函数阻止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设置常数项b
                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
 
                # 说明alpha已经发生改变
                alphaPairsChanged += 1
                print "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
 
        #如果没有更新，那么继续迭代；如果有更新，那么迭代次数归0，继续优化
        if (alphaPairsChanged == 0): iter += 1
        else: iter = 0
        print "iteration number: %d" % iter
 
    # 只有当某次优化更新达到了最大迭代次数，这个时候才返回优化之后的alpha和b
    return 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 SVC

logging.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等相关的学习资源、经典和最新的深度学习相关的论文研读，算法和其他互联网技能的学习，概率论、线性代数等高等数学知识的回顾]