上一章我們所定義的“線性可分支援向量機”要求訓練資料是線性可分的。然而在實際中，訓練資料往往包括異常值(outlier)，故而常是線性不可分的。這就要求我們要對上一章的演算法做出一定的修改，即放寬條件，將原始的硬間隔最大化轉換為軟間隔最大化。該問題最終可以寫成無約束優化的形式，目標函式由合頁損失函式和正則項構成。合頁損失函式處處連續，此時可以採用基於梯度的數值優化演算法求解（梯度下降法、牛頓法等）不過，此時的目標函式非凸，不一定保證收斂到最優解。

學習策略

軟間隔最大化

上一章我們所定義的“線性可分支援向量機”要求訓練資料是線性可分的。然而在實際中，訓練資料往往包括異常值(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\)儘量小）,又要使誤分類點個數儘量少。這稱之為軟間隔化。

線性支援向量機

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

\[\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)\)轉換為對偶問題進行求解，原始問題的對偶問題如下（推導和上一章一樣）：

\[\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\)可能不唯一。

基於合頁損失函式的的無約束優化演算法

其實，問題\((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損失函式、合頁損失函式、感知機損失函式歸納如下：

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

這三個函式的影象對比如下圖所示：

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

演算法

上一章我們嘗試過在不借助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
from scipy.optimize import minimize
import numpy as np
import math
from copy import deepcopy
import torch
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=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):
        loss_vec = 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() 
            # b有梯度，w沒有梯度？
            self.sgd([self.w, self.b])
            # print(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.5, 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)

可以看到，效果不太好，哈哈，目前還在debug

The accuracy is: 0.6

數學是符號的藝術，音樂是上界的語言。