統計學習：線性支持向量機(Pytorch實現)

1 學習策略

1.1 軟間隔最大化

上一章我們所定義的“線性可分支持向量機”要求訓練數據是線性可分的。然而在實際中，訓練數據往往包括異常值(outlier)，故而常是線性不可分的。這就要求我們要對上一章的算法做出一定的修改，即放寬條件，將原始的硬間隔最大化轉換為軟間隔最大化。
給定訓練集

\[\begin{aligned} D = \{\{\bm{x}^{(1)}, y^{(1)}\}, \{\bm{x}^{(2)}, y^{(2)}\},..., \{\bm{x}^{(m)}, y^{(m)}\}\} \end{aligned}\tag{1} \]

其中\(\bm{x}^{(i)} \in \mathcal{X} \subseteq \mathbb{R}^n\)，\(y^{(i)} \in \mathcal{Y} = \{+1, -1\}\)。
如果訓練集是線性可分的，則線性可分支持向量機等價於求解以下凸優化問題：

\[\begin{aligned} \underset{\bm{w}, b}{\max} \quad \frac{1}{2} || \bm{w}||^2\\ \text{s.t.} \quad y^{(i)} (\bm{w}^T \bm{x}^{(i)} + b) \geqslant 1 \\ \quad (i = 1, 2, ..., m) \end{aligned} \tag{2} \]

其中\(y^{(i)} (\bm{w}^T \bm{x}^{(i)} + b) -1 \geq 0\)表示樣本點\((\bm{x}^{(i)}, y^{(i)})\)滿足函數間隔大於等於1。現在我們對每個樣本點\((\bm{x}^{(i)}, y^{(i)})\)放寬條件，引入一個松弛變量\(\xi_{i} \geqslant 0\),使約束條件變為\(y^{(i)} (\bm{w}^T \bm{x}^{(i)} + b) \geq 1-\xi_{i}\)。並對每個松弛變量進行一個大小為\(\xi_{i}\)的代價懲罰，目標函數轉變為：\(\frac{1}{2} || \bm{w}||^2+C\sum_{i=1}^{m}\xi_{i}\)，此處\(C>0\)稱為懲罰系數。此時優化函數即要使間隔盡量大（使\(\frac{1}{2} || \bm{w}||^2\)盡量小）,又要使誤分類點個數盡量少。這稱之為軟間隔化。

1.2 線性支持向量機

就這樣，線性支持向量機變為如下凸二次規划問題(原始問題)：

\[\begin{aligned} \underset{\bm{w}, b}{\max} \quad \frac{1}{2} || \bm{w}||^2 + C\sum_{i=1}^{m}\xi_{i}\\ \text{s.t.} \quad y^{(i)} (\bm{w}^T \bm{x}^{(i)} + b) \geqslant 1-\xi_{i} \\ \xi_{i} \geqslant 0 \\ \quad (i = 1, 2, ..., m) \end{aligned} \tag{3} \]

因為是凸二次規划，因此關於\((\bm{w}, b, \bm{\xi})\)的解一定存在，可以證明\(\bm{w}\)的解唯一，但\(b\)的解可能不唯一，而是存在於一個區間。
設\((2)\)的解為\(\bm{w}^{*}, b^*\)，這樣可得到分離超平面\(\{\bm{x} | \bm{w}^{*T}\bm{x}+b=0\}\)和分類決策函數\(f(\bm{x})=\text{sign}(\bm{w}^{*T}\bm{x}+b^*)\)

2 算法

2.1 常規的帶約束優化算法

和上一章一樣，我們將原始問題\((2)\)轉換為對偶問題進行求解，原始問題的對偶問題如下（推導和上一章一樣）：

\[\begin{aligned} \underset{\bm{\alpha}}{\max} -\frac{1}{2}\sum_{i=1}^{m}\sum_{j=1}^{m}\alpha_i\alpha_jy^{(i)}y^{(j)}\langle \bm{x}^{(i)}, \bm{x}^{(j)}\rangle + \sum_{i=1}^{m}\alpha_i \\ s.t. \sum_{i=1}^{m}\alpha_iy^{(i)} = 0 \\ C - \alpha_i - \mu_i = 0\\ \alpha_i \geq 0 \\ \mu_i \geq 0 \\ i=1, 2,...,m \end{aligned} \tag{4} \]

接下來我們消去\(\mu_i\)，只留下\(\alpha_i\)，將約束條件轉換為：\(0 \leqslant \alpha_i \leqslant C\) ，並對目標函數求極大轉換為求極小，這樣對\((4)\)進行進一步變換得到：

\[\begin{aligned} \underset{\bm{\alpha}}{\max} -\frac{1}{2}\sum_{i=1}^{m}\sum_{j=1}^{m}\alpha_i\alpha_jy^{(i)}y^{(j)}\langle \bm{x}^{(i)}, \bm{x}^{(j)}\rangle + \sum_{i=1}^{m}\alpha_i \\ s.t. \sum_{i=1}^{m}\alpha_iy^{(i)} = 0 \\ 0 \leqslant \alpha_i \leqslant C \\ i=1, 2,...,m \end{aligned} \tag{5} \]

這就是原始優化問題\((3)\)的對偶形式。

這樣，我們得到對偶問題的解為\(\bm{\alpha}^*=(\alpha_1^*, \alpha_2^*, ..., \alpha_m^*)^T\)。
此時情況要比完全線性可分時要復雜，比如正例不再僅僅分布於軟間隔平面上和正例那一側，還可能分布於負例對應的那一側（由於誤分）；對於負例亦然。如下圖所示：

圖中實線為超平面，虛線為間隔邊界，“\(\circ\)”為正例點，“\(\times\)”為負例點，\(\frac{\xi_i}{||\bm{w}||}\)為實例\(\bm{x}^{(i)}\)到間隔邊界的距離。
此時類比上一章，我們將\(\alpha_i^* > 0\)（注意，不要求一定有\(\alpha_i^*<C\)）對應樣本點\((\bm{x}^{(i)}, y^{(i)})\)的實例\(\bm{x}^{(i)}\)稱為支持向量（軟間隔支持向量）。它的分布比硬間隔情況下要復雜得多：可以在間隔邊界上，也可以在間隔邊界和分離超平面之間，甚至可以在間隔超平面誤分類的一側。若\(0<\alpha_i^*<C\)，則\(\xi_i =0\)，支持向量恰好落在間隔邊界上；若\(\alpha_i = C, 0<\xi_i<1\)，則分類正確且\(\bm{x}^{(i)}\)在間隔邊界和分離超平面之間；若\(\alpha_i^*=C, \xi_i=1\)，則\(\bm{x}^i\)在分離超平面上；若\(\alpha_i^*=C, \xi_i>1\)則分類錯誤，\(\bm{x}^{(i)}\)分離超平面誤分類的那一側。

接下來我們需要將對偶問題的解轉換為原始問題的解。

對於對偶問題的解\(\bm{\alpha}^*=(\alpha_1^*, \alpha_2^*, ..., \alpha_m^*)^T\),若存在\(\alpha^*\)的一個分量\(0 <\alpha_s^* < C\)，則原始問題的解\(w^*\)和\(b^*\)可以按下式求得：

\[\begin{aligned} \bm{w}^{*} = \sum_{i=1}^{m}\alpha_i^{*}y^{(i)}\bm{x}^{(i)} \\ b^{*} = y^{(s)} - \sum_{i=1}^{m}\alpha_i^*y^{(i)}\langle \bm{x}^{(s)}, \bm{x}^{(i)} \rangle \end{aligned} \tag{6} \]

式 \((11)\) 的推導方法和上一章類似，由原始問題(凸二次規划問題)的解滿足KKT條件導出。


這樣，我們就可以得到分離超平面如下：

\[\{\bm{x} | \sum_{i=1}^{m}\alpha_i^*y^{(i)}\langle \bm{x}^{(i)}, \bm{x} \rangle + b^* = 0\} \tag{7} \]

分類決策函數如下:

\[f(\bm{x}) = \text{sign}(\sum_{i=1}^{m}\alpha_i^*y^{(i)}\langle \bm{x}^{(i)}, \bm{x} \rangle + b^*) \tag{8} \]

(同樣地，之所以寫成\(\langle \bm{x}^{(i)}, \bm{x} \rangle\)的內積形式是為了方便我們后面引入核函數)

綜上，按照式\((3)\)的策略求解線性可分支持向量機的算法如下：　

可以看到在算法的步驟\((2)\)中，是對隨機采的一個\(\alpha_s^*\)計算出\(b^*\)，故按照式\((3)\)的策略（原始問題）求解出的\(b\)可能不唯一。

2.2 基於合頁損失函數的的無約束優化算法

其實，問題\((3)\)還可以寫成無約束優化的形式，目標函數如下：

\[ \underset{\bm{w}, b}{\min} \space \frac{1}{m}\sum_{i=1}^{m}\max(0, 1-y^{(i)}(\bm{w}^{T}\bm{x}^{(i)} + b))+\lambda||\bm{w}||^2 \tag{9} \]

式\((8)\)的第一項為經驗風險或經驗損失（加上正則項整體作為結構風險）。函數

\[ L(y(\bm{w}^{T}\bm{x} + b))=\max(0, 1-y(\bm{w}^{T}\bm{x} + b)) \tag{10} \]

稱為合頁損失函數（hinge loss function）。s
合頁損失函數意味着：如果樣本\((\bm{x}^{(i)}, y^{(i)})\)被正確分類且函數間隔（確信度）大於1（即\(y^{(i)}(\bm{w}^{T}\bm{x}^{(i)} + b)>0\)），則損失是0，否則損失是\(1-y^{(i)}(\bm{w}^{T}\bm{x}^{(i)} + b)\)。像上面提到的分類圖中，\(\bm{x}^{(4)}\)被正確分類，但函數間隔不大於1，損失不是0。
0-1損失函數、合頁損失函數、感知機損失函數歸納如下：

\[\left\{ \begin{aligned} &\quad \text{I}(y(\bm{w}^{T}\bm{x}+b)<0) \quad \quad (0-1損失函數)\\ &\quad \max(0, 1-y(\bm{w}^{T}\bm{x} + b)) \quad \quad (合頁損失函數) \\ &\quad \max(0, -y(\bm{w}^{T}\bm{x} + b)) \quad \quad (感知機損失函數) \end{aligned} \right . \tag{11} \]

這三個函數的圖像對比如下圖所示：

可以看到合頁損失函數形似合頁，故而得名。而0-1損失函數雖然是二分類問題真正的損失函數，但它不是連續可導的，對其進行優化是NP困難的，所以我們轉而優化其上界（合頁損失函數）構成的目標函數，這時上界損失函數又稱為代理損失函數（surrogate loss function）。而對於感知機損失函數，直觀理解為：當\((\bm{x}^{(i)}, y^{(i)})\)被正確分類時，損失是0，否則是\(-y(\bm{w}^{T}\bm{x} + b)\)。相比之下，合頁損失函數不僅要求要正確分類，而且要確信度足夠大時損失才是0，也就是說合頁損失函數對學習效果有更高的要求。
合頁損失函數處處連續，此時可以采用基於梯度的數值優化算法求解（梯度下降法、牛頓法等），在此不再贅述。不過，此時的目標函數非凸，不一定保證收斂到最優解。

3 編程實現

上一章我們嘗試過在不借助SMO算法的情況下采用凸優化算法求解問題\((3)\)這一凸二次規划問題，結果發現收斂速度過慢。現在我們試試采用基於合頁損失函數的梯度下降算法來直接求解參數\(\bm{w}\)和\(b\)。

from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
import numpy as np
import math
import torch

# 數據預處理
def preprocess(X, y):
    # 對X進行min max標准化
    for i in range(X.shape[1]):
        field = X[:, i]
        X[:, i] = (field - field.min()) / (field.max() - field.min()) 
    # 標簽統一轉化為1和-1
    y = np.where(y==1, y, -1)
    return X, y

class SVM():
    def __init__(self, lamb=0.01, input_dim=256):
        self.lamb = lamb
        self.w, self.b = np.random.rand(
            input_dim), 0.0

    # 注意一個批量的loss可並行化計算
    def objective_func(self, X, y):
        # 注意此處的Hinge loss函數函數不是在每個點都可導
        loss_vec = 1 - torch.mul(y, torch.matmul(X, self.w) + self.b)
        loss_vec = torch.where(loss_vec > 0, loss_vec, 0.)
        return 1/self.m * loss_vec.sum() + self.lamb * torch.norm(self.w)**2 

    def train_data_loader(self, X_train, y_train):
        # 每輪迭代隨機采一個batch
        data = np.concatenate((X_train, y_train.reshape(-1, 1)), axis=1)                                                    
        np.random.shuffle(data)
        X_train, y_train = data[:, :-1], data[:, -1]
        for ep in range(self.epoch):
            for bz in range(math.ceil(self.m/self.batch_size)):
                start = bz * self.batch_size
                yield (
                    torch.tensor(X_train[start: start+self.batch_size], dtype=torch.float64), 
                torch.tensor(y_train[start: start+self.batch_size],dtype=torch.float64))

    def test_data_loader(self, X_test):
        # 每輪迭代隨機采一個batch
        for bz in range(math.ceil(self.m_t/self.test_batch_size)):
            start = bz * self.test_batch_size
            yield X_test[start: start+self.test_batch_size]

    def compile(self, **kwargs):
        self.batch_size = kwargs['batch_size']
        self.test_batch_size = kwargs['test_batch_size']
        self.eta = kwargs['eta'] # 學習率
        self.epoch = kwargs['epoch'] # 遍歷多少次訓練集

    def sgd(self, params):
        with torch.no_grad(): 
            for param in params:
                param -= self.eta*param.grad
                param.grad.zero_()

    def fit(self, X_train, y_train):
        self.m = X_train.shape[0] #樣本個數
        # 主義w初始化為隨機數，不能初始化為0
        self.w, self.b = torch.tensor(
            self.w, dtype=torch.float64, requires_grad=True), torch.tensor(self.b, requires_grad=True)
        for X, y in self.train_data_loader(X_train, y_train):
            loss_v = self.objective_func(X, y)
            loss_v.backward() 
            self.sgd([self.w, self.b])

        self.w = self.w.detach().numpy()
        self.b = self.b.detach().numpy()

    def pred(self, X_test):
        # 遍歷測試集中的每一個樣本
        self.m_t = X_test.shape[0]
        pred_list = []
        for x in self.test_data_loader(X_test):
            pred_list.append(np.sign(np.matmul(x, self.w) + self.b))
        return np.concatenate(pred_list, axis=0)
                
if __name__ == "__main__":
    X, y = load_breast_cancer(return_X_y=True)
    X, y = preprocess(X, y)
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
    clf = SVM(lamb=0.01, input_dim=X_train.shape[1])
    clf.compile(batch_size=256, test_batch_size=1024, eta=0.001, epoch=1000) #定義訓練參數
    clf.fit(X_train, y_train)
    y_pred = clf.pred(X_test)
    acc_score = accuracy_score(y_test, y_pred)
    print("The accuracy is: %.1f" % acc_score)

可以看到，收斂效果不太好：

The accuracy is: 0.6

經過分析，我們發現主要原因還是因為SVM的損失函數

\[ \max(0, 1-y\bm{w}^{T}\bm{x}) \]

(注意，我們采用簡法表示合並了偏置\(b\)，即\(\bm{w}\)最后一維是偏置且\(\bm{x}\)擴展一維)並不是在每個點都可導，直接用Pytorch來求解難以收斂。我們可以考慮使用其次梯度(subgradient)公式

\[\left\{ \begin{aligned} & - y \bm{x} \quad \text{if} \space y \bm{w}^T\bm{x} < 1 \\ & 0 \quad \quad \quad\text{otherwise} \end{aligned} \right. \]

使用次梯度進行求解的代碼如下：

from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
import numpy as np
from random import choice

# 數據預處理
def preprocess(X, y):
    # 對X進行min max標准化
    for i in range(X.shape[1]):
        field = X[:, i]
        X[:, i] = (field - field.min()) / (field.max() - field.min()) 
    # 標簽統一轉化為1和-1
    y = np.where(y==1, y, -1)
    return X, y

class SVM():
    def __init__(self, lamb, input_dim):
        self.lamb = lamb
        self.w = np.random.ranf(
            input_dim)

    def compile(self, **kwargs):
        self.eta = kwargs['eta'] # 學習率
        self.epoch = kwargs['epoch'] # 遍歷多少次訓練集

    def gradient(self, x, y, w):
        return - y * x if y * w.dot(x) < 1 else 0

    def fit(self, X_train, y_train):
        # 主義w初始化為隨機數，不能初始化為0
        for x, y in zip(X_train, y_train):
            grad = self.gradient(x, y, self.w)
            self.w -= self.eta*grad

    def pred(self, X_test):
        # 遍歷測試集中的每一個樣本
        pred_list = []
        for x in X_test:
            pred_list.append(np.sign(np.matmul(x, self.w)))
        return pred_list
                
if __name__ == "__main__":
    X, y = load_breast_cancer(return_X_y=True)
    X, y = preprocess(X, y)
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
    D = X_train.shape[1]
    
    # 擴充1維
    X_train = np.concatenate([X_train, np.ones([X_train.shape[0], 1])], axis=1)
    X_test = np.concatenate([X_test, np.ones([X_test.shape[0], 1])], axis=1)

    clf = SVM(lamb=0.01, input_dim=(D + 1))
    clf.compile(eta=0.01, epoch=1000) #定義訓練參數
    clf.fit(X_train, y_train)
    y_pred = clf.pred(X_test)

    acc_score = accuracy_score(y_test, y_pred)
    print("The accuracy is: %.1f" % acc_score)

可以看到，還是難以收斂：

The accuracy is: 0.6

看來用梯度求解Hinge Loss不是最佳的方法，而我們上一篇博客《統計學習：線性可分支持向量機(Scipy實現)》 也證明了用優化工具包直接求解對偶問題效率也很低。那么路在何方？接下來我們會介紹求解SVM的SMO算法，這個才是求解SVM的殺手鐧。

參考

• [1] 李航. 統計學習方法(第2版)[M]. 清華大學出版社, 2019.
• [2] 周志華. 機器學習[M]. 清華大學出版社, 2016.
• [3] Boyd S, Boyd S P, Vandenberghe L. Convex optimization[M]. Cambridge university press, 2004.