0

机器学习算法系列(十五)-软间隔支持向量机算法(Soft-margin Support Vector Machin...

 2 years ago
source link: https://segmentfault.com/a/1190000041358058
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
neoserver,ios ssh client

阅读本文需要的背景知识点:硬间隔支持向量机、松弛变量、一丢丢编程知识

  前面一节我们介绍了一种最基础的支持向量机模型——硬间隔支持向量机,该模型能对线性可分的数据集进行分类,但现实中的数据集往往是线性不可分的,那么这一节我们来介绍支持向量机中的第二种——软间隔支持向量机1(Soft-margin Support Vector Machine),来解决上面说的数据集线性不可分的问题。

二、模型介绍

  我们先来看下硬间隔支持向量机的原始模型如下:

max⁡w,b12wTw s.t. yi(wTxi+b)≥1i=1,2,⋯ ,N \begin{array}{} \underset{w, b}{\max} \frac{1}{2} w^Tw \\ \text { s.t. } \quad y_{i}\left(w^{T} x_{i}+b\right) \geq 1 \quad i=1,2, \cdots, N \end{array} w,bmax​21​wTw s.t. yi​(wTxi​+b)≥1i=1,2,⋯,N​

<center>式2-1</center>

  约束条件 y(wx + b) 大于等于 1 是为了使得所有样本点都在正确的分类下,这也是为什么称为硬间隔的原因。而现在数据集无法用一个超平面完全分开,这时就需要允许部分数据集不满足上述约束条件。

修改上面的代价函数,加上分类错误的样本点的惩罚项,其中 1(x) 在前面章节中提到过,即指示函数2(indicator function),当满足下面不等式时函数返回 1,不满足时函数返回 0。同时通过常数 C 来控制惩罚力度,其中 C > 0,可以看到当 C 取无限大时,会迫使每一个样本点分类正确,等同于硬间隔支持向量机。

min⁡w,b12wTw+C∑i=1N1xi(yi(wTxi+b)<1) s.t. C>0 \begin{array}{} \underset{w, b}{\min} \frac{1}{2} w^{T} w+C \sum_{i=1}^{N} 1_{x_{i}}\left(y_{i}\left(w^{T} x_{i}+b\right) \lt 1\right) \\ \text { s.t. } \quad C>0 \end{array} w,bmin​21​wTw+C∑i=1N​1xi​​(yi​(wTxi​+b)<1) s.t. C>0​

<center>式2-2</center>

  式2-2中的指示函数既不是凸函数又不是连续函数,处理起来比较麻烦,这时可以用max函数来替换,形式如下:

min⁡w,b12wTw+C∑i=1Nmax⁡(0,1−yi(wTxi+b)) s.t. C>0 \begin{array}{} \underset{w, b}{\min} \frac{1}{2} w^{T} w+C \sum_{i=1}^{N} \max\left(0, 1 - y_{i}\left(w^{T} x_{i}+b\right)\right) \\ \text { s.t. } \quad C>0 \end{array} w,bmin​21​wTw+C∑i=1N​max(0,1−yi​(wTxi​+b)) s.t. C>0​

<center>式2-3</center>

  这时再引入松弛变量3(slack variable)ξ,同时加上如下的条件,进一步将 max 函数去掉,式2-3就是软间隔支持向量机的原始模型:

min⁡w,b,ξ12wTw+C∑i=1Nξi s.t. C>0ξi≥0ξi≥1−yi(wTxi+b)i=1,2,⋯ ,N \begin{array}{} \underset{w, b, \xi}{\min} \frac{1}{2} w^{T} w+C \sum_{i=1}^{N} \xi_{i} \\ \text { s.t. } \quad C>0 \quad \xi_{i} \geq 0 \quad \xi_{i} \geq 1-y_{i}\left(w^{T} x_{i}+b\right) \quad i=1,2, \cdots, N \end{array} w,b,ξmin​21​wTw+C∑i=1N​ξi​ s.t. C>0ξi​≥0ξi​≥1−yi​(wTxi​+b)i=1,2,⋯,N​

<center>式2-4</center>

  与硬间隔支持向量机的处理方式一样,同样对原始模型用拉格朗日乘子法进行转换:

(1)使用拉格朗日乘子法,引入两类拉格朗日参数 λ、μ,得到拉格朗日函数

(2)拉格朗日函数对 w 求偏导数并令其为零

(3)同硬间隔支持向量机一样得到 w 的解析解

(4)拉格朗日函数对 b 求偏导数并令其为零

(5)拉格朗日函数对 ξ 求偏导数并令其为零

(6)得到 C = λ + μ

L(w,b,ξ,λ,μ)=12wTw+C∑i=1Nξi+∑i=1Nλi(1−ξi−yi(wTxi+b))−∑i=1Nμiξi(1)∂L∂w=w−∑i=1Nλiyixi=0(2)w=∑i=1Nλiyixi(3)∂L∂b=∑i=1Nλiyi=0(4)∂L∂ξi=C−λi−μi=0(5)C=λi+μi(6) \begin{aligned} L(w, b, \xi, \lambda, \mu) &=\frac{1}{2} w^{T} w+C \sum_{i=1}^{N} \xi_{i}+\sum_{i=1}^{N} \lambda_{i}\left(1-\xi_{i}-y_{i}\left(w^{T} x_{i}+b\right)\right)-\sum_{i=1}^{N} \mu_{i} \xi_{i} & (1)\\ \frac{\partial L}{\partial w} &=w-\sum_{i=1}^{N} \lambda_{i} y_{i} x_{i}=0 & (2)\\ w &=\sum_{i=1}^{N} \lambda_{i} y_{i} x_{i} & (3)\\ \frac{\partial L}{\partial b} &=\sum_{i=1}^{N} \lambda_{i} y_{i}=0 & (4)\\ \frac{\partial L}{\partial \xi_{i}} &=C-\lambda_{i}-\mu_{i}=0 & (5)\\ C &=\lambda_{i}+\mu_{i} & (6) \end{aligned} L(w,b,ξ,λ,μ)∂w∂L​w∂b∂L​∂ξi​∂L​C​=21​wTw+Ci=1∑N​ξi​+i=1∑N​λi​(1−ξi​−yi​(wTxi​+b))−i=1∑N​μi​ξi​=w−i=1∑N​λi​yi​xi​=0=i=1∑N​λi​yi​xi​=i=1∑N​λi​yi​=0=C−λi​−μi​=0=λi​+μi​​(1)(2)(3)(4)(5)(6)​

<center>式2-5</center>

  将式2-5得到的结果带回拉格朗日函数中:

(1)拉格朗日函数

(2)带入后展开括号

(3)可以看到(2)式中第2、5项相互抵消,第3、8项相互抵消,第7项为零,合并后得

L(w,b,ξ,λ,μ)=12wTw+C∑i=1Nξi+∑i=1Nλi(1−ξi−yi(wTxi+b))−∑i=1Nμiξi(1)=12∑i=1N∑j=1NλiλjyiyjxiTxi+∑i=1Nλiξi+∑i=1Nμiξi+∑i=1Nλi−∑i=1Nλiξi−∑i=1N∑j=1NλiλjyiyjxiTxi−∑i=1Nλiyib−∑i=1Nμiξi(2)=∑i=1Nλi−12∑i=1N∑j=1NλiλjyiyjxiTxi(3) \begin{aligned} L(w, b, \xi, \lambda, \mu) &=\frac{1}{2} w^{T} w+C \sum_{i=1}^{N} \xi_{i}+\sum_{i=1}^{N} \lambda_{i}\left(1-\xi_{i}-y_{i}\left(w^{T} x_{i}+b\right)\right)-\sum_{i=1}^{N} \mu_{i} \xi_{i} & (1)\\ &=\frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \lambda_{i} \lambda_{j} y_{i} y_{j} x_{i}^{T} x_{i}+\sum_{i=1}^{N} \lambda_{i} \xi_{i}+\sum_{i=1}^{N} \mu_{i} \xi_{i}+\sum_{i=1}^{N} \lambda_{i}-\sum_{i=1}^{N} \lambda_{i} \xi_{i}-\sum_{i=1}^{N} \sum_{j=1}^{N} \lambda_{i} \lambda_{j} y_{i} y_{j} x_{i}^{T} x_{i}-\sum_{i=1}^{N} \lambda_{i} y_{i} b-\sum_{i=1}^{N} \mu_{i} \xi_{i} & (2)\\ &=\sum_{i=1}^{N} \lambda_{i}-\frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \lambda_{i} \lambda_{j} y_{i} y_{j} x_{i}^{T} x_{i} & (3) \end{aligned} L(w,b,ξ,λ,μ)​=21​wTw+Ci=1∑N​ξi​+i=1∑N​λi​(1−ξi​−yi​(wTxi​+b))−i=1∑N​μi​ξi​=21​i=1∑N​j=1∑N​λi​λj​yi​yj​xiT​xi​+i=1∑N​λi​ξi​+i=1∑N​μi​ξi​+i=1∑N​λi​−i=1∑N​λi​ξi​−i=1∑N​j=1∑N​λi​λj​yi​yj​xiT​xi​−i=1∑N​λi​yi​b−i=1∑N​μi​ξi​=i=1∑N​λi​−21​i=1∑N​j=1∑N​λi​λj​yi​yj​xiT​xi​​(1)(2)(3)​

<center>式2-6</center>

  在来看下拉格朗日乘子参数的条件:

(1)拉格朗日乘子 μ 的限制条件

(2)两边同时加上拉格朗日乘子 λ

(3)结合式2-5(6)的结果得到

μi≥0(1)λi+μi≥λi(2)C≥λi(3) \begin{aligned} \mu_i &\ge 0 & (1) \\ \lambda_i + \mu_i &\ge \lambda_i & (2) \\ C &\ge \lambda_i & (3) \end{aligned} μi​λi​+μi​C​≥0≥λi​≥λi​​(1)(2)(3)​

<center>式2-7</center>

  同硬间隔支持向量机一样,引入 KKT 条件如下:

{∇w,b,ξL(w,b,ξ,λ,μ)=0λi≥0yi(wTxi+b)−1+ξi≥0λi(yi(wTxi+b)−1+ξi)=0μi≥0ξi≥0μiξi=0 \left\{\begin{aligned} \nabla_{w, b, \xi} L(w, b, \xi, \lambda, \mu) &=0 \\ \lambda_{i} & \geq 0 \\ y_{i}\left(w^{T} x_{i}+b\right)-1+\xi_{i} & \geq 0 \\ \lambda_{i}\left(y_{i}\left(w^{T} x_{i}+b\right)-1+\xi_{i}\right) &=0 \\ \mu_{i} & \geq 0 \\ \xi_{i} & \geq 0 \\ \mu_{i} \xi_{i} &=0 \end{aligned}\right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧​∇w,b,ξ​L(w,b,ξ,λ,μ)λi​yi​(wTxi​+b)−1+ξi​λi​(yi​(wTxi​+b)−1+ξi​)μi​ξi​μi​ξi​​=0≥0≥0=0≥0≥0=0​

<center>式2-8</center>

  经过运用拉格朗日乘子法后,得到了原模型的拉格朗日对偶模型并且需要满足如上所示的 KKT 条件:

max⁡λ∑i=1Nλi−12∑i=1N∑j=1NλiλjyiyjxiTxj s.t. ∑i=1Nλiyi=00≤λi≤Ci=1,2,⋯ ,N \begin{array}{} \underset{\lambda}{\max} \sum_{i=1}^{N} \lambda_{i}-\frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \lambda_{i} \lambda_{j} y_{i} y_{j} x_{i}^{T} x_{j} \\ \text { s.t. } \quad \sum_{i=1}^{N} \lambda_{i} y_{i}=0 \quad 0 \leq \lambda_{i} \leq C \quad i=1,2, \cdots, N \end{array} λmax​∑i=1N​λi​−21​∑i=1N​∑j=1N​λi​λj​yi​yj​xiT​xj​ s.t. ∑i=1N​λi​yi​=00≤λi​≤Ci=1,2,⋯,N​

<center>式2-9</center>

三、算法步骤

  将上面式2-9的对偶模型与硬间隔支持向量机的对偶模型相互比较后,会发现两者的差别就在于拉格朗日乘子的约束条件不同,前者只多了一个惩罚因子 C 的上界,所以其求解的算法与硬间隔支持向量机的大致相同,只有如下两个地方不一样。

  算法中一个有改动的地方是新 λ_b 的上下界的计算:

(1)、(2)所有 λ 都需要大于等于零并且小于等于 C

(3)分情况讨论 λ_b 的限制条件

(4)综合(1)式得到最终变量 λ_b 的限制条件

0≤λbnew ≤C(1)0≤λaold +yayb(λbold −λbnew )≤C(2)⇒{λaold +λbold −C≤λbnew ≤λaold +λbold ,yayb=+1λbold −λaold ≤λbnew ≤C−λaold +λbold ,yayb=−1(3)⇒{max⁡(0,λaold +λbold −C)≤λbnew ≤min⁡(C,λaold +λbold ),yayb=+1max⁡(0,λbold −λaold )≤λbnew ≤min⁡(C,C−λaold +λbold ),yayb=−1(4) \begin{aligned} &0 \leq \lambda_{b}^{\text {new }} \leq C & (1)\\ &0 \leq \lambda_{a}^{\text {old }}+y_{a} y_{b}\left(\lambda_{b}^{\text {old }}-\lambda_{b}^{\text {new }}\right) \leq C & (2)\\ \Rightarrow&\left\{\begin{array}{} \lambda_{a}^{\text {old }}+\lambda_{b}^{\text {old }}-C \leq \lambda_{b}^{\text {new }} \leq \lambda_{a}^{\text {old }}+\lambda_{b}^{\text {old }}, \quad y_{a} y_{b}=+1 \\ \lambda_{b}^{\text {old }}-\lambda_{a}^{\text {old }} \leq \lambda_{b}^{\text {new }} \leq C-\lambda_{a}^{\text {old }}+\lambda_{b}^{\text {old }}, \quad y_{a} y_{b}=-1 \end{array}\right. & (3)\\ \Rightarrow&\left\{\begin{array}{} \max \left(0, \lambda_{a}^{\text {old }}+\lambda_{b}^{\text {old }}-C\right) \leq \lambda_{b}^{\text {new }} \leq \min \left(C, \lambda_{a}^{\text {old }}+\lambda_{b}^{\text {old }}\right), \quad y_{a} y_{b}=+1 \\ \max \left(0, \lambda_{b}^{\text {old }}-\lambda_{a}^{\text {old }}\right) \leq \lambda_{b}^{\text {new }} \leq \min \left(C, C-\lambda_{a}^{\text {old }}+\lambda_{b}^{\text {old }}\right), \quad y_{a} y_{b}=-1 \end{array}\right. & (4) \end{aligned} ⇒⇒​0≤λbnew ​≤C0≤λaold ​+ya​yb​(λbold ​−λbnew ​)≤C{λaold ​+λbold ​−C≤λbnew ​≤λaold ​+λbold ​,ya​yb​=+1λbold ​−λaold ​≤λbnew ​≤C−λaold ​+λbold ​,ya​yb​=−1​{max(0,λaold ​+λbold ​−C)≤λbnew ​≤min(C,λaold ​+λbold ​),ya​yb​=+1max(0,λbold ​−λaold ​)≤λbnew ​≤min(C,C−λaold ​+λbold ​),ya​yb​=−1​​(1)(2)(3)(4)​

<center>式3-1</center>

  算法中另一个有改动的地方是 KKT 条件的判断:

(1)考虑 λ 等于 0 的情况

(2)根据式2-5 中(6)的结果可知 μ 等于 C

(3)由于 μ 等于 C,根据式2-8的 KKT 条件可知 ξ 等于 0

(4)将 ξ 等于 0 带入 KKT 条件得

λ=0(1)μ=C(2)ξ=0(3)yi(wTxi+b)−1≥0(4) \begin{aligned} \lambda &=0 & (1)\\ \mu &=C & (2)\\ \xi &=0 & (3)\\ y_{i}\left(w^{T} x_{i}+b\right)-1 & \geq 0 & (4) \end{aligned} λμξyi​(wTxi​+b)−1​=0=C=0≥0​(1)(2)(3)(4)​

<center>式3-2</center>

(1)考虑 λ 等于 C 的情况

(2)根据式2-5 中(6)的结果可知 μ 等于 0

(3)由于 λ 不等于 0,根据式2-8的 KKT 条件可知

(4)由于 ξ 大于等于 0 结合(3)式可得

λ=C(1)μ=0(2)yi(wTxi+b)−1+ξ=0(3)yi(wTxi+b)−1≤0(4) \begin{aligned} \lambda &=C & (1)\\ \mu &=0 & (2)\\ y_{i}\left(w^{T} x_{i}+b\right)-1 + \xi &=0 & (3)\\ y_{i}\left(w^{T} x_{i}+b\right)-1 & \leq 0 & (4) \end{aligned} λμyi​(wTxi​+b)−1+ξyi​(wTxi​+b)−1​=C=0=0≤0​(1)(2)(3)(4)​

<center>式3-3</center>

(1)考虑 λ 介于 0 与 C 之间的情况

(2)根据式2-5 中(6)的结果可知 μ 也介于 0 与 C 之间

(3)根据式2-8的 KKT 条件可知

(4)由于 λ 不等于 0,根据式2-8的 KKT 条件并结合(3)式可得

0<λ<C(1)0<μ<C(2)ξ=0(3)yi(wTxi+b)−1=0(4) \begin{aligned} 0 \lt \lambda &\lt C & (1)\\ 0 \lt \mu &\lt C & (2)\\ \xi &=0 & (3)\\ y_{i}\left(w^{T} x_{i}+b\right)-1 & = 0 & (4) \end{aligned} 0<λ0<μξyi​(wTxi​+b)−1​<C<C=0=0​(1)(2)(3)(4)​

<center>式3-4</center>

  综合上面式3-2、式3-3、式3-4可得如下的新限制条件:

{yi(wTxi+b)−1≥0,λ=0yi(wTxi+b)−1=0,0<λ<Cyi(wTxi+b)−1≤0,λ=C \left\{\begin{array}{} y_{i}\left(w^{T} x_{i}+b\right)-1 \geq 0, & \lambda=0 \\ y_{i}\left(w^{T} x_{i}+b\right)-1=0, & 0 \lt \lambda \lt C \\ y_{i}\left(w^{T} x_{i}+b\right)-1 \leq 0, & \lambda=C \end{array}\right. ⎩⎪⎨⎪⎧​yi​(wTxi​+b)−1≥0,yi​(wTxi​+b)−1=0,yi​(wTxi​+b)−1≤0,​λ=00<λ<Cλ=C​

<center>式3-5</center>

  算法步骤请参考硬间隔支持向量机的第三大节和下面的代码实现。

四、代码实现

使用 Python 实现

import numpy as np

class SMO:
    """
    软间隔支持向量机
    序列最小优化算法实现(Sequential minimal optimization/SMO)
    """

    def __init__(self, X, y):
        # 训练样本特征矩阵(N * p)
        self.X = X
        # 训练样本标签向量(N * 1)
        self.y = y
        # 拉格朗日乘子向量(N * 1)
        self.alpha = np.zeros(X.shape[0])
        # 误差向量,默认为负的标签向量(N * 1)
        self.errors = -y
        # 偏移量 
        self.b = 0
        # 权重向量(p * 1)
        self.w = np.zeros(X.shape[1])
        # 代价值
        self.cost = -np.inf

    def fit(self, C = 1, tol = 1e-4):
        """
        算法来自 John C. Platt 的论文
        https://www.microsoft.com/en-us/research/uploads/prod/1998/04/sequential-minimal-optimization.pdf
        """
        # 更新变化次数
        numChanged = 0
        # 是否检查全部
        examineAll = True
        while numChanged > 0 or examineAll:
            numChanged = 0
            if examineAll:
                for idx in range(X.shape[0]):
                    numChanged += self.update(idx, C)
            else:
                for idx in range(X.shape[0]):
                    if self.alpha[idx] <= 0:
                        continue
                    numChanged += self.update(idx, C)
            if examineAll:
                examineAll = False
            elif numChanged == 0:
                examineAll = True
            # 计算代价值
            cost = self.calcCost()
            if self.cost > 0:
                # 当代价值的变化小于容许的范围时结束算法
                rate = (cost - self.cost) / self.cost
                if rate <= tol:
                    break
            self.cost = cost

    def update(self, idx, C = 1):
        """
        对下标为 idx 的拉格朗日乘子进行更新
        """
        X = self.X
        y = self.y
        alpha = self.alpha
        # 检查当前拉格朗日乘子是否满足KKT条件,满足条件则直接返回 0
        if self.checkKKT(idx, C):
            return 0
        if len(alpha[(alpha != 0)]) > 1:
            # 按照|E1 - E2|最大的原则寻找第二个待优化的拉格朗日乘子的下标
            jdx = self.selectJdx(idx)
            # 对下标为 idx、jdx 的拉格朗日乘子进行更新,当成功更新时直接返回 1
            if self.updateAlpha(idx, jdx, C):
                return 1
        # 当未更新成功时,遍历不为零的拉格朗日乘子进行更新
        for jdx in range(X.shape[0]):
            if alpha[jdx] == 0:
                continue
            # 对下标为 idx、jdx 的拉格朗日乘子进行更新,当成功更新时直接返回 1
            if self.updateAlpha(idx, jdx, C):
                return 1
        # 当依然没有没有更新成功时,遍历为零的拉格朗日乘子进行更新
        for jdx in range(X.shape[0]):
            if alpha[jdx] != 0:
                continue
            # 对下标为 idx、jdx 的拉格朗日乘子进行更新,当成功更新时直接返回 1
            if self.updateAlpha(idx, jdx, C):
                return 1
        # 依然没有更新时返回 0
        return 0

    def selectJdx(self, idx):
        """
        寻找第二个待优化的拉格朗日乘子的下标
        """
        errors = self.errors
        if errors[idx] > 0:
            # 当误差项大于零时,选择误差向量中最小值的下标
            return np.argmin(errors)
        elif errors[idx] < 0:
            # 当误差项小于零时,选择误差向量中最大值的下标
            return np.argmax(errors)
        else:
            # 当误差项等于零时,选择误差向量中最大值和最小值的绝对值最大的下标
            minJdx = np.argmin(errors)
            maxJdx = np.argmax(errors)
            if max(np.abs(errors[minJdx]), np.abs(errors[maxJdx])) - errors[minJdx]:
                return minJdx
            else:
                return maxJdx

    def calcB(self):
        """
        计算偏移量
        分别计算每一个拉格朗日乘子不为零对应的偏移量后取其平均值
        """
        X = self.X
        y = self.y
        alpha = self.alpha
        alpha_gt = alpha[alpha > 0]
        # 拉格朗日乘子向量中不为零的数量
        alpha_gt_len = len(alpha_gt)
        # 全部为零时直接返回 0
        if alpha_gt_len == 0:
            return 0
        # b = y - Wx,具体算法请参考文章中的说明
        X_gt = X[alpha > 0]
        y_gt = y[alpha > 0]
        alpha_gt_y = np.array(np.multiply(alpha_gt, y_gt)).reshape(-1, 1)
        s = np.sum(np.multiply(alpha_gt_y, X_gt), axis=0)
        return np.sum(y_gt - X_gt.dot(s)) / alpha_gt_len

    def calcCost(self):
        """
        计算代价值
        按照文章中的算法计算即可
        """
        X = self.X
        y = self.y
        alpha = self.alpha
        cost = 0
        for idx in range(X.shape[0]):
            for jdx in range(X.shape[0]):
                cost = cost + (y[idx] * y[jdx] * X[idx].dot(X[jdx]) * alpha[idx] * alpha[jdx])
        return np.sum(alpha) - cost / 2

    def checkKKT(self, idx, C = 1):
        """
        检查下标为 idx 的拉格朗日乘子是否满足 KKT 条件
        1. alpha >= 0
        2. alpha <= C
        3. y * f(x) - 1 >= 0
        4. alpha * (y * f(x) - 1) = 0
        """
        y = self.y
        errors = self.errors
        alpha = self.alpha
        r = errors[idx] * y[idx]
        if (alpha[idx] > 0 and alpha[idx] < C and r == 0) or (alpha[idx] == 0 and r > 0) or (alpha[idx] == C and r < 0):
            return True
        return False

    def calcE(self):
        """
        计算误差向量
        E = f(x) - y
        """
        X = self.X
        y = self.y
        alpha = self.alpha
        alpha_y = np.array(np.multiply(alpha, y)).reshape(-1, 1)
        errors = X.dot(X.T).dot(alpha_y).T[0] + self.b - y
        return errors

    def calcU(self, idx, jdx, C = 1):
        """
        计算拉格朗日乘子上界,使两个待优化的拉格朗日乘子同时大于等于0
        按照文章中的算法计算即可
        """
        y = self.y
        alpha = self.alpha
        if y[idx] * y[jdx] == 1:
            return max(0, alpha[jdx] + alpha[idx] - C)
        else:
            return max(0.0, alpha[jdx] - alpha[idx])

    def calcV(self, idx, jdx, C = 1):
        """
        计算拉格朗日乘子下界,使两个待优化的拉格朗日乘子同时大于等于0
        按照文章中的算法计算即可
        """
        y = self.y
        alpha = self.alpha
        if y[idx] * y[jdx] == 1:
            return min(C, alpha[jdx] + alpha[idx])
        else:
            return min(C, C + alpha[jdx] - alpha[idx])

    def updateAlpha(self, idx, jdx, C = 1):
        """
        对下标为 idx、jdx 的拉格朗日乘子进行更新
        按照文章中的算法计算即可
        """
        if idx == jdx:
            return False
        X = self.X
        y = self.y
        alpha = self.alpha
        errors = self.errors
        # idx 的误差项
        Ei = errors[idx]
        # jdx 的误差项
        Ej = errors[jdx]
        Kii = X[idx].dot(X[idx])
        Kjj = X[jdx].dot(X[jdx])
        Kij = X[idx].dot(X[jdx])
        # 计算 K
        K = Kii + Kjj - 2 * Kij
        oldAlphaIdx = alpha[idx]
        oldAlphaJdx = alpha[jdx]
        # 计算 jdx 的新拉格朗日乘子
        newAlphaJdx = oldAlphaJdx + y[jdx] * (Ei - Ej) / K
        U = self.calcU(idx, jdx, C)
        V = self.calcV(idx, jdx, C)
        if newAlphaJdx < U:
            # 当新值超过上界时,修改其为上界
            newAlphaJdx = U
        if newAlphaJdx > V:
            # 当新值低于下界时,修改其为下界
            newAlphaJdx = V
        if oldAlphaJdx == newAlphaJdx:
            # 当新值与旧值相等时,判断为未更新,直接返回
            return False
        # 计算 idx 的新拉格朗日乘子
        newAlphaIdx = oldAlphaIdx + y[idx] * y[jdx] * (oldAlphaJdx - newAlphaJdx)
        # 更新权重向量
        self.w = self.w + y[idx] * (newAlphaIdx - oldAlphaIdx) * X[idx] + y[jdx] * (newAlphaJdx - oldAlphaJdx) * X[jdx]
        # 更新拉格朗日乘子向量
        alpha[idx] = newAlphaIdx
        alpha[jdx] = newAlphaJdx
        oldB = self.b
        # 重新计算偏移量
        self.b = self.calcB()
        # 重新计算误差向量
        eps = np.finfo(np.float32).eps
        newErrors = []
        for i in range(X.shape[0]):
            newError = errors[i] + y[idx] * (newAlphaIdx - oldAlphaIdx) * X[idx].dot(X[i]) + y[jdx] * (newAlphaJdx - oldAlphaJdx) * X[jdx].dot(X[i]) - oldB + self.b
            if np.abs(newError) < eps:
                newError = 0
            newErrors.append(newError)
        self.errors = newErrors
        return True

五、第三方库实现

scikit-learn4 实现

from sklearn.svm import SVC

svc = SVC(kernel = "linear", C = 1)
# 拟合数据
svc.fit(X, y)
# 权重系数
w = svc.coef_
# 截距
b = svc.intercept_
print("w", w, "b", b)

六、动画演示

  下图展示了一个线性不可分的数据集,其中红色表示标签值为 -1 的样本点,蓝色代表标签值为 1 的样本点:

13.png
<center>图6-1</center>

  下图展示了通过软间隔支持向量机拟合数据后的结果,其中浅红色的区域为预测值为 -1 的部分,浅蓝色的区域则为预测值为 1 的部分:

14.png
<center>图6-2</center>

七、思维导图

15.jpeg
<center>图7-1</center>

八、参考文献

完整演示请点击这里

注:本文力求准确并通俗易懂,但由于笔者也是初学者,水平有限,如文中存在错误或遗漏之处,恳请读者通过留言的方式批评指正

本文首发于——AI导图,欢迎关注


About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK