SVM代码实现

Posted 2022-11-21 redo19990701

Cvxopt解决二次规划

标准二次规划形式:

\(\beginequation\beginsplit\min\quad&\frac12\mathttx^TPx+q^Tx\\s.t\quad&\mathttGx\le h\\&\mathttAx=0\endsplit\endequation\\\)

解决步骤

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

硬间隔支持向量机

优化问题

\(\beginequation\beginsplit\min\quad&\frac12\sum_n=1^N\sum_m=1^Ny_ny_m\alpha_n\alpha_m\mathttx_n^Tx_m-\sum_n=1^N\alpha_n\\s.t\quad&\sum_n=1^N\alpha_ny_n=0\\&\alpha_n\ge0\endsplit\endequation\\\)

计算\(\alpha\)

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

计算\(\mathttw\)

由对偶问题中的\(\nabla_\mathttw\mathcalL=\mathttw-\sum_n=1^N\alpha_ny_n\mathttx_n=0\)可以得出

\(\mathttw=\sum_n=1^N\alpha_ny_n\mathttx_n\)

计算\(b\)

从\(\alpha\)中选一个非零的\(\alpha_n\),其对应的\(\mathttx'\)即为支持向量,即满足

\(y_n\mathtt(w^Tx'+b)=1\)

\(b=\frac1y_n-\mathttw^Tx'\)

计算目标函数

\(g(\mathttx)=sign(\mathttw^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))

软间隔支持向量机

优化问题

\(\beginequation\beginsplit\min\quad&\frac12\sum_n=1^N\sum_m=1^Ny_ny_m\alpha_n\alpha_m\mathttx_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\endsplit\endequation\\\)

计算\(\alpha\)

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

计算\(\mathttw\)

由对偶问题中的\(\nabla_\mathttw\mathcalL=\mathttw-\sum_n=1^N\alpha_ny_n\mathttx_n=0\)可以得出

\(\mathttw=\sum_n=1^N\alpha_ny_n\mathttx_n\)

计算\(b\)

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

\(y_n\mathtt(w^Tx'+b)=1\)

\(b=\frac1y_n-\mathttw^Tx'\)

计算目标函数

\(g(\mathttx)=sign(\mathttw^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映射到高维上,\(\mathttw\)很难计算出

计算\(b\)

从\(\alpha\)中选一个大于0且小于\(C\)的\(\alpha_n\),其对应的\(\mathttx_m\)即为支持向量,即满足

\(b=y_m-\sum\limits_\alpha>0\alpha_ny_nK(\mathttx_n,x_m)\)

计算目标函数

\(g(\mathttx)=sign(\sum\limits_\alpha_n>0\alpha_ny_nK(\mathttx_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))