標準二次規劃形式: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))
優化問題
\(\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))
本質上和軟間隔支持向量機相同,只是將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))