SVM代码实现

深碍√TFBOYSˉ_ 2023-06-03 05:56 7阅读 0赞

## Cvxopt解决二次规划 ##

标准二次规划形式:

\\(\\begin\{equation\}\\begin\{split\}\\min\\quad&\\frac\{1\}\{2\}\\mathtt\{x^TPx+q^Tx\}\\\\s.t\\quad&\\mathtt\{Gx\\le h\}\\\\&\\mathtt\{Ax\}=0\\end\{split\}\\end\{equation\}\\\\\\)

解决步骤

from cvxopt import matrix,solvers
    #组装出标准二次规划的各个矩阵参数(P,q,G,h,A,b)
    sol=solvers.qp(P,q,G,h,A,b)
    #sol['x']即为二次规划的x最优值

## 硬间隔支持向量机 ##

**优化问题**

\\(\\begin\{equation\}\\begin\{split\}\\min\\quad&\\frac\{1\}\{2\}\\sum\_\{n=1\}^N\\sum\_\{m=1\}^Ny\_ny\_m\\alpha\_n\\alpha\_m\\mathtt\{x\_n^Tx\_m\}-\\sum\_\{n=1\}^N\\alpha\_n\\\\s.t\\quad&\\sum\_\{n=1\}^N\\alpha\_ny\_n=0\\\\&\\alpha\_n\\ge0\\end\{split\}\\end\{equation\}\\\\\\)

**计算\\(\\alpha\\)**

优化问题转化为二次规划的标准形式带入即可得到\\(\\alpha\\)

**计算\\(\\mathtt\{w\}\\)**

由对偶问题中的\\(\\nabla\_\\mathtt\{w\}\\mathcal\{L\}=\\mathtt\{w\}-\\sum\_\{n=1\}^N\\alpha\_ny\_n\\mathtt\{x\_n\}=0\\)可以得出

\\(\\mathtt\{w\}=\\sum\_\{n=1\}^N\\alpha\_ny\_n\\mathtt\{x\_n\}\\)

**计算\\(b\\)**

从\\(\\alpha\\)中选一个非零的\\(\\alpha\_n\\),其对应的\\(\\mathtt\{x'\}\\)即为支持向量,即满足

\\(y\_n\\mathtt\{(w^Tx'+b)\}=1\\)

\\(b=\\frac\{1\}\{y\_n\}-\\mathtt\{w^Tx'\}\\)

**计算目标函数**

\\(g(\\mathtt\{x\})=sign(\\mathtt\{w^Tx+b\})\\)

代码实现

import numpy as np
    import matplotlib.pyplot as plt
    from cvxopt import matrix,solvers
    class Svm_h:
        def __init__(self,x,y):
            self.x=x
            self.y=y
            m=len(y)
            #下面的1.0,0.0都是为了转类型为浮点型
            P=matrix((x@x.T)*(y@y.T)*1.0)
            q=matrix(np.ones((m,1))*(-1))
            G=matrix(np.identity(m)*-1)
            h=matrix(np.zeros((m,1)))
            A=matrix(y.T*1.0)
            b=matrix(0.0)
            sol=solvers.qp(P,q,G,h,A,b)
            self.alpha=np.array(sol['x'],ndmin=2)
            w=self.alpha*y
            w=np.hstack((w,w))
            self.w=(w*x).sum(axis=0).T
            index=0
            while(self.alpha[index]<0.00001):
                #浮点数表示数字有损失,本文假设精度为0.00001,不太严谨请不必纠结
                #大于0的α对应的点为支持向量
                index+=1
            self.b=y[index]-self.w@x[index,:]
        def signal(self,x):
            #返回未通过符号函数的结果
            return self.w.T@x+self.b
        def g(self,x):
            #返回+1或-1
            return np.sign(self.signal(x))

## 软间隔支持向量机 ##

**优化问题**

\\(\\begin\{equation\}\\begin\{split\}\\min\\quad&\\frac\{1\}\{2\}\\sum\_\{n=1\}^N\\sum\_\{m=1\}^Ny\_ny\_m\\alpha\_n\\alpha\_m\\mathtt\{x\_n^Tx\_m\}-\\sum\_\{n=1\}^N\\alpha\_n\\\\s.t\\quad&\\sum\_\{n=1\}^N\\alpha\_ny\_n=0\\\\&\\ C\\ge \\alpha\_n\\ge0\\end\{split\}\\end\{equation\}\\\\\\)

**计算\\(\\alpha\\)**

优化问题转化为二次规划的标准形式带入即可得到\\(\\alpha\\)

**计算\\(\\mathtt\{w\}\\)**

由对偶问题中的\\(\\nabla\_\\mathtt\{w\}\\mathcal\{L\}=\\mathtt\{w\}-\\sum\_\{n=1\}^N\\alpha\_ny\_n\\mathtt\{x\_n\}=0\\)可以得出

\\(\\mathtt\{w\}=\\sum\_\{n=1\}^N\\alpha\_ny\_n\\mathtt\{x\_n\}\\)

**计算\\(b\\)**

从\\(\\alpha\\)中选一个大于0且小于\\(C\\)的\\(\\alpha\_n\\),其对应的\\(\\mathtt\{x'\}\\)即为支持向量,即满足

\\(y\_n\\mathtt\{(w^Tx'+b)\}=1\\)

\\(b=\\frac\{1\}\{y\_n\}-\\mathtt\{w^Tx'\}\\)

**计算目标函数**

\\(g(\\mathtt\{x\})=sign(\\mathtt\{w^Tx+b\})\\)

import numpy as np
    import matplotlib.pyplot as plt
    from cvxopt import matrix,solvers
    class Svm_s:
        def __init__(self,x,y,c):
            self.x=x
            self.y=y
            m=len(y)
            #下面的1.0,0.0都是为了转类型为浮点型
            P=matrix((x@x.T)*(y@y.T)*1.0)
            q=matrix(np.ones((m,1))*(-1))
            G=matrix(np.vstack((np.identity(m)*(-1.0),np.identity(m)*(1.0))))
            h=matrix(np.vstack((np.zeros((m,1))*1.0,np.ones((m,1))*c*1.0)))
            A=matrix(y.T*1.0)
            b=matrix(0.0)
            sol=solvers.qp(P,q,G,h,A,b)
            self.alpha=np.array(sol['x'],ndmin=2)
            w=self.alpha*y
            w=np.hstack((w,w))
            self.w=(w*x).sum(axis=0).T
            index=0
            while(0.00001>self.alpha[index] or self.alpha[index]>c-0.00001):
                #浮点数表示数字有损失,本文假设精度为0.00001,不太严谨请不必纠结
                #大于0且小于α的α对应的点为支持向量
                index+=1
            self.b=y[index]-self.w@x[index,:]
        def signal(self,x):
            #返回未通过符号函数的结果
            return self.w.T@x+self.b
        def g(self,x):
            #返回+1或-1
            return np.sign(self.signal(x))

## 软间隔RBF核函数支持向量机 ##

本质上和软间隔支持向量机相同,只是将x映射到高维上,\\(\\mathtt\{w\}\\)很难计算出

**计算\\(b\\)**

从\\(\\alpha\\)中选一个大于0且小于\\(C\\)的\\(\\alpha\_n\\),其对应的\\(\\mathtt\{x\_m\}\\)即为支持向量,即满足

\\(b=y\_m-\\sum\\limits\_\{\\alpha>0\}\\alpha\_ny\_nK(\\mathtt\{x\_n,x\_m\})\\)

**计算目标函数**

\\(g(\\mathtt\{x\})=sign(\\sum\\limits\_\{\\alpha\_n>0\}\\alpha\_ny\_nK(\\mathtt\{x\_n,x\})+b)\\)

import numpy as np
    import matplotlib.pyplot as plt
    from cvxopt import matrix,solvers
    class Svm_s:
        def __init__(self,x,y,c,r):
            self.x=x
            self.y=y
            self.r=r
            self.c=c
            m=len(y)
            #下面的1.0,0.0都是为了转类型为浮点型
            p=y@y.T*1.0
            for i in range(m):
                for j in range(m):
                    p[i][j]*=self.rbf_kernel(x[i,:],x[j,:])
            P=matrix(p)
            q=matrix(np.ones((m,1))*(-1))
            G=matrix(np.vstack((np.identity(m)*(-1.0),np.identity(m)*(1.0))))
            h=matrix(np.vstack((np.zeros((m,1))*1.0,np.ones((m,1))*c*1.0)))
            A=matrix(y.T*1.0)
            b=matrix(0.0)
            sol=solvers.qp(P,q,G,h,A,b)
            self.alpha=np.array(sol['x'],ndmin=2)
            self.index=0
            while(0.00001>self.alpha[self.index] or self.alpha[self.index]>c-0.00001):
                #浮点数表示数字有损失,本文假设精度为0.00001,不太严谨请不必纠结
                #大于0且小于α的α对应的点为支持向量
                self.index+=1
            self.b=y[self.index]-sum([y[i]*self.alpha[self.index]*self.rbf_kernel(self.x[i,:,None],self.x[self.index,:,None]) for i in range(m)])
        def rbf_kernel(self,x1,x2):
            return np.exp(-1*self.r*(x1-x2).T@(x1-x2))
        def signal(self,x):
            #返回未通过符号函数的结果
            sum=0
            for n in range(len(self.y)):
                sum+=self.y[n]*self.alpha[self.index]*self.rbf_kernel(self.x[n,:,None],x)+self.b
            return sum
        def g(self,x):
            #返回+1或-1
            return np.sign(self.signal(x))

转载于:https://www.cnblogs.com/redo19990701/p/11565361.html