学习支持向量机时参考了很多大神的博客,和经典著作,从公式推导到代码实现,亲历亲为。
遇到很多疑惑,也是各种百度+Google,虽然最后也差不多都解决了,但终归是因为数学基础不扎实。
在这里做一个总结,以便以后复习。
有些地方加入了我自己的理解,仅作参考。
当然,如有疑问,请各位不吝赐教。
线性不可分的情况和核函数的部分本文没有涉及。(主要是作者本人也不是很懂。。。
参考文献及完整代码在文末。
码字不易,不如点赞支持一下这个推公式推到头秃的博主?
文章目录
1. SVM解决的问题
SVM解决的问题是经典的二元分类问题。给出一个分类标准使得样本集可以被最好地分类。在样本特征是二维的情况下,可以用下图表示
上图是我用Python实现的SVM分类器的最终结果图。
中间的实线是我们最终需要的分割线,在三维及以上的情况下叫做划分超平面。
画圈的样本是距离分割线最近的样本,叫做支持向量,这也是支持向量机名字的由来。两条黄色的虚线之间的距离叫做间隔,显然,这个间隔越大,也就代表两边的样本离分割线越远,我们得到的分割线就越鲁棒。
SVM的最终目的就是求出分割线的所有参数,在这个过程中还要确定样本中的支持向量是哪些。
2. 目标函数
我们的目的是找到使间隔最大的支持向量和分割平面,那么下面就要找到间隔的数学表达,也就是目标函数。
首先,样本集表示为:
D
=
{
(
x
1
,
y
i
)
,
(
x
2
,
y
2
)
.
.
.
(
x
m
,
y
m
)
}
D=\{(\boldsymbol x_1,y_i), (\boldsymbol x_2, y_2)...(\boldsymbol x_m,y_m)\}
D={(x1,yi),(x2,y2)...(xm,ym)}
其中,
x
i
∈
R
n
\boldsymbol x_i \in {\rm {R}}^n
xi∈Rn,
y
i
∈
{
+
1
,
−
1
}
,
i
=
1
,
2
,
.
.
.
,
m
y_i\in\{+1, -1\}, i=1,2,...,m
yi∈{+1,−1},i=1,2,...,m.
x
i
\boldsymbol x_i
xi 是n维向量,为了便于可视化,我们以二维为例。
划分超平面表示为:
f
(
x
)
=
w
T
x
+
b
f(x) = \boldsymbol w^T\boldsymbol x + b
f(x)=wTx+b
我们将这个超平面标记为
(
w
,
b
)
(\boldsymbol w,b)
(w,b)
样本点到超平面的距离表示为:
γ
=
∣
w
T
x
+
b
∣
∣
∣
w
∣
∣
\gamma=\frac{|\boldsymbol w^T\boldsymbol x+b|}{||\boldsymbol w||}
γ=∣∣w∣∣∣wTx+b∣
假设超平面
(
w
,
b
)
(\boldsymbol w, b)
(w,b) 可以将样本点正确分类,那么每个样本集中的点到超平面的距离
γ
i
\gamma_i
γi 应该都大于等于支持向量(假设表示为
x
0
\boldsymbol x_0
x0)到超平面的距离
γ
^
\hat \gamma
γ^:
∣
w
T
x
i
+
b
∣
∣
∣
w
∣
∣
≥
∣
w
T
x
0
+
b
∣
∣
∣
w
∣
∣
,
i
=
1
,
2
,
.
.
.
m
\frac{|\boldsymbol w^T\boldsymbol x_i+b|}{||\boldsymbol w||} \geq \frac{|\boldsymbol w^T\boldsymbol x_0+b|}{||\boldsymbol w||}, i=1,2,...m
∣∣w∣∣∣wTxi+b∣≥∣∣w∣∣∣wTx0+b∣,i=1,2,...m
那么支持向量与超平面的间隔就是不等式右边的部分。但是它过于复杂,所以这里对
(
w
,
b
)
(\boldsymbol w,b)
(w,b) 进行缩放变换,使得
∣
w
T
x
0
+
b
∣
=
1
|\boldsymbol w^T\boldsymbol x_0+b| = 1
∣wTx0+b∣=1。这样处理对结果不会有影响,原因如下:
假设 w T x 0 + b = c , c ∈ R \boldsymbol w^T \boldsymbol x_0 + b = c, c \in R wTx0+b=c,c∈R,因为支持向量不在分割超平面上,所以 c ≠ 0 c \neq 0 c=0。因此只要两边同除 c c c 就可以了。
变换之后的不等式为:
∣
w
T
x
i
+
b
∣
∣
∣
w
∣
∣
≥
1
∣
∣
w
∣
∣
,
i
=
1
,
2
,
.
.
.
m
\frac{|\boldsymbol w^T\boldsymbol x_i+b|}{||\boldsymbol w||} \geq \frac{1}{||\boldsymbol w||}, i=1,2,...m
∣∣w∣∣∣wTxi+b∣≥∣∣w∣∣1,i=1,2,...m
两个异类支持向量之间的间隔表示为:
2
/
∣
∣
w
∣
∣
{2}/{||\boldsymbol w||}
2/∣∣w∣∣,就是目标函数。对上式两边同时乘以
∣
∣
w
∣
∣
||\boldsymbol w||
∣∣w∣∣ 可得,对样本集
D
D
D 的所有点满足:
{
w
T
x
i
+
b
≥
+
1
,
y
=
+
1
w
T
x
i
+
b
≤
−
1
,
y
=
−
1
\begin{cases} \boldsymbol w^T \boldsymbol x_i + b \geq+1, y=+1 \\ \boldsymbol w^T \boldsymbol x_i + b \leq -1, y=-1 \end{cases}
{wTxi+b≥+1,y=+1wTxi+b≤−1,y=−1
合并一下就是:
y
i
f
(
x
i
)
=
y
i
(
w
T
x
i
+
b
)
≥
1
y_if(x_i)=y_i(\boldsymbol w^T \boldsymbol x_i + b) \geq 1
yif(xi)=yi(wTxi+b)≥1
《机器学习》书中对这部分的描述如下:
目标函数的数学表达:
max
w
,
b
2
∣
∣
w
∣
∣
\max_{\boldsymbol w,b} \frac{2}{||\boldsymbol w||}
w,bmax∣∣w∣∣2
为了方便计算,将上式取倒数,并对
∣
∣
w
∣
∣
||\boldsymbol w||
∣∣w∣∣ 取平方,得到最终需要优化的目标函数:
min
w
,
b
1
2
∣
∣
w
∣
∣
2
s.t.
y
i
(
w
T
x
i
+
b
)
≥
1
,
i
=
1
,
2
,
.
.
.
m
\begin{aligned} & \min_{\boldsymbol w,b} \frac{1}{2}||\boldsymbol w||^2 \\[2ex] & \text{s.t. } y_i(\boldsymbol w^T \boldsymbol x_i + b) \geq 1, i=1,2,...m \end{aligned}
w,bmin21∣∣w∣∣2s.t. yi(wTxi+b)≥1,i=1,2,...m
3. 目标函数的优化
如果没有约束条件,这就是一个简单的多项式求极值的问题,只需令一阶导数等于零,二阶导数小于零求出
w
\boldsymbol w
w的所有元素即可。
但是有了约束条件就要使用些特殊方法,比如拉格朗日乘子法。
为什么要用拉格朗日乘子法?
在特征维度低的线性问题中是可以不使用的,具体参考:https://www.zhihu.com/question/36694952
3.1 拉格朗日乘子法(Lagrange Multiplier)
这里简单介绍一下拉格朗日乘子法及KKT条件的原理和使用。
对有一个等式约束的凸函数优化问题:
min
x
f
(
x
)
s.t.
h
(
x
)
=
0
\begin{aligned} & \min_{\boldsymbol x} f(\boldsymbol x) \\[2ex] & \text{s.t. } h(\boldsymbol x)=0 \end{aligned}
xminf(x)s.t. h(x)=0
在三维的情况下,
h
(
x
)
=
0
h(\boldsymbol x)=0
h(x)=0 是一个曲面,这个曲面在坐标轴 x-y 平面上的投影是一条曲线,也就是,看起来和 x-y 平面 “垂直”。
f
(
x
)
f(\boldsymbol x)
f(x)是一个凸函数,其与
h
(
x
)
h(\boldsymbol x)
h(x)相交。
上图中红色代表
f
(
x
)
f(\boldsymbol x)
f(x)的等值线,蓝色线代表
h
(
x
)
h(\boldsymbol x)
h(x)在x-y平面上的投影。
而两个曲面相交的部分是一条曲线,这条曲线在 f ( x ) f(\boldsymbol x) f(x)曲面上的最低点(就是图中黄色的点)就是 f ( x ) f(\boldsymbol x) f(x)取最小值的点。
这个点一定与
f
(
x
)
f(\boldsymbol x)
f(x)的某一条等值线相切。且在切点处两函数的梯度共线(可以同向或反向)。即存在
λ
≠
0
\lambda \neq 0
λ=0,使得在最小值点处:
∇
f
(
x
∗
)
+
λ
∇
h
(
x
∗
)
=
0
\nabla f(\boldsymbol x^*) + \lambda \nabla h(\boldsymbol x^*) = 0
∇f(x∗)+λ∇h(x∗)=0
其中,
x
∗
\boldsymbol x^*
x∗是最小值点。此时,令拉格朗日函数为:
L
(
x
,
λ
)
=
f
(
x
)
+
λ
h
(
x
)
L(\boldsymbol x,\lambda)= f(\boldsymbol x) + \lambda h(\boldsymbol x)
L(x,λ)=f(x)+λh(x)
对该函数求导得:
∇
L
(
x
,
λ
)
=
∇
f
(
x
)
+
λ
∇
h
(
x
)
\nabla L(\boldsymbol x,\lambda)=\nabla f(\boldsymbol x) + \lambda \nabla h\boldsymbol (x)
∇L(x,λ)=∇f(x)+λ∇h(x)
显然,当
x
=
x
∗
\boldsymbol x=\boldsymbol x^*
x=x∗时,
∇
L
(
x
,
λ
)
=
0
\nabla L(\boldsymbol x,\lambda) = 0
∇L(x,λ)=0。也就是说,
L
(
x
,
λ
)
L(\boldsymbol x,\lambda)
L(x,λ)与
f
(
x
)
f(\boldsymbol x)
f(x)同时取得极小值。
也就是说,带有等式约束的优化问题转化为了不带约束条件的优化问题。
令 L ( x , λ ) L(\boldsymbol x, \lambda) L(x,λ) 对 x , λ \boldsymbol x,\lambda x,λ 的导数都为0就可以得到最小值。
3.2 KKT条件
以上方法适用于在约束为等式的情况。当约束为不等式时,即:
min
x
f
(
x
)
s.t.
g
(
x
)
≤
0
\begin{aligned} & \min_{\boldsymbol x} f(\boldsymbol x) \\[2ex] & \text{s.t. } g(\boldsymbol x) \leq 0 \end{aligned}
xminf(x)s.t. g(x)≤0
极小值点的位置存在以下两种情况(以一个不等式约束为例):
- 极小值点在 g ( x ) < 0 g(\boldsymbol x)<0 g(x)<0 区域内,此时可以忽略约束条件直接对 f ( x ) f(\boldsymbol x) f(x)求极小值即可,相当于把上节中的 λ \lambda λ 置零;
- 极小值点在边界 g ( x ) = 0 g(\boldsymbol x)=0 g(x)=0 上,此时与3.1节所述情况一样。唯一不同的是,在上图所示的极值点处, f ( x ) f(\boldsymbol x) f(x) 与 g ( x ) g(\boldsymbol x) g(x) 的梯度方向一定相反。(如果它们梯度方向相同,就可以继续向两者减小的方向优化)
将以上两种情况结合可得,存在
μ
≥
0
\mu \geq 0
μ≥0 使得:
μ
g
(
x
)
=
0
\mu g(\boldsymbol x)=0
μg(x)=0
由此引出
f
(
x
)
f(\boldsymbol x)
f(x) 取得极小值的充分必要条件,即KKT(Karush-Kuhn-Tucker)条件:
{
g
(
x
)
≤
0
μ
≥
0
μ
g
(
x
)
=
0
\begin{cases} g(\boldsymbol x) \leq 0 \\ \mu \geq 0 \\ \mu g(\boldsymbol x)=0 \end{cases}
⎩⎪⎨⎪⎧g(x)≤0μ≥0μg(x)=0
原函数
f
(
x
)
f(\boldsymbol x)
f(x) 转化为与上节形式相同的拉格朗日函数。
将这几种情况结合起来,推广到具有多个不等式或等式约束的优化问题上:
对有有限个约束的最小化问题:
min
x
f
(
x
)
s.t.
h
i
(
x
)
=
0
,
i
=
1
,
,
2
,
.
.
.
m
g
j
(
x
)
≤
0
,
j
=
1
,
,
2
,
.
.
.
n
\begin{aligned} & \min_{\boldsymbol x} f(\boldsymbol x) \\[2ex] \text{s.t. }& h_i(\boldsymbol x)=0, i=1,,2,...m \\[2ex] & g_j(\boldsymbol x)\leq 0,j=1,,2,...n \\[2ex] \end{aligned}
s.t. xminf(x)hi(x)=0,i=1,,2,...mgj(x)≤0,j=1,,2,...n
其拉格朗日函数为:
L
(
x
,
λ
i
)
=
f
(
x
)
+
∑
i
=
1
m
λ
i
h
i
(
x
)
+
∑
j
=
1
n
μ
j
g
j
(
x
)
L(\boldsymbol x,\lambda_i)= f(\boldsymbol x) + \sum_{i=1}^m \lambda_i h_i (\boldsymbol x)+ \sum_{j=1}^n \mu_j g_j (\boldsymbol x)
L(x,λi)=f(x)+i=1∑mλihi(x)+j=1∑nμjgj(x)
KKT条件为:
{
h
i
(
x
)
=
0
g
j
(
x
)
≤
0
λ
i
,
μ
j
≥
0
μ
j
g
j
(
x
)
=
0
\begin{cases} h_i(\boldsymbol x) =0 \\[2ex] g_j(\boldsymbol x) \leq 0 \\[2ex] \lambda_i, \mu_j \geq 0 \\[2ex] \mu_j g_j(\boldsymbol x)=0 \end{cases}
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧hi(x)=0gj(x)≤0λi,μj≥0μjgj(x)=0
3.3 对偶问题
现在我们回到第2节末尾,再看这个优化问题:
min
w
,
b
1
2
∣
∣
w
∣
∣
2
s.t.
y
i
(
w
T
x
i
+
b
)
≥
1
,
i
=
1
,
2
,
.
.
.
m
\begin{aligned} & \min_{\boldsymbol w,b} \frac{1}{2}||\boldsymbol w||^2 \\[2ex] & \text{s.t. } y_i(\boldsymbol w^T \boldsymbol x_i + b) \geq 1, i=1,2,...m \end{aligned}
w,bmin21∣∣w∣∣2s.t. yi(wTxi+b)≥1,i=1,2,...m
这就是个适用于不等式约束的优化。拉格朗日函数如下:
L
(
w
,
b
,
α
)
=
1
2
∣
∣
w
∣
∣
2
+
∑
i
=
1
m
α
i
(
1
−
y
i
(
w
T
x
i
+
b
)
)
,
α
i
≥
0
L(\boldsymbol w, b, \alpha) = \frac{1}{2}||\boldsymbol w||^2 + \sum_{i=1}^m{\alpha_i(1-y_i(\boldsymbol w^T \boldsymbol x_i + b))}, \alpha_i \geq 0
L(w,b,α)=21∣∣w∣∣2+i=1∑mαi(1−yi(wTxi+b)),αi≥0
若要求
L
(
w
,
b
,
α
)
L(\boldsymbol w, b, \alpha)
L(w,b,α) 的最小值则首先把
w
,
b
\boldsymbol w, b
w,b 看作常量,对
L
(
w
,
b
,
α
)
L(\boldsymbol w, b, \alpha)
L(w,b,α)求关于
α
\alpha
α的最大值。因为
1
−
y
i
(
w
T
x
i
+
b
)
1-y_i(\boldsymbol w^T \boldsymbol x_i + b)
1−yi(wTxi+b) 是小于等于0的,在
α
i
≥
0
\alpha_i \geq 0
αi≥0的前提下,其最大值就是:
1
2
∣
∣
w
∣
∣
2
\frac{1}{2}||\boldsymbol w||^2
21∣∣w∣∣2。
因此,原问题可写为:
min
w
,
b
max
α
i
L
(
w
,
b
,
α
)
\min_{\boldsymbol w, b}\max_{\alpha_i}L(\boldsymbol w, b, \alpha)
w,bminαimaxL(w,b,α)
直接对该问题求解比较困难,因此这里引入对偶问题。
使用对偶问题求解的原因不止是为了计算方便,还是为线性不可分时核函数的引入作铺垫(本文不涉及这部分)。
先上结论:
min
w
,
b
max
α
i
L
(
w
,
b
,
α
)
=
max
α
i
min
w
,
b
L
(
w
,
b
,
α
)
\min_{\boldsymbol w, b}\max_{\alpha_i}L(\boldsymbol w, b, \alpha)=\max_{\alpha_i}\min_{\boldsymbol w, b}L(\boldsymbol w, b, \alpha)
w,bminαimaxL(w,b,α)=αimaxw,bminL(w,b,α)
也就是求取最大最小值的顺序发生了变化。从先对 α \alpha α求最大值变为先对 w , b \boldsymbol w, b w,b 求最小值。
等式两边的问题互为对偶问题。
当然两个对偶问题同时取得最优解是有条件的,具体证明可以参考:https://www.cnblogs.com/90zeng/p/Lagrange_duality.html
转化为对偶问题后,令
w
,
b
\boldsymbol w, b
w,b 的导数为0:
{
∂
L
∂
w
=
w
−
∑
i
=
1
m
α
i
y
i
x
i
=
0
∂
L
∂
b
=
−
∑
i
=
1
m
α
i
y
i
=
0
\begin{cases} \cfrac{\partial L}{\partial \boldsymbol w}=\boldsymbol w-\sum_{i=1}^m\alpha_iy_i\boldsymbol x_i=0 \\[2ex] \cfrac{\partial L}{\partial b}=-\sum_{i=1}^m\alpha_iy_i=0 \end{cases}
⎩⎪⎪⎨⎪⎪⎧∂w∂L=w−∑i=1mαiyixi=0∂b∂L=−∑i=1mαiyi=0
将上式代入
L
(
w
,
b
,
α
)
L(\boldsymbol w, b, \alpha)
L(w,b,α) 并化简可得:
L
(
w
,
b
,
α
)
=
∑
i
=
1
m
α
i
−
1
2
∑
i
=
1
m
∑
j
=
1
m
α
i
α
j
y
i
y
j
x
i
T
x
j
L(\boldsymbol w, b, \alpha)=\sum_{i=1}^m\alpha_i-\frac12\sum_{i=1}^m\sum_{j=1}^m\alpha_i \alpha_jy_iy_j\boldsymbol x_i^T\boldsymbol x_j
L(w,b,α)=i=1∑mαi−21i=1∑mj=1∑mαiαjyiyjxiTxj
化简过程如下:(公式太多还是手写吧emmm…)
化简后就只用求
α
\alpha
α 即可。这样问题的复杂度就从与样本的维度有关变为只与样本数量有关了。
周志华的《机器学习》一书对这个最优解的几何意义有一段描述可以参考:
3.4 SMO(Sequential Minimal Optimazation)方法
上节最后得到的关于 α \alpha α 的式子是一个二次规划问题,使用通用算法开销比较大。SMO方法是1998年提出的,用于求取SVM中的拉格朗日乘子,比通用算法更加高效。
其主要思想为,选取两个变量 α i , α j \alpha_i, \alpha_j αi,αj,固定其他参数,对这两个参数进行优化,然后重复这个过程。
Note:选取
α
i
,
α
j
\alpha_i, \alpha_j
αi,αj的标准:
不过博主有点懒,所以直接随机选取了。。。。需要深入了解可以参考《机器学习》(文末有下载方法~)
每次选取两个的原因是这些参数之间有约束 ∑ i = 1 m α i y i = 0 \sum_{i=1}^m\alpha_iy_i=0 ∑i=1mαiyi=0,只改变一个约束会被打破。
为了表示方便,记选取的两个参数为
α
1
,
α
2
\alpha_1, \alpha_2
α1,α2,固定除
α
1
,
α
2
\alpha_1, \alpha_2
α1,α2以外的参数:
α
1
y
1
+
α
2
y
2
=
−
∑
i
≠
1
,
2
α
i
y
i
=
ξ
α
1
=
y
1
(
ξ
−
α
2
y
2
)
\alpha_1y_1+ \alpha_2y_2 =- \sum_{i\neq 1,2}\alpha_i y_i=\xi \\[2ex] \alpha_1 = y_1(\xi-\alpha_2y_2)
α1y1+α2y2=−i=1,2∑αiyi=ξα1=y1(ξ−α2y2)
带入 L ( w , b , α ) L(\boldsymbol w, b, \alpha) L(w,b,α)得:
为简便表示,下式中, K i j = x i T x j K_{ij} =\boldsymbol x_i^T\boldsymbol x_j Kij=xiTxj
L
(
w
,
b
,
α
)
=
∑
i
=
1
m
α
i
−
1
2
∑
i
=
1
m
∑
j
=
1
m
α
i
α
j
y
i
y
j
K
i
j
=
∑
i
≥
3
m
α
i
+
y
1
ξ
−
α
2
y
1
y
2
+
α
2
−
1
2
∑
i
≥
3
m
∑
j
≥
3
m
α
i
α
j
y
i
y
j
K
i
j
−
1
2
(
y
1
ξ
−
α
2
y
1
y
2
)
2
K
11
−
y
2
α
2
(
ξ
−
α
2
y
2
)
K
12
−
1
2
α
2
2
K
22
−
∑
i
≥
3
m
(
ξ
−
α
2
y
2
)
α
i
y
i
K
1
i
−
∑
i
≥
3
m
α
2
y
2
α
i
y
i
K
2
i
=
−
α
2
y
1
y
2
+
α
2
−
1
2
(
y
1
ξ
−
α
2
y
1
y
2
)
2
K
11
−
y
2
α
2
(
−
α
2
y
2
)
K
12
−
1
2
α
2
2
K
22
−
∑
i
≥
3
m
(
ξ
−
α
2
y
2
)
α
i
y
i
K
1
i
−
∑
i
≥
3
m
α
2
y
2
α
i
y
i
K
2
i
+
c
o
n
s
t
a
n
t
=
−
α
2
y
1
y
2
+
α
2
−
1
2
α
2
2
K
11
+
α
2
y
2
ξ
K
11
−
α
2
y
2
ξ
K
12
+
α
2
2
K
12
−
1
2
α
2
2
K
22
+
∑
i
≥
3
m
α
2
y
2
α
i
y
i
K
1
i
−
∑
i
≥
3
m
α
2
y
2
α
i
y
i
K
2
i
+
c
o
n
s
t
a
n
t
\begin{aligned} &L(\boldsymbol w, b, \alpha) \\[2ex] =& \sum_{i=1}^m\alpha_i-\frac12\sum_{i=1}^m\sum_{j=1}^m\alpha_i \alpha_jy_iy_jK_{ij} \\[2ex] =& \sum_{i\geq 3}^m\alpha_i + y_1\xi-\alpha_2y_1y_2 + \alpha_2 - \cfrac12 \sum_{i\geq 3}^m\sum_{j\geq 3}^m\alpha_i \alpha_jy_iy_jK_{ij} - \cfrac12(y_1\xi-\alpha_2y_1y_2)^2K_{11}-y_2\alpha_2(\xi-\alpha_2y_2)K_{12}-\cfrac12 \alpha_2^2 K_{22} \\[2ex] & -\sum_{i\geq 3}^m(\xi-\alpha_2y_2)\alpha_iy_iK_{1i} - \sum_{i\geq 3}^m\alpha_2y_2\alpha_iy_iK_{2i} \\[2ex] =& -\alpha_2y_1y_2 + \alpha_2 - \cfrac12(y_1\xi-\alpha_2y_1y_2)^2K_{11}-y_2\alpha_2(-\alpha_2y_2)K_{12}-\cfrac12 \alpha_2^2 K_{22} \\[2ex] & -\sum_{i\geq 3}^m(\xi-\alpha_2y_2)\alpha_iy_iK_{1i} - \sum_{i\geq 3}^m\alpha_2y_2\alpha_iy_iK_{2i}+ constant \\[2ex] =& -\alpha_2y_1y_2 + \alpha_2 - \cfrac12\alpha_2^2K_{11} + \alpha_2y_2\xi K_{11} - \alpha_2y_2\xi K_{12} + \alpha_2^2K_{12}-\cfrac12 \alpha_2^2 K_{22} \\[2ex] & +\sum_{i\geq 3}^m\alpha_2y_2\alpha_iy_iK_{1i} - \sum_{i\geq 3}^m\alpha_2y_2\alpha_iy_iK_{2i} + constant \\[2ex] \end{aligned}
====L(w,b,α)i=1∑mαi−21i=1∑mj=1∑mαiαjyiyjKiji≥3∑mαi+y1ξ−α2y1y2+α2−21i≥3∑mj≥3∑mαiαjyiyjKij−21(y1ξ−α2y1y2)2K11−y2α2(ξ−α2y2)K12−21α22K22−i≥3∑m(ξ−α2y2)αiyiK1i−i≥3∑mα2y2αiyiK2i−α2y1y2+α2−21(y1ξ−α2y1y2)2K11−y2α2(−α2y2)K12−21α22K22−i≥3∑m(ξ−α2y2)αiyiK1i−i≥3∑mα2y2αiyiK2i+constant−α2y1y2+α2−21α22K11+α2y2ξK11−α2y2ξK12+α22K12−21α22K22+i≥3∑mα2y2αiyiK1i−i≥3∑mα2y2αiyiK2i+constant
因为我们的最终目的是让
L
(
w
,
b
,
α
)
L(\boldsymbol w, b, \alpha)
L(w,b,α) 对
α
2
\alpha_2
α2 求导,而其他
α
\alpha
α 已被固定,视为常数。且
x
,
y
x, y
x,y 也都是已知常数,所以将上式中不含
α
2
\alpha_2
α2 的常数项合并为
c
o
n
s
t
a
n
t
constant
constant。
下一步,令其 α 2 \alpha_2 α2 的导数为0:
∂ L ( w , b , α ) ∂ α 2 = − y 1 y 2 + 1 − α 2 K 11 + y 2 ξ K 11 − y 2 ξ K 12 + 2 α 2 K 12 − α 2 K 22 + ∑ i ≥ 3 m y 2 α i y i K 1 i − ∑ i ≥ 3 m y 2 α i y i K 2 i = 0 \cfrac{\partial L(\boldsymbol w, b, \alpha)}{\partial \alpha_2} = -y_1y_2 + 1 - \alpha_2K_{11} + y_2\xi K_{11} - y_2\xi K_{12} + 2\alpha_2K_{12}- \alpha_2 K_{22} + \sum_{i\geq 3}^my_2\alpha_iy_iK_{1i} - \sum_{i\geq 3}^my_2\alpha_iy_iK_{2i}=0 ∂α2∂L(w,b,α)=−y1y2+1−α2K11+y2ξK11−y2ξK12+2α2K12−α2K22+i≥3∑my2αiyiK1i−i≥3∑my2αiyiK2i=0
由上节可知:
w
=
∑
i
=
1
m
α
i
y
i
x
i
f
(
x
1
)
=
w
T
x
1
+
b
=
∑
i
=
1
m
α
i
y
i
K
1
i
+
b
f
(
x
2
)
=
w
T
x
2
+
b
=
∑
i
=
1
m
α
i
y
i
K
2
i
+
b
\boldsymbol w=\sum_{i=1}^m\alpha_iy_i\boldsymbol x_i \\[2ex] f(\boldsymbol x_1)=\boldsymbol w^Tx_1+b=\sum_{i=1}^m\alpha_iy_iK_{1i} + b \\[2ex] f(\boldsymbol x_2)=\boldsymbol w^Tx_2+b=\sum_{i=1}^m\alpha_iy_iK_{2i} + b \\[2ex]
w=i=1∑mαiyixif(x1)=wTx1+b=i=1∑mαiyiK1i+bf(x2)=wTx2+b=i=1∑mαiyiK2i+b
注意到这个形式和上面的最后两项很相似,所以可以带进去计算。
但是需要注意的是这里的 f ( x ) f(\boldsymbol x) f(x)是常量(否则带入一个变量并不利于求导结果的计算),所以这里的 f ( x ) f(\boldsymbol x) f(x)中包含的 α i ( i = 1 , 2 , . . . , m ) \alpha_i(i=1,2,...,m) αi(i=1,2,...,m) 是初始化时赋予的值,将其记作 α 2 o l d \alpha_2^{old} α2old;把其他的 α 2 \alpha_2 α2 看作待更新的变量,记作 α 2 n e w \alpha_2^{new} α2new。
∂ L ( w , b , α ) ∂ α 2 = − y 1 y 2 + 1 − α 2 n e w K 11 + y 2 ξ K 11 − y 2 ξ K 12 + 2 α 2 n e w K 12 − α 2 n e w K 22 + y 2 [ f ( x 1 ) − ( ξ − α 2 o l d y 2 ) K 11 − α 2 o l d y 2 K 12 − b ] − y 2 [ f ( x 2 ) − ( ξ − α 2 o l d y 2 ) K 12 − α 2 o l d y 2 K 22 − b ] = ( − K 11 − K 22 + 2 K 12 ) α 2 n e w + ( K 11 + K 22 − 2 K 12 ) α 2 o l d + 1 + y 1 y 2 + y 2 [ f ( x 1 ) − f ( x 2 ) ] = 0 \begin{aligned} \cfrac{\partial L(\boldsymbol w, b, \alpha)}{\partial \alpha_2} = & -y_1y_2 + 1 - \alpha_2^{new}K_{11} + y_2\xi K_{11} - y_2\xi K_{12} + 2\alpha_2^{new}K_{12}- \alpha_2^{new} K_{22} \\[2ex] & +y_2[f(x_1)-(\xi - \alpha_2^{old}y_2)K_{11}- \alpha_2^{old}y_2K_{12} - b] - y_2[f(x_2)-(\xi - \alpha_2^{old}y_2)K_{12}- \alpha_2^{old}y_2K_{22} - b] \\[2ex] = & (-K_{11}-K_{22}+2K_{12})\alpha_2^{new} + (K_{11}+K_{22}-2K_{12})\alpha_2^{old}+1+y_1y_2+y_2[f(x_1)-f(x_2)] \\[2ex] = & 0 \end{aligned} ∂α2∂L(w,b,α)===−y1y2+1−α2newK11+y2ξK11−y2ξK12+2α2newK12−α2newK22+y2[f(x1)−(ξ−α2oldy2)K11−α2oldy2K12−b]−y2[f(x2)−(ξ−α2oldy2)K12−α2oldy2K22−b](−K11−K22+2K12)α2new+(K11+K22−2K12)α2old+1+y1y2+y2[f(x1)−f(x2)]0
令
η
=
K
11
+
K
22
−
2
K
12
\eta=K_{11}+K_{22}-2K_{12}
η=K11+K22−2K12,得:
∂
L
(
w
,
b
,
α
)
∂
α
2
=
−
η
α
2
n
e
w
+
η
α
2
o
l
d
+
1
+
y
1
y
2
+
y
2
[
f
(
x
1
)
−
f
(
x
2
)
]
=
0
α
2
n
e
w
=
α
2
o
l
d
+
1
+
y
1
y
2
+
y
2
[
f
(
x
1
)
−
f
(
x
2
)
]
η
\cfrac{\partial L(\boldsymbol w, b, \alpha)}{\partial \alpha_2} = -\eta\alpha_2^{new}+\eta\alpha_2^{old} +1+y_1y_2+y_2[f(\boldsymbol x_1)-f(\boldsymbol x_2)]=0 \\[2ex] \alpha_2^{new} = \alpha_2^{old}+\cfrac{1+y_1y_2+y_2[f(\boldsymbol x_1)-f(\boldsymbol x_2)]}{\eta}
∂α2∂L(w,b,α)=−ηα2new+ηα2old+1+y1y2+y2[f(x1)−f(x2)]=0α2new=α2old+η1+y1y2+y2[f(x1)−f(x2)]
同时,因为需要求最大值,所以还要求对
α
2
\alpha_2
α2 的二阶导小于零:
∂ 2 L ( w , b , α ) ∂ 2 α 2 = − K 11 − K 22 + 2 K 12 ≤ 0 η = K 11 + K 22 − 2 K 12 ≥ 0 \cfrac{\partial^2 L(\boldsymbol w, b, \alpha)}{\partial^2 \alpha_2} = -K_{11}-K_{22}+2K_{12} \leq 0 \\[2ex] \eta = K_{11}+K_{22}-2K_{12} \geq 0 ∂2α2∂2L(w,b,α)=−K11−K22+2K12≤0η=K11+K22−2K12≥0
到这里为止,
α
2
\alpha_2
α2 的更新已经完成,
α
1
\alpha_1
α1 可以通过两者之间的约束求得。
α
1
n
e
w
y
1
+
α
2
n
e
w
y
2
=
α
1
o
l
d
y
1
+
α
2
o
l
d
y
2
=
ξ
\alpha_1^{new}y_1 + \alpha_2^{new}y_2 = \alpha_1^{old}y_1 + \alpha_2^{old}y_2=\xi \\[2ex]
α1newy1+α2newy2=α1oldy1+α2oldy2=ξ
最后一步,就是求解参数
b
b
b。
本文只采用了前一种方法,感兴趣的童鞋可以用(6.18)的方法试试,欢迎在评论区讨论~
具体操作如下:
假设更新过后的 α 2 n e w > 0 \alpha_2^{new} >0 α2new>0,即该 α 2 n e w \alpha_2^{new} α2new 对应的样本为支持向量,则满足 y 2 f n e w ( x 2 ) = 1 y_2f^{new}(\boldsymbol x_2)=1 y2fnew(x2)=1,结合更新之前的 f ( x 2 ) f(\boldsymbol x_2) f(x2):
{
f
(
x
2
)
=
∑
i
≥
3
m
α
i
y
i
K
2
i
+
α
1
o
l
d
y
1
K
11
+
α
2
o
l
d
y
2
K
12
+
b
o
l
d
y
2
(
∑
i
≥
3
m
α
i
y
i
K
2
i
+
α
1
n
e
w
y
1
K
11
+
α
2
n
e
w
y
2
K
12
+
b
n
e
w
)
=
1
\begin{cases} f(\boldsymbol x_2) = \sum_{i \geq 3}^{m}\alpha_iy_iK_{2i} + \alpha_1^{old}y_1K_{11}+\alpha_2^{old}y_2K_{12} + b^{old} \\[2ex] y_2(\sum_{i \geq 3}^{m}\alpha_iy_iK_{2i} + \alpha_1^{new}y_1K_{11}+\alpha_2^{new}y_2K_{12} + b^{new}) = 1 \end{cases}
⎩⎨⎧f(x2)=∑i≥3mαiyiK2i+α1oldy1K11+α2oldy2K12+boldy2(∑i≥3mαiyiK2i+α1newy1K11+α2newy2K12+bnew)=1
解得:
b
n
e
w
=
y
2
−
f
(
x
2
)
−
(
α
1
n
e
w
−
α
1
o
l
d
)
y
1
K
11
+
(
α
2
n
e
w
−
α
2
o
l
d
)
y
2
K
12
b^{new}= y_2 - f(\boldsymbol x_2)-(\alpha_1^{new} - \alpha_1^{old})y_1K_{11} + (\alpha_2^{new} - \alpha_2^{old})y_2K_{12}
bnew=y2−f(x2)−(α1new−α1old)y1K11+(α2new−α2old)y2K12
4 软间隔
由于不是所有的数据集都可以完美的用一条线分割开,在两个类别中间可能会有个别样本点超过分割平面,使得数据集不能被完美分开。即便可以被分开,也有可能是过拟合造成的。
这时我们的算法需要对这些异常点进行容错,使其不影响数据集正常划分。
异常点的表示就是其不满足约束条件:
y i ( w T x i + b ) ≥ 1 y_i(\boldsymbol w^T \boldsymbol x_i + b) \geq 1 yi(wTxi+b)≥1
不过我们希望不满足约束的样本越少越好,因此将原问题改为:
min w , b L ( w , b ) = min w , b 1 2 ∣ ∣ w ∣ ∣ 2 + C ∑ i = 1 m l 0 / 1 ( y i ( w T x i + b ) − 1 ) \min_{\boldsymbol w, b}L(\boldsymbol w, b) =\min_{\boldsymbol w, b} \frac{1}{2} || \boldsymbol w ||^2 +C \sum_{i=1}^m{l_{0/1}(y_i(\boldsymbol w^T \boldsymbol x_i + b ) - 1)} w,bminL(w,b)=w,bmin21∣∣w∣∣2+Ci=1∑ml0/1(yi(wTxi+b)−1)
C > 0 C>0 C>0 是一个常数。当其取较大的值时,就迫使更多的样本趋于满足约束条件,取无穷大时就和原问题等价。
l 0 / 1 l_{0/1} l0/1叫做“0/1损失函数”,数学表示为:
l 0 / 1 ( z ) = { 1 , if z < 0 0 , else l_{0/1}(z) = \begin{cases} 1, & \text{if $z<0$ } \\[2ex] 0, & \text{else} \end{cases} l0/1(z)=⎩⎨⎧1,0,if z<0 else
但是它数学性质不好,不连续,所以我们用合页损失函数(Hinge loss)代替:
l h i n g e ( z ) = max ( 0 , 1 − z ) l_{hinge}(z)=\max(0, 1-z) lhinge(z)=max(0,1−z)
如下图(图片来自《机器学习》):
此时,优化目标变为:
min w , b 1 2 ∣ ∣ w ∣ ∣ 2 + C ∑ i = 1 m max ( 0 , 1 − y i ( w T x i + b ) ) \min_{\boldsymbol w, b} \frac{1}{2}||\boldsymbol w||^2 +C \sum_{i=1}^m{\max(0, 1 - y_i (\boldsymbol w^T \boldsymbol x_i + b))} w,bmin21∣∣w∣∣2+Ci=1∑mmax(0,1−yi(wTxi+b))
这里有一个问题,按照合页损失函数的定义,把前面式子中的求和项带入,得到的结果应该是: 1 − z = 1 − ( y i ( w T x i + b ) − 1 ) = 2 − y i ( w T x i + b ) 1-z=1-(y_i(\boldsymbol w^T \boldsymbol x_i + b ) - 1) = 2-y_i(\boldsymbol w^T \boldsymbol x_i+ b ) 1−z=1−(yi(wTxi+b)−1)=2−yi(wTxi+b)
个人理解是,根据我们要解决的问题对hinge损失函数做了改动:
l h i n g e ( z ) = max ( 0 , − z ) l_{hinge}(z)=\max(0, -z) lhinge(z)=max(0,−z)
优化目标里含有max(),显然不便于计算,所以这里引入“松弛变量(slack variables)”
ξ
i
≥
0
\xi_i \geq 0
ξi≥0,表示每个样本不满足约束的程度。若样本
x
i
x_i
xi 满足约束,则对应的
ξ
i
=
0
\xi_i = 0
ξi=0;若不满足,则
1
−
y
i
(
w
T
x
i
+
b
)
≤
ξ
i
1 - y_i (\boldsymbol w^T \boldsymbol x_i + b) \leq \xi_i
1−yi(wTxi+b)≤ξi。优化目标可以重新写为:
min
w
,
b
1
2
∣
∣
w
∣
∣
2
+
C
∑
i
=
1
m
ξ
i
s.t.
1
−
y
i
(
w
T
x
i
+
b
)
≤
ξ
i
ξ
i
≥
0
\begin{aligned} &\min_{\boldsymbol w, b} \frac{1}{2}||\boldsymbol w||^2 +C \sum_{i=1}^m{\xi_i} \\[2ex] \text{s.t. } & 1 - y_i (\boldsymbol w^T \boldsymbol x_i + b) \leq \xi_i \\[2ex] & \xi_i \geq 0 \end{aligned}
s.t. w,bmin21∣∣w∣∣2+Ci=1∑mξi1−yi(wTxi+b)≤ξiξi≥0
和前面一样,重新使用拉格朗日乘子法。拉格朗日函数:
L ( w , b , α , μ ) = 1 2 ∣ ∣ w ∣ ∣ 2 + C ∑ i = 1 m ξ i + ∑ i = 1 m α i ( 1 − y i ( w T x i + b ) − ξ i ) − ∑ i = 1 m μ i ξ i , α i , μ i ≥ 0 L(\boldsymbol w, b, \alpha, \mu) = \frac{1}{2}||\boldsymbol w||^2 +C \sum_{i=1}^m{\xi_i} + \sum_{i=1}^m{\alpha_i(1-y_i(\boldsymbol w^T \boldsymbol x_i + b) - \xi_i)} - \sum_{i=1}^m\mu_i\xi_i,\\[2ex] \alpha_i,\mu_i \geq 0 L(w,b,α,μ)=21∣∣w∣∣2+Ci=1∑mξi+i=1∑mαi(1−yi(wTxi+b)−ξi)−i=1∑mμiξi,αi,μi≥0
求解过程也和前面一样,贴个图:
和前面得到的结果对比一下,唯一不同的地方在于
α
i
\alpha_i
αi 的取值范围有了上界。这个上界来源于图片上的式(6.39),意义就是表示了算法对于样本点的容错程度。
5 代码
下面,我们将把以上的算法用Python程序实现。
首先我们对上面的算法步骤作一个总结,帮助我们理清思路:
整个SMO算法分为以下几个阶段:
- 选择 α 1 , α 2 \alpha_1, \alpha_2 α1,α2
- 判断 η = K 11 + K 22 − 2 K 12 ≥ 0 \eta = K_{11}+K_{22}-2K_{12} \geq 0 η=K11+K22−2K12≥0
- 计算 α 2 n e w = α 2 o l d + 1 + y 1 y 2 + y 2 [ f ( x 1 ) − f ( x 2 ) ] η \alpha_2^{new} = \alpha_2^{old}+\cfrac{1+y_1y_2+y_2[f(x_1)-f(x_2)]}{\eta} α2new=α2old+η1+y1y2+y2[f(x1)−f(x2)]
- 对 α 1 , α 2 \alpha_1, \alpha_2 α1,α2 做上下界限定,范围为 [ 0 , C ] [0, C] [0,C]
- 判断 α 2 \alpha_2 α2 的改变量是否足够大,太小视为不变
- 计算 α 1 n e w = α 1 o l d + α 2 o l d y 1 y 2 − α 2 n e w y 1 y 2 \alpha_1^{new} = \alpha_1^{old} + \alpha_2^{old}y_1y_2 - \alpha_2^{new}y_1y_2 α1new=α1old+α2oldy1y2−α2newy1y2
- 计算 b b b
- 重复以上步骤
5.1 数据集生成及一些辅助函数
写了一个简单的函数用来随机生成数据集。
指定一条曲线 f ( x ) f(x) f(x),本例中是 f ( x ) = x f(x)=x f(x)=x,随机均匀生成100对数据,取值为[0, 10],标为两个类别并画图表示。(为了突出两数据集的间隔,这里舍弃与分割线垂直距离小于1的点)
import matplotlib.pyplot as plt
import numpy as np
def generate_dataset():
# np.random.seed(3)
class1 = []; class2 = []
label1 = []; label2 =[]
# define decision line
def f(x):
return 1*x
for _ in range(100):
x = np.random.rand() * 10
y = np.random.rand() * 10
if y-f(x) > 1:
class1.append([x, y])
label1.append(1)
elif y-f(x) < -1:
class2.append([x, y])
label2.append(-1)
x1, y1 = zip(*[data for data in class1])
x2, y2 = zip(*[data for data in class2])
plt.plot(x1, y1, 'ro')
plt.plot(x2, y2, 'bo')
plt.axis([0,10,0,10])
plt.show()
return class1+class2, label1+label2
随机选择第2个参数的下标:
def select_j(i, m):
j = i
while(j == i):
j = np.random.randint(0, m)
return j
限制 α \alpha α 的上下界:
def clip_alpha(aj, H: float, L:float):
return min(max(L, aj), H)
已知 α \alpha α 求 w \boldsymbol w w:
def get_W(alphas, dataset, label):
a, x, y = map(np.array, [alphas, dataset, label])
W = np.dot(a * y, x)
# transform W form np.array to Python List
return W.tolist()
5.2 主函数
简化版的SMO算法:
def smo_simple(dataset, labels, C, max_iter):
data = np.array(dataset, dtype=np.float)
label = np.array(labels, dtype=np.float)
b = 0
m, n = np.shape(data)
alphas = np.zeros(m, dtype=np.float)
iter = 0
while iter < max_iter:
alpha_pair_changed = 0
for i in range(m):
x_i, y_i = data[i], label[i]
fx_i = np.dot(alphas*label, np.dot(data, x_i)) + b
e_i = fx_i - y_i
j = select_j(iter, m)
x_j, y_j = data[j], label[j]
fx_j = np.dot(alphas*label, np.dot(data, x_j)) + b
e_j = fx_j - y_j
# calculate a_j_new
eta = np.dot(x_i, x_i) + np.dot(x_j, x_j) - 2*np.dot(x_i, x_j)
if eta <= 0:
print("eta <= 0, continue")
continue
a_i, a_j = alphas[i], alphas[j]
a_j_new = a_j + y_j * (e_i - e_j) / eta
# limit a_j_new to [0, C]
if y_i == y_j:
L = max(0., a_i + a_j - C)
H = min(a_i + a_j, C)
else:
L = max(0., a_j - a_i)
H = min(C - a_i + a_j, C)
if L < H:
a_j_new = clip_alpha(a_j_new, H, L)
else:
print("L >= H. (L, H) =", (L, H))
continue
# judge if a_j moves an enough distance
if abs(a_j_new - a_j) < 0.00001:
print("a_j has not moved enough, a_j_new - a_j = %f" % (a_j_new - a_j))
continue
# calculate a_i_new and update a_i, a_j
a_i_new = (a_j - a_j_new)*y_i*y_j + a_i
alphas[i], alphas[j] = a_i_new, a_j_new
alpha_pair_changed += 1
#calculate b
b_i = -e_i + (a_i - a_i_new)*y_i*np.dot(x_i, x_i) + (a_j - a_j_new)*y_j*np.dot(x_i, x_j) + b
b_j = -e_j + (a_i - a_i_new)*y_i*np.dot(x_i, x_j) + (a_j - a_j_new)*y_j*np.dot(x_j, x_j) + b
# b = b_i if b_i == b_j else (b_i + b_j)/2
if 0 < a_i_new < C:
b = b_i
elif 0 < a_j_new < C:
b = b_j
else:
b = (b_i + b_j)/2
print("(a_i, a_j) moved from (%f, %f) to (%f, %f)." % (a_i, a_j, a_i_new, a_j_new))
if alpha_pair_changed == 0:
print("Iteration %d of max_iter %d" % (iter+1, max_iter))
iter += 1
else:
iter = 0
return alphas, b
main函数:
dataset, label = generate_dataset()
# print(label)
alphas, b = smo_simple(dataset, label, C=6, max_iter=40)
print(alphas)
W = get_W(alphas, dataset, label)
print(W, b)
*5.3 使用matplotlib绘图
如果想要可视化计算结果,可以在main函数后添加以下代码:
class1, class2 = [], []
for data in zip(dataset, label):
if data[1] == 1.0:
class1.append(data[0])
elif data[1] == -1.0:
class2.append(data[0])
x11, x12 = zip(*class1)
x21, x22 = zip(*class2)
plt.plot(x11, x12, 'ro')
plt.plot(x21, x22, 'bo')
x = np.linspace(0, 10, 50)
y = -(W[0]*x + b)/W[1]
plt.plot(x, y)
for i in range(len(dataset)):
if alphas[i] > 1e-3:
xi_1, xi_2 = dataset[i][0], dataset[i][1]
plt.scatter(xi_1, xi_2, s=150, c='none', linewidths=1.5, edgecolors='#1f77b4')
x = np.linspace(0, 10, 50)
y = -W[0]/W[1] * x + (xi_2 + W[0]/W[1] * xi_1)
plt.plot(x, y, 'y--')
plt.show()
6 github 地址
https://github.com/chenyr0021/machine_learning_exercises
7 参考文献
- 《机器学习》,周志华
- 《机器学习实战》,Peter Harrington
- 简易解说拉格朗日对偶(Lagrange duality)
- 支持向量机系列(5)——SMO算法解对偶问题
- SVM为什么能够求解对偶问题,求解对偶问题为什么和原问题一样?为什么要求解对偶问题?svm的公式是什么?如果线性不可分怎么办?
- 为什么支持向量机要用拉格朗日对偶算法来解最大化间隔问题?
- 机器学习算法实践-SVM中的SMO算法
- 支持向量机通俗导论(理解SVM的三层境界)