SVM求解_SMO 机器学习

最新推荐文章于 2024-04-20 23:16:41 发布

明朝百晓生

最新推荐文章于 2024-04-20 23:16:41 发布

阅读量624

点赞数

分类专栏：人工智能

4AM_明朝百晓生

本文链接：https://blog.csdn.net/chengxf2/article/details/103118280

版权

人工智能专栏收录该内容

102 篇文章 3 订阅

订阅专栏

前言：

SMO 是求解SVM的一种算法(顺序最小化算法)sequential minimal optimization ，由1998年Platt提出。

主要解如下凸二次规划对偶问题：

$min_{\alpha} \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_i\alpha_jy_iy_jK(x_i,x_j)-\sum_{i=1}^{N}\alpha_i$

$s.t : \left\{\begin{matrix} \sum_{i=1}^{N} \alpha_i y_i=0\\ 0 \leq \alpha_i \leq C, i =1,2,..N \end{matrix}\right.$

求解过程中里面关于二次规划算法收敛性证明，看了很多文档一直没找到。

SMO简介
二次规划解析方法
变量选择方法
SMO 算法
例子

一 SMO简介

SMO是一种启发式算法：

如果所有的变量满足KKT条件，最优化问题的解就得到了。

否则选择两个变量，固定其它变量，针对这两个变量构建一个二次规划问题。

整个SMO算法包括两个部分：

a 求解二次规划的解析方法

b 选择变量的启发方法

二二次规划解析

SMO 求解的子问题为：

$min_{\alpha_1,\alpha_2} W(\alpha_1, \alpha_2)= \frac{1}{2}K_{11}\alpha_1^2+\frac{1}{2}K_{22}\alpha_2^2+y_1y_2K_{12}\alpha_1\alpha_2-(\alpha_1+\alpha_2)+y_1\alpha_1\sum_{i=3}^{N}y_i\alpha_iK_{i1}+y_2\alpha_2\sum_{i=3}^{N}y_i\alpha_iK_{i2}$

$s.t : \left\{\begin{matrix} \alpha_1y_1+\alpha_2y_2=-\sum_{i=3}^{N}y_i\alpha_i=\zeta \\ 0 \leq \alpha_i \leq C , i=1,2 \end{matrix}\right.$

首先分析约束条件，在满足约束条件下，求解极小值

假设初始可行解为 $\alpha_1^{old},\alpha_2^{old}$ ,最优解为 $\alpha_1^{new},\alpha_2^{new}$ ,其中

$L \leq \alpha_2^{new} \leq H$

因为：

$\alpha_1^{old}y_1+\alpha_2^{old}y_2=\alpha_1^{new}y_1+\alpha_2^{new}y_2$

$\alpha_2^{new}=\alpha_1^{old}y_1y_2+\alpha_2^{old}-\alpha_1^{new}y_1y_2$

2.1 当 $y_1 \neq y_2 => y_1y_2=-1$ ，根据 $\alpha_1^{new} \in [0,C]$

则： $\alpha_2^{new}=\alpha_2^{old}-\alpha_1^{old}+\alpha_1^{new}$

$L= max(0, \alpha_2^{old}-\alpha_1^{old})$

$H= min(C,C+ \alpha_2^{old}-\alpha_1^{old})$

2.2 当 $y_1= y_2 =>y_1y_2=1$

则 $\alpha_2^{new}=\alpha_1^{old}+\alpha_2^{old}-\alpha_1^{new}$

$L=max(0, \alpha_1^{old}+\alpha_2^{old}-C)$

$H=min(C,\alpha_1^{old}+\alpha_2^{old})$

为了计算方便，记

$g(x)=\sum_{i=1}^{N}\alpha_iy_iK(x_i,x)+b$

$E_i=g(x_i)-y_i=(\sum_{j=1}^{N}\alpha_jy_jK(x_j,x_i)+b)-y_i,i=1,2$

定理1：

沿着约束方向，未经剪辑的解是

$\alpha_2^{new,unc}=\alpha_2^{old}+\frac{y_2(E_1-E_2)}{\eta}$

其中：

$\eta = K_{11}+K_{22}-2K_{12}$

经过剪辑后 $\alpha_2$ 的解是：

$\alpha_2^{new}=\left\{\begin{matrix} H, \alpha_2^{new,unc}>H\\ \alpha_2^{new,unc} , L \leq \alpha_2^{new,unc} \leq H \\ L, \alpha_2^{new,unc}<L \end{matrix}\right.$

接着求 $\alpha_1^{new}$

$\alpha_1^{new}=\alpha_1^{old}+y_1y_2(\alpha_2^{old}-\alpha_2^{new})$

证明：

引入记号

$v_i=\sum_{j=3}^{N}\alpha_jy_jK(x_i,x_j)=g(x_i)-\sum_{j=1}^{2}\alpha_jy_jK(x_i,x_j)-b$

目标函数可写成：

$W(\alpha_1,\alpha_2)=\frac{1}{2}K_{11}\alpha_1^2+\frac{1}{2}K_{22}\alpha_2^2+y_1y_2K_{12}\alpha_1\alpha_2-(\alpha_1+\alpha_2)+y_1v_1\alpha_1+y_2v_2\alpha_2$

因为:

$\alpha_1y_1=\zeta -\alpha_2y_2$ $=>\alpha_1=(\zeta-y_2\alpha_2)y_1$

得到关于 $\alpha_2$ 的目标函数：

$W(\alpha_2)=\frac{1}{2}K_{11}(\zeta-\alpha_2y_2)^2+\frac{1}{2}K_{22}\alpha_2^2+y_2K_{12}(\zeta-\alpha_2y_2)\alpha_2-(\zeta-\alpha_2y_2)y_1-\alpha_2+v_1(\zeta-\alpha_2y_2)+y_2v_2\alpha_2$

对 $\alpha_2$ 求导：

$K_{11}\alpha_2+K_{22}\alpha_2-2K_{12}\alpha_2-K_{11} \zeta y_2+K_{12}\zeta y_2+y_1y_2-1-v_1y_2+y_2v_2=0$

$\alpha_2=\frac{y_2(y_2-y_1+\zeta K_{11}-\zeta K_{12}+v_1-v_2)}{K_{11}+K_{22}-2K_{12}}$

$=\frac{y_2(y_2-y_1+\zeta K_{11}-\zeta K_{12}+v_1-v_2)}{K_{11}+K_{22}-2K_{12}}$

$=\frac{y_2(y_2-y_1+(\alpha_1^{old}y_1+\alpha_2^{old} y_2)+(g(x_1)-\sum_{j=1}^{2}y_j\alpha_jK_{1j}-b) -(g(x_2)-\sum_{j=1}^{2}y_j\alpha_jK_{2j}-b))}{k_{11}+k_{22}-2k_{12}}$

带入 $\eta=K_{11}+K_{22}-2K_{12}$ 得到

$\alpha_2^{new,unc}=\alpha_2^{old}+\frac{y_2(E_1-E_2)}{\eta}$

三变量的选择

3.1 第一个变量选择（外层循环）

选取违反KKT最严重的样本点，并将其对应的变量作为第一个变量

$\alpha_i=0 =>y_ig(x_i) \geq 1$

$0<\alpha_i<C =>y_ig(x_i)=1$

$\alpha_i=C => y_ig(x_i) \leq 1$

其中：

$g(x_i)=\sum_{j=1}^{N}\alpha_jy_jK(x_i,x_j)+b$

3.2 第二个变量选择（内层循环）

$\alpha_2$ 选择 $|E_1-E_2|$ 最大的

3.3 计算阀值b 和差值 $E_i$

$\sum_{i=1}^{N} \alpha_iy_iK_{i1}+b=y_1$

$b_1^{new}=y_1-\sum_{i=3}^{N}\alpha_iy_iK_{i1}-\alpha_1^{new}y_1K_{11}-\alpha_2^{new}y_2K_{21}$ ...................1

因为：

$E_1=\sum_{i=3}^{N}\alpha_iy_iK_{i1}+\alpha_1^{old}y_1K_{11}+\alpha_2^{old}y_2K_{21}+b^{old}-y_1$

$y_1-\sum_{i=3}^{N}\alpha_iy_i K_{i1}=-E_1+\alpha_1^{old}y_1K_{11}+\alpha_2^{old}y_2K_{21}+b^{old}$

带入1：

$b_1^{new}=-E_1-y_1K_{11}(\alpha_1^{new}-\alpha_1^{old})-y_2K_{21}(\alpha_2^{new}-\alpha_2^{old})+b^{old}$

同理

$b_2^{new}=-E_2-y_1K_{12}(\alpha_1^{new}-\alpha_2^{new})-y_2K_{22}(\alpha_2^{new}-\alpha_2^{old})+b^{old}$

如果 $\alpha_1==\alpha_2 \in (0,C)$

$b_1^{new}=b_2^{new}$

如果是0 或者C, 选择它们的中点作为 $b^{new}$

四 SMO算法

输入：

训练数据集

$T={(x_1,y_1),(x_2,y_2),....(x_N,y_N)}$

$y_i \in{-1,+1}$

精度 $\varepsilon$

C: 参数上限

输出： $\alpha$

；；；；；；；；；；；；；流程；；；；；；；；；；；；；；；；；；；；

取初始值 $\alpha^{0}=0,k=0$ , k 代表迭代次数
选取优化变量 $\alpha_1^{k},\alpha_2^{k}$ , 解析求解两个变量的最优化问题, 得到新的 $\alpha_1^{k+1},\alpha_2^{k+1}$

$\alpha_2^{new,unc}=\alpha_2^{old}+\frac{y_2(E_1-E_2)}{\eta}$

$\alpha_1^{new}=\alpha_1^{old}+y_1y_2(\alpha_2^{old}-\alpha_2^{new})$

$b_1^{new}=y_1-\sum_{i=3}^{N}\alpha_iy_iK_{i1}-\alpha_1^{new}y_1K_{11}-\alpha_2^{new}y_2K_{21}$

$b_2^{new}=-E_2-y_1K_{12}(\alpha_1^{new}-\alpha_2^{new})-y_2K_{22}(\alpha_2^{new}-\alpha_2^{old})+b^{old}$

3 若在精度 $\varepsilon$ 满足停止条件

$\sum \alpha_iy_i=0$

$0 \leq \alpha_i \leq C,i =1,2,...N$

$y_i \cdot g(x_i) =\left\{\begin{matrix} \geq 1 ,\alpha_i=0\\ =1, \alpha_i \in (0,C) \\ \leq 1, \alpha_i= C \end{matrix}\right.$

$g(x_i)=\sum_{j=1}^{N}\alpha_iy_iK(x_i,x_j)+b$

转4，否则令k=k+1,转2

4 $\alpha=\alpha^{k+1}$

五例子：

高斯核

线性核

CODE

# -*- coding: utf-8 -*-
"""
Created on Wed Nov 20 17:03:44 2019

@author: chengxf2
"""
import numpy as np
import os
import matplotlib.pyplot  as plt


class SVM:
    
    
    def DrawPot(self):
        
        fig = plt.figure()
        ax = fig.add_subplot(111)
        
        for i in range(self.m):
            
            x1 = self.dataMat[i,0]
            x2 = self.dataMat[i,1]
            
            label = self.labelMat[i,0]
            
            
            if i in self.svIndex:
                if label == -1:
                    ax.scatter(x1, x2, c = 'red' ,marker= 'o',s=80)
                else:
                    ax.scatter(x1, x2, c = 'green' , marker= 's',s=80)

            else:
                
                    ax.scatter(x1, x2, c ='black', marker= 'x',s=10)
        plt.xlabel("x1")
        plt.ylabel("x2")
        plt.show()

                
            
            
    
    """
    核函数
    Args
         dataMat: 数据集
         z: 样本
    return
         k: 核函数
    """
    def Kernel(self,  dataMat, z):
        
        m,n = np.shape(dataMat)
        
        k = np.mat(np.zeros((m, 1)))
        
        key = self.opt[0]
        para = self.opt[1]
        
        #print("\n key : ",key)
  
        if "linear"==key:
            k = dataMat*z.T
        
        elif "Gauss"==key:
            
            for j in range(m):
                x = dataMat[j,:]
               # print("\n x: ",x)
                delta = x-z
                k[j]=delta*delta.T
            sigma = -para**2
            
            k = np.exp(k/sigma)
        else:
            raise NameError("Error")
        
        return k
                
        
        
    
    """
    加载数据集
    Args
       None
    return
       None
    """
    def LoadData(self):
        
        
        file = open(self.path)
        lines = file.readlines()
        data  = []
        label =[]
        for line in lines:
            
            arr = line.strip().split('\t')
            data.append([float(arr[0]),float(arr[1])])
            label.append(float(arr[2]))
        
        self.labelMat =np.mat(label).T  # 标签集 #100 *1
        self.dataMat = np.mat(data)  ##数据集
        self.m, self.n = np.shape(self.dataMat)
        self.alpha = np.mat(np.zeros((self.m,1)))  #拉格朗日乘子
        
        self.eCache = np.mat(np.zeros((self.m,2))) ##[0,Ei]， 第一位是标志位（0,1),第二位是Ei: g(xi)-yi
        self.KMat = np.mat(np.zeros((self.m, self.m)))  ##Gram矩阵,核函数
        
        for i in range(self.m):
            self.KMat[:,i] = self.Kernel(self.dataMat, self.dataMat[i,:])
            
        #print("\n   核函数：  \n",self.k)  
    
    """
    计算当前预测值和实际值的差别
    Args
      i: 样本点
    return 
      Ei: 预测值和实际值的差
    """
    def  CalcE(self, i):
          
         para = np.multiply(self.alpha, self.labelMat).T*self.KMat[:,i] ##(m,1)
         
         #print("para ",para)
         gx = float(para+self.b)
         
         ei = gx - self.labelMat[i,0]
         #print("\n ",ei)
         return ei
     
    """
    产生一个不同于i的值
    Args
      i: 当前的i
    return 
      j
    """
    def RandJ(self, i):
        j = i
        
        while(j==i):
            
             j = np.random.uniform(0,self.m)
        return int(j)
        
     
    """
    选择J
    Args
       i: 样本点
       Ei: 差值
    return
       j:  对应的alphaJ
       Ej:  对应的差值
    """
    def ChooseJ(self, i, Ei):
        
        maxJ = -1; maxDiff = 0.0; Ej = 0.0;
        
        self.eCache[i]=[1,Ei]
        
        eList = np.nonzero(self.eCache[:,0].A)[0] 
        n = len(eList)
        
        if n>1: ###循环1
            for e in eList:
                if e == i:
                    continue
                Ek = self.CalcE(e)
                diffE = abs(Ei-Ek)
                
                if diffE>maxDiff:
                    maxDiff = diffE; maxJ = e; Ej = Ek
            return maxJ, Ej  ######***********************************逻辑差一点
        else:
            j = self.RandJ(i)
            Ej = self.CalcE(j)
        return j, Ej
    
    """
    编辑参数
    Args
        aj: 参数
        H: 上限
        L : 下限
    return
        aj: 剪辑过的参数
    """      
    def ClipAlpha(self, aj, H, L):
        
        if aj>H:
            aj = H
        elif aj<L:
            aj = L
        
        return aj
    
    """
    当前的点是否违反KKT
    Args
        i: 当前的样本点
        Ei： g-lable 预测值减去实际值
        alpha+u = C
        u*tol = 0
        
        1: 当alpha (0,C) => u=(0,c),tol = 0, 则xi落在间隔边界上,yg-1 = 0=> yEi = 0
        2: 当alpha = C=>u=0, yg-1+tol=0 =>yEi+tol = 0，只能落在分离平面和间隔边界间，
        3: 当 alpha =0=> u=C, tol =0 :  yg>=1 => yEi>=0
    """
    def VolidataKKT(self, i, Ei):
        
        ZERO = 1e-3
        """
        bKKT = ((self.labelMat[i,0]*Ei<-self.tol)  and (self.alpha[i,0]<self.c)) or \
        ((self.labelMat[i,0]*Ei>self.tol) and (self.alpha[i,0]>0))
        """
        label = self.labelMat[i,0]
        
        if self.alpha[i,0]>0 and self.alpha[i,0]<self.c :## 落在支撑平面上  #yg=1
             if abs(label*Ei)>ZERO:
                 return True
             else:
                 return False
                
        if  abs(self.alpha[i,0]-self.c)<=ZERO:  ## alpha =c,ui=0,yg-1+tol=0 落在支撑平面内 yg <=1
            
            if abs(label*Ei+self.tol)>ZERO:
                return True
            else:
                return False
            
        if abs(self.alpha[i,0])<ZERO:  ## ui=C, tol =0 ,当前的点不起作用， yg>=1 就可以了
            
            if label*Ei<ZERO:
                return True
            else:
                return False
        

        
        return False
                 
        
            
        
        
      
    
    """
    保存差值矩阵
    Args:
        i
    return
        None
    """
    def SaveE(self, i):
        Ei = self.CalcE(i)
        
        self.eCache[i]=[1,Ei]
    """
    进入内循环
    Args
       i: 当前样本点
       先检测alphaI是否违反KKT，如果违反，再找到alphaJ,更新对应的参数
    return
        1: 进行了优化
        0：未优化
    """
    def InnerLoop(self, i):
        
        Ei = self.CalcE(i)
        bKKT = self.VolidataKKT(i, Ei)
        if bKKT:
            
            j, Ej = self.ChooseJ(i, Ei) ##选择合适的alphaJ,用来更新
            
            IOld = self.alpha[i,0].copy()
            JOld = self.alpha[j,0].copy()
            
            ###确定上下界
            if (self.labelMat[i,0]!= self.labelMat[j,0]):
                L = max (0, self.alpha[j,0]-self.alpha[i,0])
                H = min(self.c, self.c+self.alpha[j,0]-self.alpha[i,0] )
            else:
                L = max(0, self.alpha[j,0]+self.alpha[i,0]-self.c)
                H = min(self.c, self.alpha[j,0]+self.alpha[i,0])
            
            
            if L == H:
                return 0
            
            eta = 2.0*self.KMat[i,j]-self.KMat[i,i]-self.KMat[j,j]
            
            if eta>=0:
                return 0
            
            self.alpha[j]-= self.labelMat[j,0]*(Ei-Ej)/eta
            self.alpha[j]= self.ClipAlpha(self.alpha[j,0],H,L)
            self.SaveE(j)
          
            
            if abs(self.alpha[j,0]-JOld)<1e-5:
                return 0
            
            self.alpha[i] +=self.labelMat[i,0]*self.labelMat[j,0]*(JOld- self.alpha[j,0])
            self.SaveE(i)
            
            b1 = -Ei- self.labelMat[i]*self.KMat[i,i]*(self.alpha[i]-IOld)-self.labelMat[j]*self.KMat[i,j]*(self.alpha[j]-JOld)+self.b
            
            b2 = -Ej- self.labelMat[i]*self.KMat[i,j]*(self.alpha[i]-IOld)-self.labelMat[j]*self.KMat[j,j]*(self.alpha[j]-JOld)+self.b
            
            if (self.alpha[i]>0 and self.alpha[i]<self.c):
                self.b = b1
            elif(self.alpha[j]>0 and self.alpha[j]<self.c):
                self.b = b2
            else:
                self.b =(b1+b2)/2.0
            
            return 1
        else:
            return 0
            

        
            
        
        
        
    """
    分类
    Args
       testMat： 测试的样本
       testLabel: 测试的标签
    """
    def classification(self, testMat, testLabel):
         
        m,n = np.shape(testMat)
        print("\n *****classification******** ")
        para = np.multiply(self.sptLabel, self.alpha[self.svIndex])  ##(1,35)
        
        print("\n para ",np.shape(para))
        
        errorNum = 0
        for i in range(m):
            z = testMat[i,:]
            kMat =self.Kernel(self.sptVects, z) ##转换为核函数
            
            pred = kMat.T*para+self.b
            
            #print("\n pred ",pred)
            
            if np.sign(pred) != np.sign(testLabel[i]):
                errorNum = errorNum+1
                
        print("错误率:::::::::::: ", errorNum/m)
            
            
            
            
        
        
    def __init__(self):
        
      
        
        self.b = 0
        self.m = 0
        self.path = os.path.abspath("testSet.txt")
        self.c =1.0 ##惩罚因子 alpha取值上限
        self.tol = 0.1 ##相当于软间隔
        self.maxIter = 10000## 迭代次数
        self.opt = ("linear", 1.0)  ##核函数 linear  Gauss
        self.LoadData()
        
    """
    训练
    Args:
        none
    return
        b:
        sptVects: 支持向量的点
    """
    def Train(self):
        print("\n **********Train************** \n")
        step = 0 ##外层迭代轮数
        entireFlag = True ##完整的遍历
        alphaPairsChanged = 0 ##是否由参数发生变更
        
       
        
        while(step <self.maxIter) and ((alphaPairsChanged>0) or entireFlag):
              
              alphaPairsChanged = 0
              if entireFlag :
                  for i in range(self.m):
                      alphaPairsChanged +=self.InnerLoop(i)  
                  step +=1
              else:
                  
                  nonBoundIs = np.nonzero((self.alpha.A>0)*(self.alpha.A<self.c))[0]
                  
                  for i in nonBoundIs:
                      alphaPairsChanged +=self.InnerLoop(i)
                  step +=1
              if entireFlag:
                   entireFlag = False ##切换标志位
              elif (alphaPairsChanged ==0):
                   entireFlag = True
        print("\n **********迭代次数************** \n",step)
        #print("\n alpha:  ",self.alpha)
        print("\n: b: ",self.b)
        
        self.svIndex = np.nonzero((self.alpha.A>0)*(self.alpha.A<self.c))[0] ##支持向量的点
        self.sptVects = self.dataMat[self.svIndex]   ##支持向量的点
        self.sptLabel = self.labelMat[self.svIndex]  ##支持向量点的标签
        #print("\n svIndex ",self.svIndex)
        #print("\n svIndex ",self.sptLabel)
        self.classification(self.dataMat, self.labelMat)
        self.DrawPot()
        
        
        
    
        

SVM = SVM()        
SVM.Train()

参考文档

《机器学习与应用》

《统计学习方法》

明朝百晓生

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
SVM求解_SMO 机器学习

前言： SMO 是求解SVM的一种算法(顺序最小化算法)sequential minimal optimization ，由1998年Platt提出。主要解如下凸二次规划对偶问题：求解过程中里面关于二次规划算法收敛性证明，看了很多文档一直没找到。目录： SMO简介二次规划解析方法变量...
复制链接

扫一扫