SVM代碼實現

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))