SVM算法原理解析,计算过程和代码实现
SVM(支持向量机)是一种二分类模型,虽然目前运用较少,但是其算法思想中所运用到的凸优化问题,向量空间和优化方法的各种思想对博主的进一步学习有很大的影响。本文是博主综合各方知识和信息,对SVM的算法过程逐步解析,过程可能会有一些错漏,欢迎大家批评指正。
算法解析
SVM的基本思想是在样本的向量空间中寻找一个超平面,使得两类样本被分割在平面的两端;这样的平面理论上有无穷多个,任一个超平面对应的判别模型为感知机。
而为了提高超平面的鲁棒性,需要寻找一个最优的超平面:两侧距离超平面最近的样本点到超平面的距离被最大化了;这种最优的超平面所对应的判别模型即为支持向量机。距离超平面最近的样本点为支持向量。
如图所示,若样本的两类子样本可以超平面所完全分割,则称该样本集线性可分。 在讨论到线性不可分前,算法的上下文都为线性可分样本集。
数学表述
- 样本点: x \textbf{x} x 为n维向量, x = [ x 1 , x 2 , . . . , x n ] \textbf{x}=[x_1,x_2,...,x_n] x=[x1,x2,...,xn]
- x i \textbf{x}_i xi表示第i个样本点的向量
- X \textbf{X} X表示样本空间
- 超平面表述: ω T x + b = 0 \omega^T\textbf{x}+b=0 ωTx+b=0
- 样本点的标签值 y = ± 1 y= \pm1 y=±1
- 样本点距超平面的空间距离为: d i = ∣ ω T x i + b ∣ ∣ ∣ ω ∣ ∣ d_i=\frac{\lvert \omega^T\textbf{x}_i+b \rvert}{\lvert\lvert \omega \rvert\rvert} di=∣∣ω∣∣∣ωTxi+b∣
- 内积的计算采用标准内积: < x i , x j > = x i T x j <x_i,x_j>=x_i^Tx_j <xi,xj>=xiTxj
设样本空间中距离超平面最近的距离为
d
d
d:
d
=
min
X
d
i
d=\min_{\textbf{X}}{d_i}
d=Xmindi
对于样本空间中的所有样本向量而言,满足
d
i
≥
d
d_i \ge d
di≥d。
让我们回到SVM的思想要求,我们想要寻找一个超平面满足:两侧 距离超平面最近 的样本点到超平面的 距离最大化;
将这句话转化为数学表述即为:
max
d
s
.
t
.
d
i
≥
d
\max{d}\\ s.t. \enspace d_i \ge d
maxds.t.di≥d
为了求解这个问题,我们需要进一步的分解和具化:
- 若样本点线性可分,则对于可分离样本的超平面满足:
{ ω T x i + b > 0 ;若 y = 1 ω T x i + b < 0 ;若 y = − 1 \left\{ \begin{aligned} &\omega^Tx_i+b \gt 0 ;若y=1\\ &\omega^Tx_i+b \lt 0 ;若y=-1 \end{aligned} \right. {ωTxi+b>0;若y=1ωTxi+b<0;若y=−1
正负方向是人为设定的,为了计算方便设置成同向;此时对于任意一个样本点有: ∣ ω T x i + b ∣ = y i ( ω T x i + b ) \lvert \omega^T\textbf{x}_i+b \rvert=y_i(\omega^T\textbf{x}_i+b) ∣ωTxi+b∣=yi(ωTxi+b)。
代入样本点距超平面的距离公式: d i = y i ( ω T x i + b ) ∣ ∣ ω ∣ ∣ d_i=\frac{y_i(\omega^T\textbf{x}_i+b)}{\lvert\lvert \omega \rvert\rvert} di=∣∣ω∣∣yi(ωTxi+b) - 由于
ω
\omega
ω的数乘变换不会影响
d
i
d_i
di的结果,所以满足要求超平面的参数
ω
\omega
ω有无数个,即
ω
=
k
∗
ω
\omega=k*\omega
ω=k∗ω; 为了保证结果的唯一性,需要规定一个合适的k值;一般按照最大的最小距离去规范
ω
\omega
ω的长度:
令
∣
∣
ω
∣
∣
∗
d
=
1
令\qquad\lvert\lvert \omega \rvert\rvert *d=1
令∣∣ω∣∣∗d=1
代入距离公式中有: d i = d ∗ y i ( ω T x i + b ) ≥ d = > { y i ( ω T x i + b ) ≥ 1 d = 1 ∣ ∣ ω ∣ ∣ d_i=d*y_i(\omega^T\textbf{x}_i+b) \ge d \\ \\ =>\left\{\begin{aligned} &y_i(\omega^T\textbf{x}_i+b) \ge 1 \\ &d=\frac{1}{\lvert\lvert \omega \rvert\rvert} \end{aligned} \right. di=d∗yi(ωTxi+b)≥d=>⎩ ⎨ ⎧yi(ωTxi+b)≥1d=∣∣ω∣∣1 - 目标问题具化成:
max
1
∣
∣
ω
∣
∣
等价
↔
min
∣
∣
ω
∣
∣
\max{\frac{1}{\lvert\lvert \omega \rvert\rvert}} \enspace \underleftrightarrow{\enspace 等价 \enspace} \enspace \min{\lvert\lvert \omega \rvert\rvert}
max∣∣ω∣∣1
等价min∣∣ω∣∣
由于对目标函数的单增变换不会改变结果的解值:
min ∣ ∣ ω ∣ ∣ 等价 ↔ min 1 2 ∣ ∣ ω ∣ ∣ 2 \min{\lvert\lvert \omega \rvert\rvert} \enspace \underleftrightarrow{\enspace 等价 \enspace} \enspace \min{\frac{1}{2}\lvert\lvert \omega \rvert\rvert^2} min∣∣ω∣∣ 等价min21∣∣ω∣∣2
最终寻找最近样本点的距离最大化的超平面问题转化为:
min
1
2
∣
∣
ω
∣
∣
2
s
.
t
.
y
i
(
ω
T
x
i
+
b
)
≥
1
\min{\frac{1}{2}\lvert\lvert \omega \rvert\rvert^2} \\ s.t. \enspace y_i(\omega^T\textbf{x}_i+b) \ge 1
min21∣∣ω∣∣2s.t.yi(ωTxi+b)≥1
在标准内积的设定下,目标函数为凸二次函数,所以问题转换成 凸二次规划问题 。
KKT条件
函数最值
求解多维函数的最值问题,一般分两种情况考虑:
假设目标函数为
f
(
x
,
y
)
f(x,y)
f(x,y),以二维为例,存在二阶偏导:
- 无条件极值:
(1) 必要条件: { ∂ f ( x , y ) ∂ x = 0 ∂ f ( x , y ) ∂ y = 0 \left\{ \begin{aligned}&\frac{\partial{f(x,y)}}{\partial{x}}=0\\&\frac{\partial{f(x,y)}}{\partial{y}}=0\end{aligned}\right. ⎩ ⎨ ⎧∂x∂f(x,y)=0∂y∂f(x,y)=0
(2) 充分条件: { ∂ 2 f ( x , y ) ∂ x 2 = A ∂ 2 f ( x , y ) ∂ x ∂ y = B ∂ 2 f ( x , y ) ∂ y 2 = C ⟶ { A C − B 2 > 0 ; A > 0 : 极小 , A < 0 : 极大 A C − B 2 < 0 ; 非极值 A C − B 2 = 0 ; 另外讨论 \left\{ \begin{aligned}&\frac{\partial^2{f(x,y)}}{\partial{x^2}}=A\\&\frac{\partial^2{f(x,y)}}{\partial{x}\partial{y}}=B\\ &\frac{\partial^2{f(x,y)}}{\partial{y^2}}=C\end{aligned}\right. \longrightarrow \left\{ \begin{aligned}&AC-B^2>0;A>0:极小,A<0:极大\\&AC-B^2<0;非极值\\&AC-B^2=0;另外讨论\end{aligned}\right. ⎩ ⎨ ⎧∂x2∂2f(x,y)=A∂x∂y∂2f(x,y)=B∂y2∂2f(x,y)=C⟶⎩ ⎨ ⎧AC−B2>0;A>0:极小,A<0:极大AC−B2<0;非极值AC−B2=0;另外讨论
一般通过必要条件寻找可能的极值点,后通过充分条件验证这些可能点;最后比较各个极值点得出到最值。 - 条件极值:
若自变量 ( x , y ) (x,y) (x,y)存在限制条件: g ( x , y ) = 0 g(x,y)=0 g(x,y)=0,则需要构造 拉格朗日函数 : L ( x , y , λ ) = f ( x , y ) + λ g ( x , y ) L(x,y,\lambda)=f(x,y)+\lambda g(x,y) L(x,y,λ)=f(x,y)+λg(x,y)
分别对三个变量求偏导: { ∂ L ( x , y , λ ) ∂ x = 0 ∂ L ( x , y , λ ) ∂ y = 0 ∂ L ( x , y , λ ) ∂ λ = 0 ⟶ { ∇ f ( x , y ) + λ ∇ g ( x , y ) = 0 g ( x , y ) = 0 \left\{ \begin{aligned}&\frac{\partial{L(x,y,\lambda)}}{\partial{x}}=0\\&\frac{\partial{L(x,y,\lambda)}}{\partial{y}}=0\\&\frac{\partial{L(x,y,\lambda)}}{\partial{\lambda}}=0\end{aligned}\right. \longrightarrow \left\{ \begin{aligned}&\nabla{f(x,y)}+\lambda\nabla{g(x,y)}=0\\&g(x,y)=0\end{aligned}\right. ⎩ ⎨ ⎧∂x∂L(x,y,λ)=0∂y∂L(x,y,λ)=0∂λ∂L(x,y,λ)=0⟶{∇f(x,y)+λ∇g(x,y)=0g(x,y)=0
求解得满足限制条件的极值点
一般优化问题的KKT条件
详细的推导过程可以参考 KKT条件,原来如此简单 | 理论+算例实践 - 科研小飞的文章 - 知乎; 写得很精彩,我就省略一些解释过程了。
-
不同情况下的拉格朗日函数统一:
若限制条件为不等式: g ( x , y ) ≤ 0 g(x,y) \le 0 g(x,y)≤0 ,其本质其实是函数值在向量空间中某个区域内的最值(含边界);所以我们可以分成两种情况来讨论: { 区域内的极值 ( 无条件极值,但极值点需要满足区域条件 ) 边界上的极值 ( 条件极值 ) \left\{ \begin{aligned}&区域内的极值(无条件极值,但极值点需要满足区域条件)\\&边界上的极值(条件极值)\end{aligned}\right. {区域内的极值(无条件极值,但极值点需要满足区域条件)边界上的极值(条件极值)
这样问题就回到上面函数最值的方法上; 为了统一形式,不用分类讨论;在区域内的情况下也引入拉格朗日算子 λ \lambda λ;构造拉格朗日函数 L ( x , y , λ ) = f ( x , y ) + λ g ( x , y ) L(x,y,\lambda)=f(x,y)+\lambda g(x,y) L(x,y,λ)=f(x,y)+λg(x,y)。
此时 { λ = 0 , 区域内 g ( x , y ) = 0 , 边界上 ⟶ λ g ( x , y ) = 0 \left\{ \begin{aligned}&\lambda=0,区域内\\&g(x,y)=0, 边界上\end{aligned}\right. \longrightarrow \lambda g(x,y)=0 {λ=0,区域内g(x,y)=0,边界上⟶λg(x,y)=0 -
关于具体优化问题 λ \lambda λ的讨论:
对于最优化问题形如: max f ( x , y ) s . t . g ( x , y ) ≥ 0 \max{f(x,y)}\\ s.t. \enspace g(x,y)\ge 0 maxf(x,y)s.t.g(x,y)≥0
为了保证拉格朗日函数所求值为极大值,需要限制 λ ≥ 0 \lambda\ge0 λ≥0;(最小化问题类似)
如图所示(计算机画图技巧还不好,正在学习中):
设极值点为: ( x ∗ , y ∗ ) (x^*,y^*) (x∗,y∗)- 由于 g ( x , y ) ≥ 0 g(x,y) \ge0 g(x,y)≥0,所以 g ( x , y ) g(x,y) g(x,y)的梯度方向( ∇ g ( x , y ) \nabla{g(x,y)} ∇g(x,y))一定是指向可行域内; 可行域内的函数值大于边界的函数值。
- 由于要最大化 f ( x , y ) f(x,y) f(x,y),所以 f ( x , y ) f(x,y) f(x,y)的梯度方向( ∇ f ( x , y ) \nabla{f(x,y)} ∇f(x,y))一定是指向可行域外; 若指向可行域内,则最大值一定不在边界上,需要修改函数和问题形式。
- 取得极值的条件是, ∇ g ( x ∗ , y ∗ ) \nabla{g(x^*,y^*)} ∇g(x∗,y∗)和 ∇ f ( x ∗ , y ∗ ) \nabla{f(x^*,y^*)} ∇f(x∗,y∗) 共线且方向相反;所以有 ∇ f ( x ∗ , y ∗ ) = − λ ∇ g ( x ∗ , y ∗ ) \nabla{f(x^*,y^*)}=-\lambda\nabla{g(x^*,y^*)} ∇f(x∗,y∗)=−λ∇g(x∗,y∗) (目标函数无法通过移动再获取增量)
所以: ∇ f ( x ∗ , y ∗ ) + λ ∇ g ( x ∗ , y ∗ ) = 0 \nabla{f(x^*,y^*)}+\lambda\nabla{g(x^*,y^*)}=0 ∇f(x∗,y∗)+λ∇g(x∗,y∗)=0
我们可以得到 λ ≥ 0 \lambda\ge0 λ≥0
同时这个条件是在设定问题形式下得到的,为了保障结果的一致性和计算的便利性; 改变问题形式, λ ≥ 0 \lambda\ge0 λ≥0是可以不满足的。 -
最优化问题的KKT条件:
max f ( x , y ) s . t . g ( x , y ) ≥ 0 \max{f(x,y)}\\ s.t. \enspace g(x,y)\ge 0 maxf(x,y)s.t.g(x,y)≥0
综上,关于该问题的KKT条件为:
{ ∇ f ( x ∗ , y ∗ ) + λ ∇ g ( x ∗ , y ∗ ) = 0 λ g ( x ∗ , y ∗ ) = 0 g ( x , y ) ≥ 0 λ ≥ 0 \left\{ \begin{aligned} &\nabla{f(x^*,y^*)}+\lambda\nabla{g(x^*,y^*)}=0\\ &\lambda g(x^*,y^*) = 0\\ &g(x,y) \ge 0\\ &\lambda \ge 0 \end{aligned} \right. ⎩ ⎨ ⎧∇f(x∗,y∗)+λ∇g(x∗,y∗)=0λg(x∗,y∗)=0g(x,y)≥0λ≥0
最大最小间隔问题的KKT条件
让我们回到原问题:
min
1
2
∣
∣
ω
∣
∣
2
s
.
t
.
y
i
(
ω
T
x
i
+
b
)
≥
1
\min{\frac{1}{2}\lvert\lvert \omega \rvert\rvert^2} \\ s.t. \enspace y_i(\omega^T\textbf{x}_i+b) \ge 1
min21∣∣ω∣∣2s.t.yi(ωTxi+b)≥1
其KKT条件为:
{
∇
(
1
2
∣
∣
ω
∣
∣
2
)
+
∑
λ
i
∇
(
1
−
y
i
(
ω
T
x
i
+
b
)
)
=
0
λ
i
(
1
−
y
i
(
ω
T
x
i
+
b
)
)
=
0
(
1
−
y
i
(
ω
T
x
i
+
b
)
)
≤
0
λ
≥
0
\left\{ \begin{aligned} &\nabla{(\frac{1}{2}\lvert\lvert \omega \rvert\rvert^2)}+\sum\lambda_i\nabla{(1-y_i(\omega^T\textbf{x}_i+b))}=0\\ &\lambda_i (1-y_i(\omega^T\textbf{x}_i+b))= 0\\ &(1-y_i(\omega^T\textbf{x}_i+b)) \le 0\\ &\lambda \ge 0 \end{aligned} \right.
⎩
⎨
⎧∇(21∣∣ω∣∣2)+∑λi∇(1−yi(ωTxi+b))=0λi(1−yi(ωTxi+b))=0(1−yi(ωTxi+b))≤0λ≥0
求解这个KKT条件方程组,便可以得到超平面的
ω
∗
x
+
b
∗
\omega^*x+b^*
ω∗x+b∗; 但是直接计算会消耗大量的空间和涉及Hermit矩阵的相关计算(博主尝试去了解原始计算方法的问题,但是能力和精力有限,不得要领,有需求的小伙伴可以参考上传的英文教程)。所以引入了SMO算法,SMO算法通过将原始问题分解为一系列最小二次规划问题,避免了计算海森矩阵的逆或伪逆,从而显著提高了训练速度。
SMO算法
拉格朗日对偶性
- 极小极大化问题:
对于一般的优化问题:
{ min f ( x ) s . t . ϕ i ( x ) ≤ 0 μ j ( x ) = 0 \left\{ \begin{aligned} &\min{f(x)}\\ s.t. \enspace & \phi_i(x)\le0 \\ &\mu_j(x) = 0 \end{aligned} \right. ⎩ ⎨ ⎧s.t.minf(x)ϕi(x)≤0μj(x)=0
引入拉格朗日算子: α , α > 0 ; β \alpha,\alpha>0;\beta α,α>0;β
其拉格朗日函数为: L ( x , α , β ) = f ( x ) + ∑ α i ϕ i ( x ) + ∑ β j μ j ( x ) L(x,\alpha,\beta)=f(x)+\sum{\alpha_i\phi_i(x)}+\sum{\beta_j\mu_j(x)} L(x,α,β)=f(x)+∑αiϕi(x)+∑βjμj(x)
同时构造关于 x x x的函数: θ P ( x ) = max α , β ; α ≥ 0 L ( x , α , β ) \theta_P(x)=\max_{\alpha,\beta;\alpha\ge0}{L(x,\alpha,\beta)} θP(x)=maxα,β;α≥0L(x,α,β)
假如 x x x不满足优化限制条件:
{ 若 ϕ i ( x ) > 0 , 则令 α = + ∞ 若 μ i ( x ) ≠ 0 , 则令 β i μ i ( x ) = + ∞ ⟶ θ P ( x ) = + ∞ \left\{ \begin{aligned} &若\phi_i{(x)}>0,则令\alpha=+\infty\\ &若\mu_i(x) \ne 0, 则令\beta_i\mu_i(x)=+\infty \end{aligned} \right. \longrightarrow\theta_P(x)=+\infty {若ϕi(x)>0,则令α=+∞若μi(x)=0,则令βiμi(x)=+∞⟶θP(x)=+∞
即 θ ( x ) \theta(x) θ(x)总能等于无穷大
假如 x x x满足优化限制条件:
θ P ( x ) = max α , β ; α ≥ 0 L ( x , α , β ) = f ( x ) \theta_P(x)=\max_{\alpha,\beta;\alpha\ge0}{L(x,\alpha,\beta)}=f(x) θP(x)=α,β;α≥0maxL(x,α,β)=f(x)
所以有 θ P ( x ) = { + ∞ , 若 x 不满足限制条件 f ( x ) , 若 x 满足限制条件 \theta_P(x)= \left\{ \begin{aligned} &+\infty,若x不满足限制条件\\ &f(x),若x满足限制条件 \end{aligned} \right. θP(x)={+∞,若x不满足限制条件f(x),若x满足限制条件
所以求 min θ ( x ) \min{\theta(x)} minθ(x)等价于求限制条件下的 min f ( x ) \min{f(x)} minf(x),也就是说原问题等价于:
min x max α , β ; α ≥ 0 L ( x , α , β ) \min_x{\max_{\alpha,\beta;\alpha\ge0}{L(x,\alpha,\beta)}} xminα,β;α≥0maxL(x,α,β)
称为广义拉格朗日函数的极小极大问题。
其最优值为: p ∗ = θ P ( x ∗ ) p^*=\theta_P(x^*) p∗=θP(x∗) - 对偶问题:
构造关于 α , β \alpha,\beta α,β的函数: θ D ( α , β ) = min x L ( x , α , β ) \theta_D(\alpha,\beta)=\min_{x}{L(x,\alpha,\beta)} θD(α,β)=minxL(x,α,β)
则极大化函数: max α , β ; α ≥ 0 θ D = max α , β ; α ≥ 0 min x L ( x , α , β ) \max_{\alpha,\beta;\alpha\ge0}{\theta_D}=\max_{\alpha,\beta;\alpha\ge0}{\min_x{L(x,\alpha,\beta)}} α,β;α≥0maxθD=α,β;α≥0maxxminL(x,α,β)
称为广义拉格朗日函数的极大极小问题。
同时也是原始问题的对偶问题,其最优值为: d ∗ = θ D ( α ∗ , β ∗ ) d^*=\theta_D(\alpha^*,\beta^*) d∗=θD(α∗,β∗)
若两个问题都有最优解,则 d ∗ ≤ p ∗ d^*\le p^* d∗≤p∗
因为 d ∗ d^* d∗的最大值为 f ( x ) f(x) f(x), p ∗ p^* p∗的最小值为 f ( x ) f(x) f(x)
这个性质便叫做弱对偶性(weak duality),对于所有优化问题都成立,即使原始问题非凸。
若原始问题及其对偶问题满足:
- f ( x ) , ϕ ( x ) f(x),\phi(x) f(x),ϕ(x)为凸函数, μ ( x ) \mu(x) μ(x)为仿射函数
- ϕ ( x ) \phi(x) ϕ(x)的不等式约束严格成立
则存在 x ∗ x^* x∗是原问题的解, α ∗ , β ∗ \alpha^*,\beta^* α∗,β∗是对偶问题的解使得(强对偶性):
d ∗ = p ∗ = L ( x ∗ , α ∗ , β ∗ ) d^*=p^*=L(x^*,\alpha^*,\beta^*) d∗=p∗=L(x∗,α∗,β∗)
解 ( x ∗ , α ∗ , β ∗ ) (x^*,\alpha^*,\beta^*) (x∗,α∗,β∗)是原始问题和对偶问题解的充要条件是 KKT条件
最大最小间隔问题的对偶问题:
原始问题转为极小极大值问题为:
min
ω
,
b
max
λ
;
λ
>
0
(
1
2
∣
∣
ω
∣
∣
2
+
∑
λ
i
(
1
−
y
i
(
ω
T
x
i
+
b
)
)
)
\min_{\omega,b}{\max_{\lambda;\lambda>0}{(\frac{1}{2}\lvert\lvert \omega \rvert\rvert^2}+\sum\lambda_i(1-y_i(\omega^T\textbf{x}_i+b)))}
ω,bminλ;λ>0max(21∣∣ω∣∣2+∑λi(1−yi(ωTxi+b)))
由于强对偶性的条件都满足,所以原问题与对偶问题的解相等;对偶问题为:
max
λ
;
λ
>
0
min
ω
,
b
(
1
2
∣
∣
ω
∣
∣
2
+
∑
λ
i
(
1
−
y
i
(
ω
T
x
i
+
b
)
)
)
\max_{\lambda;\lambda>0}{\min_{\omega,b}{(\frac{1}{2}\lvert\lvert \omega \rvert\rvert^2}+\sum\lambda_i(1-y_i(\omega^T\textbf{x}_i+b)))}
λ;λ>0maxω,bmin(21∣∣ω∣∣2+∑λi(1−yi(ωTxi+b)))
创建
λ
\lambda
λ的函数:
θ
(
λ
)
=
min
ω
,
b
(
1
2
∣
∣
ω
∣
∣
2
+
∑
λ
i
(
1
−
y
i
(
ω
T
x
i
+
b
)
)
)
\theta(\lambda)=\min_{\omega,b}{(\frac{1}{2}\lvert\lvert \omega \rvert\rvert^2}+\sum\lambda_i(1-y_i(\omega^T\textbf{x}_i+b)))
θ(λ)=minω,b(21∣∣ω∣∣2+∑λi(1−yi(ωTxi+b)))
现在先求问题:
min
ω
,
b
(
1
2
∣
∣
ω
∣
∣
2
+
∑
λ
i
(
1
−
y
i
(
ω
T
x
i
+
b
)
)
)
;
λ
≥
0
\min_{\omega,b}{(\frac{1}{2}\lvert\lvert \omega \rvert\rvert^2}+\sum\lambda_i(1-y_i(\omega^T\textbf{x}_i+b))); \lambda \ge0
minω,b(21∣∣ω∣∣2+∑λi(1−yi(ωTxi+b)));λ≥0
∂
L
(
ω
,
b
,
λ
)
∂
(
ω
)
=
ω
−
∑
i
=
1
n
λ
i
y
i
x
i
=
0
∂
L
(
ω
,
b
,
λ
)
∂
(
b
)
=
∑
i
=
1
n
λ
i
y
i
=
0
\begin{aligned} &\frac{\partial{L(\omega,b,\lambda)}}{\partial(\omega)}=\omega-\sum_{i=1}^n{\lambda_i y_i \textbf{x}_i}=0\\ &\frac{\partial{L(\omega,b,\lambda)}}{\partial(b)}=\sum_{i=1}^n{\lambda_i y_i}=0 \end{aligned}
∂(ω)∂L(ω,b,λ)=ω−i=1∑nλiyixi=0∂(b)∂L(ω,b,λ)=i=1∑nλiyi=0
得:
ω
∗
=
∑
i
=
1
n
λ
i
y
i
x
i
;
b
\omega^*=\sum_{i=1}^n{\lambda_i y_i \textbf{x}_i};b
ω∗=∑i=1nλiyixi;b 随意,后面通过KKT条件求出(可以向量表示,后期再补充了。)
带入得:
θ
(
λ
)
=
∑
j
=
1
n
λ
j
−
1
2
∑
i
=
1
n
∑
j
=
1
n
λ
i
λ
j
y
i
y
j
(
x
i
T
x
j
)
\theta(\lambda)=\sum_{j=1}^n{\lambda_j}-\frac{1}{2}\sum_{i=1}^n{\sum_{j=1}^{n}{\lambda_i\lambda_j y_i y_j(\textbf{x}_i^T\textbf{x}_j)}}
θ(λ)=j=1∑nλj−21i=1∑nj=1∑nλiλjyiyj(xiTxj)
最大化问题可以表述为:
max
θ
(
λ
)
s
.
t
.
∑
i
=
1
n
λ
i
y
i
=
0
λ
i
≥
0
\begin{aligned} &\max{\theta(\lambda)}\\ s.t. & \sum_{i=1}^n{\lambda_i y_i}=0 \\ & \lambda_i \ge 0 \end{aligned}
s.t.maxθ(λ)i=1∑nλiyi=0λi≥0
SMO算法
具体的求解过程请参考: 机器学习算法实践-SVM中的SMO算法 - 邵正将的文章 - 知乎
在大数据量样本的情况下,上述问题仍旧是个计算庞大的问题。SMO的基本思想类似于坐标上升算法,每次迭代中选取多元函数中的一维,固定其他维度,请当前维度下的极值;经过多次迭代收敛达到优化函数的目的。
但是由于所求问题存在约束:
∑
i
=
1
n
λ
i
y
i
=
0
\sum_{i=1}^n{\lambda_i y_i}=0
∑i=1nλiyi=0,函数自由度仅有(n-1);所以在SMO算法中每次选取两个参数
λ
i
,
λ
j
\lambda_i,\lambda_j
λi,λj,固定其他参数;代入限制条件有(设
i
=
1
,
j
=
2
i=1,j=2
i=1,j=2):
λ
1
y
1
+
λ
2
y
2
=
−
∑
i
=
3
n
λ
i
y
i
=
C
两边
×
y
2
:
λ
2
=
C
y
2
−
λ
1
y
1
y
2
\begin{aligned} &\lambda_1 y_1 + \lambda_2 y_2 = -\sum_{i=3}^{n}\lambda_i y_i = C\\ &两边×y_2:\\ & \lambda_2=Cy_2-\lambda_1 y_1 y_2 \end{aligned}
λ1y1+λ2y2=−i=3∑nλiyi=C两边×y2:λ2=Cy2−λ1y1y2
这样就可以得到第k轮时优化的目标函数:
θ
(
λ
1
,
λ
2
)
=
λ
1
+
λ
2
−
1
2
λ
1
2
y
1
2
(
x
1
T
x
1
)
−
1
2
λ
2
2
y
2
2
(
x
2
T
x
2
)
−
λ
1
λ
2
y
1
y
2
(
x
1
T
x
2
)
−
∑
i
=
3
n
λ
1
λ
i
y
1
y
i
(
x
1
T
x
i
)
−
∑
i
=
3
n
λ
2
λ
i
y
2
y
i
(
x
2
T
x
i
)
+
C
0
\theta(\lambda_1,\lambda_2)=\lambda_1+\lambda_2-\frac{1}{2}\lambda_1^2y_1^2(\textbf{x}_1^T\textbf{x}_1)-\frac{1}{2}\lambda_2^2y_2^2(\textbf{x}_2^T\textbf{x}_2)-\lambda_1\lambda_2y_1y_2(\textbf{x}_1^T\textbf{x}_2)-\sum_{i=3}^n{\lambda_1\lambda_iy_1y_i(\textbf{x}_1^T\textbf{x}_i)}-\sum_{i=3}^n{\lambda_2\lambda_iy_2y_i(\textbf{x}_2^T\textbf{x}_i)}+C_0
θ(λ1,λ2)=λ1+λ2−21λ12y12(x1Tx1)−21λ22y22(x2Tx2)−λ1λ2y1y2(x1Tx2)−i=3∑nλ1λiy1yi(x1Txi)−i=3∑nλ2λiy2yi(x2Txi)+C0
为了表述方便,对式子中的常数做命名规定以简化算式:
K i j = x i T x j ,即两向量的标准内积;有 K j i = K j i K_{ij}=\textbf{x}_i^T\textbf{x}_j,即两向量的标准内积;有K_{ji}=K_{ji} Kij=xiTxj,即两向量的标准内积;有Kji=Kji
v i = ∑ k ! = i , j n λ k y k ( x i T x k ) ; i , j 为每轮所选的参数角标 v_i=\sum_{k!=i,j}^n{\lambda_ky_k(\textbf{x}_i^T\textbf{x}_k)};i,j为每轮所选的参数角标 vi=∑k!=i,jnλkyk(xiTxk);i,j为每轮所选的参数角标
e i = f ( x i ) − y i , 为第 i 个样本的预测误差 e_i=f(\textbf{x}_i)-y_i,为第i个样本的预测误差 ei=f(xi)−yi,为第i个样本的预测误差
λ i o l d 为上一轮参数值 ; λ i n e w 为本轮参数值 \lambda_i^{old}为上一轮参数值;\lambda_i^{new}为本轮参数值 λiold为上一轮参数值;λinew为本轮参数值
代入限制条件,去掉无关常数后:
θ
(
λ
1
)
=
λ
1
−
y
1
y
2
λ
1
−
1
2
K
11
λ
1
2
−
1
2
K
22
λ
1
2
+
K
12
λ
1
2
+
K
22
C
y
1
λ
1
−
K
12
C
y
1
λ
1
−
v
1
y
1
λ
1
+
v
2
y
1
λ
2
\theta(\lambda_1)=\lambda_1-y_1y_2\lambda_1-\frac{1}{2}K_{11}\lambda_1^2-\frac{1}{2}K_{22}\lambda_1^2+K_{12}\lambda_1^2+K_{22}Cy_1\lambda_1-K_{12}Cy_1\lambda_1-v_1y_1\lambda_1+v_2y_1\lambda_2
θ(λ1)=λ1−y1y2λ1−21K11λ12−21K22λ12+K12λ12+K22Cy1λ1−K12Cy1λ1−v1y1λ1+v2y1λ2
对函数求导:
∂
θ
(
λ
1
)
∂
λ
1
=
−
(
K
11
+
K
22
−
2
K
12
)
λ
1
+
K
22
C
y
1
−
K
12
C
y
1
−
(
v
1
−
v
2
)
y
1
=
0
\frac{\partial\theta(\lambda_1)}{\partial{\lambda_1}}= -(K_{11}+K_{22}-2K_{12})\lambda_1+K_{22}Cy_1-K_{12}Cy_1-(v_1-v_2)y_1=0
∂λ1∂θ(λ1)=−(K11+K22−2K12)λ1+K22Cy1−K12Cy1−(v1−v2)y1=0
引入上一轮的预测值
f
(
x
1
)
f(x_1)
f(x1),创建迭代过程:
f
(
x
1
)
=
ω
T
x
1
+
b
;
w
=
∑
j
=
1
n
λ
j
o
l
d
y
j
x
j
=
∑
j
=
1
n
λ
j
o
l
d
y
j
x
j
T
x
1
+
b
=
∑
j
=
1
n
λ
j
o
l
d
y
j
K
j
1
+
b
=
λ
1
y
1
K
11
+
λ
2
y
2
K
21
+
v
1
+
b
f
(
x
2
)
−
f
(
x
1
)
=
v
2
−
v
1
+
K
22
C
−
K
12
C
−
(
K
11
+
K
22
−
2
K
12
)
y
1
λ
1
o
l
d
\begin{aligned} f(x_1)&=\omega^T\textbf{x}_1+b; w=\sum_{j=1}^n{\lambda_j^{old} y_j \textbf{x}_j}\\ &=\sum_{j=1}^n{\lambda_j^{old} y_j \textbf{x}_j^T\textbf{x}_1+b}\\ &=\sum_{j=1}^n{\lambda_j^{old} y_j K_{j1}+b}\\ &=\lambda_1y_1K_{11}+\lambda_2y_2K_{21}+v_1+b\\ f(x_2)-f(x_1)&=v_2-v_1+K_{22}C-K_{12}C-(K_{11}+K_{22}-2K_{12})y_1\lambda_1^{old} \end{aligned}
f(x1)f(x2)−f(x1)=ωTx1+b;w=j=1∑nλjoldyjxj=j=1∑nλjoldyjxjTx1+b=j=1∑nλjoldyjKj1+b=λ1y1K11+λ2y2K21+v1+b=v2−v1+K22C−K12C−(K11+K22−2K12)y1λ1old
将
v
2
−
v
1
v_2-v_1
v2−v1代入极值条件式:
λ
1
n
e
w
=
λ
1
o
l
d
+
y
1
(
e
2
−
e
1
)
K
11
+
K
22
−
2
K
12
\lambda_1^{new}= \lambda_1^{old}+\frac{y_1(e_2-e_1)}{K_{11}+K_{22}-2K_{12}}
λ1new=λ1old+K11+K22−2K12y1(e2−e1)
再由
λ
1
n
e
w
y
1
+
λ
2
n
e
w
y
2
=
λ
1
o
l
d
y
1
+
λ
2
o
l
d
y
2
\lambda_1^{new}y_1+\lambda_2^{new}y_2=\lambda_1^{old}y_1+\lambda_2^{old}y_2
λ1newy1+λ2newy2=λ1oldy1+λ2oldy2得到:
λ
2
n
e
w
=
λ
2
o
l
d
+
y
2
(
e
1
−
e
2
)
K
11
+
K
22
−
2
K
12
\lambda_2^{new}=\lambda_2^{old}+\frac{y_2(e_1-e_2)}{K_{11}+K_{22}-2K_{12}}
λ2new=λ2old+K11+K22−2K12y2(e1−e2)
求出之后注意
λ
\lambda
λ的限制条件:
λ
>
0
\lambda>0
λ>0;
所以当
y
1
=
y
2
y_1=y_2
y1=y2时:
λ
1
n
e
w
=
λ
1
o
l
d
+
λ
2
o
l
d
−
λ
2
n
e
w
≤
λ
1
o
l
d
+
λ
2
o
l
d
\lambda_1^{new}=\lambda_1^{old}+\lambda_2^{old}-\lambda_2^{new}\le\lambda_1^{old}+\lambda_2^{old}
λ1new=λ1old+λ2old−λ2new≤λ1old+λ2old
当
y
1
=
−
y
2
y_1=-y_2
y1=−y2时:
λ
1
n
e
w
=
λ
1
o
l
d
−
λ
2
o
l
d
+
λ
2
n
e
w
≥
λ
1
o
l
d
−
λ
2
o
l
d
\lambda_1^{new}=\lambda_1^{old}-\lambda_2^{old}+\lambda_2^{new}\ge\lambda_1^{old}-\lambda_2^{old}
λ1new=λ1old−λ2old+λ2new≥λ1old−λ2old
Python实现
代码
线性可分,硬间隔
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import random
from sklearn.datasets import make_blobs
from pylab import mpl
#防止乱码
mpl.rcParams["font.sans-serif"] = ["SimHei"]
mpl.rcParams["axes.unicode_minus"] = False
# 创建数据,目前先考虑线性可分的情况
def create_data():
X, y = make_blobs(n_samples=100, centers=2, random_state=0, cluster_std=0.6)
return X,y
# 初始化参数
def init_para(n:int,l:int):
omega = np.zeros([n,1])
b = 0
lm = np.zeros([l,1])
return omega,b,lm
# 预测计算
def f(X,omega,b):
f = np.dot(omega.T,X)+b
return f
# 对于参数的选取采用随机的方式
def choose_j(i,l):
j_l = list(range(l))
j_l = j_l[0:i]+j_l[i+1:-1]
return random.choice(j_l)
# 超平面的计算
def plane(lm,X,y):
tao = np.multiply(lm,y)
omega = np.dot(tao.T,X).T
return omega
X,y = create_data()
y = (y.reshape(len(y),1)-0.5)*2
omega,b,lm = init_para(len(X[0]),len(X))
l,n = X.shape
# 定义最大循环次数
mt = 40
it = 0
# 定义移动阈值,若变化量小于移动阈值,则不改变
ct = 0.00001
while it < mt:
lm_changed = 0
for i in range(l):
j = choose_j(i,l)
K_ii = np.dot(X[i],X[i])
K_jj = np.dot(X[j],X[j])
K_ij = np.dot(X[i],X[j])
z = K_ii+K_jj-2*K_ij
e_i = f(X[i],omega,b)-y[i]
e_j = f(X[j],omega,b)-y[j]
# 大于0的限制条件
if y[i]==y[j]:
lm_i = max(0,min(lm[i]+lm[j],lm[i]+y[i]*(e_j-e_i)/z))
else:
lm_i = max(lm[i]-lm[j],lm[i]+y[i]*(e_j-e_i)/z,0)
lm_j = lm[j]+y[i]*y[j]*(lm[i]-lm_i)
if lm_i < 0 or lm_j <0 :
continue
if abs(lm_i-lm[i])<ct:
continue
# 更新超平面参数
b_i = -e_i-y[i]*K_ii*(lm_i-lm[i])-y[j]*K_ij*(lm_j-lm[j])+b
b_j = -e_j-y[j]*K_jj*(lm_j-lm[j])-y[i]*K_ij*(lm_i-lm[i])+b
if lm_i>0 and lm_j>0:
b = (b_i+b_j)/2
elif lm_i >0 :
b = b_i
elif lm_j >0 :
b = b_j
lm[i] = lm_i
lm[j] = lm_j
omega = plane(lm,X,y)
lm_changed +=1
# 如果没有更新 那么迭代—+1
if lm_changed == 0:
it+=1
else:
it = 0
lmd = pd.DataFrame(lm,columns=['lambda'])
support_vectors = X[lmd[lmd['lambda']>ct].index]
# 画图
plt.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.Paired)
# 画出支持向量
plt.scatter(support_vectors[:, 0], support_vectors[:, 1], s=100,
facecolors='none', edgecolors='k',label='支持向量')
# 画出决策边界和超平面
ax = plt.gca()
xlim = ax.get_xlim()
ylim = ax.get_ylim()
# 创建网格
xx = np.linspace(xlim[0], xlim[1], 30)
yy = np.linspace(ylim[0], ylim[1], 30)
YY, XX = np.meshgrid(yy, xx)
Z = omega[0]*XX+omega[1]*YY+b
# 绘制等高线图来表示决策边界和间隔
ax.contour(XX, YY, Z, colors='k', levels=[-1, 0, 1], alpha=0.5,
linestyles=['--', '-', '--'])
plt.xlabel('X')
plt.ylabel('Y')
plt.legend(loc='upper left')
plt.title('SVM示意图')
plt.show()
print('***')
运行结果
(因为数据是随意生成的,所以每次结果不同):
软间隔
软间隔引入
若干样本中存在少数异常点,导致样本线性不可分;
那么可以允许异常点不满足限制条件:
1
−
y
i
(
ω
T
x
i
+
b
)
≤
0
1-y_i(\omega^T\textbf{x}_i+b) \le0
1−yi(ωTxi+b)≤0;
加入一个松弛变量
ε
i
≥
0
→
1
−
y
i
(
ω
T
x
i
+
b
)
−
ε
i
≤
0
\varepsilon_i\ge0 \rightarrow 1-y_i(\omega^T\textbf{x}_i+b)-\varepsilon_i \le0
εi≥0→1−yi(ωTxi+b)−εi≤0;
同时在目标函数中加入惩罚项:
C
∑
ε
i
→
1
2
∣
∣
ω
∣
∣
2
+
C
∑
ε
i
C\sum{\varepsilon_i} \rightarrow \frac{1}{2}\lvert\lvert \omega \rvert\rvert^2+C\sum{\varepsilon_i}
C∑εi→21∣∣ω∣∣2+C∑εi;
其中C为大于0的常数,为异常样本的惩罚程度;有C越大,允许的松弛变量越小
软间隔条件下的KKT条件为:
{
∇
(
1
2
∣
∣
ω
∣
∣
2
+
C
∑
ε
i
)
+
∑
λ
i
∇
(
1
−
y
i
(
ω
T
x
i
+
b
)
−
ε
i
)
+
∑
μ
i
∇
(
−
ε
i
)
=
0
λ
i
(
1
−
y
i
(
ω
T
x
i
+
b
)
−
ε
i
)
=
0
μ
i
ε
i
=
0
(
1
−
y
i
(
ω
T
x
i
+
b
)
−
ε
i
)
≤
0
−
ε
i
≤
0
λ
≥
0
;
u
≥
0
\left\{ \begin{aligned} &\nabla{(\frac{1}{2}\lvert\lvert \omega \rvert\rvert^2+C\sum{\varepsilon_i})}+\sum\lambda_i\nabla{(1-y_i(\omega^T\textbf{x}_i+b)-\varepsilon_i)}+\sum{\mu_i\nabla{(-\varepsilon_i)}}=0\\ &\lambda_i (1-y_i(\omega^T\textbf{x}_i+b)-\varepsilon_i)= 0\\ &\mu_i \varepsilon_i= 0\\ &(1-y_i(\omega^T\textbf{x}_i+b)-\varepsilon_i) \le 0\\ &-\varepsilon_i \le 0\\ &\lambda \ge 0; \,u \ge 0 \\ \end{aligned} \right.
⎩
⎨
⎧∇(21∣∣ω∣∣2+C∑εi)+∑λi∇(1−yi(ωTxi+b)−εi)+∑μi∇(−εi)=0λi(1−yi(ωTxi+b)−εi)=0μiεi=0(1−yi(ωTxi+b)−εi)≤0−εi≤0λ≥0;u≥0
原始的极小极大值问题为:
min
ω
,
b
,
ε
max
λ
;
λ
>
0
(
1
2
∣
∣
ω
∣
∣
2
+
C
∑
ε
i
+
∑
λ
i
(
1
−
y
i
(
ω
T
x
i
+
b
)
−
ε
i
)
+
∑
μ
i
(
−
ε
i
)
)
\min_{\omega,b,\varepsilon}{\max_{\lambda;\lambda>0}{(\frac{1}{2}\lvert\lvert \omega \rvert\rvert^2+C\sum{\varepsilon_i}}+\sum\lambda_i(1-y_i(\omega^T\textbf{x}_i+b)-\varepsilon_i)+\sum\mu_i(-\varepsilon_i))}
ω,b,εminλ;λ>0max(21∣∣ω∣∣2+C∑εi+∑λi(1−yi(ωTxi+b)−εi)+∑μi(−εi))
转化为对偶问题为:
max
λ
;
λ
>
0
min
ω
,
b
,
ε
(
1
2
∣
∣
ω
∣
∣
2
+
C
∑
ε
i
+
∑
λ
i
(
1
−
y
i
(
ω
T
x
i
+
b
)
−
ε
i
)
+
∑
μ
i
(
−
ε
i
)
)
\max_{\lambda;\lambda>0}{\min_{\omega,b,\varepsilon}{(\frac{1}{2}\lvert\lvert \omega \rvert\rvert^2+C\sum{\varepsilon_i}}+\sum\lambda_i(1-y_i(\omega^T\textbf{x}_i+b)-\varepsilon_i)+\sum\mu_i(-\varepsilon_i))}
λ;λ>0maxω,b,εmin(21∣∣ω∣∣2+C∑εi+∑λi(1−yi(ωTxi+b)−εi)+∑μi(−εi))
设函数
θ
(
λ
)
=
min
ω
,
b
,
ε
(
1
2
∣
∣
ω
∣
∣
2
+
C
∑
ε
i
+
∑
λ
i
(
1
−
y
i
(
ω
T
x
i
+
b
)
−
ε
i
)
+
∑
μ
i
(
−
ε
i
)
)
\theta(\lambda)=\min_{\omega,b,\varepsilon}{(\frac{1}{2}\lvert\lvert \omega \rvert\rvert^2+C\sum{\varepsilon_i}}+\sum\lambda_i(1-y_i(\omega^T\textbf{x}_i+b)-\varepsilon_i)+\sum\mu_i(-\varepsilon_i))
θ(λ)=minω,b,ε(21∣∣ω∣∣2+C∑εi+∑λi(1−yi(ωTxi+b)−εi)+∑μi(−εi))
求解:
∂
L
(
ω
,
b
,
λ
,
ε
)
∂
(
ω
)
=
ω
−
∑
i
=
1
n
λ
i
y
i
x
i
=
0
∂
L
(
ω
,
b
,
λ
,
ε
)
∂
(
b
)
=
∑
i
=
1
n
λ
i
y
i
=
0
∂
L
(
ω
,
b
,
λ
,
ε
)
∂
(
ε
i
)
=
C
−
λ
i
−
μ
i
=
0
\begin{aligned} &\frac{\partial{L(\omega,b,\lambda,\varepsilon)}}{\partial(\omega)}=\omega-\sum_{i=1}^n{\lambda_i y_i \textbf{x}_i}=0\\ &\frac{\partial{L(\omega,b,\lambda,\varepsilon)}}{\partial(b)}=\sum_{i=1}^n{\lambda_i y_i}=0\\ &\frac{\partial{L(\omega,b,\lambda,\varepsilon)}}{\partial(\varepsilon_i)}=C-\lambda_i-\mu_i=0 \end{aligned}
∂(ω)∂L(ω,b,λ,ε)=ω−i=1∑nλiyixi=0∂(b)∂L(ω,b,λ,ε)=i=1∑nλiyi=0∂(εi)∂L(ω,b,λ,ε)=C−λi−μi=0
代入得:
θ
(
λ
)
=
∑
j
=
1
n
λ
j
−
1
2
∑
i
=
1
n
∑
j
=
1
n
λ
i
λ
j
y
i
y
j
(
x
i
T
x
j
)
s
.
t
.
λ
≥
0
;
∑
i
=
1
n
λ
i
y
i
=
0
;
C
−
λ
i
−
μ
i
=
0
\theta(\lambda)=\sum_{j=1}^n{\lambda_j}-\frac{1}{2}\sum_{i=1}^n{\sum_{j=1}^{n}{\lambda_i\lambda_j y_i y_j(\textbf{x}_i^T\textbf{x}_j)}} \\ s.t. \enspace \lambda \ge0;\sum_{i=1}^n{\lambda_i y_i}=0;C-\lambda_i-\mu_i=0
θ(λ)=j=1∑nλj−21i=1∑nj=1∑nλiλjyiyj(xiTxj)s.t.λ≥0;i=1∑nλiyi=0;C−λi−μi=0
在软间隔的情况下,优化目标不变,但是多了一个限制条件;该限制条件限制了
λ
\lambda
λ的上限,即
λ
≤
C
\lambda \le C
λ≤C;
所以解的形式和硬间隔类似,但是根据限制条件:
λ
>
0
\lambda>0
λ>0和
λ
≤
C
\lambda \le C
λ≤C;
当
y
1
=
y
2
y_1=y_2
y1=y2时:
λ
1
o
l
d
+
λ
2
o
l
d
−
C
≤
λ
1
n
e
w
=
λ
1
o
l
d
+
λ
2
o
l
d
−
λ
2
n
e
w
≤
λ
1
o
l
d
+
λ
2
o
l
d
\lambda_1^{old}+\lambda_2^{old}-C\le\lambda_1^{new}=\lambda_1^{old}+\lambda_2^{old}-\lambda_2^{new}\le\lambda_1^{old}+\lambda_2^{old}
λ1old+λ2old−C≤λ1new=λ1old+λ2old−λ2new≤λ1old+λ2old
当
y
1
=
−
y
2
y_1=-y_2
y1=−y2时:
λ
1
o
l
d
−
λ
2
o
l
d
≤
λ
1
n
e
w
=
λ
1
o
l
d
−
λ
2
o
l
d
+
λ
2
n
e
w
≤
λ
1
o
l
d
−
λ
2
o
l
d
+
C
\lambda_1^{old}-\lambda_2^{old}\le\lambda_1^{new}=\lambda_1^{old}-\lambda_2^{old}+\lambda_2^{new}\le\lambda_1^{old}-\lambda_2^{old}+C
λ1old−λ2old≤λ1new=λ1old−λ2old+λ2new≤λ1old−λ2old+C
在python实现时,也只用增加
λ
\lambda
λ的可行域判断即可
Python实现
引入惩罚程度C,C越大对边界内的点容忍越小,约接近硬间隔的情况;
在硬间隔的代码中修改
λ
\lambda
λ的取值限制即可;
# 大于0的限制条件
if y[i]==y[j]:
lm_i = max(0,min(lm[i]+lm[j],lm[i]+y[i]*(e_j-e_i)/z,C),lm[i]+lm[j]-C)
else:
lm_i = max(lm[i]-lm[j],min(lm[i]+y[i]*(e_j-e_i)/z,lm[i]-lm[j]+C,C),0)
lm_j = lm[j]+y[i]*y[j]*(lm[i]-lm_i)
决策边界内的向量仍对超平面的生成起作用,所以仍是支持向量:
线性不可分
基本思想
关于向量空间和其映射的特征空间以及核应该是一套完整的关于非线性映射的逻辑,博主也还是一知半解,大家可以参考一下这个文章的讲解:支持向量机的通俗导论。
若整个数据集分割线是个完全的曲面,在样本集所在的向量空间中,我们无法使用平面来分割样本集;很自然的就有拟合一个非线性关系的想法,这等价于应用一个固定的非线性映射,将数据映射到特征空间,在特征空间中使用线性学习器,即映射
x
→
ϕ
(
x
)
x\rightarrow\phi(\textbf{x})
x→ϕ(x)。
那么分类器的函数就可以表示为:
f
(
x
)
=
ω
ϕ
(
x
)
+
b
f(x)=\omega\phi(\textbf{x})+b
f(x)=ωϕ(x)+b
根据之前的求解过程,在求对偶问题的时候
ϕ
(
x
)
\phi(x)
ϕ(x)这个位置是没动的(因为是对
ω
\omega
ω的偏导),所以非线性问题的对偶问题可以表示成:
θ
(
λ
)
=
∑
j
=
1
n
λ
j
−
1
2
∑
i
=
1
n
∑
j
=
1
n
λ
i
λ
j
y
i
y
j
<
ϕ
(
x
i
)
,
ϕ
(
x
j
)
>
\theta(\lambda)=\sum_{j=1}^n{\lambda_j}-\frac{1}{2}\sum_{i=1}^n{\sum_{j=1}^{n}{\lambda_i\lambda_j y_i y_j<\phi(\textbf{x}_i),\phi(\textbf{x}_j)>}}
θ(λ)=j=1∑nλj−21i=1∑nj=1∑nλiλjyiyj<ϕ(xi),ϕ(xj)>
其中,
<
ϕ
(
x
i
)
,
ϕ
(
x
j
)
>
<\phi(\textbf{x}_i),\phi(\textbf{x}_j)>
<ϕ(xi),ϕ(xj)>是特征空间中向量的内积。
这样,如果我们选择一个正确的映射
ϕ
(
x
)
\phi(x)
ϕ(x),那么我们可以依然可以计算出对应的曲面。
如:
如图所示的圆集,其理想分割面应该是个二次曲线:
a
1
x
1
2
+
a
2
x
2
2
+
a
3
x
1
+
a
4
x
2
+
a
5
x
1
x
2
+
a
6
=
0
a_1x_1^2+a_2x_2^2+a_3x_1+a_4x_2+a_5x_1x_2+a_6=0
a1x12+a2x22+a3x1+a4x2+a5x1x2+a6=0
这样其实是将原本的二维向量
x
\textbf{x}
x上升成五维向量
z
\textbf{z}
z:
x
=
{
x
1
x
2
→
z
=
{
x
1
2
x
2
2
x
1
x
2
x
1
x
2
\textbf{x}=\left\{\begin{aligned} &x_1\\ &x_2 \end{aligned}\right. \rightarrow \textbf{z}=\left\{\begin{aligned} &x_1^2\\ &x_2^2\\ &x_1\\ &x_2\\ &x_1x_2 \end{aligned}\right.
x={x1x2→z=⎩
⎨
⎧x12x22x1x2x1x2
这样问题就解决了?
但是
- 不同数据集的理想分割曲面是不同的,每次选择映射也不同,同时新增的数据也可能改变映射关系;
- 单纯从圆集角度来看,都会从二维上升成五维;如果是三维球体,那么需要上升的维度就是十九维,依次类推,很容易引起维度爆炸,大大增加训练难度;
核函数的引入
如果我们有个函数
K
(
x
1
,
x
2
)
=
<
ϕ
(
x
1
)
,
ϕ
(
x
2
)
>
K(\textbf{x}_1,\textbf{x}_2)=<\phi(\textbf{x}_1),\phi(\textbf{x}_2)>
K(x1,x2)=<ϕ(x1),ϕ(x2)>,那么在对偶问题中我们就不需要显性的映射关系
ϕ
(
x
)
\phi(\textbf{x})
ϕ(x),而可以 在当前维度下计算核函数的值 ,进而求解对偶问题的解:
θ
(
λ
)
=
∑
j
=
1
n
λ
j
−
1
2
∑
i
=
1
n
∑
j
=
1
n
λ
i
λ
j
y
i
y
j
K
(
x
1
,
x
2
)
\theta(\lambda)=\sum_{j=1}^n{\lambda_j}-\frac{1}{2}\sum_{i=1}^n{\sum_{j=1}^{n}{\lambda_i\lambda_j y_i y_jK(\textbf{x}_1,\textbf{x}_2)}}
θ(λ)=j=1∑nλj−21i=1∑nj=1∑nλiλjyiyjK(x1,x2)
但是这里有个问题:虽然我们求出了
λ
\lambda
λ的值,但是我们怎么使用这个结果进行分类呢?
在线性可分的情况下,可以显性的求出
f
(
x
)
=
ω
x
+
b
f(\textbf{x})=\omega\textbf{x}+b
f(x)=ωx+b ,但是由于引入显性映射的复杂性,所以很难求出线性的曲面;
这时候需要回顾对偶问题的
ω
\omega
ω值:
ω
∗
=
∑
i
=
1
n
λ
i
y
i
x
i
\omega^*=\sum_{i=1}^n{\lambda_i y_i \textbf{x}_i}
ω∗=∑i=1nλiyixi
将
x
\textbf{x}
x换为映射后:
ω
∗
=
∑
i
=
1
n
λ
i
y
i
ϕ
(
x
i
)
→
f
(
x
)
=
∑
i
=
1
n
λ
i
y
i
<
ϕ
(
x
i
)
,
ϕ
(
x
)
>
+
b
→
f
(
x
)
=
∑
i
=
1
n
λ
i
y
i
K
(
x
i
,
x
)
+
b
\begin{aligned} &\omega^*=\sum_{i=1}^n{\lambda_i y_i \phi(\textbf{x}_i}) \\ &\rightarrow f(\textbf{x})=\sum_{i=1}^n{\lambda_i y_i <\phi(\textbf{x}_i}),\phi(\textbf{x})>+b\\ &\rightarrow f(\textbf{x})=\sum_{i=1}^n{\lambda_i y_i K(\textbf{x}_i,\textbf{x})}+b \end{aligned}
ω∗=i=1∑nλiyiϕ(xi)→f(x)=i=1∑nλiyi<ϕ(xi),ϕ(x)>+b→f(x)=i=1∑nλiyiK(xi,x)+b
所以,通过核函数我们不用显性求出映射,同时分割平面是高维度的平面,也不用显性的在当前空间算出。
至此,其实我们可以把线性可分情况也纳入统一的表述中,其用的核函数是一种线性核,即
K
(
x
1
,
x
2
)
=
<
x
1
,
x
2
>
K(\textbf{x}_1,\textbf{x}_2)=<\textbf{x}_1,\textbf{x}_2>
K(x1,x2)=<x1,x2>
在选用其他核的情况下,我们只需要将线性可分情况下的核计算换成对应核即可。
关于核的一些原理和推导,博主还是管中窥豹,大家有兴趣的可以参考我上传的英语教程。
常用的核函数:
- 多项式核: K ( x 1 , x 2 ) = ( x 1 T x 2 ) d K(\textbf{x}_1,\textbf{x}_2)=(\textbf{x}_1^T\textbf{x}_2)^d K(x1,x2)=(x1Tx2)d
- 高斯核: K ( x 1 , x 2 ) = e − ∣ ∣ x 1 − x 2 ∣ ∣ 2 2 σ 2 K(\textbf{x}_1,\textbf{x}_2)=e^{-\frac{||\textbf{x}_1-\textbf{x}_2||^2}{2\sigma^2}} K(x1,x2)=e−2σ2∣∣x1−x2∣∣2
高斯核将原始空间映射为无限维空间,合理的选择参数 σ \sigma σ则可以将任意的数据映射为线性可分,所以应用较为广泛。
Python实现
根据上述线性可分硬间隔代码中调整核值计算即可,同时将 ω \omega ω的输出删掉,改成直接预测值的输出。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import random
from sklearn.datasets import make_blobs,make_circles
from matplotlib.patches import Wedge
from pylab import mpl
mpl.rcParams["font.sans-serif"] = ["SimHei"]
mpl.rcParams["axes.unicode_minus"] = False
# 创建数据,目前先考虑线性可分的情况
def create_data():
# 改成圆集
X,y=make_circles(n_samples=100,factor=0.1,noise=0.1)
return X,y
# 核函数的计算
def K(X1,X2):
theta = 2
K=np.exp(-np.sum((X1-X2)**2)/theta)
return K
# 初始化参数
def init_para(n:int,l:int):
omega = np.zeros([n,1])
b = 0
lm = np.zeros([l,1])
return omega,b,lm
def f(lm,X,y,x,b):
tao = np.multiply(lm,y)
k = np.apply_along_axis(K,1,X,x)
f = np.dot(tao.T,k)+b
return f[0]
# 对于参数的选取采用随机的方式
def choose_j(i,l):
j_l = list(range(l))
j_l = j_l[0:i]+j_l[i+1:-1]
return random.choice(j_l)
X,y = create_data()
y = (y.reshape(len(y),1)-0.5)*2
# 利用高斯升维
omega,b,lm = init_para(len(X[0]),len(X))
l,n = X.shape
# 定义最大循环次数
mt = 40
it = 0
C=0.6
# 定义移动阈值,若变化量小于移动阈值,则不改变
ct = 0.00001
while it < mt:
lm_changed = 0
for i in range(l):
j = choose_j(i,l)
K_ii = K(X[i],X[i])
K_jj = K(X[j],X[j])
K_ij = K(X[i],X[j])
z = K_ii+K_jj-2*K_ij
e_i = f(lm,X,y,X[i],b)-y[i]
e_j = f(lm,X,y,X[j],b)-y[j]
# 大于0的限制条件
if y[i]==y[j]:
lm_i = max(0,min(lm[i]+lm[j],lm[i]+y[i]*(e_j-e_i)/z,C),lm[i]+lm[j]-C)
else:
lm_i = max(lm[i]-lm[j],min(lm[i]+y[i]*(e_j-e_i)/z,lm[i]-lm[j]+C,C),0)
lm_j = lm[j]+y[i]*y[j]*(lm[i]-lm_i)
if lm_i < 0 or lm_j <0 :
continue
if abs(lm_i-lm[i])<ct:
continue
# 更新超平面参数
b_i = -e_i-y[i]*K_ii*(lm_i-lm[i])-y[j]*K_ij*(lm_j-lm[j])+b
b_j = -e_j-y[j]*K_jj*(lm_j-lm[j])-y[i]*K_ij*(lm_i-lm[i])+b
if lm_i>0 and lm_j>0:
b = (b_i+b_j)/2
elif lm_i >0 :
b = b_i
elif lm_j >0 :
b = b_j
lm[i] = lm_i
lm[j] = lm_j
lm_changed +=1
# 如果没有更新 那么迭代—+1
if lm_changed == 0:
it+=1
else:
it = 0
lmd = pd.DataFrame(lm,columns=['lambda'])
support_vectors = X[lmd[lmd['lambda']>ct].index]
# 画图
plt.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.Paired)
plt.scatter(support_vectors[:, 0], support_vectors[:, 1], s=100,
facecolors='none', edgecolors='k',label='支持向量')
ax = plt.gca()
xlim = ax.get_xlim()
ylim = ax.get_ylim()
# 创建网格
xx = np.linspace(xlim[0], xlim[1], 30)
yy = np.linspace(ylim[0], ylim[1], 30)
YY, XX = np.meshgrid(yy, xx)
# 创建预测Z值
Z = []
for i in range(len(XX)):
_Z = []
for j in range(len(XX[i])):
Point = np.array([XX[i][j],YY[i][j]])
_Z.append(f(lm,X,y,Point,b))
Z.append(_Z)
ax.contour(XX, YY, Z, colors='k', levels=[-1, 0,1], alpha=0.5,
linestyles=['--', '-', '--'])
plt.xlabel('X')
plt.ylabel('Y')
plt.legend(loc='upper left')
plt.title('SVM示意图')
plt.show()
运行结果如下:
目前代码运行相对较慢,但是结果可收敛,等博主有空再回来优化下 T.T
感谢大家花时间阅读!
文章总体结构参考: 【机器学习】支持向量机 SVM(非常详细) - 阿泽的文章 - 知乎