怎樣使用Python從零實(shí)現(xiàn)多分類(lèi)SVM呢？

本文將首先簡(jiǎn)要概述支持向量機(jī)及其訓(xùn)練和推理方程，然后將其轉(zhuǎn)換為代碼以開(kāi)發(fā)支持向量機(jī)模型。之后然后將其擴(kuò)展成多分類(lèi)的場(chǎng)景，并通過(guò)使用Sci-kit Learn測(cè)試我們的模型來(lái)結(jié)束。

SVM概述

支持向量機(jī)的目標(biāo)是擬合獲得最大邊緣的超平面(兩個(gè)類(lèi)中最近點(diǎn)的距離)。可以直觀地表明，這樣的超平面(A)比沒(méi)有最大化邊際的超平面(B)具有更好的泛化特性和對(duì)噪聲的魯棒性。

為了實(shí)現(xiàn)這一點(diǎn)，SVM通過(guò)求解以下優(yōu)化問(wèn)題找到超平面的W和b:

它試圖找到W,b，使最近點(diǎn)的距離最大化，并正確分類(lèi)所有內(nèi)容(如y取±1的約束)。這可以被證明相當(dāng)于以下優(yōu)化問(wèn)題:

可以寫(xiě)出等價(jià)的對(duì)偶優(yōu)化問(wèn)題

這個(gè)問(wèn)題的解決方案產(chǎn)生了一個(gè)拉格朗日乘數(shù)，我們假設(shè)數(shù)據(jù)集中的每個(gè)點(diǎn)的大小為m:(α 1， α 2，…，α _n)。目標(biāo)函數(shù)在α中明顯是二次的，約束是線(xiàn)性的，這意味著它可以很容易地用二次規(guī)劃求解。一旦找到解，由對(duì)偶的推導(dǎo)可知:

注意，只有具有α>0的點(diǎn)才定義超平面(對(duì)和有貢獻(xiàn))。這些被稱(chēng)為支持向量。因此當(dāng)給定一個(gè)新例子x時(shí)，返回其預(yù)測(cè)y=±1的預(yù)測(cè)方程為:

這種支持向量機(jī)的基本形式被稱(chēng)為硬邊界支持向量機(jī)（hard margin SVM），因?yàn)樗鉀Q的優(yōu)化問(wèn)題(如上所述)強(qiáng)制要求訓(xùn)練中的所有點(diǎn)必須被正確分類(lèi)。但在實(shí)際場(chǎng)景中，可能存在一些噪聲，阻止或限制了完美分離數(shù)據(jù)的超平面，在這種情況下，優(yōu)化問(wèn)題將不返回或返回一個(gè)糟糕的解決方案。

軟邊界支持向量機(jī)（soft margin SVM）通過(guò)引入C常數(shù)(用戶(hù)給定的超參數(shù))來(lái)適應(yīng)優(yōu)化問(wèn)題，該常數(shù)控制它應(yīng)該有多“硬”。特別地，它將原優(yōu)化問(wèn)題修改為:

它允許每個(gè)點(diǎn)產(chǎn)生一些錯(cuò)誤λ(例如，在超平面的錯(cuò)誤一側(cè))，并且通過(guò)將它們?cè)谀繕?biāo)函數(shù)中的總和加權(quán)C來(lái)減少它們。當(dāng)C趨于無(wú)窮時(shí)(一般情況下肯定不會(huì))，它就等于硬邊界。與此同時(shí)，較小的C將允許更多的“違規(guī)行為”(以換取更大的支持;例如，更小的w (w)。

可以證明，等價(jià)對(duì)偶問(wèn)題只有在約束每個(gè)點(diǎn)的α≤C時(shí)才會(huì)發(fā)生變化。

由于允許違例，支持向量(帶有α>0的點(diǎn))不再都在邊界的邊緣。任何錯(cuò)誤的支持向量都具有α=C，而非支持向量(α=0)不能發(fā)生錯(cuò)誤。我們稱(chēng)潛在錯(cuò)誤(α=C)的支持向量為“非錯(cuò)誤編劇支持向量”和其他純粹的支持向量(沒(méi)有違規(guī);“邊界支持向量”(0<α

這樣推理方程不變:

現(xiàn)在(x?，y?)必須是一個(gè)沒(méi)有違規(guī)的支持向量，因?yàn)榉匠碳僭O(shè)它在邊界的邊緣。

軟邊界支持向量機(jī)擴(kuò)展了硬邊界支持向量機(jī)來(lái)處理噪聲，但通常由于噪聲以外的因素，例如自然非線(xiàn)性，數(shù)據(jù)不能被超平面分離。軟邊界支持向量機(jī)可以用于這樣的情況，但是最優(yōu)解決方案的超平面，它允許的誤差遠(yuǎn)遠(yuǎn)超過(guò)現(xiàn)實(shí)中可以容忍的誤差。

例如，在左邊的例子中，無(wú)論C的設(shè)置如何，軟邊界支持向量機(jī)都找不到線(xiàn)性超平面。但是可以通過(guò)某種轉(zhuǎn)換函數(shù)z=Φ(x)將數(shù)據(jù)集中的每個(gè)點(diǎn)x映射到更高的維度，從而使數(shù)據(jù)在新的高維空間中更加線(xiàn)性(或完全線(xiàn)性)。這相當(dāng)于用z替換x得到:

在現(xiàn)實(shí)中，特別是當(dāng)Φ轉(zhuǎn)換為非常高維的空間時(shí)，計(jì)算z可能需要很長(zhǎng)時(shí)間。所以就出現(xiàn)了核函數(shù)。它用一個(gè)數(shù)學(xué)函數(shù)(稱(chēng)為核函數(shù))的等效計(jì)算來(lái)取代z，并且更快(例如，對(duì)z進(jìn)行代數(shù)簡(jiǎn)化)。例如，這里有一些流行的核函數(shù)(每個(gè)都對(duì)應(yīng)于一些轉(zhuǎn)換Φ到更高維度空間):

這樣，對(duì)偶優(yōu)化問(wèn)題就變成:

直觀地，推理方程(經(jīng)過(guò)代數(shù)處理后)為:

上面所有方程的完整推導(dǎo)，有很多相關(guān)的文章了，我們就不詳細(xì)介紹了。

Python實(shí)現(xiàn)

對(duì)于實(shí)現(xiàn)，我們將使用下面這些庫(kù)：

import numpy as np                  # for basic operations over arrays
 from scipy.spatial import distance  # to compute the Gaussian kernel
 import cvxopt                       # to solve the dual opt. problem
 import copy                         # to copy numpy arrays

定義核和SVM超參數(shù)，我們將實(shí)現(xiàn)常見(jiàn)的三個(gè)核函數(shù)：

class SVM:
     linear = lambda x, x? , c=0: x @ x?.T
     polynomial = lambda x, x? , Q=5: (1 + x @ x?.T)**Q
     rbf = lambda x, x?, γ=10: np.exp(-γ*distance.cdist(x, x?,'sqeuclidean'))
     kernel_funs = {'linear': linear, 'polynomial': polynomial, 'rbf': rbf}

為了與其他核保持一致，線(xiàn)性核采用了一個(gè)額外的無(wú)用的超參數(shù)。kernel_funs接受核函數(shù)名稱(chēng)的字符串，并返回相應(yīng)的內(nèi)核函數(shù)。

繼續(xù)定義構(gòu)造函數(shù):

class SVM:
     linear = lambda x, x? , c=0: x @ x?.T
     polynomial = lambda x, x? , Q=5: (1 + x @ x?.T)**Q
     rbf = lambda x, x?, γ=10: np.exp(-γ*distance.cdist(x, x?,'sqeuclidean'))
     kernel_funs = {'linear': linear, 'polynomial': polynomial, 'rbf': rbf}
     
     def __init__(self, kernel='rbf', C=1, k=2):
         # set the hyperparameters
         self.kernel_str = kernel
         self.kernel = SVM.kernel_funs[kernel]
         self.C = C                  # regularization parameter
         self.k = k                  # kernel parameter
         
         # training data and support vectors (set later)
         self.X, y = None, None
         self.αs = None
         
         # for multi-class classification (set later)
         self.multiclass = False
         self.clfs = []

SVM有三個(gè)主要的超參數(shù)，核(我們存儲(chǔ)給定的字符串和相應(yīng)的核函數(shù))，正則化參數(shù)C和核超參數(shù)(傳遞給核函數(shù));它表示多項(xiàng)式核的Q和RBF核的γ。

為了兼容sklearn的形式，我們需要使用fit和predict函數(shù)來(lái)擴(kuò)展這個(gè)類(lèi)，定義以下函數(shù)，并在稍后將其用作裝飾器:

SVMClass = lambda func: setattr(SVM, func.__name__, func) or func

擬合SVM對(duì)應(yīng)于通過(guò)求解對(duì)偶優(yōu)化問(wèn)題找到每個(gè)點(diǎn)的支持向量α:

設(shè)α為可變列向量(α?α?…α _n);y為標(biāo)簽(y?α?…y_N)常數(shù)列向量;K為常數(shù)矩陣，其中K[n,m]計(jì)算核在(x, x)處的值。點(diǎn)積、外積和二次型分別基于索引的等價(jià)表達(dá)式:

可以將對(duì)偶優(yōu)化問(wèn)題寫(xiě)成矩陣形式如下:

這是一個(gè)二次規(guī)劃，CVXOPT的文檔中解釋如下:

可以只使用(P,q)或(P,q,G,h)或(P,q,G,h, A, b)等等來(lái)調(diào)用它(任何未給出的都將由默認(rèn)值設(shè)置，例如1)。

對(duì)于(P, q, G, h, A, b)的值，我們的例子可以做以下比較:

為了便于比較，將第一個(gè)重寫(xiě)如下:

現(xiàn)在很明顯(0≤α等價(jià)于-α≤0):

我們就可以寫(xiě)出如下的fit函數(shù):

@SVMClass
 def fit(self, X, y, eval_train=False):
     # if more than two unique labels, call the multiclass version
     if len(np.unique(y)) > 2:
         self.multiclass = True
         return self.multi_fit(X, y, eval_train)
     
     # if labels given in {0,1} change it to {-1,1}
     if set(np.unique(y)) == {0, 1}: y[y == 0] = -1
 
     # ensure y is a Nx1 column vector (needed by CVXOPT)
     self.y = y.reshape(-1, 1).astype(np.double) # Has to be a column vector
     self.X = X
     N = X.shape[0]  # Number of points
     
     # compute the kernel over all possible pairs of (x, x') in the data
     # by Numpy's vectorization this yields the matrix K
     self.K = self.kernel(X, X, self.k)
     
     ### Set up optimization parameters
     # For 1/2 x^T P x + q^T x
     P = cvxopt.matrix(self.y @ self.y.T * self.K)
     q = cvxopt.matrix(-np.ones((N, 1)))
     
     # For Ax = b
     A = cvxopt.matrix(self.y.T)
     b = cvxopt.matrix(np.zeros(1))
 
     # For Gx <= h
     G = cvxopt.matrix(np.vstack((-np.identity(N),
                                  np.identity(N))))
     h = cvxopt.matrix(np.vstack((np.zeros((N,1)),
                                  np.ones((N,1)) * self.C)))
 
     # Solve    
     cvxopt.solvers.options['show_progress'] = False
     sol = cvxopt.solvers.qp(P, q, G, h, A, b)
     self.αs = np.array(sol["x"])            # our solution
         
     # a Boolean array that flags points which are support vectors
     self.is_sv = ((self.αs-1e-3 > 0)&(self.αs <= self.C)).squeeze()
     # an index of some margin support vector
     self.margin_sv = np.argmax((0 < self.αs-1e-3)&(self.αs < self.C-1e-3))
     
     if eval_train:  
       print(f"Finished training with accuracy{self.evaluate(X, y)}")

我們確保這是一個(gè)二進(jìn)制問(wèn)題，并且二進(jìn)制標(biāo)簽按照支持向量機(jī)(±1)的假設(shè)設(shè)置，并且y是一個(gè)維數(shù)為(N,1)的列向量。然后求解求解(α?α?…α _n) 的優(yōu)化問(wèn)題。

使用(α?α?…α _n) _來(lái)獲得在與支持向量對(duì)應(yīng)的任何索引處為1的標(biāo)志數(shù)組，然后可以通過(guò)僅對(duì)支持向量和(x?，y?)的邊界支持向量的索引求和來(lái)應(yīng)用預(yù)測(cè)方程。我們確實(shí)假設(shè)非支持向量可能不完全具有α=0，如果它的α≤1e-3，那么這是近似為零(CVXOPT結(jié)果可能不是最終精確的)。同樣假設(shè)非邊際支持向量可能不完全具有α=C。

下面就是預(yù)測(cè)的方法，預(yù)測(cè)方程為:

@SVMClass
 def predict(self, X_t):
     if self.multiclass: return self.multi_predict(X_t)
     # compute (x?, y?)
     x?, y? = self.X[self.margin_sv, np.newaxis], self.y[self.margin_sv]
     # find support vectors
     αs, y, X= self.αs[self.is_sv], self.y[self.is_sv], self.X[self.is_sv]
     # compute the second term
     b = y? - np.sum(αs * y * self.kernel(X, x?, self.k), axis=0)
     # compute the score
     score = np.sum(αs * y * self.kernel(X, X_t, self.k), axis=0) + b
     return np.sign(score).astype(int), score

我們還可以實(shí)現(xiàn)一個(gè)評(píng)估方法來(lái)計(jì)算精度(在上面的fit中使用)。

@SVMClass
 def evaluate(self, X,y):  
     outputs, _ = self.predict(X)
     accuracy = np.sum(outputs == y) / len(y)
     return round(accuracy, 2)

最后測(cè)試我們的完整代碼：

from sklearn.datasets import make_classification
 import numpy as np
 
 # Load the dataset
 np.random.seed(1)
 X, y = make_classification(n_samples=2500, n_features=5, 
                            n_redundant=0, n_informative=5, 
                            n_classes=2,  class_sep=0.3)
 
 # Test Implemented SVM
 svm = SVM(kernel='rbf', k=1)
 svm.fit(X, y, eval_train=True)
 
 y_pred, _ = svm.predict(X)
 print(f"Accuracy: {np.sum(y==y_pred)/y.shape[0]}")  #0.9108
 
 # Test with Scikit
 from sklearn.svm import SVC
 clf = SVC(kernel='rbf', C=1, gamma=1)
 clf.fit(X, y)
 y_pred = clf.predict(X)
 print(f"Accuracy: {sum(y==y_pred)/y.shape[0]}")    #0.9108

多分類(lèi)SVM

我們都知道SVM的目標(biāo)是二元分類(lèi)，如果要將模型推廣到多類(lèi)則需要為每個(gè)類(lèi)訓(xùn)練一個(gè)二元SVM分類(lèi)器，然后對(duì)每個(gè)類(lèi)進(jìn)行循環(huán)，并將屬于它的點(diǎn)重新標(biāo)記為+1，并將所有其他類(lèi)的點(diǎn)重新標(biāo)記為-1。

當(dāng)給定k個(gè)類(lèi)時(shí)，訓(xùn)練的結(jié)果是k個(gè)分類(lèi)器，其中第i個(gè)分類(lèi)器在數(shù)據(jù)上進(jìn)行訓(xùn)練，第i個(gè)分類(lèi)器被標(biāo)記為+1，所有其他分類(lèi)器被標(biāo)記為-1。

@SVMClass
 def multi_fit(self, X, y, eval_train=False):
     self.k = len(np.unique(y))      # number of classes
     # for each pair of classes
     for i in range(self.k):
         # get the data for the pair
         Xs, Ys = X, copy.copy(y)
         # change the labels to -1 and 1
         Ys[Ys!=i], Ys[Ys==i] = -1, +1
         # fit the classifier
         clf = SVM(kernel=self.kernel_str, C=self.C, k=self.k)
         clf.fit(Xs, Ys)
         # save the classifier
         self.clfs.append(clf)
     if eval_train:  
         print(f"Finished training with accuracy {self.evaluate(X, y)}")

然后，為了對(duì)新示例執(zhí)行預(yù)測(cè)，我們選擇相應(yīng)分類(lèi)器最自信(得分最高)的類(lèi)。

@SVMClass
 def multi_predict(self, X):
     # get the predictions from all classifiers
     N = X.shape[0]
     preds = np.zeros((N, self.k))
     for i, clf in enumerate(self.clfs):
         _, preds[:, i] = clf.predict(X)
     
     # get the argmax and the corresponding score
     return np.argmax(preds, axis=1), np.max(preds, axis=1)

完整測(cè)試代碼：

from sklearn.datasets import make_classification
 import numpy as np
 
 # Load the dataset
 np.random.seed(1)
 X, y = make_classification(n_samples=500, n_features=2, 
                            n_redundant=0, n_informative=2, 
                            n_classes=4, n_clusters_per_class=1,  
                            class_sep=0.3)
 
 # Test SVM
 svm = SVM(kernel='rbf', k=4)
 svm.fit(X, y, eval_train=True)
 
 y_pred = svm.predict(X)
 print(f"Accuracy: {np.sum(y==y_pred)/y.shape[0]}") # 0.65
 
 # Test with Scikit
 from sklearn.multiclass import OneVsRestClassifier
 from sklearn.svm import SVC
 
 clf = OneVsRestClassifier(SVC(kernel='rbf', C=1, gamma=4)).fit(X, y)
 y_pred = clf.predict(X)
 print(f"Accuracy: {sum(y==y_pred)/y.shape[0]}")    # 0.65

繪制每個(gè)決策區(qū)域的圖示，得到以下圖:

可以看到，我們的實(shí)現(xiàn)與Sci-kit Learn結(jié)果相當(dāng)，說(shuō)明在算法實(shí)現(xiàn)上沒(méi)有問(wèn)題。注意：SVM默認(rèn)支持OVR(沒(méi)有如上所示的顯式調(diào)用)，它是特定于SVM的進(jìn)一步優(yōu)化。

總結(jié)

我們使用Python實(shí)現(xiàn)了支持向量機(jī)(SVM)學(xué)習(xí)算法，并且包括了軟邊界和常用的三個(gè)核函數(shù)。我們還將SVM擴(kuò)展到多分類(lèi)的場(chǎng)景，并使用Sci-kit Learn驗(yàn)證了我們的實(shí)現(xiàn)。希望通過(guò)本文你可以更好的了解SVM。