《Machine Learning in Action》—— 剖析支持向量机，单手狂撕线性SVM

手撕机器学习系列文章就暂时更新到此吧，目前已经完成了支持向量机SVM、决策树、KNN、贝叶斯、线性回归、Logistic回归，其他算法还请允许Taoye在这里先赊个账，后期有机会有时间再给大家补上。

更新至此，也是收到了部分读者的好评。虽然不多，但还是非常感谢大家的支持，希望每一位阅读过的读者都能够有所收获。

该系列文章的全部内容都是Taoye纯手打，也是参考了不少书籍以及公开资源，系列总字数在15W左右（含源码），后期会再慢慢填补，更多的技术文章可以来访Taoye的公众号：玩世不恭的Coder。文档可以随意传播，但注意不可修改其中的内容。

如果文章中有任何不懂的问题，都可以直接提出，Taoye看见后会第一时间回复，同时欢迎大家来此私密地向Taoye催更：玩世不恭的Coder，公众号也有Taoye的私人联系方式。有些话，Taoye只能在那里和你偷偷地说 (#`O′)

为了提高大家的阅读体验，手撕机器学习系列文章Taoye已经整理成PDF和和HTML，阅读效果都很不错，在公众号【玩世不恭的Coder】下回复【666】即可免费获取。

前面在写NumPy文章的结尾处也有提到，本来是打算按照《机器学习实战 / Machine Learning in Action》这本书来手撕其中代码的，但由于实际原因，可能需要先手撕SVM了，这个算法感觉还是挺让人头疼，其中内部太复杂了，涉及到的数学公式太多了，也涉及到了许多陌声的名词，如：非线性约束条件下的最优化、KKT条件、拉格朗日对偶、最大间隔、最优下界、核函数等等，天书或许、可能、大概就是这样的吧。

记得与SVM初次邂逅是在17年，那个时候的自己年少轻狂，视图想着站在巨人的肩膀上，把所有自己感兴趣的内容都搞懂，深入骨髓的那种。但后来残酷的现实让我明白一个道理：**你知道的越多，你不知道的也越多。**而且那个时候自己也没有能力、资源和机会去深入SVM内部，完全无法理解SVM的内部原理。所以，当时自己对SVM的收获只有一个：SVM主要是用来做分类任务的，仅此而已。

第二次接触SVM是在准备考研复试吧，当时复试并没有给出具体内容和范围，而且自己也还是个初出茅庐的小子，对这种所谓的复试有种莫名的恐惧感。也只有从上届学长学姐的口中，得知复试的时候老师会考究学生是否有科研的潜力，所以最好把机器学习熟知一下。那个时候也是处于新冠疫情的紧张时期嘛，就疯狂补习机器学习的内容，其中就包括支持向量机——SVM，主要的学习渠道是吴恩达老师的机器学习课程，感觉讲的的确不错，非常适合我这种菜鸟级选手学习。当时也算是对SVM有了一定的认识吧，也大致了解了SVM的工作原理，当然了，也只是对SVM有了个的浅显的认识，没有手撕SVM的过程，也没有完全把它整明白。尽管如此，复试的过程依然被面试导师锤的体无完肤，除了问了机器学习相关内容之外，编译原理等一些专业知识对于我这个贸易专业的学生来讲可太痛苦了，之前也没有接触过，全程阿巴阿巴。想到这，眼角又又。。。

第三次面对SVM也就是现在了，想着无论如何也要打通我的任督二脉，一定要搞清楚SVM的来龙去脉，也要像面试老师捶我那样，把SVM往死里锤。于是有了下文学习SVM之后的总结，一方面算是重新梳理一遍SVM，另一方面也希望来访的读者能够有所收获。

对于刚刚接触SVM的读者，Taoye主要有以下几条建议，也是我学习SVM过程中的一个小总结吧：

SVM内部的数学公式很多，但请不要未战先怯，犯下兵家大忌。无论是阅读该篇文章也好，学习其他相关SVM资源也罢，还请诸君耐心、认真且完整的看完。
SVM的原理过程会涉及到很多的符号或记号，一定要梳理清楚他们所代表的含义。另外，推导过程中会存在很多的向量或矩阵，一定要明白其中shape，这一点可能在不同的资料中会有不一样的处理方式。
在阅读的同时，一定要拿出稿纸手动推演SVM的过程，尽可能明白整个过程的来龙去脉，有不明白的地方可以留言或查找其他相关资料来辅助自己的理解。
阅读一遍或许有不少知识不理解，但多阅读几遍相信一定会有不少收获

本文参考了不少书籍资料以及许多大佬的技术文章，行文风格尽可能做到通俗易懂，但其中涉及到的数学公式在所难免，还请诸读者静下心来慢慢品尝。由于个人水平有限，才疏学浅，对于SVM也只是略知皮毛，可能文中有不少表述稍有欠妥、有不少错误或不当之处，还请诸君批评指正。

我是Taoye，爱专研，爱分享，热衷于各种技术，学习之余喜欢下象棋、听音乐、聊动漫，希望借此一亩三分地记录自己的成长过程以及生活点滴，也希望能结实更多志同道合的圈内朋友，更多内容欢迎来访微信公主号：玩世不恭的Coder

符号说明

符号	说明
xix_ixi	表示单个样本，其中包含多个属性，显示为一个行向量
x(i)x^{(i)}x(i)	表示单个样本中的某个属性特征
yiy_iyi	表示单个样本所对应的标签（具体分类），为整数，非1即-1
wTw^TwT	表示的是权值行向量，其中w=(w1,w2,...,wm)Tw=(w_1,w_2,...,w_m)^Tw=(w1,w2,...,wm)T，也是所需要训练的参数
bbb	表示的是决策面中的bbb，也是所需要训练的参数
γ^\hat{\gamma}γ^	表示函数间隔，具体解释见下文
γ\gammaγ	表示几何间隔，具体解释见下文
α\alphaα	α=(α1,...,αN)\alpha=(\alpha_1,...,\alpha_N)α=(α1,...,αN)表示拉格朗日乘子
KijK_{ij}Kij	在这篇文章中表示线性核函数

关于上述的符号说明，仅仅只是本篇文章的一部分，其他符号可通过正文了解。上述符号可能部分暂时不懂，但没关系，读者可在阅读的过程中，随时返回来查看，即可理解每个符号所代表的意义。

一、SVM是什么

关于SVM是什么，之前在Byte Size Biology上看到有篇文章很好的解释了SVM，在知乎上也有一位名叫“简之”的用户通过故事的形式来将其进行转述，通俗易懂，很好的向首次接触SVM的读者介绍了SVM能干嘛。

Support Vector Machines explained well（可能需要翻墙）：http://bytesizebio.net/2014/02/05/support-vector-machines-explained-well/
支持向量机（SVM）是什么意思？：https://www.zhihu.com/question/21094489/answer/86273196

油管上也有更为直观认识SVM的短视频（翻墙）：https://www.youtube.com/watch?v=3liCbRZPrZA

总结一哈：

对于一个二分类问题，样本数据是线性可分的，我们则需要通过一根直线（也就是上述例子当中枝条）将两个不同种类的样本进行分离。按道理来讲，我们已经实现了需求，但是，这根枝条的具体摆放位置可能有无数多个，而我们的最终目的是将枝条摆放到一个最好的位置，从而当我们引入了一些新样本的时候，依然能最好很好将两类的数据分离开来，也就是说我们需要将模型的“泛化性能”最大化。之前也看到过一个例子（来源忘了），这里分享下，大概就是讲：我们在通过悬崖吊桥的时候，会不自觉的尽可能往中间走，因为这样的话，当左右起风的时候，虽然我们的位置会左右稍微移动，但还不至于跌落悬崖。而越靠近边缘，风险就越大，就是这么个道理。而寻找最大“泛化性能”的过程，就是将枝条摆放在距离小球最远的位置，而小球相对于枝条的位置就是**“间隔”，而我们要做的就是将这间隔最大化**。

上述仅仅是对于线性可分数据分类的一种处理方式，但有的时候，理想是美好的，现实却是残酷的。在实际样本数据中，大多都是线性不可分的，也就是说我们无法找到合适位置的枝条将不同分类的数据分离开来，更别提“间隔最大化”了。这个时候，我们就需要将桌上的小球排起，然后用一个平面（纸）将不同分类小球分离开来。也就是说，我们将低维度映射到了高纬度，这样对于小球的分类就更加容易了。

再之后，无聊的大人们，把这些球叫做「data」，把棍子叫做「classifier」, 最大间隙trick 叫做「optimization」，拍桌子叫做「kernelling」, 那张纸叫做「hyperplane」。

二、线性可分SVM与间隔最大化

我们先来看具体看看线性可分的二分类问题。

假如说，我们这里有一堆样本，也就是我们常说的训练集S={(x1,y1),(x1,y1),(x3,y3),...,(xn,yn)}S=\{(x_1, y_1), (x_1, y_1), (x_3, y_3), ..., (x_n, y_n)\}S={(x1,y1),(x1,y1),(x3,y3),...,(xn,yn)}，且xix_ixi表示的是样本iii的属性特征向量，其内部有多个不同属性，这里我们不妨指定每个样本含有两个属性特征，也就是说xi=(xi1,xi2)Tx_i=(x_{i1}, x_{i2})^Txi=(xi1,xi2)T（之所以用列向量表示，主要是方便后面超平面的构建）。而yiy_iyi表示的是每个样本中所有属性特征所对应的标签，由于这里的问题属性二分类，所以yiy_iyi的值只存在两种。为此，我们可以通过Matplotlib在平面直角坐标系中绘制其图像，初步观察其分布规律，具体代码如下：

import numpy as np
import pylab as pl
from sklearn import svm%matplotlib inline
print(np.__version__)"""Author: Taoye微信公众号：玩世不恭的Coder
"""
if __name__ == "__main__":# np.random.seed(100)        # 可自行设置随机种子每次随机产生相同的数据X_data = np.concatenate((np.add(np.random.randn(20, 2), [3, 3]),       np.subtract(np.random.randn(20, 2), [3, 3])),axis = 0)      # random随机生成数据，+ -3达到不同类别数据分隔的目的 Y_label = np.concatenate((np.zeros([20]), np.ones([20])), axis = 0)svc = svm.SVC(kernel = "linear")svc.fit(X_data, Y_label)w = svc.coef_[0]      ww = -w[0] / w[1]                # 获取权重wxx = np.linspace(-6, 6)yy = ww * xx - (svc.intercept_[0]) / w[1]   # intercept_获取结截距b_down = svc.support_vectors_[0]         # 得到对应的支持向量yy_down = ww * xx + (b_down[1] - ww * b_down[0])b_up = svc.support_vectors_[-1]yy_up = ww * xx + (b_up[1] - ww * b_up[0])pl.plot(xx, yy, "k-"); pl.plot(xx, yy_down, "k--"); pl.plot(xx, yy_up, "k--")pl.scatter(X_data[:, 0], X_data[:, 1], c = Y_label, cmap = pl.cm.Paired)pl.show()

执行代码，可以绘制如下所示图片，注意：以上代码每次运行都会随机产生不同的二分类数据集，如想每次随机产生相同的数据集，可自行配置np.random.seed随机种子；另外，还有一点需要需要说明的是，上述代码使用到了NumPy，关于NumPy的使用，可自行参考之前写的一篇文章：print( “Hello，NumPy！” )

如上图所示，我们可以发现，棕色代表一类数据集，此时标签yi=1y_i=1yi=1，蓝色代表另一类数据集，标签yi=0y_i=0yi=0，而要想将上图两类数据集分离开来，显然不止一条直线，与上图两条虚线平行且居其之间的任意一条直线都能达到此目的。在这无数条直线中，要数上图中的三条最为特殊，中间实线居于两条虚线中间，我们一般称其为**“决策面”或“超平面”，而其所表示的方程，我们一般称作“决策方程”或“超平面方程”**，在这里可以表示为wTx+b=0w^Tx+b=0wTx+b=0。（下面会推导）

从上图我们还可以观察得到，在所有样本数据集中，虚线上的样本距离决策面最近，我们把这种比较特殊的样本数据一般称之为**“支持向量”，而支持向量到决策面之间的距离称为“间隔”**。我们不难发现，决策面的位置主要取决于支持向量，而与支持向量除外的数据样本没有关系。（因为支持向量的确定就已经确定了最大间隔）

关于上述提到的一些关于SVM的名词概念，在正式推演之前，还是有必要理解清楚的。

前面我们也有提到，关于能将两类不同数据集相互分隔开来的直线有无数种，而我们要想在这无数种直线找到一条最合适的，也就是达到一个间隔最大化的目的，这就是一个**“最优化”**问题。而最优化问题，我们需要了解两个因素，分别是目标函数和优化对象。既然我们是想要达到间隔最大化的目标，那么目标函数自然就是间隔，而优化对象就是我们的决策面方程。所以，我们首先需要用数学来明确间隔和决策面方程：

我们知道，在平面直角坐标系中，一条直线可以用其一般式方程来来表示：

Ax+By+C=0Ax+By+C=0 Ax+By+C=0

而根据上述图像，我们可以知道，横纵坐标代表的意义是一个样本的不同属性特征，而标签则是通过颜色（棕色和蓝色）来进行表示。所以上述的直线的一般式方程中的x、yx、yx、y表示的就是一个样本的两种属性特征，为了方便理解，我们不妨将其修改为x(1)、x(2)x^{(1)}、x^{(2)}x(1)、x(2)，并将A、B、CA、B、CA、B、C替换为w1、w2、bw_1、w_2、bw1、w2、b，对此，我们可以将上述方程向量化，得到：

KaTeX parse error: Undefined control sequence: \ at position 151: …rix} \right) \̲ ̲ +b=0 \tag{2…

令w=(w1,w2)T,x=(x(1),x(2))w=(w_1, w_2)^T, x=(x^{(1)}, x^{(2)})w=(w1,w2)T,x=(x(1),x(2))，则上述指向方程最终可以表示为:

wTx+b=0(2-2)w^Tx+b=0 \tag{2-2} wTx+b=0(2-2)

式（1-2）表示的就是我们优化问题的优化对象，也就是决策面方程。我们知道在平面直角坐标系中，一条直线可以通过其斜率和截距来确定，而在决策面方程里，我们不难得到：www确定了决策面的方向（斜率），而bbb确定了决策面的截距。既然我们已经得到了优化问题的优化对象——决策面方程，那么接下来就需要得到目标函数——间隔的表达式了。

在此，我们需要引入函数间隔和几何间隔的概念了。

一般来讲，我们可以通过样本点到决策面的距离来表示分类预测的正确程度，距离决策面越远，一般分类就越正确（可根据图像自行理解），而在超平面wTx+b=0w^Tx+b=0wTx+b=0确定的情况下，我们可以通过∣wTxi+b∣|w^Tx_i+b|∣wTxi+b∣的值来描述样本xix_ixi距超平面的远近（注意：这里是描述远近。而不是确切的距离）。我们知道，样本有不同的分类，所以有的时候∣wTxi+b∣|w^Tx_i+b|∣wTxi+b∣的符号具有不确定性，所以我们可以通过wTxi+bw^Tx_i+bwTxi+b和yiy_iyi的符号来判断分类结果的正确性。也就是说可以通过yi(wTxi+b)y_i(w^Tx_i+b)yi(wTxi+b)值的大小和符号来判断一个样本分类的正确性和正确程度，这个就是我们的函数间隔（这个概念务必要理解清楚），我们不妨用γi^\hat{\gamma_i}γi^来表示：

γi^=yi(wTxi+b)(2-3)\hat{\gamma_i}=y_i(w^Tx_i + b) \tag{2-3} γi^=yi(wTxi+b)(2-3)

而我们知道，上述的i=1,2,...,Ni=1,2,...,Ni=1,2,...,N，表示的是有NNN个样本，而在NNN个样本中，固然存在一个样本使得该值达到最小，该样本也就是我们前面所说的支持向量，我们把支持向量到超平面的函数间隔定义为γ^\hat{\gamma}γ^，则：

γ^=mini=1,...,Nγi^(2-4)\hat{\gamma} = \mathop{min}\limits_{i=1,...,N}\hat{\gamma_i} \tag{2-4} γ^=i=1,...,Nminγi^(2-4)

我们只有函数间隔还不够，函数间隔描述的仅仅是样本分类的正确性和正确程度，而不是确切的间隔。当我们成比例的更改www和bbb的时候，决策面并没有发生任何改变，而函数间隔却发生了改变，所以我们需要对www进行约束，从而当无论www和bbb怎么成比例变动的时候，我们的“间隔”都不会发生改变，也就是进行将www规范化，由函数间隔规范后的间隔我们称之为几何间隔，我们不妨用γ\gammaγ来表示，则某个样本到超平面的几何间隔为γi\gamma_iγi：

γi=∣wTxi+b∣∣∣w∣∣(2-5)\gamma_i=\frac{|w^Tx_i + b|}{||w||} \tag{2-5} γi=∣∣w∣∣∣wTxi+b∣(2-5)

我们可以将式子（1-5）理解成点到直线的距离公式（这个在中学时期就学过的）。对于这个二分类问题，我们不妨将二分类的标签定为 yi={1,−1}y_i=\{1, -1\}yi={1,−1}，则γi\gamma_iγi可以表示为（之所以乘以yiy_iyi，主要是把绝对值符号去掉，这样就能描述所有样本了，而省去了分类讨论的麻烦）：

γi=yiwTxi+b∣∣w∣∣(2-6)\gamma_i=y_i\frac{w^Tx_i + b}{||w||} \tag{2-6} γi=yi∣∣w∣∣wTxi+b(2-6)

而我们知道，上述的i=1,2,...,Ni=1,2,...,Ni=1,2,...,N，表示的是有NNN个样本，而在NNN个样本中，固然存在一个样本使得该值达到最小，该样本也就是我们前面所说的支持向量，我们把支持向量到超平面的几何间隔定义为γ^\hat{\gamma}γ^，则：

γ=mini=1,...,Nγi(2-7)\gamma=\mathop{min}\limits_{i=1,...,N}\gamma_i \tag{2-7} γ=i=1,...,Nminγi(2-7)

通过上述对函数间隔和几何间隔的分析，我们可以得到他们之间的关系：

γi=γi^∣∣w∣∣，即γ=γ^∣∣w∣∣(2-8)\gamma_i=\frac{\hat{\gamma_i}}{||w||}，即\gamma=\frac{\hat{\gamma}}{||w||} \tag{2-8} γi=∣∣w∣∣γi^，即γ=∣∣w∣∣γ^(2-8)

自此，我们已经分析得到了该优化问题的优化对象——决策面方程和目标函数——间隔（几何间隔）。在之前，我们提到了支持向量的概念，那么支持向量具有什么特性呢？细想一下不难发现，支持向量到决策平面的间隔是最近的，即满足如下式子：

γi=yiwTxi+b∣∣w∣∣≥γ，即yi(wTxi+b)≥γ∣∣w∣∣=γ^(2-9)\gamma_i=y_i\frac{w^Tx_i + b}{||w||}\geq\gamma ，即y_i(w^Tx_i + b)\geq\gamma||w||=\hat{\gamma} \tag{2-9} γi=yi∣∣w∣∣wTxi+b≥γ，即yi(wTxi+b)≥γ∣∣w∣∣=γ^(2-9)

对此，我们就可以通过数学来表达该优化问题：

maxw,bγ^∣∣w∣∣s.t.yi(wTxi+b)≥γ^，i=1,2,3,...,N(2-10)\begin{aligned} & \mathop{max}\limits_{w,b} \quad \frac{\hat{\gamma}}{||w||} \tag{2-10} \\ & s.t. \quad\ y_i(w^Tx_i + b)\geq\hat{\gamma}，i=1,2,3,...,N \end{aligned} w,bmax∣∣w∣∣γ^s.t. yi(wTxi+b)≥γ^，i=1,2,3,...,N(2-10)

前面，我们提到了，函数间隔描述的仅仅是样本分类的正确性和正确程度，而不是确切的间隔。当我们成比例的更改www和bbb的时候，函数间隔虽然发生了改变，但并不会影响我们的几何间隔——目标函数，也就是说，此时产生了一个等价的最优化问题。为了方便描述问题，我们不妨取γ^=1\hat{\gamma}=1γ^=1。另外，我们注意到目标函数maxw,bγ^∣∣w∣∣\mathop{max}\limits_{w,b} \frac{\hat{\gamma}}{||w||}w,bmax∣∣w∣∣γ^可以等价于minw,b12∣∣w∣∣2\mathop{min}\limits_{w,b}\frac{1}{2}||w||^2w,bmin21∣∣w∣∣2，（注意，仅仅是等价。而之所以等价为二分之一、平方的形式，主要是方便后期的求导操作），对此，我们可以将数学表达的优化问题转化为如下形式：
minw,b12∣∣w∣∣2s.t.yi(wTxi+b)−1≥0，i=1,2,3,...,N(2-11)\begin{aligned} & \mathop{min}\limits_{w,b} \quad \frac{1}{2}||w||^2 \\ & s.t. \quad\ y_i(w^Tx_i + b)-1\geq0，i=1,2,3,...,N \tag{2-11} \end{aligned} w,bmin21∣∣w∣∣2s.t. yi(wTxi+b)−1≥0，i=1,2,3,...,N(2-11)

关于如上提到的决策平面和间隔等概念，我们可以通过下图来直观感受下（理解清楚）：

至此，我们已经得到了该优化问题的数学表达，我们不妨通过一个小例子来检测下：

例子来源：李航-《统计学习方法》第七章内容

例子1：已知一个如图所示的训练数据集，其正例点是x1=(3,3)T,x2=(4,3)Tx_1=(3, 3)^T,x_2=(4, 3)^Tx1=(3,3)T,x2=(4,3)T，负例点为x3=(1,1)Tx_3=(1, 1)^Tx3=(1,1)T，试求最大间隔分离的决策面？

根据前面的优化问题表达，我们可以得到如下表示：

minw,b12(w12+w22)s.t.3w1+3w2+b≥14w1+3w2+b≥1−w1−w2−b≥1\begin{aligned} & \mathop{min}\limits_{w,b} \quad \frac{1}{2}(w_1^2+w_2^2) \\ & s.t. \quad\ \ 3w_1+3w_2+b \geq1 \\ & \qquad \quad\ 4w_1+3w_2+b \geq1 \\ & \qquad \quad -w_1-w_2-b \geq1 \end{aligned} w,bmin21(w12+w22)s.t. 3w1+3w2+b≥1 4w1+3w2+b≥1−w1−w2−b≥1

求解可以得到目标函数的最小时，w1=w2=12,b=−2w_1=w_2=\frac{1}{2},b=-2w1=w2=21,b=−2，于是最大间隔分离的决策面为：

12x(1)+12x(2)−2=0\frac{1}{2}x^{(1)}+\frac{1}{2}x^{(2)}-2=0 21x(1)+21x(2)−2=0

其中，x1=(3,3)Tx_1=(3,3)^Tx1=(3,3)T与x3=(1,1)Tx_3=(1,1)^Tx3=(1,1)T为支持向量。

三、拉格朗日乘数法和对偶问题

首先，我们有必要先了解下什么是拉格朗日乘数法。

关于对偶问题，《统计学习方法》一书中是如此描述的：

为了求解线性可分的支持向量机的最优化问题，将它作为原始最优化问题，应用拉格朗日对偶性，通过求解对偶问题（dual problem）得到原始问题（primary problem）的最优解，这就是线性可分支持向量机的对偶算法（dual algorithm）。这样做的优点，一是对偶问题往往更容易求解；二是自然引入核函数，进而推广到非线性分类问题。

按照前面对优化问题的求解思路，构造拉格朗日方程的目的是**将约束条件放到目标函数中，从而将有约束优化问题转换为无约束优化问题。**我们仍然秉承这一思路去解决不等式约束条件下的优化问题，那么如何针对不等式约束条件下的优化问题构建拉格朗日函数呢？

因为我们要求解的是最小化问题，所以一个直观的想法是如果我能够构造一个函数，使得该函数在可行解区域内与原目标函数完全一致，而在可行解区域外的数值非常大，甚至是无穷大，那么这个没有约束条件的新目标函数的优化问题就与原来有约束条件的原始目标函数的优化是等价的问题。

对此，我们重新回顾下原优化问题的数学表达：

关于拉个朗日乘数法的解题思路，其实早在大学《高等数学》这门课程中就已经提到过，我们不妨通过一个小例子来了解下其解题的一般形式：

例子来源：李亿民-《微积分方法》第二章内容

例子2：求函数f(x,y)=4x+xy2+y2f(x, y)=4x+xy^2+y^2f(x,y)=4x+xy2+y2在D:x2+y2≤1D:x^2+y^2\leq1D:x2+y2≤1上的最值。

做拉格朗日函数L(x,y,λ)L(x,y,\lambda)L(x,y,λ)：

L(x,y,λ)=4x+xy2+y2+λ(x2+y2−1)L(x, y, \lambda)=4x+xy^2+y^2+\lambda(x^2+y^2-1) L(x,y,λ)=4x+xy2+y2+λ(x2+y2−1)

上面的拉格朗日函数，我们分别对x,y,λx,y,\lambdax,y,λ求偏导，得到：

{Lx′(x,y,λ)=4+y2+2λx=0Ly′(x,y,λ)=2xy+2y+2λy=0Lλ′(x,y,λ)=x2+y2−1=0\begin{cases} L_x^{'}(x,y,\lambda)=4+y^2+2\lambda x=0 \\ L_y^{'}(x,y,\lambda)=2xy+2y+2\lambda y=0 \\ L_\lambda^{'}(x,y,\lambda)=x^2+y^2-1=0\\ \end{cases} ⎩⎪⎨⎪⎧Lx′(x,y,λ)=4+y2+2λx=0Ly′(x,y,λ)=2xy+2y+2λy=0Lλ′(x,y,λ)=x2+y2−1=0

求解得到x=1,y=0x=1,y=0x=1,y=0或x=−1,y=0x=-1,y=0x=−1,y=0，所以fmin(x,y)=f(−1,0)=−4，fmax(x,y)=f(1,0)=4f_{min}(x,y)=f(-1,0)=-4，f_{max}(x,y)=f(1,0)=4fmin(x,y)=f(−1,0)=−4，fmax(x,y)=f(1,0)=4

上例就是利用拉格朗日求解最优问题的一般思路，先构造拉格朗日函数，然后分别对参数求偏导，最终求解得到最优解。而我们要想得到最大间隔，同样可以根据该思路进行，只不过式子更加复杂，而我们还需要利用拉格朗日的对偶性来简化优化问题。

在此之前，我们先来回顾下该优化问题下的数学表达：

minw,b12∣∣w∣∣2s.t.yi(wTxi+b)−1≥0，i=1,2,3,...,N(3-1)\begin{aligned} & \mathop{min}\limits_{w,b} \quad \frac{1}{2}||w||^2 \\ & s.t. \quad\ y_i(w^Tx_i + b)-1\geq0，i=1,2,3,...,N \tag{3-1} \end{aligned} w,bmin21∣∣w∣∣2s.t. yi(wTxi+b)−1≥0，i=1,2,3,...,N(3-1)

按照我们例2的思路，我们首先构造出该优化问题的拉格朗日函数，注意：（2-1）式的限制条件是对于样本 iii 的，而这里的i=1,2,3,...,Ni=1,2,3,...,Ni=1,2,3,...,N，所以想这样的约束条件有N个，而每个约束条件都对应一个拉格朗日乘子（例2的λ\lambdaλ），我们不妨将该拉格朗日乘子定义为αi，其中i=1,2,3,...N\alpha_i，其中i=1,2,3,...Nαi，其中i=1,2,3,...N。据此，我们构建出的该优化问题的拉格朗日函数如下：

L(w,b,α)=12∣∣w∣∣2−∑i=1Nαiyi(wTxi+b)+∑i=1Nαi(3-2)L(w,b,\alpha)=\frac{1}{2}||w||^2-\sum^N_{i=1}\alpha_iy_i(w^Tx_i+b) + \sum^N_{i=1}\alpha_i \tag{3-2} L(w,b,α)=21∣∣w∣∣2−i=1∑Nαiyi(wTxi+b)+i=1∑Nαi(3-2)

其中，α=(α1,α2,α3,...,αN)T\alpha=(\alpha_1,\alpha_2,\alpha_3,...,\alpha_N)^Tα=(α1,α2,α3,...,αN)T为拉格朗日乘子向量。

令θ(w,b)=maxα>0L(w,b,α)\theta(w,b)=\mathop{max}\limits_{\alpha>0} \ L(w,b,\alpha)θ(w,b)=α>0max L(w,b,α)，原始目标函数即可转化为：

minw,b12∣∣w∣∣2=minw,bmaxαL(w,b,α)=p∗(3-3)\mathop{min}\limits_{w,b} \frac{1}{2}||w||^2= \mathop{min}\limits_{w,b}\ \mathop{max}\limits_{\alpha}L(w,b,\alpha) =p^* \tag{3-3} w,bmin21∣∣w∣∣2=w,bmin αmaxL(w,b,α)=p∗(3-3)

根据上述目标函数，我们可以发现是先求最大值，再求最小值，而内层最小值是关于α\alphaα，而外层最大值是关于w,bw,bw,b，而我们的α\alphaα又是不等式约束，这样对于求解来讲过于麻烦。由拉格朗日的对偶性，原始问题的对偶问题就是极大极小值，其对偶问题如下：

maxαminw,bL(w,b,α)=d∗(3-4)\mathop{max}\limits_{\alpha}\ \mathop{min}\limits_{w,b}L(w,b,\alpha) =d^* \tag{3-4} αmax w,bminL(w,b,α)=d∗(3-4)

且原问题和对偶问题满足如关系：d∗≤p∗d^*\leq p^*d∗≤p∗，而我们要找的是d∗=p∗d^*= p^*d∗=p∗。读到这里，我们应该会存在两个疑问：

疑问1：为什么前面我们令θ(w,b)=maxα>0L(w,b,α)\theta(w,b)=\mathop{max}\limits_{\alpha>0} \ L(w,b,\alpha)θ(w,b)=α>0max L(w,b,α)，而不是其他的表达形式呢？

主要是因为，这样替换之后，我们能使得该函数在可行解区域内与原目标函数完全一致，而在可行解区域外的数值非常大，甚至是无穷大，那么这个没有约束条件的新目标函数的优化问题就与原来有约束条件的原始目标函数的优化是等价的问题。（这句话要重点理解）

令：

θ(w,b)=maxα>0L(w,b,α)=maxα>0[12∣∣w∣∣2−∑i=1Nαiyi(wTxi+b)+∑i=1Nαi](3-5)\theta(w,b)=\mathop{max}\limits_{\alpha>0} \ L(w,b,\alpha)=\mathop{max}\limits_{\alpha>0} \ [\frac{1}{2}||w||^2-\sum^N_{i=1}\alpha_iy_i(w^Tx_i+b) + \sum^N_{i=1}\alpha_i \tag{3-5}] θ(w,b)=α>0max L(w,b,α)=α>0max [21∣∣w∣∣2−i=1∑Nαiyi(wTxi+b)+i=1∑Nαi](3-5)

之后，θ(w,b)\theta(w,b)θ(w,b)在可行解区域之内yi(wTxi+b)−1≥0\quad\ y_i(w^Tx_i + b)-1\geq0 yi(wTxi+b)−1≥0和可行解区域之外yi(wTxi+b)−1<0\quad\ y_i(w^Tx_i + b)-1<0 yi(wTxi+b)−1<0，我们通过分析得到：在可行解区域之内，

∑i=1Nαiyi(wTxi+b)−∑i=1Nαi≥0\sum^N_{i=1}\alpha_iy_i(w^Tx_i+b) - \sum^N_{i=1}\alpha_i \geq0 i=1∑Nαiyi(wTxi+b)−i=1∑Nαi≥0

此时，我们要求的是 maxα>0L(w,b,α)\mathop{max}\limits_{\alpha>0} \ L(w,b,\alpha)α>0max L(w,b,α)，而12∣∣w∣∣2\frac{1}{2}||w||^221∣∣w∣∣2减去一个大于等于0的最大值是多少？不就是等于12∣∣w∣∣2\frac{1}{2}||w||^221∣∣w∣∣2么，而12∣∣w∣∣2\frac{1}{2}||w||^221∣∣w∣∣2同样也是我们可行解区域之内的目标函数。也就是说在可行解区域之内，12∣∣w∣∣2\frac{1}{2}||w||^221∣∣w∣∣2就等价于maxα>0L(w,b,α)\mathop{max}\limits_{\alpha>0} \ L(w,b,\alpha)α>0max L(w,b,α)。同理，我们可以分析得到，在的可行解区域之外，此时maxα>0L(w,b,α)\mathop{max}\limits_{\alpha>0} \ L(w,b,\alpha)α>0max L(w,b,α)可区域无穷大。所以没有约束的新目标函数的优化问题就与原来有约束条件的原始目标函数的优化是等价的问题，即：

12∣∣w∣∣2=>maxα>0L(w,b,α)(3-6)\frac{1}{2}||w||^2 \quad => \quad \mathop{max}\limits_{\alpha>0} \ L(w,b,\alpha) \tag{3-6} 21∣∣w∣∣2=>α>0max L(w,b,α)(3-6)

疑问2：为什么将原问题转化为对偶问题之后，在什么样的条件下，才会满足d∗=p∗d^*= p^*d∗=p∗？

转化为对偶问题之后，要使得d∗=p∗d^*= p^*d∗=p∗等式成立，则需要满足如下条件，也就是该问题成立时候的KKT条件：

{αi,i=1,2,3,...,Nyi(wTxi+b)≥0,i=1,2,3,...,Nαi(yi(wTxi+b)−1)=0,i=1,2,3,...,N(3-7)\begin{cases} \alpha_i, & i=1,2,3,...,N\\ y_i(w^Tx_i+b)\geq0, & i=1,2,3,...,N\\ \alpha_i(y_i(w^Tx_i+b)-1)=0, & i=1,2,3,...,N\\ \end{cases} \tag{3-7} ⎩⎪⎨⎪⎧αi,yi(wTxi+b)≥0,αi(yi(wTxi+b)−1)=0,i=1,2,3,...,Ni=1,2,3,...,Ni=1,2,3,...,N(3-7)

前两个条件我们不难得到，前面我们也有过对其进行分析，那么第三个条件为什么需要满足呢？

在介绍支持向量和决策平面的时候，我们有提到，最终的决策平面只会与支持向量相关，而与其他大多数样本数据（非支持向量）无关。我们不妨来对它们分别介绍，当样本点为支持向量的时候，此时yi(wTxi+b)−1=0y_i(w^Tx_i+b)-1=0yi(wTxi+b)−1=0，即当所取的样本点为非支持向量的时候，要使得αi(yi(wTxi+b)−1)=0\alpha_i(y_i(w^Tx_i+b)-1)=0αi(yi(wTxi+b)−1)=0，则此时的αi=0\alpha_i=0αi=0。所以从整体样本域来讲，只有满足αi(yi(wTxi+b)−1)=0\alpha_i(y_i(w^Tx_i+b)-1)=0αi(yi(wTxi+b)−1)=0的时候，才能到达目标决策平面与支持向量有关，而与其他大多数非支持向量无关。

综上，只有满足KKT条件的时候，才能到达d∗=p∗d^*= p^*d∗=p∗，即在满足KKT条件的时候，才能将原始目标问题，转化为对偶问题，此时两者是等价的，且此时的目标函数为内层关于w,bw,bw,b的最小值，而外层是关于α\alphaα的最大值，这样一来，大大方便了我们对目标函数的求解。所以优化问题，我们转化如下：

maxαminw,b[12∣∣w∣∣2−∑i=1Nαiyi(wTxi+b)+∑i=1Nαi](3-8)\mathop{max}\limits_{\alpha}\ \mathop{min}\limits_{w,b}[\frac{1}{2}||w||^2-\sum^N_{i=1}\alpha_iy_i(w^Tx_i+b) + \sum^N_{i=1}\alpha_i] \tag{3-8} αmax w,bmin[21∣∣w∣∣2−i=1∑Nαiyi(wTxi+b)+i=1∑Nαi](3-8)

而要解决这个优化问题，我们可以分为两步来进行求解：

先求内层关于w,bw,bw,b的最小值，此时我们应该将α\alphaα视为常数
再求外层关于α\alphaα的最大值，由于经过了第一步关于w,bw,bw,b的最小值，这两个参数已经被消掉了，所以第二步只会存在关于α\alphaα的求解

经过上述对目标函数的问题分析，我们下面根据上述的两个步骤来手握手式的进行求解。

四、模型的求解

关于拉格朗日的内层求解

首先，我们需要对内层的最小值问题进行求解，即求：

minw,bL(w,b)=minw,b[12∣∣w∣∣2−∑i=1Nαiyi(wTxi+b)+∑i=1Nαi](4-1)\mathop{min}\limits_{w,b}L(w,b)=\mathop{min}\limits_{w,b}[\frac{1}{2}||w||^2-\sum^N_{i=1}\alpha_iy_i(w^Tx_i+b) + \sum^N_{i=1}\alpha_i] \tag{4-1} w,bminL(w,b)=w,bmin[21∣∣w∣∣2−i=1∑Nαiyi(wTxi+b)+i=1∑Nαi](4-1)

注意，此时L(w,b)L(w,b)L(w,b)仅仅只是关于w,bw,bw,b，而α\alphaα是由外层最大值问题进行求解的，在这里当做常数处理即可。根据例2，我们需要求出L(w,b)L(w,b)L(w,b)关于w,bw,bw,b的偏导，我们在这假设每个样本含有mmm个属性，则xi=(xi(1),xi(2),...,xi(m))Tx_i=(x_i^{(1)},x_i^{(2)},...,x_i^{(m)})^Txi=(xi(1),xi(2),...,xi(m))T，应各位“老婆们”的要求，具体求偏导的详细过程如下：

为了让读者彻底明白上述过程，所以步骤有点多，这里就不采用Latex语法来编辑上述公式的推导过程了，当然了，Taoye会尽可能地将过程写的足够详细。上述关于拉格朗日函数求偏导的过程，自认为已经写的很详细了，最主要的是要区分x、wx、wx、w的shape问题，以及各自代表的意义是什么。对于上述过程有不清楚的，随时欢迎联系Taoye或在下方留言，也欢迎更多读者来访微信公众号：玩世不恭的Coder。

对此，我们已经通过上面过程得到L(w,b)L(w,b)L(w,b)关于w,bw,bw,b偏导所要满足的式子：

w=∑i=1Nαiyixi0=∑i=1Nαiyi(4-2)w=\sum_{i=1}^N\alpha_iy_ix_i \\ 0=\sum_{i=1}^N\alpha_iy_i \tag{4-2} w=i=1∑Nαiyixi0=i=1∑Nαiyi(4-2)

之后，我们将式子（4-2）重新带回到（4-1）中的最小值问题中，即可消掉w,bw,bw,b参数，也就是说得到的式子仅仅只是关于α\alphaα，具体代回过程如下图所示：

上图就是将对w,bw,bw,b求导之后所得到的式子待会L(w,b)L(w,b)L(w,b)的完整过程，务必要好好理解清楚。

关于这个代回的过程，有两点还是有必要说一下的，这也是前几天实验室的同学存在疑问的地方。（参照上述图片来解疑）

疑问1：这里前一项难道不是和后一项同样为0么？因为不是说了∑i=1Nαiyi=0\sum_{i=1}^N\alpha_iy_i=0∑i=1Nαiyi=0啊？？？

其实这个地方的后一项是为0，而前一项并不一定为0。因为后一项中的bbb其实就相当于一个固定的常数，也就是αiyi\alpha_iy_iαiyi中的每一项所乘的数都是bbb，这样的话，固定的常数乘以0，结果当然依然等于0咯。

而我们再看看前一项，可以发现前一项中除了αiyi\alpha_iy_iαiyi之外，还有wTxiw^Tx_iwTxi的存在，而我们的iii是会随着样本的变化而改变的，所以每次乘的数可能并不一定相同。举个例子理解一下吧：2+3−5=02+3-5=02+3−5=0，这个我们都应该知道。当我们在2、3、52、3、52、3、5的前面同时乘以一个相同的数，这个等式是依然成立的，然而当我在2、3、52、3、52、3、5前面分别乘以一个并不相同的数，那么这个等式就不成立了，比如 2∗2+2∗3−2∗5=02*2+2*3-2*5=02∗2+2∗3−2∗5=0依然成立，而 2∗2+3∗3−4∗5=02*2+3*3-4*5=02∗2+3∗3−4∗5=0就不成立了。

就是这个道理，务必要好好理解清楚。

疑问2：为什么这一步可以推导至下面那一步呢？？？

其实这个问题很好理解，因为这两个成绩形式的式子是相互独立的，也就是说∑\sum∑虽然都有iii，但是这两个iii并不一定是同时取同一个值的。这就像一种笛卡尔积的形式，所以将前一个∑\sum∑式子中的iii换成jjj
依然成立，所以能推导至下面那一步。。。

通过上述过程，我们已经得到了代回之后的式子，如下：

minw,bL(w,b)=∑i=1Nαi−12∑i=1N∑j=1NαiαjyiyjxiTxi(4-3)\mathop{min}\limits_{w,b}L(w,b)=\sum_{i=1}^N\alpha_i-\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_jx_i^Tx_i \tag{4-3} w,bminL(w,b)=i=1∑Nαi−21i=1∑Nj=1∑NαiαjyiyjxiTxi(4-3)

并且我们观察可以发现，式子此时仅仅只存在α\alphaα，而w,bw,bw,b已经成功被消掉了。注意：上式中的xi,yjx_i,y_jxi,yj表示的样本，也就说这些样本的各个属性特征以及标签都是已知的，所以上式只有α\alphaα是未知的。

至此，我们已经解决了对偶问题的内层最小值问题，接下来我们就要求解外层的最大值问题了，将最小值的式子代回原对偶问题，我们更新下对偶问题，得到如下：

maxαminw,bL(w,b,α)=maxαminw,b[12∣∣w∣∣2−∑i=1Nαiyi(wTxi+b)+∑i=1Nαi]=maxα[∑i=1Nαi−12∑i=1N∑j=1NαiαjyiyjxiTxj]=minα[12∑i=1N∑j=1NαiαjyiyjxiTxj−∑i=1Nαi]s.t.αi≥0,i=1,2,...,N∑i=1Nαiyi=0(4-4)\begin{aligned} & \mathop{max}\limits_{\alpha}\ \mathop{min}\limits_{w,b} L(w,b,\alpha) \\ = & \mathop{max}\limits_{\alpha}\ \mathop{min}\limits_{w,b}[\frac{1}{2}||w||^2-\sum^N_{i=1}\alpha_iy_i(w^Tx_i+b) + \sum^N_{i=1}\alpha_i] \tag{4-4} \\ = & \mathop{max}\limits_{\alpha}[\sum_{i=1}^N\alpha_i-\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_jx_i^Tx_j] \\ = & \mathop{min}\limits_{\alpha}[\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_jx_i^Tx_j-\sum_{i=1}^N\alpha_i] \\ & s.t. \quad \alpha_i\geq0, \quad i=1,2,...,N \\ & \qquad \quad \sum_{i=1}^N\alpha_iy_i=0 \end{aligned} ===αmax w,bminL(w,b,α)αmax w,bmin[21∣∣w∣∣2−i=1∑Nαiyi(wTxi+b)+i=1∑Nαi]αmax[i=1∑Nαi−21i=1∑Nj=1∑NαiαjyiyjxiTxj]αmin[21i=1∑Nj=1∑NαiαjyiyjxiTxj−i=1∑Nαi]s.t.αi≥0,i=1,2,...,Ni=1∑Nαiyi=0(4-4)

如上，已经将原对偶转换为了上式的样子，下面我们据此再来看之前的例1

例子来源：李航-《统计学习方法》第七章内容

例子3：已知一个如图所示的训练数据集，其正例点是x1=(3,3)T,x2=(4,3)Tx_1=(3, 3)^T,x_2=(4, 3)^Tx1=(3,3)T,x2=(4,3)T，负例点为x3=(1,1)Tx_3=(1, 1)^Tx3=(1,1)T，试求最大间隔分离的决策面？

根据所给数据，其对偶问题是：

minα[12∑j=1NαiαjyiyjxiTxj−∑i=1Nαi]=12(18α12+25α22+2α32+42α1α2−12α1α3−14α2α3)−α1−α2−α3s.t.α1+α2−α3=0α1≥0,i=1,2,3\begin{aligned} \mathop{min}\limits_{\alpha} & [\frac{1}{2}\sum_{j=1}^N\alpha_i\alpha_jy_iy_jx_i^Tx_j-\sum_{i=1}^N\alpha_i] \\ & = \frac{1}{2}(18\alpha_1^2+25\alpha_2^2+2\alpha_3^2+42\alpha_1\alpha_2-12\alpha_1\alpha_3-14\alpha_2\alpha_3)-\alpha_1-\alpha_2-\alpha_3 \\ s.t. & \alpha_1+ \alpha_2-\alpha_3=0 \\ & \alpha_1 \geq0, \quad i=1,2,3 \end{aligned} αmins.t.[21j=1∑NαiαjyiyjxiTxj−i=1∑Nαi]=21(18α12+25α22+2α32+42α1α2−12α1α3−14α2α3)−α1−α2−α3α1+α2−α3=0α1≥0,i=1,2,3

我们将α3=α1+α2\alpha_3=\alpha_1+\alpha_2α3=α1+α2代入到目标函数并记为：

s(α1,α2)=4α12+132α22+10α1α2−2α1−2α2s(\alpha_1,\alpha_2)=4\alpha_1^2+\frac{13}{2}\alpha_2^2+10\alpha_1\alpha_2-2\alpha_1-2\alpha_2 s(α1,α2)=4α12+213α22+10α1α2−2α1−2α2
对α1,α2\alpha_1,\alpha_2α1,α2求偏导数并令其为0，易知s(α1,α2)s(\alpha_1,\alpha_2)s(α1,α2)在点(32,−1)T(\frac{3}{2},-1)^T(23,−1)T取极值，但该点不满足约束条件α2≥0\alpha_2\geq0α2≥0，所以最小值应在边界上达到。

当α1=0\alpha_1=0α1=0时，最小值为s(0,213)=−213s(0,\frac{2}{13})=-\frac{2}{13}s(0,132)=−132；当α2=0\alpha_2=0α2=0时，最小值s(14,0)=−14s(\frac{1}{4},0)=-\frac{1}{4}s(41,0)=−41。所以，当α1=14,α2=0,α3=α1+α2=14\alpha_1=\frac{1}{4},\alpha_2=0,\alpha_3=\alpha_1+\alpha_2=\frac{1}{4}α1=41,α2=0,α3=α1+α2=41时，此时s(α1,α2)s(\alpha_1,\alpha_2)s(α1,α2)达到最小。

这样的话，就说明α1=α3=14\alpha_1=\alpha_3=\frac{1}{4}α1=α3=41所对应的点x1,x3x_1,x_3x1,x3就是支持向量了，根据w=∑i=1Nαiyixiw=\sum_{i=1}^N\alpha_iy_ix_iw=∑i=1Nαiyixi，我们可以求得此时w1=w2=12w_1=w_2=\frac{1}{2}w1=w2=21，再由b=y−∑i=1NαiyixiTxb=y-\sum_{i=1}^N\alpha_iy_ix_i^Txb=y−∑i=1NαiyixiTx，可以得到b=−2b=-2b=−2

综上，我们可以得到决策面最终为：

12x(1)+12x(2)−2=0\frac{1}{2}x^{(1)}+\frac{1}{2}x^{(2)}-2=0 21x(1)+21x(2)−2=0

其中，x1=(3,3)Tx_1=(3,3)^Tx1=(3,3)T与x3=(1,1)Tx_3=(1,1)^Tx3=(1,1)T为支持向量。

历经九九八十一难，终于来到了这最后一步了，只要我们能求得上式的最小值，那么模型的求解也就该到一段落了。

那么，问题来了，关于上式，我们应当如何求解呢？？？

下面就应该是我们的重头戏闪亮登场了——SMO算法，各位看官掌声欢迎一下，有请SMO大佬登台表演。。。

基于SMO算法的外层求解

关于SMO算法，李航老师的《统计学习方法》一书中是这么描述的，

关于外层问题的求解，我们有许多方法可供使用。当我们的数据样本比较少的时候，可以将支持向量机的学习问题转换为凸二次规划问题，这样的凸二次规划问题具有全局最优解，并且有许多最优化算法可以用于这一问题的求解。但是当训练样本容量很大时，这些算法往往变得非常低效，以致无法使用。所以，如何高效地实现支持向量机学习就成为一个重要的问题。目前人们已提出许多快速实现算法。而在这些算法中，要数Platt的序列最小最优化（sequential minimal optimization SMO）算法最为人所知。

**SMO算法是一种启发式算法，其基本思想是：**如果所有变量的解都满足此最优化问题的KKT条件，那么这个最优化问题的解就得到了。因为KKT条件是该最优化问题的充分必要条件。否则，选择两个变量，固定其他变量，针对这两个变量构建一个二次规划问题。这个二次规划问题关于这两个变量的解就应该更接近原始二次规划问题的解，因为这会使得原始二次规划问题的目标函数值变得更小。更重要的是，这时子问题可以通过解析方法求解，这样就可以大大提高整个算法的计算速度。子问题有两个变量，一个是违反KKT条件最严重的那一个，另一个由约束条件自动确定。如此，SMO算法将原问题不断分解为子问题并对子问题求解，进而达到求解原问题的目的。

上述内容来自李航——《统计学习方法》第二版。

不知道大伙读完上述关于内容是什么感受，这里简单总结一下李航老师所表达的意思吧。

在Taoye的印象里，小学时期上语文课的时候学习过一篇文章叫做**《走一步，再走一步》**。（具体几年级就记不清楚了）

嘿！您还别说，刚刚去搜索了下这篇课文，还真就叫这个名儿。第一次读李航老师《统计学习方法》中关于SMO的内容之后，就让我想起这篇文章。我还专门重新读了一下这篇文章，主要讲的内容是这样的：

文章中主人公名叫小亨特，他不是天生体弱怯懦嘛，在一次和小伙伴攀登悬崖的时候，由于内心的恐惧、害怕，在攀登途中上不去，也下不来。然后呢，他的小伙伴杰里就把小亨特的父亲找来了，父亲对小亨特说：“不要想有多远，有多困难，你需要想的是迈一小步。这个你能做到。看着手电光指的地方。看到那块石头没有？”，最终通过父亲的鼓励，小亨特成功脱险。文末作者还总结道：在我生命中有很多时刻，每当我遇到一个遥不可及、令人害怕的情境，并感到惊慌失措时，我都能够应付——因为我回想起了很久以前自己上过的那一课。我提醒自己不要看下面遥远的岩石，而是注意相对轻松、容易的第一小步，迈出一小步、再一小步，就这样体会每一步带来的成就感，直到完成了自己想要完成的，达到了自己的目标，然后再回头看时，不禁对自己走过的这段漫漫长路感到惊讶和自豪。

把《走一步，再走一步》这篇文章搬出来，真的不是在凑字数从而给大家阅读带来压力，只是觉得SMO算法描述的就是这么个理儿。算了，不多说了，说多了还真的会有凑字数的嫌疑。(ノへ￣、)

下面我们开始进入到SMO吧。。。

在这之前，我们把外层的最小值问题再搬出来：

minα[12∑i=1N∑j=1NαiαjyiyjxiTxj−∑i=1Nαi]s.t.αi≥0,i=1,2,...,N∑i=1Nαiyi=0(4-5)\begin{aligned} & \mathop{min}\limits_{\alpha}[\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_jx_i^Tx_j-\sum_{i=1}^N\alpha_i] \\ & s.t. \quad \alpha_i\geq0, \quad i=1,2,...,N \\ & \qquad \quad \sum_{i=1}^N\alpha_iy_i=0 \end{aligned} \tag{4-5} αmin[21i=1∑Nj=1∑NαiαjyiyjxiTxj−i=1∑Nαi]s.t.αi≥0,i=1,2,...,Ni=1∑Nαiyi=0(4-5)

在这里，我们是假设对于所有的样本数据都是100%线性可分的。

对于该优化问题的SMO算法，我们可以这样理解：因为在我们的数据集中，xix_ixi属于每个样本的属性特征，yiy_iyi为样本所对应的标签，而这些都是已知的，上述优化问题的目标函数只存在αi\alpha_iαi为未知变量，且未知变量有NNN个。而根据SMO算法的思想，我们每次只将其中两个α\alphaα看做变量，而其他α\alphaα仅仅只是常数。在却确定其中一个变量的时候，另一个变量就可以通过∑i=1Nαiyi=0\sum_{i=1}^N\alpha_iy_i=0∑i=1Nαiyi=0约束得到。我们不妨将该两个变量定为α1,α2\alpha_1,\alpha_2α1,α2，则SMO不断执行如下两个步骤直至收敛：

选取一对需要更新的变量α1和α2\alpha_1和\alpha_2α1和α2
固定α1和α2\alpha_1和\alpha_2α1和α2以外的参数，根据求解式（4-5）获得更新后α1和α2\alpha_1和\alpha_2α1和α2
更新好α1和α2\alpha_1和\alpha_2α1和α2之后，重新选取两个α\alphaα进行不断更新迭代（重复1、2步骤）

讲到这里，SMO算法是不是和《走一步，再走一步》中主人公类似呢？将一个大的、复杂的问题转换成多个小问题，然后不断的迭代更新。

为什么我们每次都同时优化两个参数，而不是一个呢？因为每次更新两个参数，才能确保∑i=1Nαiyi=0\sum_{i=1}^N\alpha_iy_i=0∑i=1Nαiyi=0约束条件成立。而当我们仅仅只是修改一个α\alphaα时，那么就将违背这个约束条件了。

据SMO的思想，我们不妨把目标函数中的α1和α2\alpha_1和\alpha_2α1和α2单独拎出来，如下：

注意：因为我们的W(α1,α2)W(\alpha_1,\alpha_2)W(α1,α2)仅仅只是对α1，α2\alpha_1，\alpha_2α1，α2作为参数，而α3,α4,...,αN\alpha_3,\alpha_4,...,\alpha_Nα3,α4,...,αN只是作为一个参数而存在。还有一点需要注意的是，因为我们是二分类问题，且样本数据的标签为非1即-1，所以**yi2=1y_i^2=1yi2=1，这个在化简过程中需要用到。**此时我们得到关于α1,α2\alpha_1,\alpha_2α1,α2的目标函数为：

W(α1,α2)=12K11α12+12K22α22+α1α2y1y2K12−(α1+α2)+y1α1∑i=3NyiαiK1i+y2α2∑i=3NyiαiK2i+H其中Kij=xiTxj,H=∑i=3N∑j=3NαiαjyiyjxiTxjT−∑i=3Nαi(4-6)\begin{aligned} W(\alpha_1, \alpha_2)= & \frac{1}{2}K_{11}\alpha_1^2+\frac{1}{2}K_{22}\alpha_2^2+\alpha_1\alpha_2y_1y_2K_{12} -(\alpha_1+\alpha_2) \\ & + y_1\alpha_1\sum_{i=3}^Ny_i\alpha_iK_{1i}+y_2\alpha_2\sum_{i=3}^Ny_i\alpha_iK_{2i} + H \\ & 其中K_{ij}=x_i^Tx_j,H=\sum_{i=3}^N\sum_{j=3}^N\alpha_i\alpha_jy_iy_jx_i^Tx_j^T-\sum_{i=3}^N\alpha_i \end{aligned} \tag{4-6} W(α1,α2)=21K11α12+21K22α22+α1α2y1y2K12−(α1+α2)+y1α1i=3∑NyiαiK1i+y2α2i=3∑NyiαiK2i+H其中Kij=xiTxj,H=i=3∑Nj=3∑NαiαjyiyjxiTxjT−i=3∑Nαi(4-6)

我们知道对于这个式子是有一个约束条件的，我们可以根据这个用α2\alpha_2α2来表示α1\alpha_1α1（注意：yi2=1y_i^2=1yi2=1），如下：

∑i=1Nαiyi=0=>α1y1+α2y2=−∑i=3Nαiyi=G=>α1=y1(G−α2y2)(4-7)\begin{aligned} & \qquad \sum_{i=1}^N\alpha_iy_i=0 \\ & => \alpha_1y_1+\alpha_2y_2=-\sum_{i=3}^N\alpha_iy_i=G \\ & => \alpha_1=y_1(G-\alpha_2y_2) \end{aligned} \tag{4-7} i=1∑Nαiyi=0=>α1y1+α2y2=−i=3∑Nαiyi=G=>α1=y1(G−α2y2)(4-7)

通过上式，用α2\alpha_2α2来表示α1\alpha_1α1，我们不妨将α1\alpha_1α1带到W(α1,α2)W(\alpha_1,\alpha_2)W(α1,α2)中，得到一个只关于α2\alpha_2α2的式子：

W(α2)=12K11(G−α2y2)2+12K22α22+α2(G−α2y2)y2K12−[y1(G−α2y2)+α2]+(G−α2y2)∑i=3NyiαiK1i+y2α2∑i=3NyiαiK2i+H(4-8)\begin{aligned} W(\alpha_2)= & \frac{1}{2}K_{11}(G-\alpha_2y_2)^2+\frac{1}{2}K_{22}\alpha_2^2+\alpha_2(G-\alpha_2y_2)y_2K_{12}\\ & -[y_1(G-\alpha_2y_2) + \alpha_2]+(G-\alpha_2y_2)\sum_{i=3}^Ny_i\alpha_iK_{1i} \\ & + y_2\alpha_2\sum_{i=3}^Ny_i\alpha_iK_{2i}+H \end{aligned} \tag{4-8} W(α2)=21K11(G−α2y2)2+21K22α22+α2(G−α2y2)y2K12−[y1(G−α2y2)+α2]+(G−α2y2)i=3∑NyiαiK1i+y2α2i=3∑NyiαiK2i+H(4-8)

此时的W(α2)W(\alpha_2)W(α2)仅仅只是关于α2\alpha_2α2的函数，我们将W(α2)W(\alpha_2)W(α2)对α2\alpha_2α2进行求导，并令导数为0：

上述就是SMO中限制其中两个变量的推到过程的推到过程（公式太多，过程有点复杂，确实不方便使用Latex语法，不过过程都已经写的很详细了，还是需要静下心来慢慢手动推导的）下面总结一下上述SMO算法的过程吧：

前面我们不是得到了仅仅关于α2\alpha_2α2为变量的W(α2)W(\alpha_2)W(α2)么，也就是说此时的未知数只有一个，我们要求W(α2)W(\alpha_2)W(α2)的最值应该怎么求呢？当然是对其进行求导咯，然后对导数为0，即可解出W(α2)W(\alpha_2)W(α2)取最值时候的α2\alpha_2α2，整理之后得到如下式子：

(K11+K22−2K12)α2=y2(y2−y1+GK11−GK12+∑i=3NαiyiK1i−∑i=3NαiyiK2i)(4-9)(K_{11}+K_{22}-2K_{12})\alpha_2=y_2(y_2-y_1+GK_{11}-GK_{12}+\sum_{i=3}^N\alpha_iy_iK_{1i}-\sum_{i=3}^N\alpha_iy_iK_{2i}) \tag{4-9} (K11+K22−2K12)α2=y2(y2−y1+GK11−GK12+i=3∑NαiyiK1i−i=3∑NαiyiK2i)(4-9)

此时，我们可以发现除了数据样本相关信息和α2\alpha_2α2之外，还有GGG的存在，而我们前面也又说到，SMO算法本身是一个不断迭代更新的过程，我们需要的是可以通过的更新之前的α2\alpha_2α2来修改α2\alpha_2α2，从而得到一个新的α2\alpha_2α2，我们不妨令新的α\alphaα为αnew\alpha^{new}αnew、旧的α\alphaα为αold\alpha^{old}αold。而我们知道旧的α\alphaα之间需要满足一个限制条件：

α1oldy1+α2oldy2=G(4-9)\alpha_1^{old}y_1 + \alpha_2^{old}y_2=G \tag{4-9} α1oldy1+α2oldy2=G(4-9)

所以，我们重新将GGG代回到(4−9)(4-9)(4−9)式，用α1old、α2old\alpha_1^{old}、\alpha_2^{old}α1old、α2old来代替GGG，(要时刻注意：yi2=1y_i^2=1yi2=1)得到：

(K11+K22−2K12)α2new=y2(y2−y1+α1oldy1K11+α2oldy2K11−α1oldy1K12−α2oldy2K12+∑i=3NαiyiK1i−∑i=3NαiyiK2i)(4-10)\begin{aligned} (K_{11}+K_{22}-2K_{12})\alpha_2^{new}= & y_2(y_2-y_1+\alpha_1^{old}y_1K_{11}+\alpha_2^{old}y_2K_{11} \\ & -\alpha_1^{old}y_1K_{12}-\alpha_2^{old}y_2K_{12} \\ & + \sum_{i=3}^N\alpha_iy_iK_{1i}-\sum_{i=3}^N\alpha_iy_iK_{2i}) \end{aligned} \tag{4-10} (K11+K22−2K12)α2new=y2(y2−y1+α1oldy1K11+α2oldy2K11−α1oldy1K12−α2oldy2K12+i=3∑NαiyiK1i−i=3∑NαiyiK2i)(4-10)

之后我们通过前面拉格朗日得到的关系式，用f(x1)、f(x2)f(x_1)、f(x_2)f(x1)、f(x2)来代替（4-10）后面的两个级数，整理最终得到：
α2new=y2(E1−E2)η+α2oldα2new=y1(G−α2newy2)(4-11)\begin{aligned} & \alpha_2^{new}=\frac{y_2(E_1-E_2)}{\eta}+\alpha_2^{old} \\ & \alpha_2^{new}=y_1(G-\alpha_2^{new}y_2) \end{aligned} \tag{4-11} α2new=ηy2(E1−E2)+α2oldα2new=y1(G−α2newy2)(4-11)

PS：关于η、E1、E2\eta、E_1、E_2η、E1、E2是什么，请见上图中的内容。

通过如上式子，我们就能求得更新之后的α1、α2\alpha_1、\alpha_2α1、α2，而SMO算法的核心在于两两不断的更新迭代，所以最终我们会得到，每个样本xi、yix_i、y_ixi、yi都会对应一个αi\alpha_iαi，而前面我们也有说过，决策面最终只会与支持向量有关，而与非支持向量的样本无关，所以大多数的αi\alpha_iαi都等于0，只有少数αi\alpha_iαi为非0。如此一来，我们就能求解得到α\alphaα向量：α=(α1,α2,...,αN)\alpha=(\alpha_1,\alpha_2,...,\alpha_N)α=(α1,α2,...,αN)，随后，我们就能通过αi、xi、yi\alpha_i、x_i、y_iαi、xi、yi就能求得参数www：

w=∑i=1Nαiyixiw=\sum_{i=1}^N\alpha_iy_ix_i w=i=1∑Nαiyixi

还有一点需要注意的是，上述过程都是默认所有样本数据都是线性可分的，也就是说没有一个样本会被误分类。但这只是理想状态下，而实际不免会有个别数据不得不被误分类，这时我们需要定义惩罚参数和容错率，而容错率是用来不断优化α\alphaα的，主要通过实际值与真实值得到。而惩罚参数我们定为CCC，而要想成功更新得到α2new\alpha_2^{new}α2new，则需要确定α2new\alpha_2^{new}α2new的范围。对此，我们不妨定义如下：

0≤αi≤C,i=1,2,3,...,NL≤α2new≤R(4-12)0\leq\alpha_i\leq C ,\qquad i=1,2,3,...,N\\ L\leq\alpha_2^{new}\leq R \tag{4-12} 0≤αi≤C,i=1,2,3,...,NL≤α2new≤R(4-12)

而我们知道：

α1newy1+α2newy2=α1oldy1+α2oldy2=G(4-13)\alpha_1^{new}y_1+\alpha_2^{new}y_2=\alpha_1^{old}y_1+\alpha_2^{old}y_2=G \tag{4-13} α1newy1+α2newy2=α1oldy1+α2oldy2=G(4-13)

综合上式，可以确定L、RL、RL、R的范围：

L=max(0,α2old−α1old),R=min(C,C+α2old−α1old)y1≠y2L=max(0,α2old+α1old−C),R=min(C,α2old+α1old)y1=y2(4-14)L=max(0, \alpha_2^{old}-\alpha_1^{old}),R=min(C,C+\alpha_2^{old}-\alpha_1^{old}) \quad y_1 \neq y_2\\ L=max(0, \alpha_2^{old}+\alpha_1^{old}-C),R=min(C,\alpha_2^{old}+\alpha_1^{old})\quad y_1=y_2 \tag{4-14} L=max(0,α2old−α1old),R=min(C,C+α2old−α1old)y1=y2L=max(0,α2old+α1old−C),R=min(C,α2old+α1old)y1=y2(4-14)

而这个在不同情况下的L,RL,RL,R的范围，我们会在下面实际编程的时候需要用到，主要是用来更新α\alphaα值。

接下来，就是更新bbb值了。前面我们已经定义了惩罚参数CCC，且0≤αinew≤C0\leq\alpha_i^{new} \leq C0≤αinew≤C，此时通过E1、E2E_1、E_2E1、E2更新得到的b1、b2b_1、b_2b1、b2固然是相等的；但假如我们更新之后的α\alphaα不在这个区间之内，则此时得到的b1、b2b_1、b_2b1、b2是不相等的，所以我们需要确定在不同情况下更新之后的bbb值：

前面，我们已经得到：

y1=∑i=1Nαiyixix1+by_1=\sum_{i=1}^N\alpha_iy_ix_ix_1+b y1=i=1∑Nαiyixix1+b

因为我们是打算通过α1new、α2new\alpha_1^{new}、\alpha_2^{new}α1new、α2new的值来得到更新之后的bnewb^{new}bnew，所以把α1new、α2new\alpha_1^{new}、\alpha_2^{new}α1new、α2new单独拎出来，得到：

b1new=y1−∑i=3NαiyixiTx1−α1newy1x1Tx1−α2newy2x2Tx1(4-15)b_1^{new}=y_1-\sum_{i=3}^N\alpha_iy_ix_i^Tx_1-\alpha_1^{new}y_1x_1^Tx_1-\alpha_2^{new}y_2x_2^Tx_1 \tag{4-15} b1new=y1−i=3∑NαiyixiTx1−α1newy1x1Tx1−α2newy2x2Tx1(4-15)

同时，因为上述中的级数形式可以使用E1和boldE_1和b^{old}E1和bold来表示，所以整理之后得：

b1new=bold−E1−y1(α1new−α1old)x1Tx1−y2(α2new−α2old)x2Tx1(4-16)b_1^{new}=b^{old}-E_1-y_1(\alpha_1^{new}-\alpha_1^{old})x_1^Tx_1-y_2(\alpha_2^{new}-\alpha_2^{old})x_2^Tx_1 \tag{4-16} b1new=bold−E1−y1(α1new−α1old)x1Tx1−y2(α2new−α2old)x2Tx1(4-16)

同理，可以得到b2newb_2^{new}b2new：

b2new=bold−E2−y1(α1new−α1old)x1Tx1−y2(α2new−α2old)x2Tx2(4-17)b_2^{new}=b^{old}-E_2-y_1(\alpha_1^{new}-\alpha_1^{old})x_1^Tx_1-y_2(\alpha_2^{new}-\alpha_2^{old})x_2^Tx_2 \tag{4-17} b2new=bold−E2−y1(α1new−α1old)x1Tx1−y2(α2new−α2old)x2Tx2(4-17)

当更新之后的α1new、α2new\alpha_1^{new}、\alpha_2^{new}α1new、α2new都有效的时候，即在[0−C][0-C][0−C]区间之内时，此时bnew=b1new=b2newb^{new}=b_1^{new}=b_2^{new}bnew=b1new=b2new，而在不满足上述条件的时候，我们更新之后的bnewb^{new}bnew取b1和b2b_1和b_2b1和b2的均值，即：

bnew={b1new,0≤α1new≤Cb2new,0≤α1new≤Cb1new+b2new2,其他(4-18)b^{new} = \begin{cases} b_1^{new}, & 0\leq\alpha_1^{new}\leq C \\ b_2^{new}, & 0\leq\alpha_1^{new}\leq C \\ \frac{b_1^{new}+b_2^{new}}{2}, & 其他 \\ \end{cases} \tag{4-18} bnew=⎩⎪⎨⎪⎧b1new,b2new,2b1new+b2new,0≤α1new≤C0≤α1new≤C其他(4-18)

如此一来，我们就已经完成了SMO算法的流程，该有的参数都已经求解出来了。

说实话，写到这，Taoye的确有点累了，脑细胞也严重不足了，但为了各位“老婆们”的正常阅读，还是得继续写下去才行。

下面，我们就通过编程来实现线性SVM算法吧！（本次手撕SVM的数据集依然采用我们前面所随机创建的）

五、编程手撕线性SVM

在前面，我们其实已经实现了线性SVM的分类，只不过那个时候使用的是sklearn内置的接口。但既然是手撕SVM，当然是需要自己来实现这一功能的。

在这里需要提前说明的是，上述代码大量使用到了NumPy操作，关于NumPy的使用，可自行参考之前写的一篇文章：print( “Hello，NumPy！” )

训练SVM模型，没数据集可不行，本次手撕SVM的数据集依然采用我们前面所随机创建，对此，我们定义一个etablish_data方法来随机创建一个SVM二分类数据集：

"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 用于生成训练数据集Parameters:data_number: 样本数据数目Return:x_data: 数据样本的属性矩阵y_label: 样本属性所对应的标签
"""
def etablish_data(data_number):x_data = np.concatenate((np.add(np.random.randn(data_number, 2), [3, 3]),       np.subtract(np.random.randn(data_number, 2), [3, 3])),axis = 0)      # random随机生成数据，+ -3达到不同类别数据分隔的目的 temp_data = np.zeros([data_number])temp_data.fill(-1)y_label = np.concatenate((temp_data, np.ones([data_number])), axis = 0)return x_data, y_label

前面，我们在讲解SMO算法的时候提到，每次都会选取随机**两个α\alphaα**来进行更新，这里我们不妨将第一个αi\alpha_iαi通过遍历的形式逐个选取，而另一个则通过np.random模块来随机选取，这里需要主要的是，第二个选取的αj\alpha_jαj不能与第一个相同。为此，我们定义一个方法random_select_alpha_j来随机选取第二个αj\alpha_jαj（第一个已经通过遍历得到）:

"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 随机选取alpha_jParameters:alpha_i_index: 第一个alpha的索引alpha_number: alpha总数目Return:alpha_j_index: 第二个alpha的索引
"""
def random_select_alpha_j(alpha_i_index, alpha_number):alpha_j_index = alpha_i_indexwhile alpha_j_index == alpha_i_index:alpha_j_index = np.random.randint(0, alpha_number)return alpha_j_index

我们知道，每一个更新之后的αjnew\alpha_j^{new}αjnew都需要满足L≤αjnew≤RL\leq\alpha_j^{new}\leq RL≤αjnew≤R。为此，我们定义一个方法modify_alpha来确定α\alphaα在区间之内：

"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 使得alpha_j在[L, R]区间之内Parameters:alpha_j: 原始alpha_jL: 左边界值R: 右边界值Return:L,R,alpha_j: 修改之后的alpha_j
"""
def modify_alpha(alpha_j, L, R):if alpha_j < L: return Lif alpha_j > R: return Rreturn alpha_j

我们模型训练是一个迭代更新的过程，而更新的前提是误差比较大，所以我们需要定义一个方法calc_E_i来计算误差，但误差又怎么计算呢？这一点其实我们在最开始就已经提到过了，误差是通过模型计算出来的值与其真实值最差得到，也就是前面提到的EjE_jEj（下面的推导务必要理解清楚，矩阵变换要十分熟悉）：

f(xj)=wTxj+b=∑i=1NαiyixiTxj+b=α1y1x1Txj+α2y2x2Txj+...+αNyNxNTxj+b=(α1y1,α2y2,...,αNyN)(x1Txjx2Txj⋮xNTxj)+b=[(α1α2⋮αN)⋅(y1y2⋮yN)]T∗xxjT+b(5-1)\begin{aligned} f(x_j) & =w^Tx_j+b \\ & = \sum_{i=1}^N\alpha_iy_ix_i^Tx_j+b \\ & = \alpha_1y_1x_1^Tx_j+\alpha_2y_2x_2^Tx_j+...+\alpha_Ny_Nx_N^Tx_j+b \\ & = (\alpha_1y_1,\alpha_2y_2,...,\alpha_Ny_N)\left( \begin{matrix} x_1^Tx_j\\ x_2^Tx_j \\ \vdots\\ x_N^Tx_j \\ \end{matrix} \right)+b \\ & = [\left( \begin{matrix} \alpha_1\\ \alpha_2 \\ \vdots\\ \alpha_N \\ \end{matrix} \right) \cdot \left( \begin{matrix} y_1\\ y_2 \\ \vdots\\ y_N \\ \end{matrix} \right)]^T*xx_j^T+b \end{aligned} \tag{5-1} f(xj)=wTxj+b=i=1∑NαiyixiTxj+b=α1y1x1Txj+α2y2x2Txj+...+αNyNxNTxj+b=(α1y1,α2y2,...,αNyN)⎝⎜⎜⎜⎛x1Txjx2Txj⋮xNTxj⎠⎟⎟⎟⎞+b=[⎝⎜⎜⎜⎛α1α2⋮αN⎠⎟⎟⎟⎞⋅⎝⎜⎜⎜⎛y1y2⋮yN⎠⎟⎟⎟⎞]T∗xxjT+b(5-1)

Ei=f(xi)−yi(5-2)E_i=f(x_i)-y_i \tag{5-2} Ei=f(xi)−yi(5-2)

根据上述误差的推导，我们现在就可以通过代码来计算误差了（上面的推导务必要理解清楚，矩阵变换要十分熟悉，才能理解下面代码所表达的含义）：

"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 计算误差并返回
"""
def calc_E(alphas, y_lable, x_data, b, i):f_x_i = float(np.multiply(alphas, y_lable).T * (x_data * x_data[i, :].T)) + breturn f_x_i - float(y_label[i])

同样的，我们把其他一些用于整体代换的单独拎出来，并通过方法进行返回，除了上述的误差EiE_iEi之外，还有η和b\eta和bη和b，相关推导过程读者可自行根据上述EiE_iEi来进行（一定要会），相关公式和代码如下：

η=Kii+Kjj−2Kijbinew=bold−Ei−yi(αinew−αiold)xiTxi−yj(αjnew−αjold)xjTxibjnew=bold−Ej−yi(αinew−αiold)xiTxi−yj(αjnew−αjold)xjTxj(5-3)\begin{aligned} & \eta = K_{ii}+K_{jj}-2K_{ij} \\ & b_i^{new}=b^{old}-E_i-y_i(\alpha_i^{new}-\alpha_i^{old})x_i^Tx_i-y_j(\alpha_j^{new}-\alpha_j^{old})x_j^Tx_i\\ & b_j^{new}=b^{old}-E_j-y_i(\alpha_i^{new}-\alpha_i^{old})x_i^Tx_i-y_j(\alpha_j^{new}-\alpha_j^{old})x_j^Tx_j \tag{5-3} \end{aligned} η=Kii+Kjj−2Kijbinew=bold−Ei−yi(αinew−αiold)xiTxi−yj(αjnew−αjold)xjTxibjnew=bold−Ej−yi(αinew−αiold)xiTxi−yj(αjnew−αjold)xjTxj(5-3)

"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 计算eta并返回
"""
def calc_eta(x_data, i, j):eta = 2.0 * x_data[i, :] * x_data[j, :].T \- x_data[i, :] * x_data[i, :].T \- x_data[j, :] * x_data[j,:].Treturn eta"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 计算b1, b2并返回
"""
def calc_b(b, x_data, y_label, alphas, alpha_i_old, alpha_j_old, E_i, E_j, i, j):b1 = b - E_i \- y_label[i] * (alphas[i] - alpha_i_old) * x_data[i, :] * x_data[i, :].T \- y_label[j] * (alphas[j] - alpha_j_old) * x_data[i, :] * x_data[j, :].Tb2 = b - E_j \- y_label[i] * (alphas[i] - alpha_i_old) * x_data[i, :] * x_data[j, :].T \- y_label[j] * (alphas[j] - alpha_j_old) * x_data[j, :] * x_data[j, :].Treturn b1, b2

OK，准备工作已经完成了，接下来是时候放出我们的核心SMO算法的代码了，大家可根据前面的SMO思想来理解，下面代码也会放出详细的注释来帮助大家理解：

"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: SMO核心算法，求得b和laphas向量Parameters:x_data: 样本属性特征矩阵y_label: 属性特征对应的标签C：惩罚参数toler：容错率max_iter：迭代次数Return:b: 决策面的参数balphas：获取决策面参数w所需要的alphas
"""
def linear_smo(x_data, y_label, C, toler, max_iter):x_data = np.mat(x_data); y_label = np.mat(y_label).T     # 将数据转换为矩阵类型m, n = x_data.shape                                      # 得到数据样本的shape信息b, alphas, iter_num = 0, np.mat(np.zeros((m, 1))), 0     # 初始化参数b和alphas和迭代次数while (iter_num < max_iter):                            # 最多迭代max_iter次alpha_optimization_num = 0   # 定义优化次数for i in range(m):          # 遍历每个样本，一次选取一个样本计算误差E_i = calc_E(alphas, y_label, x_data, b, i)      # 样本i的误差计算if ((y_label[i] * E_i < -toler) and (alphas[i] < C)) or ((y_label[i] * E_i > toler) and (alphas[i] > 0)):j = random_select_alpha_j(i, m)              # 随机选取一个不与i重复jE_j = calc_E(alphas, y_label, x_data, b, j)  # 计算样本j的误差alpha_i_old = alphas[i].copy(); alpha_j_old = alphas[j].copy();if (y_label[i] != y_label[j]):               # 重新规范alphas的左右区间L, R = max(0, alphas[j] - alphas[i]), min(C, C + alphas[j] - alphas[i])else:L, R = max(0, alphas[j] + alphas[i] - C), min(C, alphas[j] + alphas[i])if L == R: print("L==R"); continue          # L==R时选取下一个样本eta = calc_eta(x_data, i, j)                  # 计算eta值if eta >= 0: print("eta>=0"); continuealphas[j] -= y_label[j] * (E_i - E_j) / etaalphas[j] = modify_alpha(alphas[j], L, R)     # 修改alpha[j]if (abs(alphas[j] - alpha_j_old) < 0.00001): print("alpha_j更改太小"); continuealphas[i] += y_label[j] * y_label[i] * (alpha_j_old - alphas[j])    # 修改alphas[i]b1, b2= calc_b(b, x_data, y_label, alphas, alpha_i_old, alpha_j_old, E_i, E_j, i, j)    # 计算b值if (0 < alphas[i]) and (C > alphas[i]): b = b1elif (0 < alphas[j]) and (C > alphas[j]): b = b2else: b = (b1 + b2)/2.0alpha_optimization_num += 1print("迭代次数:%d，样本:%d，alphas向量的优化次数:%d" % (iter_num, i+1, alpha_optimization_num))if (alpha_optimization_num == 0): iter_num += 1else: iter_num = 0print("迭代次数:%d" % iter_num)return b, alphas

上述SMO核心方法，我们可以通过定义输入样本的属性特征、标签以及迭代次数等来得到α\alphaα和bbb。随后，我们可以通过www与α\alphaα之间的关系来的计算出www，关系和相关代码如下所示：

w=∑i=1Nαiyixi(5-4)w=\sum_{i=1}^N\alpha_iy_ix_i \tag{5-4} w=i=1∑Nαiyixi(5-4)

"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 根据公式计算出w权值向量Parameters:x_data: 样本属性特征矩阵y_label: 属性特征对应的标签alphas：linear_smo方法所返回的alphas向量Return:w: 决策面的参数w
"""
def calc_w(x_data, y_label, alphas):x_data, y_label, alphas = np.array(x_data), np.array(y_label), np.array(alphas)return np.dot((np.tile(y_label.reshape(1, -1).T, (1, 2)) * x_data).T, alphas).tolist()

好的，bbb有了，www也有了，该有的都有了，接下来就是验证模型效果了，这里我们使用Matplotlib来绘制，定义一个plot_result方法来展示模型分类结果：

import pylab as pl"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 绘制出分类结果Parameters:x_data: 样本属性特征矩阵y_label: 属性特征对应的标签w：决策面的w参数b：决策面的参数b
"""
def plot_result(x_data, y_label, w, b):data_number, _ = x_data.shape; middle = int(data_number / 2)plt.scatter(x_data[:, 0], x_data[:, 1], c = y_label, cmap = pl.cm.Paired)x1, x2 = np.max(x_data), np.min(x_data)w1, w2 = w[0][0], w[1][0]y1, y2 = (-b - w1 * x1) / w2, (-b - w1 * x2) / w2plt.plot([float(x1), float(x2)], [float(y1), float(y2)])    # 绘制决策面for index, alpha in enumerate(alphas):if alpha > 0:b_temp = - w1 * x_data[index][0] - w2 * x_data[index][1]y1_temp, y2_temp = (-b_temp - w1 * x1) / w2, (-b_temp - w1 * x2) / w2plt.plot([float(x1), float(x2)], [float(y1_temp), float(y2_temp)], "k--")    # 绘制支持向量plt.scatter(x_data[index][0], x_data[index][1], s=150, c='none', alpha=0.7, linewidth=2, edgecolor='red')   # 圈出支持向量plt.show()if __name__ == "__main__":x_data, y_label = etablish_data(50)b, alphas = linear_smo(x_data, y_label, 0.8, 0.0001, 40)w = calc_w(x_data, y_label, alphas)plot_result(x_data, y_label, w, b)

完整代码：

import numpy as np
import pylab as pl
from matplotlib import pyplot as pltclass TearLinearSVM:def __init__(self):pass"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 用于生成训练数据集Parameters:data_number: 样本数据数目Return:x_data: 数据样本的属性矩阵y_label: 样本属性所对应的标签"""def etablish_data(self, data_number):x_data = np.concatenate((np.add(np.random.randn(data_number, 2), [3, 3]),       np.subtract(np.random.randn(data_number, 2), [3, 3])),axis = 0)      # random随机生成数据，+ -3达到不同类别数据分隔的目的 temp_data = np.zeros([data_number])temp_data.fill(-1)y_label = np.concatenate((temp_data, np.ones([data_number])), axis = 0)return x_data, y_label"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 随机选取alpha_jParameters:alpha_i_index: 第一个alpha的索引alpha_number: alpha总数目Return:alpha_j_index: 第二个alpha的索引"""def random_select_alpha_j(self, alpha_i_index, alpha_number):alpha_j_index = alpha_i_indexwhile alpha_j_index == alpha_i_index:alpha_j_index = np.random.randint(0, alpha_number)return alpha_j_index"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 使得alpha_j在[L, R]区间之内Parameters:alpha_j: 原始alpha_jL: 左边界值R: 右边界值Return:L,R,alpha_j: 修改之后的alpha_j"""def modify_alpha(self, alpha_j, L, R):if alpha_j < L: return Lif alpha_j > R: return Rreturn alpha_j"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 计算误差并返回"""def calc_E(self, alphas, y_lable, x_data, b, i):f_x_i = float(np.dot(np.multiply(alphas, y_lable).T, x_data * x_data[i, :].T)) + breturn f_x_i - float(y_label[i])"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 计算eta并返回"""def calc_eta(self, x_data, i, j):eta = 2.0 * x_data[i, :] * x_data[j, :].T \- x_data[i, :] * x_data[i, :].T \- x_data[j, :] * x_data[j,:].Treturn eta"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 计算b1, b2并返回"""def calc_b(self, b, x_data, y_label, alphas, alpha_i_old, alpha_j_old, E_i, E_j, i, j):b1 = b - E_i \- y_label[i] * (alphas[i] - alpha_i_old) * x_data[i, :] * x_data[i, :].T \- y_label[j] * (alphas[j] - alpha_j_old) * x_data[i, :] * x_data[j, :].Tb2 = b - E_j \- y_label[i] * (alphas[i] - alpha_i_old) * x_data[i, :] * x_data[j, :].T \- y_label[j] * (alphas[j] - alpha_j_old) * x_data[j, :] * x_data[j, :].Treturn b1, b2"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: SMO核心算法，求得b和laphas向量Parameters:x_data: 样本属性特征矩阵y_label: 属性特征对应的标签C：惩罚参数toler：容错率max_iter：迭代次数Return:b: 决策面的参数balphas：获取决策面参数w所需要的alphas"""def linear_smo(self, x_data, y_label, C, toler, max_iter):x_data = np.mat(x_data); y_label = np.mat(y_label).T     # 将数据转换为矩阵类型m, n = x_data.shape                                      # 得到数据样本的shape信息b, alphas, iter_num = 0, np.mat(np.zeros((m, 1))), 0     # 初始化参数b和alphas和迭代次数while (iter_num < max_iter):                            # 最多迭代max_iter次alpha_optimization_num = 0   # 定义优化次数for i in range(m):          # 遍历每个样本，一次选取一个样本计算误差E_i = self.calc_E(alphas, y_label, x_data, b, i)      # 样本i的误差计算if ((y_label[i] * E_i < -toler) and (alphas[i] < C)) or ((y_label[i] * E_i > toler) and (alphas[i] > 0)):j = self.random_select_alpha_j(i, m)              # 随机选取一个不与i重复jE_j = self.calc_E(alphas, y_label, x_data, b, j)  # 计算样本j的误差alpha_i_old = alphas[i].copy(); alpha_j_old = alphas[j].copy();if (y_label[i] != y_label[j]):               # 重新规范alphas的左右区间L, R = max(0, alphas[j] - alphas[i]), min(C, C + alphas[j] - alphas[i])else:L, R = max(0, alphas[j] + alphas[i] - C), min(C, alphas[j] + alphas[i])if L == R: print("L==R"); continue          # L==R时选取下一个样本eta = self.calc_eta(x_data, i, j)                  # 计算eta值if eta >= 0: print("eta>=0"); continuealphas[j] -= y_label[j] * (E_i - E_j) / etaalphas[j] = self.modify_alpha(alphas[j], L, R)     # 修改alpha[j]if (abs(alphas[j] - alpha_j_old) < 0.00001): print("alpha_j更改太小"); continuealphas[i] += y_label[j] * y_label[i] * (alpha_j_old - alphas[j])    # 修改alphas[i]b1, b2= self.calc_b(b, x_data, y_label, alphas, alpha_i_old, alpha_j_old, E_i, E_j, i, j)    # 计算b值if (0 < alphas[i]) and (C > alphas[i]): b = b1elif (0 < alphas[j]) and (C > alphas[j]): b = b2else: b = (b1 + b2)/2.0alpha_optimization_num += 1print("迭代次数:%d，样本:%d，alphas向量的优化次数:%d" % (iter_num, i+1, alpha_optimization_num))if (alpha_optimization_num == 0): iter_num += 1else: iter_num = 0print("迭代次数:%d" % iter_num)return b, alphas"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 根据公式计算出w权值向量Parameters:x_data: 样本属性特征矩阵y_label: 属性特征对应的标签alphas：linear_smo方法所返回的alphas向量Return:w: 决策面的参数w"""def calc_w(self, x_data, y_label, alphas):x_data, y_label, alphas = np.array(x_data), np.array(y_label), np.array(alphas)return np.dot((np.tile(y_label.reshape(1, -1).T, (1, 2)) * x_data).T, alphas).tolist()"""Author: Taoye微信公众号: 玩世不恭的CoderExplain: 绘制出分类结果Parameters:x_data: 样本属性特征矩阵y_label: 属性特征对应的标签w：决策面的w参数b：决策面的参数b"""def plot_result(self, x_data, y_label, w, b):data_number, _ = x_data.shape; middle = int(data_number / 2)plt.scatter(x_data[:, 0], x_data[:, 1], c = y_label, cmap = pl.cm.Paired)x1, x2 = np.max(x_data), np.min(x_data)w1, w2 = w[0][0], w[1][0]y1, y2 = (-b - w1 * x1) / w2, (-b - w1 * x2) / w2plt.plot([float(x1), float(x2)], [float(y1), float(y2)])    # 绘制决策面for index, alpha in enumerate(alphas):if alpha > 0:b_temp = - w1 * x_data[index][0] - w2 * x_data[index][1]y1_temp, y2_temp = (-b_temp - w1 * x1) / w2, (-b_temp - w1 * x2) / w2plt.plot([float(x1), float(x2)], [float(y1_temp), float(y2_temp)], "k--")    # 绘制支持向量plt.scatter(x_data[index][0], x_data[index][1], s=150, c='none', alpha=0.7, linewidth=2, edgecolor='red')   # 圈出支持向量plt.show()if __name__ == '__main__':linear_svm = TearLinearSVM()x_data, y_label = linear_svm.etablish_data(50)b, alphas = linear_svm.linear_smo(x_data, y_label, 0.8, 0.0001, 40)w = linear_svm.calc_w(x_data, y_label, alphas)linear_svm.plot_result(x_data, y_label, w, b)

呼呼呼！可算是结束了，做个小总结吧。

SVM是学习机器学习必然接触到的一个重要算法，所以一定要对其内在原理了解清楚，并不是说一定要手撕SVM的完整代码，但最起码使用框架的时候要了解内部都做了什么“小动作”，不要为了用而用。

本文介绍了线性SVM的算法原理，主要分为了五个部分的内容。一、首先通过参考比较权威的书籍以及优秀资料对SVM做了一个比较“良心”的介绍，让读者对SVM有一个比较宏观的概念，这小子（SVM）究竟是谁？竟如此骚气，让不少研究者拜倒其石榴裙下。二、其次向读者介绍了线性SVM以及最大间隔，这部分也是手撕SVM必须要掌握的一些基本概念，并且最终得到了SVM最初的优化问题。三、利用拉格朗日乘数法构建最值问题，将优化问题中的约束问题集成到了目标函数本身，之后利用拉格朗日的对偶性，将最初的优化问题转化成了内层关于w、bw、bw、b的最小值，外层关于α\alphaα的对偶问题。四、对偶问题的求解，这也是SVM算法的最核心内容，先是对内层关于w、bw、bw、b的函数求导，然后代回原式，从而消掉w,bw,bw,b参数，只留下未知的α\alphaα，随后利用SMO算法求得迭代更新之后的α\alphaα。五、手撕线性SVM的代码实现，结果证明，分类的效果还不错，这是个好家伙！！！

说实话，这篇文章有点肝，也是挪用了不少其他任务的时间。

这篇文章仅仅只是手撕线性SVM，也就是说大多数据样本都可以被正确分类，但在实际中，许多的数据集都是线性不可分的，这个时候可能就要引入核函数的概念了。关于非线性SVM，我们留在之后再来肝。。。

本文参考了不少书籍资料以及许多大佬的技术文章，行文风格尽可能做到了通俗易懂，但其中涉及到的数学公式在所难免，还请诸读者静下心来慢慢品尝。由于个人水平有限，才疏学浅，对于SVM也只是略知皮毛，可能文中有不少表述稍有欠妥、有不少错误或不当之处，还请诸君批评指正，有任何问题欢迎在下方留言。

参考资料：

[1] 《机器学习实战》：Peter Harrington 人民邮电出版社
[2] 《统计学习方法》：李航第二版清华大学出版社
[3] 《机器学习》：周志华清华大学出版社
[4] 《微积分方法》：李亿民中国海洋出版社
[5] Support Vector Machines explained well（翻墙）：http://bytesizebio.net/2014/02/05/support-vector-machines-explained-well/
[6] 关于更为直观认识SVM的video（翻墙）：https://www.youtube.com/watch?v=3liCbRZPrZA
[7] 支持向量机（SVM）是什么意思？：https://www.zhihu.com/question/21094489/answer/86273196
[8] 看了这篇文章你还不懂SVM你就来打我：https://zhuanlan.zhihu.com/p/49331510
[9] 拉格朗日乘数法：https://www.cnblogs.com/maybe2030/p/4946256.html

推荐阅读

print( "Hello，NumPy！" )
干啥啥不行，吃饭第一名
Taoye渗透到一家黑平台总部，背后的真相细思极恐
《大话数据库》-SQL语句执行时，底层究竟做了什么小动作？
那些年，我们玩过的Git，真香
基于Ubuntu+Python+Tensorflow+Jupyter notebook搭建深度学习环境
网络爬虫之页面花式解析
手握手带你了解Docker容器技术
一文详解Hexo+Github小白建站
打开ElasticSearch、kibana、logstash的正确方式