December 19, 2017

#machine-learning
#notes
#support-vector-machine

Notes on Support Vector Machine

这篇文章是我在学习支持向量机时写的综述性质的笔记，简要介绍了支持向量机的理论推导以及如何实现（二次规划和 SMO 算法）。这篇文章综合了许多资料中的结果，在这里向这些资料的作者表示感谢，我在文章最后的参考一节列出来这些资料，并且注明了每篇资料所讲述的侧重点。如果看不懂我的文章，这些资料也是很好的入门点。

问题表述

支持向量机所要解决的问题表述起来非常简单：有一簇点被分为了互不相交的两类，我们要找一个界把这两类点恰好分开。在这一节，我们从二维形式中开始，把问题讲清楚。

二维空间中的表述

直线可以表示为 $y=kx+b$ ，然而这样是不够的，我们需要使用更一般的形式 $Ax+By+C=0$ ，其实就是把向量形式 $(A,B)\cdot(x,y)+C=0$ 拆开写。有了这个直线后，任意给空间中的一个点 $(x_0,y_0)$ ，把它带入直线方程，可以得到三种结果：

$Ax_0+By_0+C>0$ ：点在直线的某一侧。
$Ax_0+By_0+C=0$ ：点在直线上。
$Ax_0+By_0+C<0$ ：点在直线的另一侧。

不考虑在直线上的点，即第二种情况，我们的直线就是一个完美的分界标准，如果设 $f=Ax+By+C$ ，这样 $f$ 就是一个完美的分类器。

高维空间中的表述

好了，现在把它拓展到 $n$ 维空间上，我们就得改改叫法了——点还是原来的点，直线变成了超平面（hyper-plane）。直线的表示方法也不一样了，我们不把点的每一个分量都写出来，而是用向量的形式来表示： $\boldsymbol w^T\boldsymbol x+b=0$ ，其中 $\boldsymbol x$ 是点， $\boldsymbol w$ 是点的每个分量所乘的系数， $b$ 是一个数，它的作用类似于二维直线中截距的东西，没有它的话我们的超平面一定经过原点。

在本节（问题表述）中，我们混用「超平面」和「线」这两个术语，举例也在二维空间中举例（大于三维的我也举不出来啊）。

最优化目标

我们肯定不会找一条线就草草了事，我们要选一条最好的线，问题来了——什么叫最好的线？我们需要一个可以量化的标准，这时就很自然地想到了线的方程 $\boldsymbol w^T\boldsymbol x+b=0$ ，我们设 $f(\boldsymbol x)=\boldsymbol w^T\boldsymbol x+b$ ，我们可以很容易地验证，点到这条线的距离越远， $f$ 的值越大。但是用 $f$ 的值作为度量并不好，因为它是与 $\boldsymbol w$ 和 $b$ 的取值相关的，如果它们两个成比例增加， $f$ 的值会在直线不变的情况下增加。这并不好。但这个值也并非没有用，我们管它叫作函数间隔（functional margin）。

这时候就会考虑一下几何意义，我们可以选择点到直线的几何距离作为度量，设点为 $\boldsymbol x$ ，直线为 $\boldsymbol w^T\boldsymbol x+b=0$ ，点到直线的距离为 $\gamma$ ，令 $\boldsymbol x_0$ 为 $\boldsymbol x$ 在直线上的投影，我们可以写出它们之间的关系

\boldsymbol x-\boldsymbol x_0=\gamma\frac{\boldsymbol w}{\|\boldsymbol w\|}

要理解这个式子，首先要明白 $\boldsymbol w$ 是直线的法向量。现在对式子的两边都应用函数 $f$ ，会有

\begin{aligned} f(\boldsymbol x-\boldsymbol x_0)&=f\left(\gamma\frac{\boldsymbol w}{\|\boldsymbol w\|}\right)\\ \boldsymbol w^T\boldsymbol x-\boldsymbol w^T\boldsymbol x_0+b&=\boldsymbol w^T\gamma\frac{\boldsymbol w}{\|\boldsymbol w\|}+b\\ \boldsymbol w^T\boldsymbol x+b-\boldsymbol w^T\boldsymbol x_0-b&=\gamma\frac{\boldsymbol w^T\boldsymbol w}{\|\boldsymbol w\|}\\ \boldsymbol w^T\boldsymbol x+b-0&=\gamma\|\boldsymbol w\|\\ \frac{\boldsymbol w^T\boldsymbol x+b}{\|\boldsymbol w\|}&=\gamma \end{aligned}

这样我们就求出了这个距离 $\gamma$ ，很巧的是，它恰好是我们的函数间隔除以 $\|\boldsymbol w\|$ 。

To be done.

对偶问题到二次规划

我们的原始优化目标是

\min_\boldsymbol w \frac12 \|\boldsymbol w\|^2 \quad \text{subject to} \quad(\forall i)y_i(\boldsymbol w^T \boldsymbol x_i-b)\ge 1

但是解这个问题的话我们的算法就不能叫做支持向量了，因为解这个问题不会告诉我们有关支持向量的任何信息。这时候就要用到一个叫做拉格朗日对偶的东西，简而言之，就是把约束条件当作优化目标的一部分，把上面的优化目标里的函数改写成下面的形式

\mathscr L(\boldsymbol w, b, \boldsymbol \alpha)=\frac 12 \|\boldsymbol w\|^2-\sum_i\alpha_i\left[y_i\left(\boldsymbol w^T\boldsymbol x_i-b\right)-1\right] \quad (\forall i)\alpha_i\ge 0

注意！上面的式子并不是我们的优化目标，我们的优化目标实际上是

\min_{\boldsymbol w,b}\max_\boldsymbol\alpha\mathscr L(\boldsymbol w, b, \boldsymbol \alpha) \quad \text{subject to} \quad (\forall i) \alpha_i\ge 0

记住这一点，然后我们回来讨论上面的这个变换是什么意思。首先介绍这个 $\boldsymbol \alpha$ ，这些 $\alpha_i$ 有一个酷炫的名字「拉格朗日乘子」，它的目的是控制住原来的约束条件，使之保持成立。举个例子，假设存在某个 $i$ ，使得 $y_i(\boldsymbol w^T\boldsymbol x_i - b)<1$ ，我们的原约束条件肯定不会成立，而且在对偶问题的约束条件里，这个 $i$ 会使得 $\mathscr L(\boldsymbol w,b,\boldsymbol \alpha)$ 变成一个很大的值，这显然不是被优化目标中的 $\min$ 所取到，故优化目标实际上和原来的优化目标等价。

这个优化目标的另一个特性就是它满足 KKT 条件，故我们可以把这个问题变形成

\max_\boldsymbol\alpha\min_{\boldsymbol w,b}\mathscr L(\boldsymbol w, b, \boldsymbol \alpha) \quad \text{subject to} \quad (\forall i) \alpha_i\ge 0

的形式，也就是说，交换 $\min$ 和 $\max$ 的位置以后，问题仍然等价（注意只有满足 KKT 条件的问题才可以这样做）。上面这个式子非常方便我们求解，为什么呢？因为 $\min$ 被换到内部了，这样我们就可以求它们的极小值

\begin{aligned} \frac{\partial L}{\partial\boldsymbol w}=0&\Rightarrow\boldsymbol w=\sum_i \alpha_i y_i \boldsymbol x_i \\ \frac{\partial L}{\partial b}=0&\Rightarrow\sum_i\alpha_i y_i=0 \end{aligned}

这样我们就有了 $\boldsymbol w$ 和 $\boldsymbol\alpha$ 的关系了，而 $b$ 已经在我们的最优化目标中消失了。此时把 $\boldsymbol w$ 代回优化目标中。

\begin{aligned} \mathscr L(\boldsymbol w,b,\boldsymbol\alpha) &=\frac12\|\boldsymbol w\|^2-\sum_i\alpha_i[y_i(\boldsymbol w^T\boldsymbol x_i+b)-1]\\ &=\frac12\boldsymbol w^T\boldsymbol w-\sum_i\alpha_iy_i\boldsymbol w^T\boldsymbol x_i-\sum_i\alpha_i y_ib+\sum_i\alpha_i\\ &=\frac12\boldsymbol w^T\boldsymbol w-\boldsymbol w^T\sum_i\alpha_i y_i\boldsymbol x_i-0+\sum_i\alpha_i\\ &=\frac12\boldsymbol w^T\boldsymbol w-\boldsymbol w^T\boldsymbol w+\sum_i\alpha_i\\ &=-\frac12\boldsymbol w^T\boldsymbol w+\sum_i\alpha_i\\ &=-\frac12\left(\sum_i \alpha_i y_i \boldsymbol x_i\right)^T\left(\sum_i \alpha_i y_i \boldsymbol x_i\right)+\sum_i\alpha_i\\ &=-\frac12\left(\sum_i \alpha_i y_i \boldsymbol x_i^T\right)\left(\sum_i \alpha_i y_i \boldsymbol x_i\right)+\sum_i\alpha_i\\ &=-\frac12\sum_i\sum_j\alpha_i\alpha_jy_iy_j\boldsymbol x_i^T\boldsymbol x+\sum_i\alpha_i \end{aligned}

这样原来的优化目标 $\max_\boldsymbol\alpha\min_{\boldsymbol w,b}\mathscr L(\boldsymbol w,b,\boldsymbol \alpha)$ 就变成了

\max_\boldsymbol\alpha\left(-\frac12\sum_i\sum_j\alpha_i\alpha_jy_iy_j\boldsymbol x_i^T\boldsymbol x_j+\sum_i\alpha_i\right)

在前面加一个符号，把它变成一个最小化问题

\min_\boldsymbol\alpha\left(\frac12\sum_i\sum_j\alpha_i\alpha_jy_iy_j\boldsymbol x_i^T\boldsymbol x_j-\sum_i\alpha_i\right)

最后，别忘了限制条件

\min_\boldsymbol\alpha\left(\frac12\sum_i\sum_j\alpha_i\alpha_jy_iy_j\boldsymbol x_i^T\boldsymbol x_j-\sum_i\alpha_i\right) \\ \text{subject to} \quad (\forall i)\alpha_i\ge 0 \quad \text{and} \quad \sum_i\alpha_iy_i=0

这个问题是一个非常标准的二次规划（quadratic programming）。二次规划的问题描述如下

\min_\boldsymbol x \left( \frac12\boldsymbol x^TP\boldsymbol x+\boldsymbol q^T\boldsymbol x \right) \\ \text{subject to} \quad G\boldsymbol x\le\boldsymbol h \quad \text{and} \quad A\boldsymbol x=\boldsymbol b

我们设

$\boldsymbol x=\boldsymbol\alpha$
$P_{ij} = y_iy_j\boldsymbol x_i^T\boldsymbol x_j$
$\boldsymbol q=-\boldsymbol 1$
$G=-I$
$\boldsymbol h=\boldsymbol 0$
$A=\mathrm{diag}(\boldsymbol y)$
$\boldsymbol b=\boldsymbol 0$

就可以把我们的优化目标完美地套在了二次规划问题的模子上。有很多现成的库可以解二次规划问题，比如 Python 中的 CVXOPT，这个库在 Linux 和 macOS 上都可以直接使用 pip 安装，在 Windows 上可以使用这里的预编译包。代码如下：

def model(x, y):
    n_samples, n_features = x.shape
    yx = y[:, np.newaxis] * x
    P = cvxopt.matrix(np.dot(yx, yx.T))
    q = cvxopt.matrix(-np.ones(n_samples, 1))
    G = cvxopt.matrix(-np.eye(n_samples))
    h = cvxopt.matrix(np.zeros(n_samples))
    A = cvxopt.matrix(y.reshape(1, -1))
    b = cvxopt.matrix(np.zeros(1))
    solution = cvxopt.solvers.qp(P, q, G, h, A, b)
    return np.aligned(solution['x'])

这个函数中的 x 代表所有点的集合，y 代表值域为 $\{-1,1\}$ 的分类集合。

如果你不想用库，自己写一个也没什么问题，这篇文章简单介绍了如何手动求解二次规划问题。

解决噪音问题

我们现在的支持向量机解决不了什么实际问题，因为它能处理的情况实在是太理想了——数据必须线性可分而且不能容忍噪音。线性可分这个问题我们下一节再讲，这一节先来解决噪音的问题。

噪音会使我们现在的支持向量机直接崩溃，因为如果出现噪音我们的约束条件就不能被满足，二次规划问题也会变得不可解。怎么办呢？我们的目标实际上是「容忍」噪音，换句话说，我们允许有一定的噪音出现在数据集里。怎么去「容忍」呢？答案当然是修改我们的约束条件了。我们加入一个变量 $\boldsymbol\xi=(\xi_i)$ ，描述第 $i$ 个点最多偏离超平面的距离。原始问题就变成了

\min_{\boldsymbol w,\boldsymbol\xi}\left(\frac12\|\boldsymbol w\|^2+C\sum_i\xi_i\right) \\ \text{subject to} \quad (\forall i)y_i(\boldsymbol w^T\boldsymbol x+b)\ge 1-\xi_i

我们不能无限度地「容忍」，所以变量 $\boldsymbol\xi$ 就被加入了优化目标中，前面的常数 $C$ 起到了类似权重的作用，叫做惩罚参数或正则化参数。惩罚参数的值越小，代表我们对噪音的容忍度越小。

这样一来，模型就被改变了，故我们也需要重新推导其对偶问题了。首先用拉格朗日对偶重写，得到

\begin{aligned} \mathscr L(\boldsymbol w, b, \boldsymbol \xi,\boldsymbol \alpha, \boldsymbol \mu) &=\frac12\|\boldsymbol w\|^2+C\sum_i\xi_i\\ &-\sum_i\alpha_i[y_i(\boldsymbol w^T\boldsymbol x - b) - 1+\xi_i]-\sum_i\mu_i\xi_i \end{aligned}

我们现在要优化的目标函数就变成了

\min_{\boldsymbol w,b,\boldsymbol\xi}\max_{\alpha_i\ge 0,\mu\ge 0}\mathscr L(\boldsymbol w,b,\boldsymbol\xi,\boldsymbol\alpha,\boldsymbol\mu)

如我们所愿，这个问题也满足 KKT 条件（证明以后补充），还是同样的方法，它可以被转换成对偶问题

\max_{\alpha_i\ge 0,\mu\ge 0}\min_{\boldsymbol w,b,\boldsymbol\xi}\mathscr L(\boldsymbol w,b,\boldsymbol\xi,\boldsymbol\alpha,\boldsymbol\mu)

而保证问题的等价性。

\begin{aligned} \frac{\partial\mathscr L}{\partial\boldsymbol w}=0 &\Rightarrow\boldsymbol w=\sum_i\alpha_iy_i\boldsymbol x_i\\ \frac{\partial\mathscr L}{\partial b}=0 &\Rightarrow \sum_i\alpha_iy_i=0\\ \frac{\partial\mathscr L}{\partial\boldsymbol\xi}=0 &\Rightarrow C-\alpha_i-\mu_i=0 \end{aligned}

把得到的三个结果代回 $\mathscr L$ 中

\begin{aligned} &\mathscr L(\boldsymbol w, b, \boldsymbol \xi,\boldsymbol \alpha, \boldsymbol \mu)\\ &=\frac12\|\boldsymbol w\|^2+C\sum_i\xi_i-\sum_i\alpha_i[y_i(\boldsymbol w^T\boldsymbol x_i + b) - 1+\xi_i]-\sum_i\mu_i\xi_i\\ &=\frac12\boldsymbol w^T\boldsymbol w+C\sum_i\xi_i-\sum_i\alpha_iy_i\boldsymbol w^T\boldsymbol x_i-\sum_i\alpha_iy_ib+\sum_i\alpha_i-\sum_i\alpha_i\xi_i-\sum_i\mu_i\xi_i\\ &=\left(\frac12\boldsymbol w^T\boldsymbol w-\boldsymbol w^T\sum_i\alpha_iy_i\boldsymbol x_i\right)+\sum_i\alpha_iy_ib+\sum_i\xi_i\left(C-\alpha_i-\mu_i\right)+\sum_i\alpha_i\\ &=\left(\frac12\boldsymbol w^T\boldsymbol w-\boldsymbol w^T\boldsymbol w\right)+0+0+\sum_i\alpha_i\\ &=-\frac12\boldsymbol w^T\boldsymbol w+\sum_i\alpha_i\\ &=-\frac12\sum_i\sum_j\alpha_i\alpha_jy_iy_j\boldsymbol x_i\boldsymbol y_i+\sum_i\alpha_i \end{aligned}

我们会发现，得到的结果和不能容忍噪音时的是一样的。这样一来，不同的就只有约束条件了

\begin{aligned} \sum_i\alpha_iy_i&=0\\ C-\alpha_i-\mu_i&=0\\ \alpha_i&\ge0\\ \mu_i&\ge0 \end{aligned}

因为 $\mu_i$ 已经从我们的优化目标中消失了，并且 $C$ 是一个常量，所以我们可以把后三个条件合成一个 $0\le\alpha_i\le C$ 。这样对偶问题最后就可以被写成

\min_{\boldsymbol\alpha}\frac12\sum_i\sum_j\alpha_i\alpha_jy_iy_j\boldsymbol x_i\boldsymbol y_i-\sum_i\alpha_i\\ \text{subject to} \quad (\forall i)0\le\alpha_i\le C \quad \text{and} \quad (\forall i)\sum_i\alpha_iy_i=0

如果要用二次规划求解这个问题，我们需要把它映射到上一节提到的二次规划问题模型。 $P$ 和 $\boldsymbol q$ 不需要改变，需要改变的只有 $G$ 和 $\boldsymbol h$ ——因为我们加入了新的限制条件，即 $\alpha_i \le C$ 。

$G=\begin{pmatrix}I \\ -I\end{pmatrix}$
$\boldsymbol h=\begin{pmatrix}\boldsymbol 0 & \boldsymbol C\end{pmatrix}$

如果使用 CVXOPT 来解这个问题，只需要在上一节的代码上做一个小修改。

def model(x, y, C):
    n_samples, n_features = x.shape
    yx = y[:, np.newaxis] * x
    P = cvxopt.matrix(np.dot(yx, yx.T))
    q = cvxopt.matrix(-np.ones(n_samples, 1))
    G = cvxopt.matrix(np.vstack((-np.eye(n_samples), np.eye(n_samples))))
    h = cvxopt.matrix(np.vstack((np.zeros(n_samples), np.ones(n_samples) * C)))
    A = cvxopt.matrix(y.reshape(1, -1))
    b = cvxopt.matrix(np.zeros(1))
    solution = cvxopt.solvers.qp(P, q, G, h, A, b)
    return np.aligned(solution['x'])

从线性到非线性

上一节的开始我们提到我们的支持向量机无法处理线性不可分的情况，比如经典的同心圆数据集。我们从同心圆数据开始，讨论如何处理非线性可分的数据。

我们知道，在不变换坐标系到极坐标的情况下，我们是无法把二维空间中的椭圆用线性的形式表达出来的，我们只能用

Ax^2+Bxy+Cy^2+Dx+Ey+F=0

这样的二次方程去表示椭圆，这很明显不是线性的，但是如果我们设 $\boldsymbol x=(x,y,xy,x^2,y^2)$ ，即用额外的三个维度表示原来两个维度的平方和积的话，原来的方程就可以写作

Ax_4+Bx_3+Cx_5+Dx_1+Ex_2+F=0

了。这是一个非常完美的线性形式，在把数据用这种方式做过变换后，我们的数据就变成线性可分的（这里先不考虑噪音）了。

然而不要沉浸在喜悦中太久，问题的本质并没有变，我们的数据上并不是五维的，用增加维度的方式解决问题并不是什么好手段（别忘了隔壁做图像处理和神经网络的都在拼命想着怎么给数据降维呢）。如果我们原本的数据变成三维的，类似的处理就会把数据变成 19 维的（即 $T_3^1+T_3^2+T_3^3$ ，其中 $T_n^k$ 为可重组合数）。换句话说，就是维度灾难。

意识到这一点是不可行的，我们就要找一个妥协的办法，回头看我们的对偶问题的优化目标，我们就会发现，我们的数据只在做内积的时候被用到了，一个大胆的想法浮现出来：我们不必真刀真枪地把数据变换到高维空间上，而是只要在做内积的时候计算它们在高维空间的内积就行了。即用一个函数 $\kappa(\boldsymbol x_i,\boldsymbol x_j)$ 来代替内积 $\boldsymbol x_i^T\boldsymbol x_j$ 。我们既不需要更改模型，也不需要变换数据，就完成了线性到非线性的转换，并且之前提到的处理噪音的方法仍然适用。这一技巧被称之为「核技巧（kernel trick）」。

更快的算法

To be done.

参考资料

支持向量机系列文章 by pluskid——介绍支持向量机理论的系列文章，深入浅出。
支持向量机原理 by 刘建平 Pinard——介绍理论的系列文章，每个小的点都讲得很透，公式推导的细节也都写了出来。
Implementing and Visualizing SVM in Python with CVXOPT by Hardik Goel——讲了如何用 CVXOPT 实现一个最简单的支持向量机