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

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

前面在写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

符号说明

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

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

一、SVM是什么

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

image

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

总结一哈：

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

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

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

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

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

假如说，我们这里有一堆样本，也就是我们常说的训练集S=\{(x_1, y_1), (x_1, y_1), (x_3, y_3), ..., (x_n, y_n)\}，且x_i表示的是样本i的属性特征向量，其内部有多个不同属性，这里我们不妨指定每个样本含有两个属性特征，也就是说x_i=(x_{i1}, x_{i2})^T（之所以用列向量表示，主要是方便后面超平面的构建）。而y_i表示的是每个样本中所有属性特征所对应的标签，由于这里的问题属性二分类，所以y_i的值只存在两种。为此，我们可以通过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]                # 获取权重w
    xx = 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！" )

image

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

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

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

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

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

Ax+By+C=0

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

\left( \begin{matrix} w_1, w_2 \\ \end{matrix} \right) \times \left( \begin{matrix} x^{(1)} \\ x^{(2)} \end{matrix} \right) \ +b=0 \tag{2-1}

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

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

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

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

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

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

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

\hat{\gamma} = \mathop{min}\limits_{i=1,...,N}\hat{\gamma_i} \tag{2-4}

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

\gamma_i=\frac{|w^Tx_i + b|}{||w||} \tag{2-5}

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

\gamma_i=y_i\frac{w^Tx_i + b}{||w||} \tag{2-6}

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

\gamma=\mathop{min}\limits_{i=1,...,N}\gamma_i \tag{2-7}

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

\gamma_i=\frac{\hat{\gamma_i}}{||w||}，即\gamma=\frac{\hat{\gamma}}{||w||} \tag{2-8}

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

\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}

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

\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和b的时候，函数间隔虽然发生了改变，但并不会影响我们的几何间隔——目标函数，也就是说，此时产生了一个等价的最优化问题。为了方便描述问题，我们不妨取\hat{\gamma}=1。另外，我们注意到目标函数\mathop{max}\limits_{w,b} \frac{\hat{\gamma}}{||w||}可以等价于\mathop{min}\limits_{w,b}\frac{1}{2}||w||^2，（注意，仅仅是等价。而之所以等价为二分之一、平方的形式，主要是方便后期的求导操作），对此，我们可以将数学表达的优化问题转化为如下形式：
\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}

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

图片来源：西瓜书

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

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

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

image

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

\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_1=w_2=\frac{1}{2},b=-2，于是最大间隔分离的决策面为：

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

其中，x_1=(3,3)^T与x_3=(1,1)^T为支持向量。

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

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

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

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

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

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

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

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

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

例子2：求函数f(x, y)=4x+xy^2+y^2在D:x^2+y^2\leq1上的最值。

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

L(x, y, \lambda)=4x+xy^2+y^2+\lambda(x^2+y^2-1)

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

\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}

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

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

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

\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}

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

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}

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

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

\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}

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

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

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

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

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

令：

\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}]

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

\sum^N_{i=1}\alpha_iy_i(w^Tx_i+b) - \sum^N_{i=1}\alpha_i \geq0

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

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

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

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

\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}

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

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

综上，只有满足KKT条件的时候，才能到达d^*= p^*，即在满足KKT条件的时候，才能将原始目标问题，转化为对偶问题，此时两者是等价的，且此时的目标函数为内层关于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{3-8}

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

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

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

四、模型的求解

关于拉格朗日的内层求解

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

\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}

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

image

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

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

w=\sum_{i=1}^N\alpha_iy_ix_i \\ 0=\sum_{i=1}^N\alpha_iy_i \tag{4-2}

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

image

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

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

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

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

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

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

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

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

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

\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}

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

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

\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}

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

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

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

image

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

\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}

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

s(\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
对\alpha_1,\alpha_2求偏导数并令其为0，易知s(\alpha_1,\alpha_2)在点(\frac{3}{2},-1)^T取极值，但该点不满足约束条件\alpha_2\geq0，所以最小值应在边界上达到。

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

这样的话，就说明\alpha_1=\alpha_3=\frac{1}{4}所对应的点x_1,x_3就是支持向量了，根据w=\sum_{i=1}^N\alpha_iy_ix_i，我们可以求得此时w_1=w_2=\frac{1}{2}，再由b=y-\sum_{i=1}^N\alpha_iy_ix_i^Tx，可以得到b=-2

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

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

其中，x_1=(3,3)^T与x_3=(1,1)^T为支持向量。

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

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

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

基于SMO算法的外层求解

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

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

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

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

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

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

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

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

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

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

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

\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}

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

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

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

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

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

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

image

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

\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}

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

\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}

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

\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(\alpha_2)仅仅只是关于\alpha_2的函数，我们将W(\alpha_2)对\alpha_2进行求导，并令导数为0：

image

image

image

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

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

(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}

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

\alpha_1^{old}y_1 + \alpha_2^{old}y_2=G \tag{4-9}

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

\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}

之后我们通过前面拉格朗日得到的关系式，用f(x_1)、f(x_2)来代替（4-10）后面的两个级数，整理最终得到：
\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}

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

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

w=\sum_{i=1}^N\alpha_iy_ix_i

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

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

而我们知道：

\alpha_1^{new}y_1+\alpha_2^{new}y_2=\alpha_1^{old}y_1+\alpha_2^{old}y_2=G \tag{4-13}

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

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,R的范围，我们会在下面实际编程的时候需要用到，主要是用来更新\alpha值。

接下来，就是更新b值了。前面我们已经定义了惩罚参数C，且0\leq\alpha_i^{new} \leq C，此时通过E_1、E_2更新得到的b_1、b_2固然是相等的；但假如我们更新之后的\alpha不在这个区间之内，则此时得到的b_1、b_2是不相等的，所以我们需要确定在不同情况下更新之后的b值：

前面，我们已经得到：

y_1=\sum_{i=1}^N\alpha_iy_ix_ix_1+b

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

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}

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

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}

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

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}

当更新之后的\alpha_1^{new}、\alpha_2^{new}都有效的时候，即在[0-C]区间之内时，此时b^{new}=b_1^{new}=b_2^{new}，而在不满足上述条件的时候，我们更新之后的b^{new}取b_1和b_2的均值，即：

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}

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

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

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

五、编程手撕线性SVM

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

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

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

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: 用于生成训练数据集
    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来进行更新，这里我们不妨将第一个\alpha_i通过遍历的形式逐个选取，而另一个则通过np.random模块来随机选取，这里需要主要的是，第二个选取的\alpha_j不能与第一个相同。为此，我们定义一个方法random_select_alpha_j来随机选取第二个\alpha_j（第一个已经通过遍历得到）:

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: 随机选取alpha_j
    Parameters:
        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_index
    while alpha_j_index == alpha_i_index:
        alpha_j_index = np.random.randint(0, alpha_number)
    return alpha_j_index

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

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

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

\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}

E_i=f(x_i)-y_i \tag{5-2}

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

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: 计算误差并返回
"""
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)) + b
    return f_x_i - float(y_label[i])

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

\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}

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: 计算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,:].T
    return eta
    
"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: 计算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, :].T
    b2 = 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, :].T
    return b1, b2

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

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: SMO核心算法，求得b和laphas向量
    Parameters:
        x_data: 样本属性特征矩阵
        y_label: 属性特征对应的标签
        C：惩罚参数
        toler：容错率
        max_iter：迭代次数
    Return:
        b: 决策面的参数b
        alphas：获取决策面参数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重复j
                E_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"); continue
                alphas[j] -= y_label[j] * (E_i - E_j) / eta
                alphas[j] = modify_alpha(alphas[j], L, R)     # 修改alpha[j]
                if (abs(alphas[j] - alpha_j_old) < 0.00001): print("alpha_j更改太小"); continue
                alphas[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 = b1
                elif (0 < alphas[j]) and (C > alphas[j]): b = b2
                else: b = (b1 + b2)/2.0
                alpha_optimization_num += 1
                print("迭代次数:%d，样本:%d，alphas向量的优化次数:%d" % (iter_num, i+1, alpha_optimization_num))
        if (alpha_optimization_num == 0): iter_num += 1
        else: iter_num = 0
        print("迭代次数:%d" % iter_num)
    return b, alphas

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

w=\sum_{i=1}^N\alpha_iy_ix_i \tag{5-4}

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: 根据公式计算出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()

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

import pylab as pl

"""
    Author: Taoye
    微信公众号: 玩世不恭的Coder
    Explain: 绘制出分类结果
    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) / w2
    plt.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) / w2
            plt.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 plt

class TearLinearSVM:
    def __init__(self):
        pass

    """
        Author: Taoye
        微信公众号: 玩世不恭的Coder
        Explain: 用于生成训练数据集
        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
        微信公众号: 玩世不恭的Coder
        Explain: 随机选取alpha_j
        Parameters:
            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_index
        while alpha_j_index == alpha_i_index:
            alpha_j_index = np.random.randint(0, alpha_number)
        return alpha_j_index

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

    """
        Author: Taoye
        微信公众号: 玩世不恭的Coder
        Explain: 计算误差并返回
    """
    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)) + b
        return f_x_i - float(y_label[i])

    """
        Author: Taoye
        微信公众号: 玩世不恭的Coder
        Explain: 计算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,:].T
        return eta

    """
        Author: Taoye
        微信公众号: 玩世不恭的Coder
        Explain: 计算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, :].T
        b2 = 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, :].T
        return b1, b2

    """
        Author: Taoye
        微信公众号: 玩世不恭的Coder
        Explain: SMO核心算法，求得b和laphas向量
        Parameters:
            x_data: 样本属性特征矩阵
            y_label: 属性特征对应的标签
            C：惩罚参数
            toler：容错率
            max_iter：迭代次数
        Return:
            b: 决策面的参数b
            alphas：获取决策面参数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重复j
                    E_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"); continue
                    alphas[j] -= y_label[j] * (E_i - E_j) / eta
                    alphas[j] = self.modify_alpha(alphas[j], L, R)     # 修改alpha[j]
                    if (abs(alphas[j] - alpha_j_old) < 0.00001): print("alpha_j更改太小"); continue
                    alphas[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 = b1
                    elif (0 < alphas[j]) and (C > alphas[j]): b = b2
                    else: b = (b1 + b2)/2.0
                    alpha_optimization_num += 1
                    print("迭代次数:%d，样本:%d，alphas向量的优化次数:%d" % (iter_num, i+1, alpha_optimization_num))
            if (alpha_optimization_num == 0): iter_num += 1
            else: iter_num = 0
            print("迭代次数:%d" % iter_num)
        return b, alphas

    """
        Author: Taoye
        微信公众号: 玩世不恭的Coder
        Explain: 根据公式计算出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
        微信公众号: 玩世不恭的Coder
        Explain: 绘制出分类结果
        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) / w2
        plt.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) / w2
                plt.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、b的最小值，外层关于\alpha的对偶问题。四、对偶问题的求解，这也是SVM算法的最核心内容，先是对内层关于w、b的函数求导，然后代回原式，从而消掉w,b参数，只留下未知的\alpha，随后利用SMO算法求得迭代更新之后的\alpha。五、手撕线性SVM的代码实现，结果证明，分类的效果还不错，这是个好家伙！！！

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

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

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

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

参考资料：

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

推荐阅读

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

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 212,080评论 6赞 493
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,422评论 3赞 385
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 157,630评论 0赞 348
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,554评论 1赞 284
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 65,662评论 6赞 386
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,856评论 1赞 290
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,014评论 3赞 408
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,752评论 0赞 268
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,212评论 1赞 303
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,541评论 2赞 327
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,687评论 1赞 341
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,347评论 4赞 331
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,973评论 3赞 315
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,777评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,006评论 1赞 266
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,406评论 2赞 360
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,576评论 2赞 349

推荐阅读更多精彩内容