支持向量机

目录

• 支持向量机简介
• 线性可分支持向量机与硬间隔最大化
• 线性可分支持向量机
• 函数间隔和几何间隔
• 间隔最大化
• 学习的对偶算法
• 线性支持向量机与软间隔最大化
• 线性支持向量机
• 学习的对偶算法
• 支持向量
• 合页损失函数
• 非线性支持向量机与核函数
• 核技巧
• 正定核
• 常用核函数
• 非线性支持向量分类机
• 序列最小最优化算法
• 两个变量二次规划的求解方法
• 变量的选择方法
• SMO算法
• SMO算法实现
• 简化版SMO算法
• 完整版SMO算法
• 非线性支持向量机SMO算法
• 附录
• 凸二次规划（convex quadratic programming)问题
• Gram矩阵
• 希尔伯特空间
• 赋范向量空间
• 核或正定核
• Reference

本文是《统计学习方法》第七章《支持向量机》的笔记，附带了少量注解和背景知识的补充；后半部分将《机器学习实战》支持向量机的Python代码加以整理注释，与公式放到一起形成对照，辅助理解。私以为，没有泛函分析基础的人是无法深刻理解支持向量机的，更不能随便写“入门指南”误导别人。所以这篇笔记既不敢删教材的内容，也不敢加以发挥，只是将李航老师的教材搬运过来，按李航老师的博客上的勘误表修正了两个错误，用括号注解了自己的一点理解，并在文末以附录的形式补充了一些背景知识。

本笔记服务于自己备忘，不做其他用途。差点烂在草稿箱里，还是发出来方便阅读一些，毕竟正在一边参考这篇笔记，一边看libsvm和liblinear的源码。

支持向量机简介

支持向量机（support vector machines，SVM)是一种二类分类模型。它的基本模型是定义在特征空间上的间隔最大的线性分类器，间隔最大使它有别于感知机；支持向量机还包括核技巧，这使它成为实质上的非线性分类器。支持向量机的学习策略就是间隔最大化，可形式化为一个求解凸二次规划（convex quadratic programming，不怕，附录有解释)的问题，也等价于正则化的合页损失函数（后面也有解释）的最小化问题。支持向量机的学习算法是求解凸二次规划的最优化算法。

支持向量机学习方法包含构建由简至繁的模型：线性可分支持向量机（linear support vector machine in linearly separable case)、线性支持向量机（linear support vector machine)及非线性支持向量机（non-linear support vector machine)。简单模型是复杂模型的基础，也是复杂模型的特殊情况。当训练数据线性可分时，通过硬间隔最大化（hard margin maximization)，学习一个线性的分类器，即线性可 分支持向量机，又称为硬间隔支持向量机；当训练数据近似线性可分时，通过软间隔最大化（soft margin maximization)，也学习一个线性的分类器，即线性支持向量机，又称为软间隔支持向量机；当训练数据线性不可分时，通过使用核技巧 (kernel trick)及软间隔最大化，学习非线性支持向量机。

本章按照上述思路介绍3类支持向量机、核函数及一种快速学习算法——序列最小最优化算法（SMO)。

线性可分支持向量机与硬间隔最大化

线性可分支持向量机

考虑一个二类分类问题。假设输入空间与特征空间为两个不同的空间。输入空间为欧氏空间或离散集合，特征空间为欧氏空间或希尔伯特空间。线性可分支持向量机、线性支持向量机假设这两个空间的元素一一对应，并将输入空间中的输入映射为特征空间中的特征向量。非线性支持向量机利用一个从输入空间到特 征空间的非线性映射将输入映射为特征向量。所以，输入都由输入空间转换到特征空间，支持向量机的学习是在特征空间进行的。

假设给定一个特征空间上的训练数据集

其中，，xi为第i个特征向量，也称为实例，yi为xi的类标记，正负时分别称xi为正例或负例；称为样本点，再假设训练数据集是线性可分的。

学习的目标是在特征空间中找到一个分离超平面，能将实例分到不同的类。分离超平面对应于方程,它由法向量w和截距b决定，可用(w,b) 来表示。分离超平面将特征空间划分为两部分，一部分是正类，一部分是负类。法向量指向的一侧为正类，另一侧为负类。

一般地，当训练数据集线性可分时，存在无穷个分离超平面可将两类数据正确分开。感知机利用误分类最小的策略，求得分离超平面，不过这时的解有无穷多个。线性可分支持向量机利用间隔最大化求最优分离超平面，这时，解是唯一的。

定义 (线性可分支持向量机）给定线性可分训练数据集，通过间隔最大化或等价地求解相应的凸二次规划问题学习得到的分离超平面为

以及相应的分类决策函数

称为线性可分支持向量机。

考虑如图所示的二维特征空间中的分类问题。图中“。”表示正例，“x ”表示负例。训练数据集线性可分，这时有许多直线能将两类数据正确划分。线性可分支持向量机对应着将两类数据正确划分并且间隔最大的直线，如图所示。

间隔最大及相应的约束最优化问题将在下面叙述。这里先介绍函数间隔和几何间隔的概念。

函数间隔和几何间隔

在上图中，有A、B、 C三个点，表示3个实例，均在分离超平面的正类一侧，预测它们的类。点A距分离超平面较远，若预测该点为正类，就比较确信预测是正确的；点C距分离超平面较近，若预测该点为正类就不那么确信；点B介于点A与C之间，预测其为正类的确信度也在A与C之间。

一般来说，一个点距离分离超平面的远近可以表示分类预测的确信程度。在超平面w*x + b = 0确定的情况下，| w*x + b|能够相对地表示点x距离超平面的远近。而w*x + b 的符号与类标记y的符号是否一致能够表示分类是否正确。所以可用量y(w*x + b)来表示分类的正确性及确信度，这就是函数间隔(functional margin)的概念。

定义(函数间隔）对于给定的训练数据集T和超平面(w,b),定义超平 面(w,b)关于样本点(xi，yi）的函数间隔为

定义超平面(w,b)关于训练数据集T的函数间隔为超平面(w,b)关于T中所有样本点的函数间隔之最小值，即

函数间隔可以表示分类预测的正确性及确信度。但是选择分离超平面时，只有函数间隔还不够。因为只要成比例地改变w和b,例如将它们改为2w和2b,超平面并没有改变，但函数间隔却成为原来的2倍。这一事实启示我们，可以对分离 超平面的法向量w加某些约束，如规范化，||w||=1,使得间隔是确定的。这时函数间隔成为几何间隔（geometric margin)。

下图给出了超平面(w，b)及其法向量w。点A表示某一实例xi，其类标记 为yi=+1。点A与超平面(w, b)的距离由线段AB给出，记作γi。

其中，||w||为w的L2范数。这是点A在超平面正的一侧的情形。如果点A在超平面负的一侧，即yi=-1,那么点与超平面的距离为

由这一事实导出几何间隔的概念。

定义(几何间隔）对于给定的训练数据集T和超平面(w，b),定义超平 面(w，b)关于样本点的几何间隔为

定义超平面(w，b)关于训练数据集T的几何间隔为超平面(w，b)关于T中所有样本点的几何间隔之最小值，即

超平面(w，b）关于样本点的几何间隔一般是实例点到超平面的带符号的距离（signed distance),当样本点被超平面正确分类时就是实例点到超平面的距离。

从函数间隔和几何间隔的定义可知，函数间隔和几何间隔有下面的关系：

如果||w||=l，那么函数间隔和几何间隔相等。如果超平面参数w和b成比例地改变（超平面没有改变)，函数间隔也按此比例改变，而几何间隔不变。

间隔最大化

支持向量机学习的基本想法是求解能够正确划分训练数据集并且几何间隔最大的分离超平面。对线性可分的训练数据集而言，线性可分分离超平面有无穷多个（等价于感知机)，但是几何间隔最大的分离超平面是唯一的。这里的间隔最大化又称为硬间隔最大化（与将要讨论的训练数据集近似线性可分时的软间隔最大化相对应)。

间隔最大化的直观解释是：对训练数据集找到几何间隔最大的超平面意味着以充分大的确信度对训练数据进行分类。也就是说，不仅将正负实例点分开，而且对最难分的实例点（离超平面最近的点）也有足够大的确信度将它们分开。这样的超平面应该对未知的新实例有很好的分类预测能力。

1、最大间隔分离超平面

下面考虑如何求得一个几何间隔最大的分离超平面，即最大间隔分离超平面。具体地，这个问题可以表示为下面的约束最优化问题：

即我们希望最大化超平面(w,b)关于训练数据集的几何间隔γ。约束条件表示的是超平面(w,b)关于每个训练样本点的几何间隔至少是γ。

考虑几何间隔和函数间隔的关系式,可将这个问题改写为

也就是把几何间隔换为函数间隔。

函数间隔的取值并不影响最优化问题的解。事实上，假设将w和b按比例改变为λw和λb，这时函数间隔成为λ。函数间隔的这一改变对上面最优化问题的不等式约束没有影响，对目标函数的优化也没有影响，也就是说，它产生一个等价的最优化问题。这样，就可以取=l。将=l代入上面的最优化问题，注意到最大化和最小化是等价的，于是就得到下面的线性可分支持向量机学习的最优化问题

这是一个凸二次规划（convex quadratic programming)问题。 

凸优化问题是指约束最优化问题

其中，目标函数f(w)和约束函数都是上的连续可微的凸函数，约束函数是上的仿射函数：

当目标函数f(w)是二次函数且约束函数g(w)是仿射函数时，上述凸最优化问题被称作凸二次规划问题。

如果求出了约束最优化问题的解,那么就可以得到最大间隔分离超平面及分类决策函数,即线性可分支持向量机模型。

综上所述，就有下面的线性可分支持向量机的学习算法——最大间隔法(maximum margin method)。

算法(线性可分支持向量机学习算法——最大间隔法）

估计是本着严密的作风，《统计学习方法》展开了两项证明：

2、最大间隔分离超平面的存在唯一性

线性可分训练数据集的最大间隔分离超平面是存在且唯一的。似乎也不难理解，“最大”从感觉上讲是个很强的约束。

定理(最大间隔分离超平面的存在唯一性）若训练数据集r线性可分，则可将训练数据集中的样本点完全正确分开的最大间隔分离超平面存在且唯一。

证明

（1)存在性

由于训练数据集线性可分，所以算法7.1中的最优化问题—定存在可行解。又由于目标函数有下界，所以这个最优化问题必有解，记作。由于训练数据集中既有正类点又有负类点，所以不是最优化的可行解，因而最优解必满足。由此得知分离超平面的存在性。

(2)唯一性

首先证明最优化问题解中的唯一性。假设问题存在两个最优解和。显然（因为min值只有一个呀）,其中c是一个常数。令，易知(w,b)是问题的可行解（这里应该指的是满足s.t.，但不一定能得到最小值），从而有

上式表明，式中的不等号可变为等号，即，从而有（两个向量相等或反向）。若λ=-1，则w=0，(w,b)不是问题的可行解，矛盾。因此必有λ=1,即

由此可以把两个最优解和分别写成和。再证。设是集合中分别对应于和使得问题的不等式等号成立的点，是集合中分别对应于和使得问题的不等式等号成立的点，则由，得

又因为(看s.t.)

所以，。同理有。因此，

由和可知，两个最优解和是相同的，解的唯一性得证。

由问题解的唯一性即得分离超平面是唯一的。

(3)分离超平面能将训练数据集中的两类点完全正确地分开。

由解满足问题的约束条件即可得知（这其实是一句废话）。

3. 支持向量和间隔边界

在线性可分情况下，训练数据集的样本点中与分离超平面距离最近的样本点的实例称为支持向量（support vector)。支持向量是使约束条件式等号成立的点，即

对的正例点，支持向量在超平面

上，对的负例点，支持向量在超平面

上。如图所示，在H1和H2上的点就是支持向量。

注意到H1和H2平行，并且没有实例点落在它们中间。在H1和H2之间形成一条长带，分离超平面与它们平行且位于它们中央。长带的宽度，即H1和H2之间的距离称为间隔（margin)。间隔依赖于分离超平面的法向量w，等于。H1和H2称为间隔边界。

在决定分离超平面时只有支持向量起作用，而其他实例点并不起作用。如果移动支持向量将改变所求的解；但是如果在间隔边界以外移动其他实例点，甚至去掉这些点，则解是不会改变的。由于支持向量在确定分离超平面中起着决定性作用，所以将这种分类模型称为支持向量机。支持向量的个数一般很少，所以支持向量机由很少的“重要的”训练样本确定。

学习的对偶算法

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

首先构建拉格朗日函数(Lagrange function)。为此，对每一个不等式约束引进拉格朗日乘子（Lagrauge multiplier),定义拉格朗日函数：

其中，为拉格朗日乘子向量。

根据拉格朗日对偶性，原始问题的对偶问题是极大极小问题：

所以，为了得到对偶问题的解，需要先求对的极小，再求对α的极大。

(1)求

将拉格朗日函数分别对求偏导数并令其等于0。

得

将式代入拉格朗日函数，并利用式,即得

式子长了点，但是并不难。是w的L2范数，所谓L2范数，指的是向量各元素的平方和然后求平方根（长度）。因为有个平方，所以上式中出现了i和j遍历相乘的现象。

那么就有

(2)求对α的极大，即是对偶问题

将上式的目标函数由求极大转换成求极小，就得到下面与之等价的对偶最优化问题： 

也就是将目标函数取个相反数而已，极大就变成极小了。

考虑原始最优化问题和对偶最优化问题,原始问题满足定理C.2的条件，所以存在使是原始问题的解，是对偶问题的解。这意味着求解原始问题可以转换为求解对偶问题。

对线性可分训练数据集，假设对偶最优化问题对α的解为，可以由求得原始最优化问题对(w,b)的解。有下面的定理。

定理7.2 设是对偶最优化问题的解，则存在下标j，使得,并可按下式求得原始最优化问题的解(w,b)：

证明 根据定理C.3,KKT条件成立，即得

原始问题中的其实就是定理C.3中的ci的相反数，这样其实就是。

由此得（也就是第一个式子）

其中至少有一个(用反证法，假设,由式可知,而不是原始最优化问题的解，产生矛盾)。既然，又因为刚才的，所以对此j有

将式代入式并注意到，即得

由此定理可知，分离超平面可以写成

分类决策函数可以写成

这就是说，分类决策函数只依赖于输入x和训练样本输入的内积。式称为线性可分支持向量机的对偶形式。

综上所述，对于给定的线性可分训练数据集，可以首先求对偶问题的解再利用式和式求得原始问题的解;从而得到分离超平面及分类决策函数。这种算法称为线性可分支持向量机的对偶学习算法，是线性可分支持向量机学习的基本算法。

算法(线性可分支持向量机学习算法）

输入:线性可分训练集，其中

输出：分离超平面和分类决策函数。

(1)构造并求解约束最优化问题

求得最优解。

(2)计算

并选择的一个正分量,计算

(3)求得分离超平面

分类决策函数：

在线性可分支持向量机中，由式、式可知，只依赖于训练数据中对应于的样本点，而其他样本点对没有影响。我们将训练数据中对应于的实例点称为支持向量。

定义(支持向量）考虑原始最优化问题及对偶最优化问题,将训练数据集中对应于的样本点的实例称为支持向量。

根据这一定义，支持向量一定在间隔边界上。由KKT互补条件之和-ci=可知，

对应于的实例，有

又因为yi=±1，所以上式等效于

即xi一定在间隔边界上。这里的支持向量的定义与前面给出的支持向量的定义是—致的。

对于线性可分问题，上述线性可分支持向量机的学习（硬间隔最大化）算法是完美的。但是，训练数据集线性可分是理想的情形。在现实问题中，训练数据集往往是线性不可分的，即在样本中出现噪声或特异点。此时，有更一般的学习算法。

线性支持向量机与软间隔最大化

线性支持向量机

线性可分问题的支持向量机学习方法，对线性不可分训练数据是不适用的，因为这时上述方法中的不等式约束并不能都成立。怎么才能将它扩展到线性不可分问题呢？这就需要修改硬间隔最大化，使其成为软间隔最大化。

假设给定一个特征空间上的训练数据集

其中，为第i个特征向童，yi为xi的类标记。再假设训练数据集不是线性可分的。通常情况是，训练数据中有一些特异点（outlier),将这些特异点除去后，剩下大部分的样本点组成的集合是线性可分的。

线性不可分意味着某些样本点不能满足函数间隔大于等于1的约束条件。为了解决这个问题，可以对每个样本点引进一个松池变量，使函数间隔加上松弛变量大于等于1。这样，约束条件变为

同时，对每个松弛变量支付一个代价目标函数由原来的变成

这里，C>0称为惩罚参数，一般由应用问题决定，C值大时对误分类的惩罚增大，C值小时对误分类的惩罚减小。最小化目标函数包含两层含义:使尽量小即间隔尽量大，同时使误分类点的个数尽量小，C是调和二者的系数。

有了上面的思路，可以和训练数据集线性可分时一样来考虑训练数据集线性不可分时的线性支持向量机学习问题。相应于硬间隔最大化，它称为软间隔最大化。线性不可分的线性支持向量机的学习问题变成如下凸二次规划（convex quadratic programming)问题（原始问题）：

上述原始问题是一个凸二次规划问题，因而关于的解是存在的。可以证明w的解是唯一的，但b的解不唯一，b的解存在于一个区间。（然而没有给出证明）

设问题的解是,于是可以得到分离超平面及分类决策函数。称这样的模型为训练样本线性不可分时的线性支持向量机，简称为线性支持向量机。显然，线性支持向量机包含线性可分支持向量机。由于现实中训练数据集往往是线性不可分的，线性支持向量机具有更广的适用性。

下面给出线性支持向量机的定义。

定义(线性支持向置机）对于给定的线性不可分的训练数据集，通过求解凸二次规划问题，即软间隔最大化问题，得到的分离超平面为

以及相应的分类决策函数

称为线性支持向量机。

学习的对偶算法

原始问题的对偶问题是

原始最优化问题的拉格朗日函数是

其中，。

对偶问题是拉格朗日函数的极大极小问题。首先求对的极小，由

得

将上面三个式子代入,得

再对求α的极大，即得对偶问题：

将上述对偶最优化问题进行变换：利用等式约束消去，从而只留下变量，并将约束写成

再将对目标函数求极大转换为求极小，于是得到对偶问题。

可以通过求解对偶问题而得到原始问题的解，进而确定分离超平面和决策函数。为此，就可以定理的形式叙述原始问题的最优解和对偶问题的最优解的关系。

定理 设是对偶问题的一个解，若存在的一个分量,则原始问题的解可按下式求得：

证明 原始问题是凸二次规划问题，解满足KKT条件。即得

由式易知式成立。再由式可知，若存在，，则（其实我觉得这么想更好理解，由和知，括号里面的只能取0。然后由和知＞0。而，所以。最后回过头来利用=0这个推论，得到）。由此即得式。

由此定理可知，分离超平面可以写成

分类决策函数可以写成

综合前面的结果，有下面的算法。

算法(线性支持向量机学习算法）

输入:线性可分训练集，其中

输出：分离超平面和分类决策函数。

(1)选择惩罚参数C>0,构造并求解凸二次规划问题

求得最优解。

(2)计算

选择的一个分量适合条件,计算

(3)求得分离超平面

分类决策函数：

步骤(2)中，对任一适合条件的，按式都可求出,但是由于原始问题对b的解并不唯一，所以实际计算时可以取在所有符合条件的样本点上的平均值。

支持向量

在线性不可分的情况下，将对偶问题的解中对应于的样本点的实例,称为支持向量（软间隔的支持向量)。如图所示，

这时的支持向量要比线性可分时的情况复杂一些。图中，分离超平面由实线表示，间隔边界由虚线表示，正例点由“。”表示，负例点由“X”表示。图中还标出了实例到间隔边界的距离。

软间隔的支持向量或者在间隔边界上，或者在间隔边界与分离超平面之间，或者在分离超平面误分一侧。若则（证明参考）,支持向量恰好落在间隔边界上：若,则分类正确，在间隔边界与分离超平面之间：若，则在分离超平面上：若,则位于分离超平面误分一侧。

合页损失函数

对于线性支持向量机学习来说，其模型为分离超平面及决策函数,其学习策略为软间隔最大化，学习算法为凸二次规划。线性支持向量机学习还有另外一种解释，就是最小化以下目标函数：

目标函数的第1项是经验损失或经验风险，函数 

称为合页损失函数（hinge loss function)。下标“+”表示以下取正值的函数。

这就是说，当样本点被正确分类且函数间隔（确信度）大于1时，损失是0,否则损失是。现在回过头来看这张图

注意到在图中的实例点x4被正确分类，但损失不是0。目标函数的第2项是系数为的w的L2范数，是正则化项。

定理 线性支持向量机原始最优化问题：

等价于最优化问题

证明 可将最优化问题写成问题。令

则。于是满足约束条件。由有，所以最优化问题可写成

若取则

与式等价。

反之，也可将最优化问题表示成问题。

合页损失函数的图形如图所示

横轴是函数间隔，纵轴是损失。由于函数形状像一个合页（我怎么就看不出来）

故名合页损失函数。

图中还画出0-1损失函数，可以认为它是二类分类问题的真正的损失函数，而合页损失函数是0-1损失函数的上界。由于0-1损失函数不是连续可导的，直接优化由其构成的目标函数比较困难，可以认为线性支持向童机是优化由0-1损失函数的上界（合页损失函数）构成的目标函数。这时的上界损失函数又称为代理损失函数（surrogate loss function)。

图中虚线显示的是感知机的损失函数。这时，当样本点被正确分类时，损失是0,否则损失是。相比之下，合页损失函数不仅要分类正确，而且确信度足够高时损失才是0。也就是说，合页损失函数对学习有更高的要求。

非线性支持向量机与核函数

对解线性分类问题，线性分类支持向量机是一种非常有效的方法。但是，有时分类问题是非线性的，这时可以使用非线性支持向量机。本节叙述非线性支持向量机，其主要特点是利用核技巧（kernel trick)。为此，先要介绍核技巧。核技巧不仅应用于支持向量机，而且应用于其他统计学习问题。

核技巧

1、非线性分类问题

非线性分类问题是指通过利用非线性模型才能很好地进行分类的问题。先看一个例子。

如左图，是一个分类问题，图中“•”表示正实例点，“x”表示负实例点。由图可见，无法用直线（线性模型）将正负实例正确分开，但可以用一条椭圆曲线（非线性模型）将它们正确分开。

一般来说，对给定的一个训练数据集，其中实例属于输入空间，，对应的标记有两类。如果能用中的一个超曲面将正负例正确分开，则称这个问题为非线性可分问题。

非线性问题往往不好求解，所以希望能用解线性分类问题的方法解决这个问题。所采取的方法是进行一个非线性变换，将非线性问题变换为线性问题，通过解变换后的线性问题的方法求解原来的非线性问题。对图7.7所示的例子，通过变换，将左图中椭圆变换成右图中的直线，将非线性分类问题变换为线性分类问题。

设原空间为，新空间为，定义从原空间到新空间的变换（映射)：

经过变换,原空间变换为新空间,原空间中的点相应地变换为新空间中的点，原空间中的椭圆

变换成为新空间中的直线

在变换后的新空间里，直线可以将变换后的正负实例点正确分开。这样，原空间的非线性可分问题就变成了新空间的线性可分问题。

上面的例子说明，用线性分类方法求解非线性分类问题分为两步：首先使用一个变换将原空间的数据映射到新空间：然后在新空间里用线性分类学习方法从训练数据中学习分类模型。核技巧就属于这样的方法。

核技巧应用到支持向量机，其基本想法就是通过一个非线性变换将输入空间(欧氏空间或离散集合）对应于一个特征空间（希尔伯特空间),使得在输入空间中的超曲面模型对应于特征空间中的超平面模型（支持向量机)。这样，分类问题的学习任务通过在特征空间中求解线性支持向量机就可以完成。

2、核函数的定义

定义(核函数）设是输入空间（欧氏空间的子集或离散集合)，又设为特征空间（希尔伯特空间)，如果存在一个从到的映射

使得对所有，函数满足条件

则称为核函数，为映射函数，式中为和的内积。

核技巧的想法是，在学习与预测中只定义核函数，而不显式地定义映射函数。通常，直接计算比较容易，而通过和计算并不容易。注意，是输入空间到特征空间的映射，特征空间—般是高维的，甚至是无穷维的。可以看到，对于给定的核,特征空间和映射函数的取法并不唯一，可以取不同的特征空间，即便是在同一特征空间里也可以取不同的映射。

下面举一个简单的例子来说明核函数和映射函数的关系。

例 假设输入空间是R2,核函数是,试找出其相关的特征空间和映射。

解 取特征空间,记，由于

所以可以取映射

容易验证。

仍取以及

同样有。

还可以取和

3.核技巧在支持向量机中的应用

我们注意到在线性支持向量机的对偶问题中，无论是目标函数还是决策函数(分离超平面）都只涉及输入实例与实例之间的内积。在对偶问题的目标函数中的内积可以用核函数来代替。此时对偶问题的目标函数成为 

同样，分类决策函数中的内积也可以用核函数代替，而分类决策函数式成为

这等价于经过映射函数将原来的输入空间变换到一个新的特征空间，将输入空间中的内积变换为特征空间中的内积,在新的特征空间里从训练样本中学习线性支持向量机。当映射函数是非线性函数时，学习到的含有核函数的支持向量机是非线性分类模型。

也就是说，在核函数给定的条件下，可以利用解线性分类问题的方法求解非线性分类问题的支持向量机。学习是隐式地在特征空间进行的，不需要显式地定义特征空间和映射函数。这样的技巧称为核技巧，它是巧妙地利用线性分类学习方法与核函数解决非线性问题的技术。在实际应用中，往往依赖领域知识直接选择核函数，核函数选择的有效性需要通过实验验证。

正定核

接下来的内容涉及到大量泛函分析的背景知识，然而我并没有系统地学习过泛函分析，虽然在附录中硬着头皮找了些基础的定义，然而还是理解得不够透彻。

已知映射函数，可以通过和的内积求得核函数。不用构造映射能否直接判断一个给定的函数是不是核函数？或者说，函数满足什么条件才能成为核函数？

本节叙述正定核的充要条件。通常所说的核函数就是正定核函数（positive definite kernel function)。为证明此定理先介绍有关的预备知识。

假设是定义在上的对称函数，并且对任意的，关于的Gram矩阵是半正定的。可以依据函数，构成一个希尔伯特空间（Hilbert space),其步骤是：首先定义映射并构成向量空间;然后在上定义内积构成内积空间；最后将完备化构成希尔伯特空间。

1.定义映射，构成向量空间先定义映射

根据这一映射，对任意，定义线性组合

考虑由线性组合为元素的集合。由于集合对加法和数乘运算是封闭的，所以构成一个向量空间。

2.在上定义内积，使其成为内积空间

在上定义一个运算*:对任意,

定义运算*

（我觉得上式的·应该是等于号的意思，不知道对不对。果然，从李航博士的勘误表中找到了正确的式子）

证明运算*是空间的内积。为此要证：

其中，(1)〜(3)由式及的对称性容易得到。现证⑷式之。由式及式可得：

由Gram矩阵的半正定性知上式右端非负，即。

书上是如是说的，但是我找到了更好懂的证明：

Reference：机器学习——核函数讲义.pdf

再证(4)之式。充分性显然。为证必要性，首先证明不等式：

设，则，于是，

其左端是的二次三项式，非负，其判别式小于等于0,即

于是式得证。现证若,则。事实上，若 

则按运算*的定义式，对任意的，有

于是，

由式和式有

由式有

此式表明，当时，对任意的x都有。

至此，证明了*为向量空间的内积。赋予内积的向量空间为内积空间。因此是一个内积空间。既然*为的内积运算，那么仍然用表示，即若

则

3.将内积空间完备化为希尔伯特空间

现在将内积空间完备化。由式定义的内积可以得到范数

因此，是一个赋范向量空间。根据泛函分析理论，对于不完备的赋范向量空间,一定可以使之完备化，得到完备的赋范向量空间。—个内积空间，当作为一个赋范向量空间是完备的时候，就是希尔伯特空间。这样，就得到了希尔伯特空间。

这一希尔伯特空间称为再生核希尔伯特空间（reproducing kernel Hilbertspace,RKHS)。这是由于核尺具有再生性，即满足 

及

称为再生核。

4.正定核的充要条件

定理(正定核的充要条件） 设是对称函数，则为正定核函数的充要条件是对任意,对应的Gram矩阵：

是半正定矩阵。

证明 必要性。由于是上的正定核，所以存在从到希尔伯特空间的映射，使得

于是，对任意,构造关于的Gram矩阵

表明关于的Gram矩阵是半正定的。

充分性。己知对称函数对任意，关于的Gram矩阵是半正定的。根据前面的结果，对给定的，可以构造从到某个希尔伯特空间的映射：

由式可知:

由式即得

表明是上的核函数。 

定理给出了正定核的充要条件，因此可以作为正定核，即核函数的另一定义。

定义(正定核的等价定义）设,是定义在上的对称函数，如果对任意,对应的Gram矩阵

是半正定矩阵，则称是正定核。

这一定义在构造核函数时很有用。但对于一个具体函数来说，检验它是否为正定核函数并不容易，因为要求对任意有限输入集验证对应的Gram矩阵是否为半正定的。在实际问题中往往应用已有的核函数。另外，由Mercer定理可以得到Mercer核（Mercer Kernel)，正定核比Mercer核更具一般性。下面介绍一些常用的核函数。

常用核函数

1、多项式核函数（polynomial kernel function）

对应的支持向量机是一个p次多项式分类器。在此情形下，分类决策函数成为

2、高斯核函数（Gaussian kernel function)

对应的支持向量机是高斯径向基函数（radial basis function)分类器。在此情形下，分类决策函数成为

3.字符串核函数（string kernel function)

核函数不仅可以定义在欧氏空间上，还可以定义在离散数据的集合上。比如，字符串核是定义在字符串集合上的核函数。字符串核函数在文本分类、信息检索、生物信息学等方面都有应用。

考虑一个有限字符表，字符串是从中取出的有限个字符的序列，包括空字符串。字符串的长度用表示，它的元素记作。两个字符串和的连接记作。所有长度为n的字符串的集合记作，所有字符串的集合记作。

考虑字符串的子串。给定一个指标序列（其实是下标序列吧），的子串定义为，其长度记作。如果i是连续的，则，否则，。

假设是长度大于或等于n字符串的集合，是的元素。现在建立字符串集合到特征空间的映射。表示定义在上的实数空间，其每一维对应一个字符串，映射将字符串对应于空间的一个向量，其在维上的取值为

这里，是一个衰减参数，表示字符串i的长度，求和在中所有与相同的子串上进行。

例如，假设为英文字符集，n为3,为长度大于或等于3的字符串的集合。考虑将字符集映射到特征空间。的一维对应于字符串asd。这时，字符串“Nasdaq”与“lass das”在这一维上的值分别是和，（□为空格)。在第1个字符串里，asd/是连续的子串。在第2个字符串里，是长度为5的不连续子串，共出现2次。

两个字符串和上的字符串核函数是基于映射的特征空间中的内积：

字符串核函数给出了字符串和中长度等于n的所有子串组成的特征向量的余弦相似度（cosine similarity)。直观上，两个字符串相同的子串越多，它们就越相似，字符串核函数的值就越大。字符串核函数可以由动态规划快速地计算。

非线性支持向量分类机

如上所述，利用核技巧，可以将线性分类的学习方法应用到非线性分类问题中去。将线性支持向量机扩展到非线性支持向量机，只需将线性支持向量机对偶形式中的内积换成核函数。

定义(非线性支持向量机）从非线性分类训练集，通过核函数与软间隔最大化，或凸二次规划,学习得到的分类决策函数 

称为非线性支持向量，是正定核函数。

下面叙述非线性支持向量机学习算法。

算法(非线性支持向置机学习算法）

输入:线性可分训练集，其中

输出：分类决策函数。

(1)选取适当的核函数和适当的参数C,构造并求解最优化问题

求得最优解。

(2)选择的一个正分量,计算

(3)构造决策函数：

当是正定核函数时，问题是凸二次规划问题，解是存在的。

序列最小最优化算法

本节讨论支持向量机学习的实现问题。我们知道，支持向量机的学习问题可以形式化为求解凸二次规划问题。这样的凸二次规划问题具有全局最优解，并且有许多最优化算法可以用于这一问题的求解。但是当训练样本容量很大时，这些算法往往变得非常低效，以致无法使用。所以，如何高效地实现支持向量机学习就成为一个重要的问题。目前人们已提出许多快速实现算法。本节讲述其中的序列最小最优化（sequential minimal optimisation,SMO)算法，这种算法1998年由Platt提出。

SMO算法要解如下凸二次规划的对偶问题：

在这个问题中，变量是拉格朗日乘子，一个变量对应于一个样本点;变量的总数等于训练样本容量N。

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

注意，子问題的两个变童中只有一个是自由变量。假设、为两个变量，固定，那么由等式约束可知

如果确定，那么也随之确定。所以子问题中同时更新两个变量。

整个SMO算法包括两个部分：求解两个变量二次规划的解析方法和选择变量的启发式方法。

两个变量二次规划的求解方法

不失一般性，假设选择的两个变量是、，其他变量是固定的。于是SMO的最优化问题的子问题可以写成：

其中，，是常数，目标函数式中省略了不含、的常数项。

为了求解两个变量的二次规划问题,首先分析约束条件，然后在此约束条件下求极小。

由于只有两个变量,约束可以用二维空间中的图形表示（如图所示)。 

不等式约束使得在盒子内，等式约束使在平行于盒子的对角线的直线上（因为y不是+1就是-1呀）。因此要求的是目标函数在一条平行于对角线的线段上的最优值。这使得两个变量的最优化问题成为实质上的单变量的最优化问题，不妨考虑为变量的最优化问题。

假设问题的初始可行解为，最优解为并且假设在沿着约束方向未经剪辑时的最优解为。

由于需满足不等式约束,所以最优值的取值范围必须满足条件

其中，与是所在的对角线段端点的界。如果(如上图左图所示)，则

如果(如上图右图所示)，则

下面，首先求沿着约束方向未经剪辑即未考虑不等式约束时的最优解；然后再求剪辑后的解。我们用定理来叙述这个结果。为了叙述简单，记

令

当i = l，2时，为函数对输入的预测值与真实输出之差，

定理 最优化问题沿着约束方向未经剪辑时的解是

其中，

是输入空间到特征空间的映射，由式给出（读作Eta，待会儿在代码里会见到，先混个脸熟，详见：希腊字母读音表.pdf）。 经剪辑后的解是

由求得是

证明 引进记号

目标函数可写成

由及，可将表示为

代入式,得到只是的函数的目标函数：

对求导数

令其为0,得到

将代入，得到

将代入，于是得到

要使其满足不等式约束必须将其限制在区间内，从而得到的表达式。由等式约束,得到的表达式。于是得到最优化问题的解。

变量的选择方法

SMO算法在每个子问题中选择两个变量优化，其中至少一个变量是违反KKT条件的。

1.第1个变量的选择

SMO称选择第1个变量的过程为外层循环。外层循环在训练样本中选取违反KKT条件最严重的样本点，并将其对应的变量作为第1个变量。具体地，检验训练样本点是否满足KKT条件，即

其中，。

该检验是在范围内进行的。在检验过程中，外层循环首先遍历所有满足条件的样本点，即在间隔边界上的支持向量点，检验它们是否满足KKT条件。如果这些样本点都满足KKT条件，那么遍历整个训练集，检验它们是否满足KKT条件。

2.第2个变量的选择

SMO称选择第2个变量的过程为内层循环。假设在外层循环中已经找到第1个变量，现在要在内层循环中找第2个变量。第2个变量选择的标准是希望能使有足够大的变化。

由式和式可知，是依赖于的，为了加快计算速度，一种简单的做法是选择，使其对应的最大。因为已定，也确定了。如果是正的，那么选择最小的作为;如果是负的，那么选择最大的作为。为了节省计算时间，将所有值保存在一个列表中。

在特殊情况下，如果内层循环通过以上方法选择的不能使目标函数有足够的下降，那么采用以下启发式规则继续选择。遍历在间隔边界上的支持向量点，依次将其对应的变量作为试用，直到目标函数有足够的下降。若找不到合适的，那么遍历训练数据集；若仍找不到合适的，则放弃第1个，再通过外层循环寻求另外的。

3.计算阈值和差值

在每次完成两个变量的优化后，都要重新计算阈值。当时，由KKT条件可知：

于是，

由的定义式有

式的前两项可写成：

代入式,可得

同样，如果，那么，

如果同时满足条件,那么。如果是0或者C,那么和以及它们之间的数都是符合KKT条件的阈值，这时选择它们的中点作为。

在每次完成两个变量的优化之后，还必须更新对应的值，并将它们保存在列表中。值的更新要用到值，以及所有支持向量对应的:

其中，S是所有支持向量的集合。

SMO算法

算法(SMO算法）

输入：训练数据集，其中精度。

输出：近似解。

(1)取初始值；

(2)选取优化变量，解析求解两个变量的最优化问题,求得最优解,更新为

(3)若在精度范围内满足停机条件

其中，

则转(4>;否则令,转(2);

(4)取。

SMO算法实现

参考《机器学习实战》中的代码，单独开一章来剖析SMO算法的具体实现。老规矩，看代码之前先看数据。有什么样的食材做什么样的菜；有什么样的数据写什么样的算法。

训练集：testSet.txt

每个实例就是一个二维向量和一个±1的label：

  
  
3.542485	1.977398	-1
3.018896	2.556416	-1
7.551510	-1.580030	1
2.114999	-0.004466	-1
8.127113	1.274372	1
7.108772	-0.986906	1
8.610639	2.046708	1
2.326297	0.265213	-1
3.634009	1.730537	-1
0.341367	-0.894998	-1

加载代码如下：

  
  
def loadDataSet(fileName):
    """
    加载数据集
    :param fileName:
    :return:
    """
    dataMat = []; labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat

顺便介绍两个辅助的函数，其作用如注释所示：

  
  
def selectJrand(i,m):
    """
    随机从0到m挑选一个不等于i的数
    :param i:
    :param m:
    :return:
    """
    j=i             #排除i
    while (j==i):
        j = int(random.uniform(0,m))
    return j
 
def clipAlpha(aj,H,L):
    """
    将aj剪裁到L(ow)和H(igh)之间
    :param aj:
    :param H:
    :param L:
    :return:
    """
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

简化版SMO算法

《机器学习实战》称完整版SMO算法的内外两层循环比较复杂，所以简化为去掉外层循环，并且去掉了内存循环的启发式规则，以达到循序渐进便于理解的目的：

  
  
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    """
    简化版SMO算法
    :param dataMatIn:       X
    :param classLabels:     Y
    :param C:               惩罚参数
    :param toler:           容错率
    :param maxIter:         最大循环次数
    :return:
    """
    dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose()
    b = 0; m,n = shape(dataMatrix)  # m:=训练实例的个数；n:=每个实例的维度
    alphas = mat(zeros((m,1)))
    iter = 0
    while (iter < maxIter):
        alphaPairsChanged = 0   #alpha是否已经进行了优化
        for i in range(m):
            #   w = alpha * y * x;  f(x_i) = w^T * x_i + b
            fXi = float(multiply(alphas,labelMat).T*dataMatrix*dataMatrix[i,:].T) + b     # 预测的类别
            Ei = fXi - float(labelMat[i])   #得到误差，如果误差太大，检查是否可能被优化
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)): #必须满足约束
                j = selectJrand(i,m)
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy()                # 教材中的α_1^old和α_2^old
                if (labelMat[i] != labelMat[j]):                                          # 两者所在的对角线段端点的界
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                if L==H: print "L==H"; continue
                # Eta = -(2 * K12 - K11 - K22)，且Eta非负，此处eta = -Eta则非正
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                if eta >= 0: print "eta>=0"; continue
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                alphas[j] = clipAlpha(alphas[j],H,L)
                #如果内层循环通过以上方法选择的α_2不能使目标函数有足够的下降，那么放弃α_1
                if (abs(alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
                b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
                if (0 < alphas[i]) and (C > alphas[i]): b = b1
                elif (0 < alphas[j]) and (C > alphas[j]): b = b2
                else: b = (b1 + b2)/2.0
                alphaPairsChanged += 1
                print "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
        if (alphaPairsChanged == 0): iter += 1
        else: iter = 0
        print "iteration number: %d" % iter
    return b,alphas

函数有点长，将其拆解开来，对照教材逐步理解：

  
  
eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T

这一句是计算：

  
  
alphas[j] -= labelMat[j]*(Ei - Ej)/eta

这句是计算新的α2：

新的α2出来之后，需要按照来更新α1：

  
  
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])

接下来按照

计算b1和b2:

  
  
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

如果同时满足条件,那么。如果是0或者C,那么和以及它们之间的数都是符合KKT条件的阈值，这时选择它们的中点作为:

  
  
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

调用方法：

  
  
dataArr, labelArr = loadDataSet('testSet.txt')
b,alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)
print b

输出：

  
  
[[-3.884627]]

由于上文提到的随机函数，最终结果未必是这个。

还可以这样观察alpha数组：

  
  
print alphas[alphas > 0]

得到：

  
  
[[ 0.12971213  0.2385676   0.36827972]]

这实际表明数据集中只有3个支持向量：

  
  
for i in range(len(alphas)):
    if alphas[i] > 0.0 : print dataArr[i], labelArr[i]

输出：

  
  
[4.658191, 3.507396] -1.0
[3.457096, -0.082216] -1.0
[6.080573, 0.418886] 1.0

根据支持向量和数据集，我们基于公式可以得到分离超平面的w参数：

  
  
def calcWs(alphas,dataArr,classLabels):
    """
    根据支持向量计算分离超平面(w,b)的w参数
    :param alphas:拉格朗日乘子向量
    :param dataArr:数据集x
    :param classLabels:数据集y
    :return: w=∑alphas_i*y_i*x_i
    """
    X = mat(dataArr); labelMat = mat(classLabels).transpose()
    m,n = shape(X)
    w = zeros((n,1))
    for i in range(m):
        w += multiply(alphas[i]*labelMat[i],X[i,:].T)
    return w

写个可视化函数看看是什么效果：

  
  
def plotSVM(dataArr, labelArr, w, b, svList):
    xcord0 = []
    ycord0 = []
    xcord1 = []
    ycord1 = []
    markers =[]
    colors =[]
    for i in range(len(labelArr)):
        xPt = dataArr[i][0]
        yPt = dataArr[i][1]
        label = labelArr[i]
        if (label == -1):
            xcord0.append(xPt)
            ycord0.append(yPt)
        else:
            xcord1.append(xPt)
            ycord1.append(yPt)
 
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord0,ycord0, marker='s', s=90)
    ax.scatter(xcord1,ycord1, marker='o', s=50, c='red')
    plt.title('Support Vectors Circled')
    for sv in svList:
        circle = Circle((dataArr[sv][0], dataArr[sv][1]), 0.5, facecolor='none', edgecolor=(0,0.8,0.8), linewidth=3, alpha=0.5)
        ax.add_patch(circle)
 
    w0=w[0][0]; w1=w[1][0]; b=float(b)
    x = arange(-2.0, 12.0, 0.1)
    y = (-w0*x - b)/w1
    ax.plot(x,y)
    ax.axis([-2,12,-8,6])
    plt.show()

画到图上去：

  
  
dataArr, labelArr = loadDataSet('testSet.txt')
b,alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)
print b
w = calcWs(alphas, dataArr, labelArr)
print w
svList = []
for i in range(len(alphas)):
    if abs(alphas[i]) > 0.0000001 :
        print dataArr[i], labelArr[i]
        svList.append(i)
 
plotSVM(dataArr, labelArr, w, b, svList)

得到：

其他的点由于alpha为0，根本无法影响分类决策函数。

完整版SMO算法

简化版虽然代码很少，但是速度慢。本着循序渐进的原则，先不考虑核函数，直接看SMO的完整实现。

在《机器学习实战》的实现中，为了避免频繁填充函数参数，定义了一个中间结构：

  
  
class optStructK:
    def __init__(self,dataMatIn, classLabels, C, toler):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2)))

所有域都是自解释的，除了eCache。这个矩阵（或者说列表）是储存每个样本点的误差的，第1列是标志位，代表第二位的误差是否是有效的。

误差是怎么算的呢？误差是这么算的：

  
  
def calcEkK(oS, k):
    """
    计算第k个样本点的误差
    :param oS:
    :param k:
    :return:
    """
    fXk = float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T)) + oS.b
    Ek = fXk - float(oS.labelMat[k])
    return Ek

当改变了模型参数后，通过如下方法更新eCache中的误差：

  
  
def updateEkK(oS, k):
    Ek = calcEkK(oS, k)
    oS.eCache[k] = [1,Ek]

由于通过简化版，已经熟悉了内层循环，那么还是先看看内层循环。先把简化版中省略掉的启发式规则补充过来：

  
  
def selectJK(i, oS, Ei):
    """
    选择第2个变量的过程，即内层循环中的启发规则。选择标准是使alpha_2有足够大的变化。
    :param i:第一个变量
    :param oS:中间结果
    :param Ei:第i个点的误差
    :return:
    """
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]                           #设为有效
    # 找寻产生最大误差变化的alpha
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]  #有效位为1的误差列表
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:                   #遍历列表找出最大误差变化
            if k == i: continue                     #第二个alpha不应该等于第一个
            Ek = calcEkK(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej
    else:                                           #找不到，只好随机选择一个
        j = selectJrand(i, oS.m)
        Ej = calcEkK(oS, j)
    return j, Ej

在此基础上看内层循环，其实除了selectJ，其他都差不多：

  
  
def innerLK(i, oS):
    Ei = calcEkK(oS, i)
    #如果误差太大，且alpha满足约束，则尝试优化它
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        j,Ej = selectJK(i, oS, Ei)   #不再是 selectJrand 那种简化的选择方法
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy() # 教材中的α_1^old和α_2^old
        if (oS.labelMat[i] != oS.labelMat[j]):                           # 两者所在的对角线段端点的界
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L==H: print "L==H"; return 0
        eta = 2.0 * oS.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - oS.X[j,:]*oS.X[j,:].T
        if eta >= 0: print "eta>=0"; return 0
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEkK(oS, j)                                                 # 更新误差缓存
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])    #更新α_1
        updateEkK(oS, i)                                                 # # 更新误差缓存
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[j,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].T
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
        else: oS.b = (b1 + b2)/2.0
        return 1
    else: return 0

对于外层循环，这里的实现并没有在训练样本中选取违反KKT条件最严重的样本点，而是优先顺序遍历间隔边界上的支持向量点，若无法优化模型则遍历整个数据集。

  
  
def smoPK(dataMatIn, classLabels, C, toler, maxIter):
    """
    完整版的Platt SMO算法
    :param dataMatIn:
    :param classLabels:
    :param C:
    :param toler:
    :param maxIter:
    :return:
    """
    oS = optStructK(mat(dataMatIn),mat(classLabels).transpose(),C,toler)
    iter = 0
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:           #遍历所有值
            for i in range(oS.m):        
                alphaPairsChanged += innerLK(i,oS)
                print "fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
            iter += 1
        else:                   #遍历间隔边界上的支持向量点
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerLK(i,oS)
                print "non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
            iter += 1
        if entireSet: entireSet = False                 #翻转entireSet
        elif (alphaPairsChanged == 0): entireSet = True  
        print "iteration number: %d" % iter
    return oS.b,oS.alphas

调用方法：

  
  
dataArr, labelArr = loadDataSet('testSet.txt')
b,alphas = smoPK(dataArr, labelArr, 0.6, 0.001, 40)
print b
w = calcWs(alphas, dataArr, labelArr)
print w
svList = []
for i in range(len(alphas)):
    if abs(alphas[i]) > 0.0000001 :
        print dataArr[i], labelArr[i]
        svList.append(i)
 
plotSVM(dataArr, labelArr, w, b, svList)

输出：

直观的感觉是支持向量变多了，速度变快了。

非线性支持向量机SMO算法

看过上面两版的代码，的确发现少了核函数的踪影，碰到下图这种线性不可分的数据还是无能为力的：

那就来写一个核函数吧：

  
  
def kernelTrans(X, A, kTup):
    """
    核函数
    :param X: 全部向量
    :param A: 某个向量
    :param kTup: 核函数名称
    :return: :raise NameError:
    """
    m,n = shape(X)
    K = mat(zeros((m,1)))
    if kTup[0]=='lin': K = X * A.T  #线性核函数
    elif kTup[0]=='rbf':
        for j in range(m):
            deltaRow = X[j,:] - A   #高斯径向基核函数
            K[j] = deltaRow*deltaRow.T
        K = exp(K/(-1*kTup[1]**2))  #K依然是一个向量
    else: raise NameError('Houston We Have a Problem -- \
    That Kernel is not recognized')
    return K

当传入lin时，实际上就相当于不使用核函数，亦即“核函数”相当于线性支持向量机分类决策函数中的原始形式。

传入rbf时，使用的是《统计学习方法》中介绍的高斯核函数（Gaussian kernel function)：

这时，如果x和z很相近（），那么核函数值为1，如果x和z相差很大（），那么核函数值约等于0。由于这个函数类似于高斯分布，因此称为高斯核函数，也叫做径向基函数(Radial Basis Function 简称RBF)。它能够把原始特征映射到无穷维。

lin和rbf的区别如下图所示：

当传入预设之外的核函数名称时，这个函数将抛出异常：

Houston We Have a Problem

这是阿波罗13号中的一句名台词。

有了核函数，在构造中间结构的时候，可以一次性地将中的计算出来存着：

  
  
class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters 
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2)))                        #第一列是有效与否的标记
        self.K = mat(zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)    #软cache

之后每次计算决策函数的时候就用这个软cache了：

  
  
def calcEk(oS, k):
    """
    计算第k个样本点的误差
    :param oS:
    :param k:
    :return:
    """
    fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek

eta和b的计算也要把内积换成核函数：

  
  
def innerL(i, oS):
    Ei = calcEkK(oS, i)
    #如果误差太大，且alpha满足约束，则尝试优化它
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        j,Ej = selectJ(i, oS, Ei)   #不再是 selectJrand 那种简化的选择方法
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy() # 教材中的α_1^old和α_2^old
        if (oS.labelMat[i] != oS.labelMat[j]):                           # 两者所在的对角线段端点的界
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L==H: print "L==H"; return 0
        eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j]
        if eta >= 0: print "eta>=0"; return 0
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS, j)                                                 # 更新误差缓存
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])    #更新α_1
        updateEk(oS, i)                                                 # # 更新误差缓存
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[i,j]
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]- oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,j]
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
        else: oS.b = (b1 + b2)/2.0
        return 1
    else: return 0

调用方法：

  
  
def testRbf(k1=1.3):
    dataArr,labelArr = loadDataSet('testSetRBF.txt')
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1))
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]
    sVs=datMat[svInd]                                                   #获取支持向量
    labelSV = labelMat[svInd]
    print "there are %d Support Vectors" % shape(sVs)[0]
    m,n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print "the training error rate is: %f" % (float(errorCount)/m)
    dataArr,labelArr = loadDataSet('testSetRBF2.txt')
    errorCount = 0
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    m,n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1    
    print "the test error rate is: %f" % (float(errorCount)/m)

用到了训练集testSetRBF.txt和测试集testSetRBF2.txt。

输出：

  
  
there are 16 Support Vectors
the training error rate is: 0.460000
the test error rate is: 0.360000

上面特地留了个参数k1，指的是高斯径向基函数的，试试将这个参数减小一点，设为0.10试试：

  
  
there are 57 Support Vectors
the training error rate is: 0.280000
the test error rate is: 0.360000

发现支持向量数目增多了，同时错误率也下降了，但是训练速度显著减慢。

当减小的时候，每个支持向量对决策函数的影响减小，因此需要更多的支持向量。当非常小的时候，支持向量非常多，相当于利用整个数据集进行分类，这种分类方法被称为k近邻。

那到底该取多少比较好呢？没有定论，事实上，这个问题被称为调参。想当年，袁隆平在试验田里折腾多少年的杂交水稻，那也是调参，是个体力活，也是个技术活。

《机器学习实战》通过手写数字的识别数据集给出了一个调参的例子，其结果如下：

最好的参数是RBF,10；这说明支持向量并不是越多越好，支持向量越多，虽然训练错误率降低，但测试错误率升高，对未知数据的分类能力降低。

另一方面，还说明了线性核函数也能取得不错的性能，有时候，朴素的东西就是好的。

附录

凸二次规划（convex quadratic programming)问题

二次规划简介

二次规划（Quadratic Programming，简称QP)是最简单的约束非线性规划问题，指的是带有二次目标函数和线性约束的最优化问题。

二次规划问题可以表述成如下标准形式：

其中为n 阶实对称矩阵，A为维矩阵，c为n 维列向量，b为m 维列向量。

特别的，当正定（设M是n阶方阵，如果对任何非零向量z，都有zTMz> 0，其中zT表示z的转置，就称M是正定矩阵；如果zTMz≥ 0，就称M是半正定矩阵）时，目标函数为凸函数，线性约束下可行域又是凸集，问题称为凸二次规划。凸二次规划是一种最简单的非线性规划，且具有如下性质：

(1) K-T条件不仅是最优解的必要条件，而且是充分条件；

(2) 局部最优解就是全局最优解。

Reference：第16讲_二次规划.ppt

Gram矩阵

Gram方阵的定义

在实数域上的欧氏空间中,我们总可以定义向量的内积.设为维欧氏空间中的任意一组向量,用这组向量的一切可能的内积作成一个方阵,即

这样的方阵定义为向量组的Gram方阵,记为简记为.并称为的Gram行列式.

Reference：http://pan.baidu.com/s/1i3lOVfZ

希尔伯特空间

Reference：http://pan.baidu.com/s/1eQGxibo

赋范向量空间

Reference：http://pan.baidu.com/s/1mgB1pzE

核或正定核

设是中的一个子集，称定义在上的函数是核函数，如果存在一个从到Hilbert空间的映射

使得对任意的，

都成立。其中表示Hilbert空间中的内积。

Reference：http://pan.baidu.com/s/1qW9Xm4s

Reference

《统计学习方法》

《机器学习实战》

http://www.cnblogs.com/jerrylead/archive/2011/03/18/1988406.html

转载：码农场 » 支持向量机  http://www.hankcs.com/ml/support-vector-machine.html