历经一个月，终于搞定了SVM(支持向量机)-附源代码解析

前言

其实整体算下来，断断续续的也得有快两个月了(原谅博主比较笨)。中间也有好几次放弃，不想写这篇总结了，但是之前立下的誓言，要将学习到的每一个机器学习算法写成博客总结，一方面呢，检验自己是否真的明白了，另一方面，也希望自己的理解能够帮助到一些人。
源码已上传到GitHub上，有需要的小伙伴自取：源码地址，如果有帮到你，不要吝啬你的小⭐⭐

随便唠唠

1.支持向量机名字由来

因为这个算法的关键点就是支持向量，那么，什么叫支持向量呢？就是离分割超平面(在二维中叫做分割线)最近的那些点。

看上面这个图，中间的实线就是分割超平面，两个虚线所经过的点就是支持向量。这时候大家可能又有疑问了，图b 和图c有什么区别呢？我们可以看到，图b分割虚线和实线之间的距离比较大，而图c分割虚线和实线之间的距离比较小。SVM要找的就是使间隔最大。

2. 为什么很多书籍和老师都说支持向量机比较难呢？

原理我们从面上很好理解，就是找到使支持向量与分割超平面距离最大的实线即可。但是具体的数学推导，就需要大量的数学方面的知识，比如拉格朗日函数、凸优化、KKT条件、对偶函数。这些知识本身不算太难，但是难点就在于你能不能完整的从头到尾将这个思路理下来，这里我指的是用笔把公式推导下来。我认为还是挺难的，最起码对我来说。之前跟着博客大佬推了一遍公式，两天没看，第三天就立刻忘记了推到思路了，以及里面的一些步骤为什么那么做。

网络上有很多真大佬，写的非常详细，推导过程，以及使用的理论讲的非常清楚，但是，我看完之后，我觉得，这种方式对我这种水平比较低的人来说，不太友好，来来回回看了不下五遍，愣是没有形成一个完整的思路。

所以，在这篇文章中，我想以我的角度，去给有同样困惑的人带来一些帮助。这篇文章不会给你解释数学原理，但是我会在文末给你提供相应知识的链接。

3. 读懂源代码对我们的理解帮助非常巨大

我之前在学习K-mean算法和逻辑回归的时候，认为只要懂算法思想就可以了，并且这些算法思想本身也不复杂。

但是在我学习SVM的时候，我发现，源代码对我的帮助太大了。我当时遇到的问题是，我能够将公式推导完毕，但是不知道如何将数学公式理论和实际代码工程相结合，我相信这不仅是我薄弱的地方，也是许多新人薄弱的地方。也许文字对于算法思路的描述能够让我们感受到原来如此，那么，源代码对于算法思路的描述，能够让我们恍然大悟。这篇文章，我就主要带大家关联源代码和数学公式。

4. 学习机器学习算法是不是只明白思想就行，需不需要动手去写代码实现一下？

很多人跟我的答案肯定是一样的，那就是需要去实现一下。

原因大概有两点：其一，我们需要更加深入的了解算法原理。我们学习算法的目的我们得搞清楚，也许我们是做工程，需要了解大量得算法，在具体项目中选择合适的算法；也许我们是做研究，需要了解算法思想以及其优缺点，对算法进行二次开发改进；再或许我们只是为了了解思想，利用这个思想去思考这个世界，无论你属于哪一个，都需要对算法有一个更加深入的了解。现在有很多现成的机器学习工具包，使用起来非常方便，如果你对实现细节不清楚的话，对于特殊问题，就无法根据现有的算法接口进行定制化修改。

其二，提高我们的工程能力。所谓的工程能力，就是从原来一个想法，到最终实现的一个过程。有人说，我是搞研究的，只需要明白原理，不用写太多。那你告诉我，为什么公司招的都是算法工程师呢？即便是你去研究所工作，没有一定的工程能力，也是会有很大劣势的。

支持向量机公式推导

一、推导

我们单独拿出来图b

我们前面也说过了，我们需要找的就是使d值最大的的分割超平面。

首先超平面公式[见参1]：

w 为法向量，b为截距，法向量指向的一侧为正类，另一侧为负类。

分类决策函数公式：

函数间隔公式[见参2]：

函数间隔体现的是分类的正确性，如果h>0,表示分类正确，如果h<0，表示分类错误。同时也体现着数据点到超平面的远近。

几何间隔公式[见参2]：

几何间距体现的点到超平面的距离，这里和图中支持向量到超平面距离d的一样的。

我们同时将w和b等比例变化，超平面并没有发生改变。比如3w和3b。此时函数间隔也会相应的扩大3倍，但是对于几何间隔来讲，并没有发生改变。

我们要获得最大间隔，就需要几何间距达到最大：

前面我们提到，等比例变换w和b，超平面并不会发生改变，但是函数间隔会发生改变。此时，我们令函数间隔h=1，此时公式变为：

等价于：

等价的原因是求解这两个式子的目标都是一样的。

2. 推导

在推导1中，我们推导出了需要求解的公式，除了公式，它还有一些限定条件：

这里限定条件为函数间隔要≥1.原因是当函数间隔=1时，当前的数据点是距离超平面最近的点，这些点就是支持向量。所以，我们的数据点要满足，所有点的函数距离要≥1，也就是我们前面给定的1.。

如何把条件和公式结合在一块呢？这里用到了拉格朗日函数。我们使用拉格朗日乘子将条件和公式写到一起：

λ为拉格朗日乘子，λ≥0.

这里拉格朗日函数有一个性质，就是对偶性，根据对偶性，原始问题L的对偶问题为其极大极小问题

即：

接下来我们对

进行计算：

将计算结果带入到式子中：

接下来计算

我们上面的计算先放在这里，我们先说一下KKT条件

前面我们使用了对偶函数来求解问题的解，但是原问题和对偶问题需要满足一定的条件才能等价，这个条件，就叫做KKT条件

上述公式用在我们这里应该是这样的：

我们将上面的式子整理一下：

接着我们重写之前的分类决策公式**【源代码公式1】**：

3.推导-引入松弛变量

为什么要引入松弛变量？

答：很多数据集都是线性不可分的数据集，这意味着我们无法通过一个线性分类器(超平面)去实现数据的分类。同时，也就不能满足函数间隔大于等于1的约束条件（这个条件的前提数据是线性可分的），于是我们引入一个变量，这个变量就叫所松弛变量，去解决这个问题，具体解决办法看下面的公式。

我们令松弛变量为ε。

我们的公式约束函数和公式可修改为：

这里C称之为惩罚参数，当C的值很大时，对误分类的惩罚就越大。

然后改进后的拉格朗日函数是：

紧接着我们进行极大极小对偶问题处理：

最终化简形式为：

4. 推导3-SMO算法

SMO算法是现在运行速度非常快的强大算法，用于训练SVM。

它的主要思路：

每次循环中选择两个λ进行优化处理，一旦找到一对满足条件的λ，那么就增大其中一个同时减小另外一个。
这里的条件指的是两个λ在间隔边界之外，而第二个条件是两个λ没有进行过处理。

根据条件，我们可写出下面的式子：

当y1≠y2,λ1+λ2=d

当y1=y2,λ1-λ2=d

由此我们可以先画出一个图：

由于我们不知道d的取值，所以，有我们画出了两条线，分别对应着d>0和d<0的情况。

我们设初始可行解为
，,最优解为

，
我们需要对

进行条件约束，根据不等式条件和图中两种不同的等式条件。

我们设

，由此可得出以下两种范围**【源代码公式2】**：

可能大家不理解为什么能写成这样，下面我用一张图来解释一下：

图中我们能够看到，由于我们不知道d的值的正负，所以会出现图中的情况。

对于y1≠y2：

三个绿圈是我们可以取的值，但是为了满足这两种情况和不等式条件的要求，我们取这两种情况的交集，所以，针对下限L，我们取绿点中最大的即

对于上限H，我们取灰色点中的最小值，即

对于y1=y2，同理，这里我就不拿图来再次说明了，直接给出结论了：

对于下限L，我们的取值范围是

对于上限H，我们的取值范围是

接下来，又到了公式推导的时候了，大家挺住！

我们前面选出了两个λ，进行修改，那么，对于相对于这两个λ来讲，其他的λ就是常数了，我们可以对下面公式进行二次推导：

接下来，我们根据λ1和λ2的等式关系，消除掉λ1：

为了方便，我们令y1d=a,y1y2=b

此时，

（这个b和前面决策函数式子中得b不是一回事，我的一个失误，所以，将前面得b在后面以b*的形式出现，大家见谅。）

将我们上面得到的关系带入到w公式中：

还记得我们前面的任务吗？

所以，此时我们对

进行求导：

紧接着，我们根据这个等式可以提出来lλ2：

接着，我们先设几个数(你要是问我为什么设，我只能告诉你是半猜的，根据式子的形式)【源代码公式3】：

然后，我们对λ2的等式进行带入化简**【源代码公式4】**：

之后呢，我们对它进行修剪，使用我们之前设定了L和H**【源代码公式5】**：

接着，我们根据来求解λ1的新值**【源代码公式6】**：

求出来了λ之后，我们接下来要求解b*。

根据KKT条件，我们能够得出以下结论：

在我们完成两个变量的更新之后，当λ满足第二种条件时，则有**【源代码公式7】**：

由于我们对于λ满足条件也是假设，所以对于b也会进行一些处理**【源代码公式8】**：

好了，到这里，整个公式全部推导完毕，终于完事了。如果大家有什么问题，欢迎博客底下留言。

支持向量机源代码分析

我们先来看一张源代码输出得图片：

中间蓝线就是分割超平面，图中带有红圈得是支持向量，也就是距离超平面最近得数据点。

注：代码中使用得alphas来替代λ

1. 随机选择另一个alphas

def selectJrand(i, m):
    j = i                                 #选择一个不等于i的j
    while (j == i):
        j = int(random.uniform(0, m))
    return j

这里我们通过随机数得方式，选取一个与i不相等得数。

2. 修剪λ

def clipAlpha(aj,H,L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

3. 主函数

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):#
    #转换为numpy的mat存储
    dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose()#转置成列向量
    #初始化b参数，统计dataMatrix的维度
    b = 0; m,n = np.shape(dataMatrix)# m=100,n=2
    #初始化alpha参数，设为0
    alphas = np.mat(np.zeros((m,1)))# 100个alpha
    #初始化迭代次数
    iter_num = 0
    #最多迭代matIter次
    while (iter_num < maxIter):
        alphaPairsChanged = 0 #用于记录alpha是否已经进行优化
        for i in range(m):
            #步骤1：计算误差Ei
            #此公式为分类决策公式
            fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            #print("fXi:",fXi)
            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_i成对优化的alpha_j
                j = selectJrand(i,m)
                #步骤1：计算误差Ej
                fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                #保存更新前的aplpha值，使用深拷贝
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();
                #步骤2：计算上下界L和H
                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

                #步骤3：计算eta，即η。
                eta = dataMatrix[i,:]*dataMatrix[i,:].T +dataMatrix[j,:]*dataMatrix[j,:].T-2.0 * dataMatrix[i,:]*dataMatrix[j,:].T
                if eta == 0: print("eta=0"); continue
                #步骤4：更新alpha_j
                alphas[j] += labelMat[j]*(Ei - Ej)/eta
                #步骤5：修剪alpha_j
                alphas[j] = clipAlpha(alphas[j],H,L)
                #abs 函数  返回数字的绝对值
                if (abs(alphas[j] - alphaJold) < 0.00001): print("alpha_j变化太小"); continue
                #步骤6：更新alpha_i
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
                #步骤7：更新b_1和b_2
                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
                #步骤8：根据b_1和b_2更新b
                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("第%d次迭代 样本:%d, alpha优化次数:%d" % (iter_num,i,alphaPairsChanged))
        #更新迭代次数
        if (alphaPairsChanged == 0): iter_num += 1
        else: iter_num = 0
        print("迭代次数: %d" % iter_num)
    return b,alphas

看着特别长，没关系，咱们慢慢说：

方法参数分别为数据集、类别标签、常数C、容错率、最大循环数
首先我们初始化alphas和b。
接下来一个while循环，循环最大次数就是参数maxIter
然后我们使用一个变量alphaPairsChanged来控制alphas是否进行优化
fXi为分类决策公式，对应得公式为【源代码公式1】
接着我们求解误差E，对应公式为【源代码公式2】
通过我们设定得条件，误差值得绝对值是否大于我们参数中得容错率和alpha的取值
随机选择一个不等于i的数j
接着往下，通过if语句，来对下限L和上限H进行赋值
计算eta，也就是我们公式中的η，对应的是公式为【源代码公式3】
更新alpha[ j ],对应的公式为【源代码公式4】
紧接着对alpha[ j ]进行修剪，对应【源代码公式5】
更新alpha[ i ]的值，对应公式为【源代码公式6】
更新b1、b2的值，对应公式为【源代码公式7】
根据b1、b2更新b，对应公式为【源代码公式8】

基本上到这里就全部结束了，非常感谢大家能够耐心的看完，如果有需要完整源代码，请在博客底下留言。第一次写这么长的文章，如哪里有不合适的地方，欢迎大家批评指正。

参考资料：

超平面公式推导
《统计学习方法第二版》 p133
支持向量机原理篇之手撕线性SVM
如何理解超平面？
《机器学习实战》

更多推荐

历经一个月,终于搞定了SVM(支持向量机)-附源代码解析

历经一个月，终于搞定了SVM(支持向量机)-附源代码解析