SVM代码实现

时间 2019-11-13

标签 svm 代码实现繁體版

原文原文链接

Cvxopt解决二次规划

标准二次规划形式:python

\(\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}\\\)ide

解决步骤函数

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}\\\)spa

计算\(\alpha\)code

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

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

由对偶问题中的\(\nabla_\mathtt{w}\mathcal{L}=\mathtt{w}-\sum_{n=1}^N\alpha_ny_n\mathtt{x_n}=0\)能够得出class

\(\mathtt{w}=\sum_{n=1}^N\alpha_ny_n\mathtt{x_n}\)import

计算\(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))

软间隔支持向量机

优化问题

计算\(\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))