机器学习（五）SVM

优化算法

像拟合，分类算法基本上都在用已知的样本数据，进行优化来得到权重和偏置值来达到想要的目的。

最优化问题一般是指对于某一个函数而言，求解在其指定作用域上的全局最小值 问题，一般分为以下三种情况(备注:以下几种方式求出来的解都有可能是局部 极小值，只有当函数是凸函数的时候，才可以得到全局最小值):

• 无约束问题:求解方式一般求解方式梯度下降法、牛顿法、坐标轴下降法等. $minf(x)$
• 等式约束条件:求解方式一般为拉格朗日乘子法.$minf(x)s.t.h_k(x)=0,k=1,2...,p$
• 不等式约束条件:求解方式一般为KKT条件$minf(x)s.t:(1)h_k(x)=0,k=1,2...p(2)g_i(x)⩽0，j=1,2,..p$

---

在SVM中我们会用到的是拉格朗日乘子法，这里简单介绍一下并不加证明：
格朗日乘子法就是当我们的优化函数存在等值约束的情况下的一种最优化求解 方式;其中参数α被称为拉格朗日乘子，要求α不等于0
$$\begin{gathered} minf(x) \\ s.t::h_i(x)=0,i=1,2,...p \\ \downarrow \\ minf(x)+\sum _{i=1}^p{\alpha }_i{h}_i(x) \end{gathered}$$
求解的方式对每个维度进行求导为0

---

在这里进行计算时，我们在求正向求解的过程中可能会出现问题，有一个非常🐮的方法是将问题变成对偶问题，这里用到的条件是KKT条件。正向会出现什么问题稍后会举出例子🌰。
这里先解释什么是对偶问题和KKT条件：

对偶问题特征：

• 对偶问题的对偶是原问题;
• 无论原始问题是否是凸的，对偶问题都是凸优化问题;
• 对偶问题可以给出原始问题的一个下界;
• 当满足一定条件的时候，原始问题和对偶问题的解是完美等价的。

KKT条件
KKT条件是泛拉格朗日乘子法的一种形式;主要应用在当我们的优化函数存在不 等值约束的情况下的一种最优化求解方式;KKT条件即满足不等式约束情况下的条件。
$$\begin{gathered} \underset{x}{\min }f(x)s.t:h_k(x)=0,k=1,2...p \\ {g}_j(x)⩽0,j=1,2...p \\ \downarrow \\ L(x,\alpha ,\beta )=f(x)+\sum _{i=1}^p{\alpha }_i{h}_i(x)+\sum _{i=1}^q{\beta }_i{g}_i(x);{\alpha }_i\ne 0,{\beta }_i\ge 0 \\ \underset{x}{\min }L(x,\alpha ,\beta ) \end{gathered}$$
其中KKT条件中的具体过程，瑞典皇家理工学院的一个ppt里面有一个完整的过程。

算法原理

在正式算法开始之前需要先考虑几个问题（这里仅讨论最简单的二分类的情况）：

• Q1为什么叫支持 向量 ，支持向量是什么❓
• Q2他跟之前的分类算法有什么区别❓
• Q3他跟优化有什么关系，处理什么样的数据，怎么样进行计算❓
• Q4为什么说SVM是将数据映射到高维，说他是一个非线形的算法❓
• Q5如何将数据变为非线形❓
• Q6SVM是如何利用拉格朗日乘子法的❓
• Q7提出SVM的核心是什么❓

---

接下来整体介绍算法：

首先，先假设数据是完全可分的。我们可以看到有非常多的线都可以将数据点完美的分开但是选择哪一个具有最强的泛化性。这里就是一个点，提高数据可分的泛化性。
我们可以到整一条粗的黄线和中间的黑线，黄线的边缘和中间的黑线是平行的。我们希望边界到中间的距离越大越好（margin），决定边缘位置的是离分界线的最近的几个点决定的。那如果我们先将最近的几个点选择出来，通过这几个点作出边界最大的分界线。这就是SVM第一个贡献，提出了最大边界的概念。这里我们想要计算出这个线的话需要线计算出距离。数据为y，x其中x是一个向量表示有多个特征。我们将整样本的y定义为1，负样本的值定义为-1

公式推导

首先我们需要求出点到直线的距离,并假设所有的样本都能够正确的分类。
那么问题就变成了距离分界线最近的点离分界线最远，将其写成数学表达式：
a r g m a x w , b { m i n i ( y i ∗ ( w T x i + b ) ) ∗ 1 ∣ ∣ w ∣ ∣ } arg\mathop{max}\limits_{w,b}\{\mathop{min}\limits_{i}(y_i*(w^Tx_i+b))*\frac{1}{||w||} \} argw,bmax​{imin​(yi​∗(wTxi​+b))∗∣∣w∣∣1​}
先求的最小距离的点，后将其平面距离最大化。

---

距离

先看大括号内部。
$$\begin{gathered} d=\frac{|W^{T}x+b|}{||w||} \\ (w^T\ast x+b)y\ge \gamma ,\gamma >0 \\ \frac{(w^T\ast x+b)y}{\gamma }\ge 1 \end{gathered}$$
其中$\gamma$在是在点到平面距离最小时成立，$w,b$同时进行缩放并不会影响分界线，故$\frac{(w^T\ast x+b)y}{\gamma }\ge 1$可以化简写成$(w^T\ast x+b)y\ge 1$其中在距离小时取到1.
min里面的部分就完成了

---

拉格朗日

原来的问题就变成
$$\begin{gathered} arg\underset{w,b}{\max }\frac{1}{||w|{|}_2} \\ s.t:y^i\ast (w^T{x}^{(i)}+b\ge 1),i=1,2...n \end{gathered}$$
然后将最大化转换为最小化,这个1/2是为了后续计算方便加的，不影响最终结果。
$$\begin{gathered} arg\underset{w,b}{\min }\frac{1}{2}||w|{|}_2 \\ s.t:y^i\ast (w^T{x}^{(i)}+b\ge 1),i=1,2...n \end{gathered}$$

这就转换成我们一开始提到的拉格朗日问题：
转换成优化公式，其中$\lambda$>0

KKT

$$L(w,b,\lambda )=\frac{1}{2}||w|{|}_2+\sum _{i=1}^n\underset{>0}{{\lambda }_i}\underset{<0}{(1-y_i(w^{T}x+b)})$$
上式可以转换为，这里在解释一下为什么拉格朗日函数的最大值的最小化
$$\begin{gathered} 1-y_i(w^T{x}_i+b)\ge 0成立\underset{\lambda }{\max }L(w,b,\lambda )=\frac{1}{2}||w|{|}_2+\infty =\infty \\ 这样子看是不是没有意义 \\ 1-y_i(w^T{x}_i+b)\le 0成立\underset{\lambda }{\max }L(w,b,\lambda )=\frac{1}{2}||w|{|}_2+0=\frac{1}{2}||w|{|}_2 \\ 这就与原要求解的内容相同只是添加了约束条件 \end{gathered}$$
$$\begin{gathered} \underset{w,b}{\min }\underset{\lambda }{\max }L(w,b,\lambda ) \\ s.t:\lambda \ge 0 \end{gathered}$$
再通过KKT条件转换为：
$$\begin{gathered} \underset{\lambda }{\max }\underset{w,b}{\min }L(w,b,\lambda ) \\ s.t:\lambda \ge 0 \end{gathered}$$

先进行对内部的min进行求解，对b和w进行求偏导
$$\begin{gathered} \frac{\partial L}{\partial b}=0\to \sum _{i=1}^n{\lambda }_i{y}_i=0 \\ \frac{\partial L}{\partial w}=0\to w=\sum _{i=1}^n{\lambda }_i{y}_i{x}_i \\ thanL(w,b,\lambda )=\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{x}_j \\ than\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{x}_j \\ \leftrightarrow \underset{\lambda }{\min }-\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{x}_j \end{gathered}$$
最后就是对$\lambda$进行求偏微分即可求出${\lambda }_i$的值
$$\begin{gathered} \exists (x_k,y_k),(1-y_k(w^T{x}_k+b))=0 \\ {w}^{\ast }=\sum _{i=1}^n{\lambda }_i{y}_i{x}_i \\ {b}^{\ast }=y_k-\sum _{i=1}^n{\lambda }_i{y}_i{x}_i^T{x}_k \\ f(x)=w^{{\ast }^T}x+b^{\ast } \end{gathered}$$
$w^{\ast },b^{\ast }$就是我们想要的结果。

---

softmargin

但是上面的内容有一个比较大的问题就是，实现的过程中默认所有的点都能够分对的情况，但是往往事与愿违，我们需要考虑存在点无法被分对，使整体保证最优。

这里就需要允许本身存在一些失误。
我们在公式中加入一个常数$C>0$
$$\begin{gathered} min\frac{1}{2}||w|{|}^2+loss \\ \underset{w,b}{\min }\frac{1}{2}||w|{|}^2+C\sum _{i=1}^n{\zeta }_i(y_i(w^T{x}_i+b)-1) \\ \zeta (z)=\{\begin{array}{lr}1& ifz<0\\ 0& otherwise\end{array} \end{gathered}$$
$\zeta$表示是否满足约束条件（该点是否正确分类）
显然，当C为无穷时，迫使样本均满足上式约束，与我们没有softmargin 的结果一样，当C为有限值时，允许一些样本不满足约束。
在计算的过程中$\zeta$是不满足可微的，因此适当的修改loss
$$\begin{gathered} \min \frac{1}{2}||w|{|}^2+C\sum _{i=1}^n{\zeta }_i \\ s.t:y_i(w^{T}x+b)\ge 1-{\zeta }_i \\ {\zeta }_i\ge 0 \end{gathered}$$
分类正确的时候loss会趋近于0，错误会有一个正值。

---

SMO

优化方程
$$\begin{gathered} L(w,b,\zeta ,\mu ,\lambda )=\frac{1}{2}||w|{|}^2+C\sum _{i=1}^m{\zeta }_i+\sum _{i=1}^m{\lambda }_i[1-{\zeta }_i-y_i(w^T{x}_i+b)]-\sum _{i=1}^m{\mu }_i{\zeta }_i \\ \underset{w,b,\zeta }{\min }\underset{\lambda ,\mu }{\max }L(w,b,\zeta ,\mu ,\lambda )\leftrightarrow \underset{\lambda ,\mu }{\max }\underset{w,b,\zeta }{\min }L(w,b,\zeta ,\mu ,\lambda ) \end{gathered}$$
然后再对每一项进行求偏导
对于SMO如何选择变量具体推到过程在链接里面，这里不再推导，过程比较复杂。核心思想就是将需要多个有优化计算的参数进行分解计算。

算法实现

from numpy import *
def loadDataSet(fileName):
    dataMat = []; labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat

def selectJrand(i,m):
    j=i #we want to select any J not equal to i
    while (j==i):
        j = int(random.uniform(0,m))
    return j

def clipAlpha(aj,H,L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

data, label = loadDataSet('SVMtestSet.txt')

计算和迭代的过程都是和上面的SMO的优化过程是一致的。

from numpy import *
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 checks if an example violates KKT conditions
            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])#update i by the same amount as j
                                                                        #the update is in the oppostie direction
                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
b,alpha = smoSimple(data,label,0.6,0.001,40)

其中若alpha不为0，则表示该向量为支持向量，在下面结果展示图中标记为蓝色

w = (label.reshape(-1)*np.array(alpha).reshape(-1)).T@data
x = np.arange(12)
y = np.array((-b-w[0]*x)/w[1]).reshape(-1)

ata=np.array(data)
label=np.array(label)
plt.scatter(data[:,0][label==1],data[:,1][label==1])
plt.scatter(data[:,0][label==-1],data[:,1][label==-1])
plt.plot(x,y,'g--')
y1 = np.array((-b-w[0]*x-1)/w[1]).reshape(-1)
plt.plot(x,y1,'m-')
y2 = np.array((-b-w[0]*x+1)/w[1]).reshape(-1)
plt.plot(x,y2,'r-')
plt.plot(data[:,0][alpha>0],data[:,1][alpha>0],'bs')

结果如下：