支持向量机SVM

手机

iphone11大小尺寸是多少？苹果iPhone11和iPhone13的区别是什么？

警方通报辅警执法直播中被撞飞：犯罪嫌疑人已投案

支持向量机SVM

2023-03-20 09:52:32 来源：博客园

模型形式

SVM用于连续特征二分类。

(相关资料图)

在逻辑回归中，我们习惯在特征中加入常数 \(1\) 来代表偏差，从而将预测写作 \(g(\theta^Tx)\)。而在 SVM 中，我们习惯把偏差单独表示，即 \(g(w^Tx+b)\)。

此外，我们还要求用 \(\pm1\) 代表分类（在后续推导中会看到这样要求的便利之处），此时取 \(g(x)=\begin{cases}1&x\ge0\\-1&x<0\end{cases}\)。

逻辑回归，函数间隔与几何间隔

可以看出，SVM 和逻辑回归都是依靠线性分界来分类，并且，为了使预测更加可信，我们希望当 \(y^{(i)}=1\) 时，\(\theta^Tx^{(i)}\gg0\)（或者 \(w^Tx^{(i)}+b\gg0\)），反之 \(\theta^Tx^{(i)}\ll0\)（或者 \(w^Tx^{(i)}+b\ll0\)）。

但是线性回归的优化目标倾向于得到整体最优：

\[\max\prod_{i}P(y^{(i)}\mid x^{(i)})\]

这样可能导致会有一些点极其靠近分界面。

而 SVM 的优化目标是最大化数据点与分界面的最小距离。那么，我们需要量化定义距离——

函数间隔 functional margin

考虑如下定义距离的方式 \(\hat\gamma^{(i)}:=y^{(i)}(w^Tx^{(i)}+b)\)，其实这一定义就是整合了 \(y^{(i)}=\pm1\) 的两种情况，我们只需要使 \(\hat\gamma^{(i)}\gg0\) 即可。

但是这样定义存在问题。考虑到分界线就是 \(w^Tx+b=0\)，我们同倍数放缩 \(w,b\) 是不影响分界线的，但是 \(\hat\gamma\) 也会同比放缩。

我们可以加入一些限制，比如 \(|w|=1\)，这样就解决了同比放缩的问题。

几何间隔 geometric margin

顾名思义，\(\gamma^{(i)}\) 就是数据点 \(i\) 到分界面的欧氏距离，但是考虑到 \(y^{(i)}=\pm1\)，这个距离是有向的（可以从表达式中的 \(y^{(i)}\) 看出有向）。

写出表达式即为

\[\gamma^{(i)}=\frac{y^{(i)}(w^Tx^{(i)}+b)}{|w|}\]

看起来很熟悉？取 \((w,b)\leadsto\left(\frac{w}{|w|},\frac{b}{|w|}\right)\)，即为限制 \(|w|=1\) 的函数间隔。或者说：

\[\gamma^{(i)}=\frac{\hat\gamma^{(i)}}{|w|}\]

实际上，由于分界面是线性的，函数间隔也反应了几何上的相对关系，仅是相差一些倍数。

SVM优化目标

取 \(\gamma=\min\{\gamma^{(i)}\}\)，\(\hat{\gamma}=\min\{\hat\gamma^{(i)}\}\)，则 SVM 要求最大化 \(\gamma\)。可以如下形式化地表达：

\[\begin{aligned}\underset{\gamma,w,b}{\operatorname{maximize}}&\ \gamma\\\mathrm{s.t.}&\forall i, \frac{y^{(i)}(w^Tx^{(i)}+b)}{|w|}\ge\gamma\end{aligned}\]

注意到，分母上的 \(|w|\) 导致约束非线性，我们利用函数间隔与几何间隔的关系作如下变形：

\[\begin{aligned}\underset{\gamma,w,b}{\operatorname{maximize}}&\ \frac{\hat\gamma}{|w|}\\\mathrm{s.t.}&\forall i, y^{(i)}(w^Tx^{(i)}+b)\ge\hat\gamma\end{aligned}\]

仅仅是形式化的表达，我们并不能从上式马上得到我们想要的 \(w,b\)，因为符合条件的 \(w,b,\hat\gamma\) 是成比的一系列取值，我们只需要其中的任意一个。

实际上也很好解决，只要限制任意一个变量即可。仍然限制 \(|w|=1\)？我们并不希望这样做，因为这不是线性约束。考虑限制 \(\hat\gamma=1\)，也即

\[\begin{aligned}\underset{w,b}{\operatorname{maximize}}&\ \frac{1}{|w|}\\\mathrm{s.t.}&\forall i, y^{(i)}(w^Tx^{(i)}+b)\ge1\end{aligned}\]

可以看到，这样的确有唯一解，并且都是线性约束。但是这样做的本质含义是什么？——通过放缩 \((w,b)\)，使得最小的函数间隔取为 \(1\)，这也是一种得到唯一解的方法。

进一步，把优化目标转化为易求的值。

\[\begin{aligned}\underset{w,b}{\operatorname{minimize}}&\ \frac{1}{2}|w|^2\\\mathrm{s.t.}&\forall i, y^{(i)}(w^Tx^{(i)}+b)\ge1\end{aligned}\]

（注：从以上构造性推导无法看出 \(\frac{1}{2}|w|^2\) 如何得来，虽然显然正确，但是应该没有解释到本质。）

这就是一个二次规划问题，可以用求解。下面我们分析广泛的 KTT 条件下的非线性规划如何求解。

拉格朗日乘子法

目标

用于解决有等式约束的最值问题，也即求（\(\max\) 同理，其中的 \(w\) 是向量，代表多元函数）

\[\min\{f(w)\mid g_1(w)=0,\dots,g_m(w)=0\}\]

此处的 \(f\) 通常是连续，或者分块连续的，并且保证所求最值存在。

前置

一些在之后的推导会用到的结论，如果不会严格证明可以凭借几何直观理解。

函数在某点的等高线的法向向量与该点的梯度平行；
设 \(f(w)\) 在 \(g(w)=0\) 限制下，在 \(w_0\) 处取得最值，则等高线 \(f(w)=f(w_0)\) 与 \(g(w)=0\) 相切。

推导

当没有约束时，我们取得最值的点 \(w_0\) 满足 \(\nabla f(w_0)=0\)，于是可以解 \(\nabla f(w)=0\) 从而求得最值。

考虑有一个约束 \(g(w)=0\) 时，设 \(f(w)\) 在 \(w_0\) 取得最值，则有 \(\nabla f(w_0)\) 平行于 \(\nabla g(w_0)\)。也即 \(\nabla f(w_0)=-\lambda\nabla g(w_0)\)，其中的 \(\lambda\) 就被称作拉格朗日算子。

如下构造 \(L(w,\lambda)\)：

\[L(w,\lambda):=f(w)-\lambda g(w)\]

\(L(w,\lambda)\) 有如下有用的性质：

若 \(g(w)=0\)，则 \(L(w,\lambda)=f(w)\)，所以我们也可以是在 \(g(w)=0\) 限制下求 \(L(w,\lambda)\) 的最值。
\(\frac{\operatorname{d}f}{\operatorname{d}\lambda}=g\)，那么 \(g(w)=0\) 的限制自然写作 \(\frac{\operatorname{d}f(w,\lambda)}{\operatorname{d}\lambda}=0\)。
\(\nabla L(w,\lambda)=\nabla f(w)+\lambda\nabla g(w)\)，则根据拉格朗日乘子的定义，在最值点 \(w_0\) 有 \(\nabla L(w_0,\lambda)=0\)。

于是我们可以发现这样构造的巧妙之处——有限制的最值问题化为了无限制的求 \(L(w,\lambda)\) 的最值问题。

推广到多个限制，同样解如下函数的无限制最值（几何上不太直观，但还是用上述性质证）：

\[L(w,\lambda_1,\dots,\lambda_m):=f(w)+\lambda_1g_1(w)+\dots+\lambda_mg_m(w)\]

KTT 条件下非线性规划求解

考虑在拉格朗日乘子法所解决的问题中引入不等式限制，形式化地：

\[\min_w\{f(w)\mid g_i(w)=0,h_j(w)\le0\}\]

仍然考虑引入拉格朗日乘子：

\[L(w,\alpha_i,\beta_j)=f(w)+\sum_i\alpha_ig_i(w)+\sum_j\beta_jh_j(w)\]

我们希望引入拉格朗日乘子能够做到以下的操作：

如果最值点在不等式的边界取得，也即 \(h_i(w)=0\)，则 \(\beta_i\) 产生与一般的拉格朗日乘子相同的效果；
如果最值点在不等式的内部取得，也即 \(h_i(w)<0\)，则 \(\beta_i=0\)。

如果以上两点满足，则可以直接求 \(L\) 的最值代替求 \(f\) 的最值。但是很明显不满足，问题出在我们的第二个要求——如果我们直接求 \(L\) 的最大值，而有 \(h_j(w)<0\)，我们显然可以让 \(\beta_j\to-\infty\)，那么 \(L\) 无最值。

上述问题的根源是 \(\beta\) 能取负数，我们想要限制 \(\beta_j\ge0\)，这样就能解决这个问题。于是作如下构造：

\[\hat L(w)=\max_{\beta\ge0,\alpha}L(w,\alpha,\beta)\]

那么无论最值点落在何处，都有 \(\hat L(w)=f(w)\)，所以我们只需要求

\[\min_w\hat{L}(w)=\min_w\max_{\beta\ge0,\alpha}L(w,\alpha,\beta)\]

利用一些现有的结论（对偶问题的结论，我还没学:(），对于这个问题，

\[\begin{aligned}\min_w\hat{L}(w)&=\min_w\max_{\beta\ge0,\alpha}L(w,\alpha,\beta)\\&=\max_{\beta\ge0,\alpha}\min_wL(w,\alpha,\beta)\end{aligned}\]

求解 SVM 最优化问题

在上文，我们推导出了 SVM 最优化目标：

\[\begin{aligned}\underset{w,b}{\operatorname{minimize}}&\ \frac{1}{2}|w|^2\\\mathrm{s.t.}&\forall i, y^{(i)}(w^Tx^{(i)}+b)\ge1\end{aligned}\]

改写为标准的二次规划问题：

\[\begin{aligned}\underset{w,b}{\operatorname{minimize}}&\ \frac{1}{2}|w|^2\\\mathrm{s.t.}&\forall i, 1-y^{(i)}(w^Tx^{(i)}+b)\le0\end{aligned}\]

引入拉格朗日乘子：

\[\begin{aligned}L(w,b,\alpha)&=\frac{1}{2}|w|^2+\sum_{i}\alpha_i\left(1-y^{(i)}(w^Tx^{(i)}+b)\right)\\&=\frac{1}{2}w^Tw-w^T\sum_{i}\alpha_iy^{(i)}x^{(i)}+\sum_i\alpha_i(1-by^{(i)})\end{aligned}\]

在 \(\alpha\) 给定且 \(\alpha_i\ge0\) 的情况下，求偏导并置零：

\[\begin{aligned}\nabla L(w,b,\alpha)&=w-\sum_{i}\alpha_iy^{(i)}x^{(i)}=0\\\frac{\partial}{\partial b}L(w,b,\alpha)&=-\sum_i\alpha_iy^{(i)}=0\end{aligned}\]

于是得到

\[w=\sum_i\alpha_iy^{(i)}x^{(x)}\]

SVM 的预测值是 \(w^Tx+b\)，于是代入 \(w\) 可得

\[h(x)=\sum_{i}\alpha_iy^{(i)}x^{(i)T}x\]

并且当且仅当下式成立，原函数有最小值：

\[\sum_i\alpha_iy^{(i)}=0\]

代入 \(L(w,b,\alpha)\) 得

\[\begin{aligned}J(\alpha)&=-\frac{1}{2}\sum_{i}\alpha_iy^{(i)}x^{(i)T}\sum_{j}\alpha_jy^{(j)}x^{(j)}+\sum_i\alpha_i(1-by^{(i)})\\&=-\frac{1}{2}\sum_{i}\alpha_iy^{(i)}x^{(i)T}\sum_{j}\alpha_jy^{(j)}x^{(j)}+\sum_i\alpha_i-b\sum_i\alpha^{(i)}y^{(i)}\\&=\sum_i\alpha_i-\frac{1}{2}\sum_{ij}\alpha_i\alpha_jy^{(i)}y^{(j)}x^{(i)T}x^{(j)}\\&=\sum_i\alpha_i-\frac{1}{2}\sum_{ij}\alpha_i\alpha_jy^{(i)}y^{(j)}\langle x^{(i)},x^{(j)}\rangle\end{aligned}\]

最后一步引入内积记号，\(\langle v,u\rangle:=v^Tu\)，只是为了下一步的推导。

同样也将 SVM 的预测值写作内积形式

\[h(x)=\sum_{i}\alpha_iy^{(i)}\langle x^{(i)},x\rangle\]

接下来就是以 \(\alpha\) 为变量求最大值，即

\[\max_{\alpha}\{J(\alpha)\mid \sum_i\alpha_iy^{(i)}=0\}\]

暂时搁置这一问题，作为 SMO 算法章节的目标。

假设已经求得 \(\alpha\)，我们可以直接算出 \(w\)：

\[w=\sum_i\alpha_iy^{(i)}x^{(x)}\]

但是如何确定 \(b\)？此时需要回归最开始的最优化目标：

\[\begin{aligned}\underset{b}{\operatorname{maximize}}\ &\hat\gamma\\\mathrm{s.t.}&y^{(i)}(w^Tx^{(i)}+b)\ge\hat\gamma\end{aligned}\]

于是可以得到

\[b=-\frac{1}{2}\left(\max_{y^{(i)}=1}\{w^Tx^{(i)}\}+\min_{y^{(i)}=-1}\{w^Tx^{(i)}\}\right)\]

也即正负分分界的中点。

核函数 Kernel

如果仅仅是完成以上的最优化目标，SVM仍然仅仅做到了线性分类，与逻辑回归的唯一区别就是最优化的思想不同。

然而结合 Kernel，SVM 能够做到非线性分类。

回顾在逻辑回归中，我们如何实现一个非常局限的非线性分类——比如我们观察到分类边界大概是个二次函数，那么我们把特征手动扩展至平方。

可以用这样一个特征变换函数表示（以 \(3\) 个特征值为例）：

\[\begin{aligned}\phi(x)&=x^2\\\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}&\mapsto\begin{pmatrix}x_1x_1\\x_1x_2\\x_1x_3\\\dots\\x_3x_3\end{pmatrix}\end{aligned}\]

进行这样的坐标变换后，原来的最优化问题也要相应修改。之前用内积表达最优化目标的作用就体现在此处——SVM的最优化目标中所有涉及特征 \(x\) 的位置，都是一个内积，我们只需要把 \(\langle x^{(i)},x^{(j)}\rangle\) 更换为 \(\langle\phi(x^{(i)}),\phi(x^{(j)})\rangle\)。

但是我们很快发现，进行上述坐标变换后，计算代价非常大，计算一次内积由原来的 \(O(n)\) 变为了 \(O(n^2)\)。是否有方法加速计算？

我们考虑用核函数 \(K\) 表示特征变换后的内积。

\[K(u,v):=\langle\phi(u),\phi(v)\rangle\]

观察以下核函数，并推导其对应的坐标变换：

\[\begin{aligned}K(u,v)&=(u^Tv)^2\\&=\left(\sum_{i=1}^{n} u_iv_i\right)^2\\&=\sum_{i=1}^{n}\sum_{j=1}^{n} u_iu_jv_iv_j\\&=(u_1u_1,u_1u_2,\dots,u_nu_n)(v_1v_1,v_1v_2,\dots,v_nv_n)^T\end{aligned}\]

发现它对应的坐标变换就是之前提到的 \(\phi(x)=x^2\)。也就是说我们可以直接计算

\[\langle\phi(x^{(i)}),\phi(x^{(j)})\rangle=K(x^{(i)},x^{(j)})=(x^{(i)T}x^{(j)})^2\]

复杂度仍为 \(O(n)\)。

常用的核函数

考虑如下核函数：

\[K_1(u,v)=(u^Tv+c)^d\]

同样的推导，可以知道它对应的变换为

\[\begin{aligned}\phi(x)&=\begin{pmatrix}x^d\\\sqrt{cd}x^{d-1}\\\vdots\\\sqrt{c^{d-1}d}x\\\sqrt{c^d}\end{pmatrix}\end{aligned}\]

合理设置 \(c\) 的取值，我们可以得到一个 \(d\) 次的坐标变换，计算代价仍然是 \(O(n)\)（如果 \(d\) 是个小常数）。

再考虑如下核函数：

\[K_2(u,v)=\exp(-\frac{|u-v|^2}{2\sigma^2})\]

这个核函数称为高斯核函数，因为它看起来很像高斯分布的概率密度。用泰勒展开为无穷级数，可知 \(K_2\) 对应的特征变换是将原特征变换为无穷维特征，而计算代价仍然较小。这很明显地体现了核函数的强大。

有效的核函数

由 Mercer 定理，对于函数 \(K(u,v)\)，存在特征变换 \(\phi(x)\) 使得 \(K\) 是 \(\phi\) 的核函数的充要条件是 \(K\) 是对称半正定的。

在这里我们可以把 \(K\) 视作一个矩阵：\(K_{ij}=K(x^{(i)},x^{(j)})\)。显然 \(K\) 是对称的。并且由于内积的线性性：

\[K(u,v)=u^TKv\]

必要性容易证明：

\[\begin{aligned}K(z,z)&=z^TKz\\&=\sum_{i}\sum_{j}z_iK_{ij}z_j\\&=\sum_{i}\sum_{j}z_i\phi^T(x^{(i)})\phi(x^{(j)})z_j\\&=\sum_{i}\sum_{j}z_iz_j\sum_{k}\phi_k(x^{(i)})\phi_k(x^{(j)})\\&=\sum_{k}\sum_{i}z_i\phi_k(x^{(i)})\sum_{j}z_j\phi_k(x^{(j)})\\&=\sum_{k}\left(\sum_{i}z_i\phi_k(x^{(i)})\right)^2\end{aligned}\]

充分性略过。

这意味着我们没有必要找到原变换 \(\phi\) 来判断 \(K\) 是否是一个有意义的核函数，而可以直接验证 \(\forall i,K(x^{(i)},x^{(i)})\ge0\) 来判断。

容错——Regularization

回顾 SVM 的优化目标：

\[\begin{aligned}\underset{w,b}{\operatorname{minimize}}&\ \frac{1}{2}|w|^2\\\mathrm{s.t.}&\forall i, y^{(i)}(w^Tx^{(i)}+b)\ge1\end{aligned}\]

然而，如果数据本身不可分（可能是噪声），也即条件不能对所有 \(i\) 满足，我们仍然希望得到一个尽可能好的分界。

引入松弛变量 \(\xi_i\)，将限制条件改为

\[\forall i,y^{(i)}(w^Tx^{(i)}+b)\ge1-\xi_i,\xi_i\ge0\]

从而允许一些点不符合分界要求，但是不能无限制的松弛，考虑修改代价函数——松弛越多代价越高。

\[\begin{aligned}\underset{w,b}{\operatorname{minimize}}&\ \frac{1}{2}|w|^2+C\sum_i\xi_i\\\mathrm{s.t.}&\forall i, 1-\xi_i - y^{(i)}(w^Tx^{(i)}+b)\le0,\\&-\xi_i\le0\end{aligned}\]

其中 \(C\) 是一个设置的参数，\(C\) 越大限制越强。

于是相应的修改拉格朗日乘子法，再引入乘子 \(\xi\)：

\[\begin{aligned}\mathcal{L}(w,b,\xi,\alpha,\beta)&=\frac{1}{2}|w|^2+C\sum_i\xi_i\\&+\sum_i\alpha_i(1-\xi_i - y^{(i)}(w^Tx^{(i)}+b))\\&+\sum_i\beta_i(-\xi_i)\end{aligned}\]

于是我们要求（需要注意 \(w,b,\xi\) 都是算法的参数，而 \(\alpha,\beta\) 是乘子）：

\[\min_{w,b,\xi}\max_{\alpha\ge0,\beta\ge0}\big\{\mathcal{L}(w,b,\alpha,\beta)\big\}=\max_{\alpha\ge0,\beta\ge0}\min_{w,b}\big\{\mathcal{L}(w,b,\alpha,\beta)\big\}\]

求偏导并置零：

\[\begin{aligned}\frac{\partial\mathcal{L}}{\partial w}&=w-\sum_{i}\alpha_iy^{(i)}x^{(i)}\\\frac{\partial\mathcal{L}}{\partial b}&=-\sum_{i}\alpha_iy^{(i)}\\\frac{\partial\mathcal{L}}{\partial \xi_i}&=C-\alpha_i-\beta_i\Rightarrow\frac{\partial\mathcal{L}}{\partial \xi}&=C-\alpha-\beta\end{aligned}\]

于是乎，\(\mathcal{L}\) 对 \(w,b,\xi\) 有最小值当且仅当

\[\alpha^Ty=0,C-\alpha-\beta=0\]

考虑到 \(\alpha_i\ge0,\beta_i\ge0\)，上述第二条限制可以写为：

\[0\le\alpha_i\le C\]

此时 \(w=\sum_{i}\alpha_iy^{(i)}x^{(i)}\)，

\[\begin{aligned}\mathcal{L}(w,b,\xi,\alpha,\beta)&=-\frac{1}{2}w^Tw+\sum_{i}\alpha_i\\&=-\frac{1}{2}\sum_{i,j}\alpha_i\alpha_jy^{(i)}y^{(j)}\langle x^{(i)},x^{(j)}\rangle+\sum_{i}\alpha_i\end{aligned}\]

总结即为：

\[\max_{0\le\alpha_i\le C,\alpha^Ty=0}\left\{\sum_{i}\alpha_i-\frac{1}{2}\sum_{i,j}\alpha_i\alpha_jy^{(i)}y^{(j)}\langle x^{(i)},x^{(j)}\rangle\right\}\]

当如下 KKT 条件满足时，原问题得解：

\[\begin{aligned}\alpha_i=0&\Leftrightarrow y_ih(x_i)\le1\\\alpha_i=C&\Leftrightarrow y_ih(x_i)\ge1\\0<\alpha_i坐标上升法 & SMO

最后我们开始着手解决上述最优化问题。需要引入梯度下降以及牛顿法以外的另一个迭代求解的算法——坐标上升法，以及该方法应用到 SVM 上的一个优化——SMO。

坐标上升法

思路同样非常简单。设 \(f(\alpha_1,\dots,\alpha_n)\) 是 \(n\) 元函数，\(\alpha_i\) 的取值在实数域上没有任何限制。我们想要求 \(f\) 的最小值（以及对应的参数取值）。

我们每次选择一个参数 \(\alpha_k\)，视其他参数为常数，则 \(f\) 成为关于 \(\alpha_k\) 的一元函数，此时我们求使得 \(f\) 取得最小值的 \(\alpha_k"\)，令 \(\alpha_k\leadsto a_k"\)。形式化地：

\[\alpha_k:=\underset{t}{\operatorname{argmin}}f(\alpha_1,\dots,\alpha_{k-1},t,\alpha_{k+1},\dots,\alpha_n)\]

如果 \(f\) 是一个凸函数（对于最大值则是凹函数），则算法正确性是显然的。SVM 的优化目标是凸函数，由这一点可以看出之前取优化目标为 \(\frac{1}{2}|w|^2\) 的明智性。

此外，对于 SVM，此方法还需修改，因为 SVM 的 \(\alpha_1, \dots, \alpha_n\) 是相关的——

\[\sum_{i}\alpha_iy^{(i)}=0\]

考虑如下修改：每次取两个坐标 \(\alpha_i,\alpha_j\)，则 \(\alpha_j\) 是关于 \(\alpha_i\) 的函数，进而 \(f\) 是关于 \(\alpha_i\) 的函数。

下面具体描述如何用坐标上升法实现 SVM。

SMO

常规的坐标上升法，我们往往会选择依次优化每一个参数，迭代若干次直至收敛。

如果 \(n\) 很大（SVM 就是这种情况），那么这个算法迭代一次无疑非常低效。我们很自然地想要选择下降最快的一个坐标。

先研究如何优化 \(\alpha_i,\alpha_j\) 两个参数。

优化参数

考虑选定 \(\alpha_i,\alpha_j\) 两个参数进行优化，优化需要满足

\(\alpha_i,\alpha_j\in[0,C]\)；
\(y^{(i)}\alpha_i+y^{(j)}\alpha_j\) 恒定，因为其满足线性约束。

于是可以针对 \(y^{(i)},y^{(j)}\) 相等、不等两种情况分别计算出 \(\alpha_j\) 的取值范围，记为 \([L,H]\)，具体参考 Plott 的论文。

我们把优化目标函数记为 \(\Psi(\alpha)\)，即

\[\Psi(\alpha)=\sum_i\alpha_i-\frac{1}{2}\sum_{ij}\alpha_i\alpha_jy^{(i)}y^{(j)}K(x^{(i)},x^{(j)})\]

既然研究 \(\alpha_i,\alpha_j\) 的优化，就把 \(\Psi(\alpha)\) 整理成只与这两个变量相关的函数（用 \(K_{pq}\) 表示 \(K(x^{(p)},x^{(q)})\)）：

\[\begin{aligned}\Psi(\alpha)&=(\alpha_i+\alpha_j)-\frac{1}{2}\alpha_i^2K_{ii}-\frac{1}{2}\alpha_j^2K_{jj}-y^{(i)}y^{(j)}\alpha_i\alpha_jK_{ij}\\&-\alpha_iy^{(i)}\sum_{k\neq i,j}\alpha_ky^{(k)}K_{ik}-\alpha_jy^{(j)}\sum_{k\neq i,j}\alpha_ky^{(k)}K_{jk}\\&+Const\end{aligned}\]

我们知道 \(\alpha_i,\alpha_j\) 之间有线性限制，具体可以写为：

\[\alpha_i+y^{(i)}y^{(j)}\alpha_j=\lambda:=-y^{(i)}\sum_{k\neq i,j}y^{(k)}\alpha_k\]

设 \(s=y^{(i)}y^{(j)}\)，则有 \(\alpha_i=\lambda-s\alpha_j\)，\(\frac{\partial \alpha_i}{\partial \alpha_j}=-s\)。

于是

\[\begin{aligned}\frac{\partial}{\partial\alpha_j}\Psi(\alpha)&=1-s-(\lambda-s\alpha_j)(-s)K_{ii}-\alpha_jK_{jj}\\&-s(\lambda-s\alpha_j)K_{ij}+\alpha_jK_{ij}\\&+y^{(j)}(v_i-v_j)\end{aligned}\]

其中 \(v_i=\sum_{k\neq i,j}\alpha_ky^{(k)}K_{ik}\)，\(v_j\) 同理。

我们将与 \(\alpha_i,\alpha_j\) 有关的变量在修改之前的值用 \(\widetilde{x}\) 标注。注意到

\[\begin{aligned}\widetilde{u}_i:=h(x^{(i)})=\sum_{k}\alpha_ky^{(k)}K_{ik}=v_i+\widetilde{\alpha}_iy^{(i)}K_{ii}+\widetilde{\alpha}_jy^{(j)}K_{ij}\end{aligned}\]

则有

\[\begin{aligned}v_i=\widetilde{u}_i-\widetilde{\alpha}_iy^{(i)}K_{ii}-\widetilde{\alpha}_jy^{(j)}K_{ij}\\v_j=\widetilde{u}_j-\widetilde{\alpha}_iy^{(i)}K_{ij}-\widetilde{\alpha}_jy^{(j)}K_{jj}\end{aligned}\]

将以上两式以及 \(\lambda=y^{(i)}\widetilde{\alpha}_i+y^{(j)}\widetilde{\alpha}_j\) 代入 \(\frac{\partial}{\partial\alpha_j}\Psi(\alpha)\)（注：笔者并未成功化简出如下式子，结论参考自 Plott 的论文）：

\[\begin{aligned}\frac{\partial}{\partial\alpha_j}\Psi(\alpha)&=(\alpha_j-\widetilde\alpha_j)(K_{ii}+K_{jj}-2K_{ij})\\&-y^{(j)}(\widetilde u_i-y^{(i)}-\widetilde u_j+y^{(j)})\end{aligned}\]

则有

\[\hat\alpha_j=\widetilde{\alpha}_j+\frac{y^{j}(\widetilde E_i-\widetilde E_j)}{\eta}\]

其中 \(E_i=u_i-y^{(i)}\)，表示样本 \(i\) 的误差（预测-实际），这一变量需要在程序中维护（是SMO算法的重要优化，笔者尚未理解优化原因）。\(\eta=K_{ii}+K_{jj}-2K_{ij}\)，仅是方便表示。

为什么写作 \(\hat\alpha_j\) 而不是 \(\alpha_j\)？因为 \(\alpha_j\) 还必须满足 \(\alpha_j\in[L,H]\) 的限制，而我们现在相当于只是在整个平面上考虑 \(\Psi(\alpha)\) 的最值。所幸 \(\Psi\) 关于 \(\alpha\) 是二次的，所以如果 \(\alpha_j\) 越过了 \([L,H]\) 边界，则最值在边界取得：

\[\alpha_j=\begin{cases}L&\hat\alpha_j\le L\\\hat\alpha_j&L<\hat\alpha_j更新其他值：

\[\begin{aligned}\alpha_i&=\lambda-s\alpha_j\\E_k&=\widetilde E_k+(\alpha_i-\widetilde\alpha_i)y^{(i)}K_{ik}+(\alpha_j-\widetilde\alpha_j)y^{(j)}K_{jk}+b-\widetilde b\end{aligned}\]

注意到偏差 \(b\) 也会随着 \(\alpha_i,\alpha_j\) 的变化而更新。下面研究怎样计算 \(b\)。我们需要调整 \(b\) 使得 \(\alpha_i\) 或 \(\alpha_j\) 对应的 KKT 条件满足：

若 \(0<\alpha_i\[\begin{aligned}u_i&=\widetilde{u}_i+(\alpha_i-\widetilde{\alpha}_i)y^{(i)}K_{ii}+(\alpha_j-\widetilde{\alpha}_j)y^{(j)}K_{ij}+b-\widetilde{b}\end{aligned}\]

则

\[\begin{aligned}b_1&=\widetilde{b}-(\alpha_i-\widetilde\alpha_i)y^{(i)}K_{ii}-(\alpha_j-\widetilde{\alpha}_j)y^{(j)}K_{ij}+y^{(i)}-\widetilde{u}_i\\&=\widetilde{b}-(\alpha_i-\widetilde\alpha_i)y^{(i)}K_{ii}-(\alpha_j-\widetilde{\alpha}_j)y^{(j)}K_{ij}-\widetilde E_i\end{aligned}\]

若 \(0<\alpha_j\[b_2=\widetilde{b}-(\alpha_i-\widetilde\alpha_i)y^{(i)}K_{ij}-(\alpha_j-\widetilde{\alpha}_j)y^{(j)}K_{jj}-\widetilde E_j\]

于是我们这样更新 \(b\)：

如果 \(\alpha_i\) 不在边界（\(0<\alpha_i
如果 \(\alpha_i\) 在边界而 \(\alpha_j\) 不在边界，则取 \(b=b_2\)；
如果 \(\alpha_i,\alpha_j\) 都不在边界，则取 \(b=\frac{1}{2}(b_1+b_2)\)。

选取优化参数

下面描述 SMO 的启发式流程，笔者尚未理解其优化原理。

外层循环（Plott 论文中的 outer loop）用于选取第一个参数 \(\alpha_j\)。

首先枚举每一个 \(\alpha\)，如果其不满足 KKT 约束，并且可以优化，则选取其为 \(\alpha_j\) 并进一步选取 \(\alpha_i\) 进行优化。完成这一轮循环后，如果没有找到可以优化的 \(\alpha_j\)，则全局优化结束，否则，下一轮循环只考虑不在边界的 \(\alpha\) 作为 \(\alpha_j\)。
枚举所有不在边界的 \(\alpha\) 为 \(\alpha_j\)，如果可以优化，则进一步寻找 \(\alpha_i\)。若这一轮循环没有找到可以优化的 \(\alpha_j\)，则下一轮循环枚举所有 \(\alpha\)（步骤 1）；否则，下一轮循环仍然只考虑不在边界的 \(\alpha\)。

内层循环用于选取第二个参数 \(\alpha_i\)。

首先尝试最大化 \(|E_i-E_j|\)——如果 \(E_i>0\)，则尝试找到最小的 \(E_j\)，否则尝试找到最大的 \(E_j\)。
如果步骤 1 找不到合法的 \(E_j\)，则枚举所有不在边界的 \(\alpha_i\) 尝试优化。注意，为了防止总是用较小的 \(i\) 进行优化而导致模型偏向靠前的数据，枚举的起点是随机选取的。
如果步骤 2 也没有找到可以优化的 \(\alpha_i\)，再枚举所有的 \(\alpha_i\) 尝试优化，同样要随机选取起点。

参考代码

**Show Code**

class SVM:    def __init__(self, c, kernel, if_linear_kernel, error_tol=0.01,                 answer_eps=0.01, max_loop=2000):        self.alpha = None        self.w = None        self.b = None        self.error = None        self.m = None        self.n = None        self.xs = None        self.ys = None        self.c = c        self.error_tol = error_tol        self.eps = answer_eps        self.kernel = kernel        self.if_linear_kernel = if_linear_kernel        self.max_loop = max_loopdef decision_function_output(self, i):    if self.if_linear_kernel:        return np.dot(self.w, self.xs[i]) + self.b    else:        return np.sum(            [self.alpha[j] * self.ys[j]             * self.kernel(self.xs[i], self.xs[j]) for j in range(self.m)])\            + self.bdef take_step(self, i1, i2):    if i1 == i2:        return False    alpha1, alpha2 = self.alpha[i1], self.alpha[i2]    y1, y2 = self.ys[i1], self.ys[i2]    e1, e2 = self.get_error(i1), self.get_error(i2)    s = y1 * y2    if y1 != y2:        l = max(0, alpha2 - alpha1)        h = min(self.c, self.c + alpha2 - alpha1)    else:        l = max(0, alpha1 + alpha2 - self.c)        h = min(self.c, alpha1 + alpha2)    if l == h:        return False    k11 = self.kernel(self.xs[i1], self.xs[i1])    k12 = self.kernel(self.xs[i1], self.xs[i2])    k22 = self.kernel(self.xs[i2], self.xs[i2])    eta = k11 + k22 - 2 * k12    if eta > 0:        new_alpha2 = alpha2 + y2 * (e1 - e2) / eta        if new_alpha2 < l:            new_alpha2 = l        elif new_alpha2 > h:            new_alpha2 = h    else:        f1 = y1 * (e1 - self.b) - alpha1 * k11 - s * alpha2 * k12        f2 = y2 * (e2 - self.b) - s * alpha1 * k12 - alpha2 * k22        l1 = alpha1 + s * (alpha2 - l)        h1 = alpha1 + s * (alpha2 - h)        lobj = l1 * f1 + l * f2 + 0.5 * (h1 * h1) * k11\               + 0.5 * (l * l) * k22 + s * l * l1 * k12        hobj = h1 * f1 + h * f2 + 0.5 * (h1 * h1) * k11\               + 0.5 * (h * h) * k22 + s * h * h1 * k12        if lobj < hobj - self.eps:            new_alpha2 = l        elif lobj > hobj + self.eps:            new_alpha2 = h        else:            new_alpha2 = alpha2    if new_alpha2 < 1e-8:        new_alpha2 = 0.0    elif new_alpha2 > self.c - 1e-8:        new_alpha2 = self.c    if np.abs(new_alpha2 - alpha2)\            < self.eps * (new_alpha2 + alpha2 + self.eps):        return False    new_alpha1 = alpha1 + s * (alpha2 - new_alpha2)    b1 = self.b - e1 - y1 * (new_alpha1 - alpha1) * k11 \         - y2 * (new_alpha2 - alpha2) * k12    b2 = self.b - e2 - y1 * (new_alpha1 - alpha1) * k12 \         - y2 * (new_alpha2 - alpha2) * k22    origin_b = self.b    if 0 < new_alpha1 < self.c:        self.b = b1    elif 0 < new_alpha2 < self.c:        self.b = b2    else:        self.b = (b1 + b2) / 2    if self.if_linear_kernel:        self.w = self.w + y1 * (new_alpha1 - alpha1) * self.xs[i1] \                 + y2 * (new_alpha2 - alpha2) * self.xs[i2]    self.alpha[i1] = new_alpha1    self.alpha[i2] = new_alpha2    self.error[i1] = 0    self.error[i2] = 0    for i in range(self.m):        if 0 < self.alpha[i] < self.c:            self.error[i] += y1 * (new_alpha1 - alpha1) * self.kernel(self.xs[i1], self.xs[i]) + \                             y2 * (new_alpha2 - alpha2) * self.kernel(self.xs[i2], self.xs[i]) + \                             self.b - origin_b    return Truedef get_error(self, i1):    if 0 < self.alpha[i1] < self.c:        return self.error[i1]    else:        return self.decision_function_output(i1) - self.ys[i1]def exam_example(self, i2):    y2 = self.ys[i2]    alpha2 = self.alpha[i2]    e2 = self.get_error(i2)    r2 = e2 * y2    if (r2 < -self.error_tol and alpha2 < self.c) or (r2 > self.error_tol and alpha2 > 0):        if len(self.alpha[(self.alpha != 0) & (self.alpha != self.c)]):            if self.error[i2] > 0:                i1 = np.argmin(self.error)            else:                i1 = np.argmax(self.error)            if self.take_step(i1, i2):                return 1        # Randomly select start point        shift = np.random.choice(np.arange(self.m))        for i1 in np.roll(np.where((self.alpha != 0) & (self.alpha != self.c))[0], shift):            if self.take_step(i1, i2):                return 1        shift = np.random.choice(np.arange(self.m))        for i1 in np.roll(np.arange(self.m), shift):            if self.take_step(i1, i2):                return 1    return 0def outer_loop(self):    exam_all = True    cnt_changed = 0    cnt_loop = 0    while cnt_changed > 0 or exam_all:        cnt_changed = 0        if cnt_loop > self.max_loop:            break        cnt_loop += 1        if exam_all:            for i in range(self.m):                cnt_changed += self.exam_example(i)            exam_all = False        else:            # alpha 是列向量，因此只取行坐标            for i in np.where((self.alpha != 0) & (self.alpha != self.c))[0]:                cnt_changed += self.exam_example(i)            if cnt_changed == 0:                exam_all = Truedef fit(self, xs, ys):    self.xs, self.ys = xs, ys    self.m, self.n = self.xs.shape    self.alpha = np.zeros(self.m)    self.b = 0    self.error = np.dot((self.alpha * ys), self.kernel(xs, xs)) + self.b    if self.if_linear_kernel:        self.w = np.zeros(self.n)    self.outer_loop()def predict(self, xs):    ys = []    for x in xs:        eva = np.sum(            [self.alpha[j] * self.ys[j] * self.kernel(x, self.xs[j])             for j in range(self.m)]        ) + self.b        if eva > 0:            ys.append(1)        else:            ys.append(-1)    return np.array(ys)

关键词：