机器学习：支持向量机(SVM)与Python实现第(四)篇

本文链接：https://blog.csdn.net/liweibin1994/article/details/77567498

本文是支持向量机(SVM)与Python实现系列的第四篇，主要介绍SMO算法，该算法是解决SVM对偶问题的有效方法。文章首先回顾坐标上升算法，然后详细阐述SMO算法的优化过程，包括如何选择参数变量、计算L和H、求解α2以及更新b的值。最后，提供了基于Python的SVM分类器实现代码片段。

摘要由CSDN通过智能技术生成

前言

最近看了Andrew Ng的机器学习视频中的支持向量机，视频的内容比较浅显，没有深入解释支持向量机中的数学原理。但是对于一个比较执着于知道为什么的人，笔者还是去网上查找了有关支持向量机原理以及实现的相关资料。在查找的过程中，笔者发现支持向量机的内容还是蛮多的，于是笔者根据自己的理解，并且参考了一些相关资料，最终写下了支持向量机的四篇博客。
机器学习：支持向量机(SVM)与Python实现第(一)篇——此篇主要介绍了分类间隔，引入SVM。
机器学习：支持向量机(SVM)与Python实现第(二)篇——此篇主要介绍了使用拉格朗日乘子来简化SVM问题的优化。
机器学习：支持向量机(SVM)与Python实现第(三)篇——此篇主要介绍非线性分类(核函数)以及松弛变量。
机器学习：支持向量机(SVM)与Python实现第(四)篇——此篇主要介绍SMO算法并用python实现了简单的SVM分类器。

前面讲了这么多，终于要用Python来真正实现SVM训练数据集了！在这里，我们使用的是SMO(swquential minimal optimization)算法，一种高效的解决SVM中的对偶问题的方法。

在我们介绍SMO算法之前，我们先来了解一下另外一个算法，它是SMO算法基础。这个算法就是坐标上升(下降)算法。

坐标上升算法(coordinate ascent)

假如我们现在要来解决下面这个无限制的优化问题：

m a x α W (α 1, α 2, . . ., α m)

$\begin{aligned} \underset{\alpha}{max}\ W(\alpha_1,\alpha_2,...,\alpha_m) \end{aligned}$ 其中，W是一个关于参数

α $\alpha$ 的函数。那么我可以使用坐标上升算法来解决，即：
这里写图片描述

在内层循环中，我们每一次都只优化一个

αi $\alpha_i$ ，对于其他的参数，都保持不变。上面的示例中，优化是按照

α1,α2,...,αm $\alpha_1,\alpha_2,...,\alpha_m$ 的顺序进行优化的，事实上，我们也可以选择那个能使W增加最快的

α $\alpha$ 。下面这张图是对这个算法的直观理解：
这里写图片描述

图中显示的是一个二维的情况，每一次都会选择一个方向进行优化，而不是像梯度下降去找下降最快的方向。其实它相当于把梯度下降分成两步来做。

SMO

现在我们回到我们的SMO算法上来，我们上一篇博文最后得到了下面的优化问题：

m a x α s . t W (α) = \sum i = 1 m α i - 1 2 \sum i, j = 1 m y (i) y (j) α i α j ⟨ x (i), x (j) ⟩ (1) 0 \leq α i \leq C, i = 1, . . ., m (2) \sum i = 1 m α i y (i) = 0 (3)

$\begin{aligned} \underset{\alpha}{max}\ &W(\alpha)=\sum_{i=1}^{m}\alpha_i-\frac{1}{2}\sum_{i,j=1}^{m}y^{(i)}y^{(j)}\alpha_i\alpha_j \left \langle x^{(i)}, x^{(j)}\right \rangle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1)\\ s.t\ \ & 0\leq \alpha_i\leq C,\ i=1,...,m \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (2)\\ &\sum_{i=1}^{m}\alpha_iy^{(i)}=0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (3) \end{aligned}$ 假设我们有一组

α $\alpha$ 满足了方程(2)和(3)，我们先固定住

α2,...,αm $\alpha_2,...,\alpha_m$ ，只单独对

α1 $\alpha_1$ 进行优化。这样有意义吗？

显然没有。原因在于方程(3)是限制条件，一旦 $\alpha_2,...,\alpha_m$ 都确定了，那么 $\alpha_1$ 也就确定了，所以这样是优化不了的。因此，如果我们想使用坐标上升算法，我们应该每次优化两个参数变量。这就是SMO算法：

这里写图片描述
一般来说，为了判断是否收敛了，可以检查KKT条件是否收敛到了一定程度，我们可以用一个收敛参数来判断，一般来说这个参数的值在0.001到0.01之间。

接下来我们就来看看这个算法最主要的过程。首先我们选择了两个参数 $\alpha_1和\alpha_2$ (其他的参数固定不变)，根据方程(3)我们可以得到：

α 1 y (1) + α 2 y (2) = - \sum i = 3 m α i y (i)

$\begin{aligned} \alpha_1y^{(1)}+\alpha_2y^{(2)}=-\sum_{i=3}^{m}\alpha_iy^{(i)} \end{aligned}$ 由于

α3,...,αm $\alpha_3,...,\alpha_m$ 都是固定的，所以方程右边其实是一个常数，我们记为

ζ $\zeta$ 也就是：

α 1 y (1) + α 2 y (2) = ζ (4)

$\begin{aligned} \alpha_1y^{(1)}+\alpha_2y^{(2)}=\zeta \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (4) \end{aligned}$
方程(4)其实就是二维平面上的一条直线，我们可以画出图像：
这里写图片描述

从不等式(2)我们可以知道

α1和α2 $\alpha_1和\alpha_2$ 都是大于0，小于C的。所以它们的取值范围是在这个矩形范围内。同时，它们又满足方程(4)的条件，所以是在这条直线上。在这种限制下，其实

α2 $\alpha_2$ 的取值范围就为

L≤α2≤H $L\leq \alpha_2\leq H$ 。在上图这个例子中，L=0。通过方程(4)，我们也可以用

α2 $\alpha_2$ 来表示

α1 $\alpha_1$ ：

α 1 = (ζ - α 2 y (2)) y (1)

$\begin{aligned} \alpha_1=(\zeta-\alpha_2y^{(2)})y^{(1)} \end{aligned}$
根据方程(1)的形式，而且现在变量只有

α2 $\alpha_2$ ，你可以发现现在

W(α) $W(\alpha)$ 是关于

α2 $\alpha_2$ 的二次函数(即可以写成

W=aα22+bα2+c $W=a\alpha_2^2+b\alpha_2+c$ 的形式)。如果我们不去考虑上面说的限制(即

L≤α2≤H $L\leq \alpha_2\leq H$ )，那么我们可以轻易地算出最大值，只需要对其求导并让导数为0解出

α2 $\alpha_2$ 即可。假设我们现在解出了

α2 $\alpha_2$ 的结果为

αnew,unclipped2 $\alpha_2^{new, unclipped}$ ，将限制条件考虑进来，我们就可以得到：

α n e w 2 = ⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ H i f α n e w, u n c l i p p e d 2 > H α n e w, u n c l i p p e d 2 i f L \leq α n e w, u n c l i p p e d 2 \leq H L i f α n e w, u n c l i p p e d 2 < L

$\alpha_2^{new} = \left\{\begin{matrix} \begin{aligned} &H\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ if\ \alpha_2^{new, unclipped}>H\\ &\alpha_2^{new, unclipped}\ \ \ if\ L\leq \alpha_2^{new, unclipped}\leq H\\ &L\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ if\ \alpha_2^{new, unclipped}<L \end{aligned} \end{matrix}\right.$ 这样获得了

αnew2 $\alpha_2^{new}$ ，就可以解出

αnew1 $\alpha_1^{new}$ 当然这里

α1和α2 $\alpha_1和\alpha_2$ 是随便选的，数据量小的时候没什么问题，数据量大了就可能会很慢。其实为了让算法更高效，可以参考Platt的paper，那里有更有效的

α $\alpha$ 选择方法。

选择了 $\alpha_1和\alpha_2$ 之后，我们要计算L和H，L和H的计算如下：
这里写图片描述

接下来我们就要去求解 $\alpha_2$ 了。再次列出我们的优化问题是：

m a x α s . t W (α) = \sum i = 1 m α i - 1 2 \sum i, j = 1 m y (i) y (j) α i α j ⟨ x (i), x (j) ⟩ (1) 0 \leq α i \leq C, i = 1, . . ., m (2) \sum i = 1 m α i y (i) = 0 (3)

W (α 1, α 2) = 1 2 α 21 y 21 K 11 + 1 2 α 22 y 22 K 22 + α 1 α 2 y 1 y 2 K 12 + α 1 y 1 v 1 + α 2 y 2 v 2 - α 1 - α 2 + W c o n s t a n t (5)

$W(\alpha_1,\alpha_2) = \frac{1}{2}\alpha_1^2y_1^2K_{11}+\frac{1}{2}\alpha_2^2y_2^2K_{22}+\alpha_1\alpha_2y_1y_2K_{12}+\alpha_1y_1v_1+\alpha_2y_2v_2-\alpha_1-\alpha_2+W_{constant}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (5)$
其中，

K i j = K (x (i)), x (j)) v i = \sum j = 3 m y j α * j K i j = u i + b * - y 1 α * 1 K 1 j - y 2 α * 2 K 2 j u i = \sum i = 1 m α i y (i) ⟨ x (i), x ⟩ + b

$\begin{aligned} &K_{ij}=K(x^{(i)}),x^{(j)})\\ &v_i=\sum_{j=3}^{m}y_j\alpha_j^*K_{ij} = u_i+b^*-y_1\alpha_1^*K_{1j}-y_2\alpha_2^*K_{2j}\\ &u_i=\sum_{i=1}^{m}\alpha_iy^{(i)}\left \langle x^{(i)},x \right \rangle+b \end{aligned}$ 这里，带”*”号表示的是某次迭代前的原始值，是常数。

Wconstant $W_{constant}$ 也是常数。利用

α1和α2 $\alpha_1和\alpha_2$ 之间的关系，可以将式(5)化为只含有

α2 $\alpha_2$ 的式子，对其求导并让导数为0。化简最终可以得到：

α 2 (K 11 + K 22 - 2 K 12) = α * 2 (K 11 + K 22 - 2 K 12) + y 2 ((u 1 - y 1) - (u 2 - y 2))

$\begin{aligned} \alpha_2(K_{11}+K_{22}-2K_{12})=\alpha_2^*(K_{11}+K_{22}-2K_{12})+y_2((u_1-y_1)-(u_2-y_2)) \end{aligned}$ 也就是：

α 2 = α * 2 + y 2 ( E 1 - E 2 ) η

$\alpha_2=\alpha_2^*+\frac{y_2(E_1-E_2)}{\eta }$ 其中，

E k = u k - y (k) η = K 11 + K 22 - 2 K 12

$E_k=u_k-y^{(k)}\\ \eta = K_{11}+K_{22}-2K_{12}$ 与之前说的一样，计算出了新的

α2 $\alpha_2$ 的值，要进行约束。

α n e w 2 = ⎧ ⎩ ⎨ ⎪ ⎪ H i f α 2 > H α 2 i f L \leq α 2 \leq H L i f α 2 < L

$\alpha_2^{new} = \left\{\begin{matrix} \begin{aligned} &H\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ if\ \alpha_2>H\\ &\alpha_2\ \ \ \ \ \ if\ L\leq \alpha_2\leq H\\ &L\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ if\ \alpha_2<L \end{aligned} \end{matrix}\right.$ 根据

α1和α2 $\alpha_1和\alpha_2$ 的关系就可以计算出新的

α1 $\alpha_1$ 。

接下来是更新b的值了。b的更新有以下几种情况：
这里写图片描述

这样我们就完成了一次更新，多次迭代之后收敛到一定程度就可以了。下面是根据上面的推导完成的代码：(代码与所用的数据也放在了github上了)

# -*- coding: utf-8 -*-
from numpy import *
import operator
import time

def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    with open(fileName) as fr:
        for line in fr.readlines():
            lineArr = line.strip().split()
            dataMat.append([float(lineArr[0]), float(lineArr[1])])
            labelMat.append(float(lineArr[2]))
    return dataMat, labelMat

# The first one, i , is the index of our first alpha,
# and m is the total number of alphas. A value is randomly chosen
# and returned as long as it’s not equal to the input i.
def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j

#剪辑大于H或小于L的α值。
#sometime alpha may larger or smaller than H or L,so we have to constrain it
def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

def smoSimple(dataMatIn,classLabels,C,toler,maxIter):
    dataMatrix=mat(dataMatIn)
    labelMat=mat(classLabels).transpose()
    b=0
    m,n=shape(dataMatrix)
    alphas = mat(zeros((m,1)))
    iter=0
    while (iter<maxIter):
        alphaPairsChanged=0
        for i in range(m):
            fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[i])
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                j = selectJrand(i,m)#随便选一个
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                if (labelMat[i] != labelMat[j]):
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                if L==H: print "L==H"; continue
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T -dataMatrix[i,:]*dataMatrix[i,:].T -dataMatrix[j,:]*dataMatrix[j,:].T
                if eta >= 0: print "eta>=0"; continue
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                alphas[j] = clipAlpha(alphas[j],H,L)
                if (abs(alphas[j] - alphaJold) < 0.00001):
                    print "j not moving enough"
                    continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
                b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T -labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T -labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
                if (0 < alphas[i]) and (C > alphas[i]): b = b1
                elif (0 < alphas[j]) and (C > alphas[j]): b = b2
                else: b = (b1 + b2)/2.0
                alphaPairsChanged += 1
                print "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
        if (alphaPairsChanged == 0): iter += 1
        else: iter = 0
        print "iteration number: %d" % iter
    return b,alphas