本文已参加「新人创造礼」活动，一起开启创造之路。

前语

博主大大小小参加过数十场数学建模竞赛，SVM经常在各种建模竞赛的优秀论文上见到该模型，一般直接运用SVM算法是比较少的，现在都是在此根底理论之上提出优化算法。可是SVM的根底理论是十分重要的思想，放眼整个分类算法中，SVM是最好的现成的分类器。这儿说的‘现成’指的是分类器不加修正即可直接运用。在神经网络没有呈现之前，SVM的优化模型能够算得上是猜测分类神器了，在机器学习中SVM依旧是最为知名的算法之一了，本篇博客将致力于将SVM算法以及原理每一个知识点都讲明白，期望没有讲明白的点咱们能够在评论区指出。

一、引论

咱们运用SVM支撑向量机一般用于分类，得到低错误率的成果。SVM能够对练习集意外的数据点做出很好的分类决议计划。那么首先咱们应该从数据层面上去看SVM到底是怎样做决议计划的，这儿来看这样一串数据集调集在二维平面坐标系上描绘的图：

现在咱们需求考虑，是否能够画出一条直线将圆形点和星星点分隔。像first第一张图片来看，圆点和星点就分的很开，很简单就能够在图中画出一条直线将两组数据分隔。而看第二张图片，圆点和星点几乎都聚合在一起，要区分的话十分困难。

咱们要划线将他们区分隔来的话，有有很多条能够画，可是咱们难以找到一条最好区分度最高的线条将它们几乎彻底区分。那么在此咱们需求了解两个关于数据集的基本概念：

二、理论铺垫

线性可分性（linear separability）

第一张图片的数据便是线性可分的，肉眼可见的线性可分。而第二张图片的数据一看便是线性不行分的。当然咱们这样直爽的观测也没有错，可是这仅是建立在二维平面数据可视化的根底之上，若是特征维度更高，例如：

超出咱们可视化的技能之外咱们怎样判别他们是否可分呢？

一切要判别线性可分性需求对线性可分性下个界说：

在分类问题中给定输入数据和学习方针：,,其间输入数据的每个样本都包含多个特征并由此构成特征空间（feature space）：,而学习方针为二元变量表明负类（negative class）和正类（positive class）。

若输入数据地点的特征空间存在作为决议计划鸿沟（decision boundary）的超平面将学习方针按正类和负类分隔，并使任意样本的点到平面距离(point to plane distance)大于等于1： decision boundary：

point to plane distance:

则称该分类问题具有线性可分性，参数,分别为超平面的法向量和截距。

到这儿肯定有很多小伙伴会问啥是决议计划鸿沟什么是超平面啊，莫慌立刻给出界说：

超平面

而对机器学习来说，触及的多是高维空间（多维度）的数据分类，高维空间的SVM，即为超平面。机器学习的最终意图便是要找到最适宜的（也即最优的）一个分类超平面（Hyper plane），从而应用这个最优分类超平面将特征数据很好地区分为两类。

咱们看图就懂了：

在数学中，超平面（Hyperplane）是n维欧式空间中余维度等于1的线性质空间。这是平面中的直线、空间中的平面之推行。设F为域（可考虑）。

n维空间中的超平面是由方程：

界说的子集，其间是不全为零的常数。

在线性代数的头绪下，F-矢量空间V中的超平面是指形如：

的子空间，其间是任一非零的线性映射。

在射影几许中，同样可界说射影空间中的超平面。在齐次坐标（)下超平面可由以下方程界说：

，其间是不全为零的常数。

总而言之：超平面H是从n维空间到n-1维空间的一个映射子空间，它有一个n维向量和一个实数界说。设d是n维欧式空间R中的一个非零向量，a是实数，则R中满意条件dX=a的点X所组成的调集称为R中的一张超平面。

决议计划鸿沟

SVM是一种优化的分类算法，其动机是寻找一个最佳的决议计划鸿沟，使得从决议计划鸿沟与各组数据之间存在margin，而且需求使各侧的margin最大化。那么这个决议计划鸿沟便是不同类之间的界限。

总而言之：在具有两个类的核算分类问题中，决议计划鸿沟或决议计划外表是超平面，其将根底向量空间划分为两个调集，一个调集。分类器将决议计划鸿沟一侧的一切点分类为归于一个类，而将另一侧的一切点分类为归于另一个类。

支撑向量（supportvector）

在了解了超平面和决议计划鸿沟咱们发现SVM的核心使命是找到一个超平面作为决议计划鸿沟。那么满意该条件的决议计划鸿沟实际上结构了2个平行的超平面作为距离鸿沟以判别样本的分类：

一切在上距离鸿沟上方的样本归于正类，鄙人距离鸿沟下方的样本归于负类。两个距离鸿沟的距离

被界说为边距（margin），坐落距离鸿沟上的正类和负类样本为支撑向量（support vector）。

丢失函数（loss function）

拿咱们的first的数据集来看：

不论咱们怎样画直线总有丢失的点没有正确分类到 。在一个分类问题不具有线性可分性时，运用超平面作为决议计划鸿沟会带来分类丢失，即部分支撑向量不再坐落距离鸿沟上，而是进入了距离鸿沟内部，或落入决议计划鸿沟的错误一侧。

丢失函数能够对分类丢失进行量化，其按数学含义能够得到的办法是0-1丢失函数：

0-1丢失函数是一个不连续的分段函数，不利于求解其最小化问题，因而在应用可结构其署理丢失（surrogate loss）。署理丢失是与原丢失函数具有相合性（consistency）的丢失函数，最小化署理丢失所得的模型参数也是最小化原丢失函数的解。也便是说咱们需求丢失函数核算最小的决议计划鸿沟。这儿给出二元分类（binary classification）中0-1丢失函数的署理丢失：

SVM运用的是铰链丢失函数:.

经历危险（empirical risk）与结构危险（structural risk）

按核算学习理论，分类器在经过学习并应用于新数据时会产生危险，危险的类型可分为经历危险和结构危险。

式中表明分类器，经历危险由丢失函数界说，描绘了分类器所给出的分类成果的精确程度；结构危险由分类器参数矩阵的范数界说，描绘了分类器自身的杂乱程度以及安稳程度，杂乱的分类器简单产生过拟合，因而是不安稳的。若一个分类器经过最小化经历危险和结构危险的线性组合以确定其模型参数：

则对该分类器的求解是一个正则化问题，常数是正则化系数。当=2时，该式被称为正则化或Tikhonov正则化（Tikhonov regularization）。SVM的结构危险表明，在线性可分问题下硬鸿沟SVM的经历危险能够归0，因而其是一个彻底最小化结构危险的分类器；在线性不行分问题中，软鸿沟SVM的经历危险不行归0，因而其是一个正则化分类器，最小化结构危险和经历危险的线性组合。

核办法

一些线性不行分的问题或许是非线性可分的，即特征空间存在超平面将正类和负类分隔。运用非线性函数能够将非线性可分问题从原始的特征空间映射至更高维的空间（希尔伯特空间（Hilbert space））从而转化为线性可分问题，此刻作为决议计划鸿沟的超平面表明如下：

式中为映射函数。因为映射函数具有杂乱的办法，难以核算其内积，因而可运用核办法（kernel method），即界说映射函数的内积为核函数：

以回避内积的显式核算。

常见的核函数

当多项式核的阶为1时，其被称为线性核，对应的非线性分类器退化为线性分类器。RBF核也被称为高斯核（Gaussian kernel），其对应的映射函数将样本空间映射至无限维空间。

三、算法流程

在了解了上述算法原理今后，咱们便对SVM算法有了个大致明晰的认知，那么SVM是怎样挖掘到丢失最小的超平面的呢？

SVM的求解能够运用二次凸优化问题的数值办法，例如内点法（Interior Point Method, IPM）和序列最小优化算法（Sequential Minimal Optimization, SMO），在拥有满足学习样本时也可运用随机梯度下降。

这儿咱们选用SMO序列最小优化算法核算：

SMO序列最小优化算法

SMO是一种坐标下降法（coordinate descent）以迭代办法求解SVM的对偶问题，其规划是在每个迭代步挑选拉格朗日乘子中的两个变量并固定其他参数，将原优化问题简化至一维子可行域，此刻约束条件有如下等价办法：

将上式右侧带入SVM的对偶问题并消去求和项中的能够得到仅关于的二次规划问题，该优化问题有闭式解能够快速核算。在此根底上，SMO有如下核算框架:

1. 初始一切化拉格朗日乘子；
2. 辨认一个不满意KKT条件的乘子，并求解其二次规划问题；
3. 重复执行上述过程直到一切乘子满意KKT条件或参数的更新量小于设定值 能够证明，在二次凸优化问题中，SMO的每步迭代都严格地优化了SVM的对偶问题，且迭代会在有限步后收敛于全局极大值。SMO算法的迭代速度与所选取乘子对KKT条件的偏离程度有关，因而SMO一般选用启发式办法选取拉格朗日乘子。

Python sklearn代码实现：

为了方便展示作用，咱们依旧用first的数据集进行：

sklearn.svm.SVC语法格式为：

class sklearn.svm.SVC(  *,
                        C=1.0, 
                        kernel='rbf',
                        degree=3, 
                        gamma='scale', 
                        coef0=0.0, 
                        shrinking=True, 
                        probability=False, 
                        tol=0.001, 
                        cache_size=200, 
                        class_weight=None, 
                        verbose=False, 
                        max_iter=- 1, 
                        decision_function_shape='ovr', 
                        break_ties=False, 
                        random_state=None)

该实现根据 libsvm，更多详细信息能够去官网上查阅：sklearn.svm.SVC.

# 导入模块
import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm, datasets
# 鸢尾花数据
iris = datasets.load_iris()
X = iris.data[:, :2] # 为便于绘图仅挑选2个特征
y = iris.target
# 测验样本（制作分类区域）
xlist1 = np.linspace(X[:, 0].min(), X[:, 0].max(), 200)
xlist2 = np.linspace(X[:, 1].min(), X[:, 1].max(), 200)
XGrid1, XGrid2 = np.meshgrid(xlist1, xlist2)
# 非线性SVM：RBF核，超参数为0.5，正则化系数为1，SMO迭代精度1e-5, 内存占用1000MB
svc = svm.SVC(kernel='rbf', C=1, gamma=0.5, tol=1e-5, cache_size=1000).fit(X, y)
# 猜测并制作成果
Z = svc.predict(np.vstack([XGrid1.ravel(), XGrid2.ravel()]).T)
Z = Z.reshape(XGrid1.shape)
plt.contourf(XGrid1, XGrid2, Z, cmap=plt.cm.hsv)
plt.contour(XGrid1, XGrid2, Z, colors=('k',))
plt.scatter(X[:, 0], X[:, 1], c=y, edgecolors='k', linewidth=1.5, cmap=plt.cm.hsv)
plt.show()

Python源代码实现+手写字辨认分类：

from numpy import *
def selectJrand(i, m):  # 在某个区间范围内随机挑选一个整数
    # i为第一个alhpa的下表，m是一切alpha的数目
    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):  # 在数值太大的时分对其进行调整
    # aj是H是下限，是L的上限
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj
def kernelTrans(X, A, kTup):  # calc the kernel or transform data to a higher dimensional space
    # X是数据，A是
    m, n = shape(X)
    K = mat(zeros((m, 1)))
    # 这儿为了简单，咱们只内置了两种核函数，可是原理是相同的，需求的话再写其他类型便是了
    # 线性核函数：k(x,x_i)=x*x_i,它不需求再传入参数。
    if kTup[0] == 'lin':
        K = X * A.T  # linear kernel，线性核函数
    elif kTup[0] == 'rbf':  # 高斯核函数，传入的参数作为detla
        for j in range(m):
            deltaRow = X[j, :] - A
            K[j] = deltaRow * deltaRow.T  # 2范数
        K = exp(K / (-1 * kTup[1] ** 2))  # divide in NumPy is element-wise not matrix like Matlab
        # numpy中，/表明对矩阵元素进行核算而不是核算逆（MATLAB）
    else:
        raise NameError('Houston We Have a Problem -- \
    That Kernel is not recognized')  
    return K
# 界说了一个类来进行SMO算法
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters
        # kTup储存核函数信息，它是一个元组，元组第一个元素是一个描绘核函数类型的字符串，其他两个元素是核函数或许需求的参数
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]  # m是样本个数，也是a的个数
        self.alphas = mat(zeros((self.m, 1)))  # 初始化a序列，都设置为0
        self.b = 0
        # 第一列给出的是eCache是否有用的标志位，而第二位是实际的E值
        self.eCache = mat(zeros((self.m, 2)))  # first column is valid flag
        # 假如第一位是1，阐明现在的这个Ek是有用的
        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 calcEk(oS, k):  # 核算第k个样本的Ek
    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
    # 选定了a_1之后挑选a_2
    # 挑选a_2
    maxK = -1;
    maxDeltaE = 0;
    Ej = 0  # 挑选|E1-E2|最大的E2并回来E2和j
    # 先将E1存进去，以便于之后的统一化进行
    oS.eCache[i] = [1, Ei]  # set valid #choose the alpha that gives the maximum delta E
    '''
    numpy函数回来非零元素的目录。
    回来值为元组， 两个值分别为两个维度， 包含了相应维度上非零元素的目录值。   
    能够经过a[nonzero(a)]来获得一切非零值。
    .A的意思是：getArray(),也便是将矩阵转换为数组
    '''
    # 获取哪些样本的Ek是有用的，ValidEcacheList里边存的是一切有用的样本行Index
    validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
    # 对每一个有用的Ecache都比较一遍
    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
            # 假如选到了和a1相同的，就持续，因为a1和a2有必要选不相同的样本
            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]
    # 核算Ek并保存在类oS中
# 内循环寻找适宜的a_2
def innerL(i, oS):
    Ei = calcEk(oS, i)  # 为什么这儿要从头算呢？因为a_1刚刚更新了，和存储的不相同
    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)):
        # 假如a1选的适宜的话，不适宜就直接结束了
        # 剩下的逻辑都相同，只不过不是运用x_ix_j，而是运用核函数
        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 smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):  # 默许核函数是线性，参数为0（那便是它本身了）
    # 这个kTup先不论，之后用
    oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup)  # 初始化这一结构
    iter = 0
    # entireSet是操控开关，一次循环对一切样本点都遍历，第二次就只遍历非鸿沟点
    entireSet = True;
    alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged: int = 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
            # 把大于0且小于C的a_i挑出来，这些是非鸿沟点，只从这些点上遍历
            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
        elif (alphaPairsChanged == 0):  # 假如alpha没有更新,核算全样本遍历
            entireSet = True
        print("iteration number: %d" % iter)
    return oS.b, oS.alphas
# 使用SVM进行分类，回来的是函数距离，大于0归于1类，小于0归于2类。
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
# 将图画转换为向量
def img2vector(filename):
    # 一共有1024个特征
    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 = []
    # 使用listdir读文件名，这儿的label写在了文件名里边
    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))
testDigits(kTup=('rbf', 20))

点重视，防走丢，如有疏忽之处，请留言指教，非常感谢

以上便是本期全部内容。我是fanstuck ，有问题咱们随时留言评论 ，咱们下期见

我正在参加技能社区创造者签约计划招募活动，点击链接报名投稿。