支持向量机SVM
模型形式
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\)）。
但是线性回归的优化目标倾向于得到整体最优：
这样可能导致会有一些点极其靠近分界面。
而 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)}\) 看出有向）。
写出表达式即为
看起来很熟悉？取 \((w,b)\leadsto\left(\frac{w}{|w|},\frac{b}{|w|}\right)\)，即为限制 \(|w|=1\) 的函数间隔。或者说：
实际上，由于分界面是线性的，函数间隔也反应了几何上的相对关系，仅是相差一些倍数。
SVM优化目标
取 \(\gamma=\min\{\gamma^{(i)}\}\)，\(\hat{\gamma}=\min\{\hat\gamma^{(i)}\}\)，则 SVM 要求最大化 \(\gamma\)。可以如下形式化地表达：
注意到，分母上的 \(|w|\) 导致约束非线性，我们利用函数间隔与几何间隔的关系作如下变形：
仅仅是形式化的表达，我们并不能从上式马上得到我们想要的 \(w,b\)，因为符合条件的 \(w,b,\hat\gamma\) 是成比的一系列取值，我们只需要其中的任意一个。
实际上也很好解决，只要限制任意一个变量即可。仍然限制 \(|w|=1\)？我们并不希望这样做，因为这不是线性约束。考虑限制 \(\hat\gamma=1\)，也即
可以看到，这样的确有唯一解，并且都是线性约束。但是这样做的本质含义是什么？——通过放缩 \((w,b)\)，使得最小的函数间隔取为 \(1\)，这也是一种得到唯一解的方法。
进一步，把优化目标转化为易求的值。
（注：从以上构造性推导无法看出 \(\frac{1}{2}|w|^2\) 如何得来，虽然显然正确，但是应该没有解释到本质。）
这就是一个二次规划问题，可以用求解。下面我们分析广泛的 KTT 条件下的非线性规划如何求解。
拉格朗日乘子法
目标
用于解决有等式约束的最值问题，也即求（\(\max\) 同理，其中的 \(w\) 是向量，代表多元函数）
此处的 \(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)\) 有如下有用的性质：
- 若 \(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)\) 的最值问题。
推广到多个限制，同样解如下函数的无限制最值（几何上不太直观，但还是用上述性质证）：
KTT 条件下非线性规划求解
考虑在拉格朗日乘子法所解决的问题中引入不等式限制，形式化地：
仍然考虑引入拉格朗日乘子：
我们希望引入拉格朗日乘子能够做到以下的操作：
- 如果最值点在不等式的边界取得，也即 \(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)=f(w)\)，所以我们只需要求
利用一些现有的结论（对偶问题的结论，我还没学:(），对于这个问题，
求解 SVM 最优化问题
在上文，我们推导出了 SVM 最优化目标：
改写为标准的二次规划问题：
引入拉格朗日乘子：
在 \(\alpha\) 给定且 \(\alpha_i\ge0\) 的情况下，求偏导并置零：
于是得到
SVM 的预测值是 \(w^Tx+b\)，于是代入 \(w\) 可得
并且当且仅当下式成立，原函数有最小值：
代入 \(L(w,b,\alpha)\) 得
最后一步引入内积记号，\(\langle v,u\rangle:=v^Tu\)，只是为了下一步的推导。
同样也将 SVM 的预测值写作内积形式
接下来就是以 \(\alpha\) 为变量求最大值，即
暂时搁置这一问题，作为 SMO 算法章节的目标。
假设已经求得 \(\alpha\)，我们可以直接算出 \(w\)：
但是如何确定 \(b\)？此时需要回归最开始的最优化目标：
于是可以得到
也即正负分分界的中点。
核函数 Kernel
如果仅仅是完成以上的最优化目标，SVM仍然仅仅做到了线性分类，与逻辑回归的唯一区别就是最优化的思想不同。
然而结合 Kernel，SVM 能够做到非线性分类。
回顾在逻辑回归中，我们如何实现一个非常局限的非线性分类——比如我们观察到分类边界大概是个二次函数，那么我们把特征手动扩展至平方。
可以用这样一个特征变换函数表示（以 \(3\) 个特征值为例）：
进行这样的坐标变换后，原来的最优化问题也要相应修改。之前用内积表达最优化目标的作用就体现在此处——SVM的最优化目标中所有涉及特征 \(x\) 的位置，都是一个内积，我们只需要把 \(\langle x^{(i)},x^{(j)}\rangle\) 更换为 \(\langle\phi(x^{(i)}),\phi(x^{(j)})\rangle\)。
但是我们很快发现，进行上述坐标变换后，计算代价非常大，计算一次内积由原来的 \(O(n)\) 变为了 \(O(n^2)\)。是否有方法加速计算？
我们考虑用核函数 \(K\) 表示特征变换后的内积。
观察以下核函数，并推导其对应的坐标变换：
发现它对应的坐标变换就是之前提到的 \(\phi(x)=x^2\)。也就是说我们可以直接计算
复杂度仍为 \(O(n)\)。
常用的核函数
考虑如下核函数：
同样的推导，可以知道它对应的变换为
合理设置 \(c\) 的取值，我们可以得到一个 \(d\) 次的坐标变换，计算代价仍然是 \(O(n)\)（如果 \(d\) 是个小常数）。
再考虑如下核函数：
这个核函数称为高斯核函数，因为它看起来很像高斯分布的概率密度。用泰勒展开为无穷级数，可知 \(K_2\) 对应的特征变换是将原特征变换为无穷维特征，而计算代价仍然较小。这很明显地体现了核函数的强大。
有效的核函数
由 Mercer 定理，对于函数 \(K(u,v)\)，存在特征变换 \(\phi(x)\) 使得 \(K\) 是 \(\phi\) 的核函数的充要条件是 \(K\) 是对称半正定的。
在这里我们可以把 \(K\) 视作一个矩阵：\(K_{ij}=K(x^{(i)},x^{(j)})\)。显然 \(K\) 是对称的。并且由于内积的线性性：
必要性容易证明：
充分性略过。
这意味着我们没有必要找到原变换 \(\phi\) 来判断 \(K\) 是否是一个有意义的核函数，而可以直接验证 \(\forall i,K(x^{(i)},x^{(i)})\ge0\) 来判断。
容错——Regularization
回顾 SVM 的优化目标：
然而，如果数据本身不可分（可能是噪声），也即条件不能对所有 \(i\) 满足，我们仍然希望得到一个尽可能好的分界。
引入松弛变量 \(\xi_i\)，将限制条件改为
从而允许一些点不符合分界要求，但是不能无限制的松弛，考虑修改代价函数——松弛越多代价越高。
其中 \(C\) 是一个设置的参数，\(C\) 越大限制越强。
于是相应的修改拉格朗日乘子法，再引入乘子 \(\xi\)：
于是我们要求（需要注意 \(w,b,\xi\) 都是算法的参数，而 \(\alpha,\beta\) 是乘子）：
求偏导并置零：
于是乎，\(\mathcal{L}\) 对 \(w,b,\xi\) 有最小值当且仅当
考虑到 \(\alpha_i\ge0,\beta_i\ge0\)，上述第二条限制可以写为：
此时 \(w=\sum_{i}\alpha_iy^{(i)}x^{(i)}\)，
总结即为：
当如下 KKT 条件满足时，原问题得解：
坐标上升法 & 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'\)。形式化地：
如果 \(f\) 是一个凸函数（对于最大值则是凹函数），则算法正确性是显然的。SVM 的优化目标是凸函数，由这一点可以看出之前取优化目标为 \(\frac{1}{2}|w|^2\) 的明智性。
此外，对于 SVM，此方法还需修改，因为 SVM 的 \(\alpha_1, \dots, \alpha_n\) 是相关的——
考虑如下修改：每次取两个坐标 \(\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)\)，即
既然研究 \(\alpha_i,\alpha_j\) 的优化，就把 \(\Psi(\alpha)\) 整理成只与这两个变量相关的函数（用 \(K_{pq}\) 表示 \(K(x^{(p)},x^{(q)})\)）：
我们知道 \(\alpha_i,\alpha_j\) 之间有线性限制，具体可以写为：
设 \(s=y^{(i)}y^{(j)}\)，则有 \(\alpha_i=\lambda-s\alpha_j\)，\(\frac{\partial \alpha_i}{\partial \alpha_j}=-s\)。
于是
其中 \(v_i=\sum_{k\neq i,j}\alpha_ky^{(k)}K_{ik}\)，\(v_j\) 同理。
我们将与 \(\alpha_i,\alpha_j\) 有关的变量在修改之前的值用 \(\widetilde{x}\) 标注。注意到
则有
将以上两式以及 \(\lambda=y^{(i)}\widetilde{\alpha}_i+y^{(j)}\widetilde{\alpha}_j\) 代入 \(\frac{\partial}{\partial\alpha_j}\Psi(\alpha)\)（注：笔者并未成功化简出如下式子，结论参考自 Plott 的论文）：
则有
其中 \(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]\) 边界，则最值在边界取得：
更新其他值：
注意到偏差 \(b\) 也会随着 \(\alpha_i,\alpha_j\) 的变化而更新。下面研究怎样计算 \(b\)。我们需要调整 \(b\) 使得 \(\alpha_i\) 或 \(\alpha_j\) 对应的 KKT 条件满足：
若 \(0<\alpha_i<C\)，要使 \(\alpha_i\) 的 KKT 条件满足，则 \(y^{(i)}u_i=1\)，即 \(u_i=y^{(i)}\)。又
则
若 \(0<\alpha_j<C\)，要使 \(\alpha_j\) 满足 KKT 限制，同理有
于是我们这样更新 \(b\)：
- 如果 \(\alpha_i\) 不在边界（\(0<\alpha_i<C\)），则取 \(b=b_1\)；
- 如果 \(\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_loop
def 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.b
def 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 True
def 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 0
def 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 = True
def 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)
