【機器學習實戰】-- Titanic 資料集（4）-- 支援向量機
1. 寫在前面:
本篇屬於實戰部分，更注重於演算法在實際專案中的應用。如需對感知機演算法本身有進一步的瞭解，可參考以下連結，在本人學習的過程中，起到了很大的幫助：
統計學習方法 李航
支援向量機原理https://www.cnblogs.com/pinard/p/6097604.html 等共5篇
空間中任意一點到超平面距離的公式推導 https://www.cnblogs.com/yanganling/p/8007050.html
拉格朗日對偶https://www.cnblogs.com/ooon/p/5723725.html
2. 資料集：
資料集地址：https://www.kaggle.com/c/titanic
Titanic資料集是Kaggle上參與人數最多的專案之一。資料本身簡單小巧，適合初學者上手，深入瞭解比較各個機器學習演算法。
資料集包含11個變數：PassengerID、Pclass、Name、Sex、Age、SibSp、Parch、Ticket、Fare、Cabin、Embarked，通過這些資料來預測乘客在Titanic事故中是否倖存下來。
3. 演算法簡介：
支援向量機原本是一種分類模型，和感知機一樣，定義了特徵空間中的一個超平面，但間隔最大化使其區別於感知機（通常是優於）。此外，和感知機一樣，利用Gram矩陣，方便地使用核函式，使其可以成為一個非線性分類器。
支援向量機的學習策略就是間隔最大化，並且等價與合頁損失函式（hinge loss）的最小化問題。
3.1 線性可分支援向量機模型：
給定一個數據集：$T=\left \{ \left ( x_{1}, y_{1} \right ), \left ( x_{2}, y_{2} \right ), ..., \left ( x_{N}, y_{N} \right ) \right \}$，其中$x_{i}\in \chi\subseteq \bf{R^{n}}$，$y_{i} \in \Upsilon= \left \{+1, -1 \right \}$，$i = 1,2,...,N$。這代表資料集共有 N 對 例項，每個例項 $x_{i}$都是n維的。
假設通過間隔最大化求解得到的分離超平面為：
$$w^{\ast } \cdot x + b^{\ast} = 0 $$
則如下分類函式即線性可分支援向量機模型：
$$ f(x) = \text{sign}(w^{\ast} \cdot x + b^{\ast}) $$
3.2 線性可分支援向量機損失函式：
首先，在定義損失函式之前，需要了解函式間隔以及幾何間隔：
超平面$(w,b)$關於樣本點$(x_{i}, y_{i})$的函式間隔（使用於感知機損失函式）：
$$\hat{\gamma_{i}} = y_{i}(w \cdot x_{i} + b) $$
超平面$(w,b)$關於訓練資料集$T$的函式間隔：
$$ \hat{\gamma} =\min_{i=1,...,N} \hat{\gamma_{i}} $$
超平面$(w,b)$關於樣本點$(x_{i}, y_{i})$的幾何間隔：
$$\gamma_{i} = \frac{y_{i} (w \cdot x_{i})}{\left \| w \right \|} $$
超平面$(w,b)$關於訓練資料集$T$的幾何間隔：
$$ \gamma = \min_{i=1,...,N}\gamma_{i} $$
回想一下，支援向量機的學習策略就是間隔最大化，即最大化超平面$(w,b)$關於訓練資料集的幾何間隔$\gamma$，可以將其表示為如下的最優化問題：
$$ \max_{w,b} \gamma \quad s.t.\frac{y_{i} (w \cdot x_{i})}{\left \| w \right \|} \geq \gamma , \quadi=1,2,...,N $$
通過函式間隔與幾何間隔的關係，將其改寫成如下形式：
$$ \max_{w,b}\frac{\hat{\gamma}}{\left \| w \right \|} \quad s.t. y_{i}(w \cdot x_{i} + b) \geq \hat{\gamma}, \quad i=1,2,...,N$$
顯而易見，函式間隔$\hat{\gamma}$的取值不會影響上述優化問題的解，我們可以令其為1。同時將最大化$\frac{1}{||w||}$問題轉化為最小化$||w||^{2}$問題，同時為了便於後續的計算，加上$\frac{1}{2}$的係數，於是優化問題轉變為：
$$ \min_{w,b} \frac{1}{2}||w||^{2} \quad s.t. y_{i}(w \cdot x_{i} + b) \geq 1, \quad i=1,2,...,N$$
為了求解上述不等式優化問題，可以應用拉格朗日乘子法，得到拉格朗日函式$L$，L即線性可分支援向量機的損失函式：
$$ L(w, b, \alpha) = \frac{1}{2}\left \| w \right \|^2 + \sum_{i=1}^{N}\alpha_{i}\left [ 1 - y_{i}(w \cdot x + b) \right ] $$
通過拉格朗日函式可以得到：
$$\min_{w,b}\frac{1}{2}\left \| w \right \|^2 = \min_{w,b} \max_{\alpha_{i} \geq 0} L(w,b,\alpha) $$
然後，將其轉化成對偶問題：
$$\min_{w,b}\frac{1}{2}\left \| w \right \|^2 = \max_{\alpha_{i} \geq 0}\min_{w,b} L(w,b,\alpha) $$
3.3 線性可分支援向量機的學習過程：
線性可分支援向量機的學習過程即求解對偶問題的過程：首先求$L(w,b,\alpha)$對$w,b$的極小，再求其對$\alpha$的極大。
第一步：對求$L(w,b,\alpha)$對$w,b$的極小：
$$\nabla_{w}L(w,b,\alpha) = w - \sum_{i=1}^{N} \alpha_{i}y_{i}x_{i} = 0 $$
$$ \nabla_{b}L(w,b,\alpha) = \sum_{i=1}^{N} \alpha_{i}y_{i} = 0 $$
可以得到如下等式：
$$ w = \sum_{i=1}^{N} \alpha_{i}y_{i}x_{i} $$
$$ \sum_{i=1}^{N} \alpha_{i}y_{i} = 0 $$
其中，第一個式子將帶回原式$L(w,b,\alpha)$，第二個式子將作為第二步求解對$\alpha$的極大時的約束條件。
將第一個式子帶回原式$L(w,b,\alpha)$可得：
$$
\begin{align*}
L(w,b,\alpha) &= \frac{1}{2}\left \| w \right \|^2 + \sum_{i=1}^{N}\alpha_{i}(1 - y_{i}w \cdot x_{i}) - y_{i}b \\
&= \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}\left ( x_{i} \cdot x_{j} \right ) + \sum_{i=1}^{N} \alpha_{i} - \sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}\left ( x_{i} \cdot x_{j} \right ) + 0 \\
&= - \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}\left ( x_{i} \cdot x_{j} \right ) + \sum_{i=1}^{N} \alpha_{i}
\end{align*}
$$
第二步：求$\min_{w,b} L(w,b,\alpha)$對$\alpha$的極大，配合約束條件$\sum_{i=1}^{N} \alpha_{i}y_{i} = 0$ 可得：
$$\max_{\alpha} - \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}\left ( x_{i} \cdot x_{j} \right ) + \sum_{i=1}^{N} \alpha_{i} \quad s.t. \quad \sum_{i=1}^{N}\alpha_{i}y_{i}=0, \quad \alpha_{i} \geq0, \quad i=1,2,...,N $$
$$ \Downarrow $$
$$\min_{\alpha} \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}\left ( x_{i} \cdot x_{j} \right ) - \sum_{i=1}^{N} \alpha_{i} \quad s.t. \quad \sum_{i=1}^{N}\alpha_{i}y_{i}=0, \quad \alpha_{i} \geq0, \quad i=1,2,...,N $$
以上問題可以通過SMO演算法(這裡不詳細解釋)求得$\alpha$的解$\alpha^{\ast}$
最終求得$w^{\ast}$ 和 $ b^{\ast} $ 為：
$$ w^{\ast} = \sum_{i=1}^{N} \alpha_{i}^{\ast}y_{i}x_{i} $$
$$ b_{s}^{\ast} = y_{s} - \sum_{i=1}^{N} \alpha_{i}^{\ast}y_{i}x_{i} \cdot x_{s}， \quad \text{ for $\alpha_{s} > 0$， $(x_{s}, y_{s})$是支援向量 } $$
3.4 線性支援向量機：
3.3中介紹的線性可分支援向量機僅適用於線性可分訓練資料集，當訓練資料集線性不可分時，約束條件 $y_{i}(w \cdot x_{i} + b) \geq \hat{\gamma}, \quad i=1,2,...,N$ 並不對所有資料點都成立。。我們將原來的約束條件稱之為硬間隔最大化，為了將支援向量機擴充套件到線性不可分資料集，我們對其進行一定的修改，修改後的約束條件稱之為軟間隔最大化。
既然有些樣本點$(x_{i}, y_{i})$不能滿足函式間隔大於等於1的約束條件，我們對每個樣本點$(x_{i}, y_{i})$引入一個鬆弛變數$\xi_{i} \geq 0$，使函式間隔加上鬆弛變數大於等於1。同時，由於加入了鬆弛變數，需要付出代價，在目標函式中加對每個鬆弛變數加入代價。
$$ \begin{align*}
& \min_{w,b,\xi} \frac{1}{2}\left \| w \right \|^{2} + C\sum_{i=1}^{N}\xi_{i} \\
\text{s.t.} \quad & y_{i}(w \cdot x_{i} + b) \geq 1 - \xi_{i}, \\
& \xi_{i} \geq 0, \quad i=1,2,...,N
\end{align*} $$
依照3.3中相似的步驟，可以得到拉格朗日函式$L$：
$$L(w,b,\xi,\alpha,\mu) = \frac{1}{2}\left \| w \right \|^2 + C\sum_{i=1}^{N}\xi_{i} + \sum_{i=1}^{N}\alpha_{i}(1-\xi_{i}-y_{i}w \cdot x_{i} - y_{i}b) - \sum_{i=1}^{N} \mu_{i}\xi_{i} \quad s.t. \quad \alpha_{i} \geq 0, \quad \mu_{i} \geq 0 $$
接著求解拉格朗日對偶問題 $ \max_{\alpha, \mu} \min_{w,b,\xi} L(w,b,\xi, \alpha, \mu) $ 可得：
$$ \begin{align*}
&\min_{\alpha} \quad \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}(x_{i}\cdot x_{j}) - \sum_{i=1}^{N}\alpha_{i} \\
s.t. \quad & \sum_{i=1}^{N}\alpha_{i}y_{i}=0, \\
& 0 \leq \alpha_{i} \leq C, \quad i=1,2,...,N
\end{align*} $$
顯而易見，其形式和線性可分支援向量機一致，只是約束條件變了， 同樣可以利用SMO演算法求解$\alpha^{\ast}$，進而解得$w^{\ast}$和$b_{s}^{\ast}$。
3.5 線性支援向量機與合頁損失函式：
如3.4節所述：線性支援向量機，其模型為分離超平面$w^{\ast} \cdot x + b^{\ast} = 0$ 以及決策函式 $f(x) = \text{sign}(w^{\ast} \cdot x + b^{\ast})$，其學習策略為軟間隔最大化。但其還有另一種解釋，就是帶有正則化項的合頁損失函式最小化：
$$ \min_{w,b} \quad \sum_{i=1}^{N}\left [ 1-y_{i}(w \cdot x_{i} +b) \right ]_{+} + \lambda \left \| w \right \|^2$$
上述式子中第一項$ [1-y(w \cdot x + b)]_{+} $為合頁損失函式(hinge loss)，第二項為引數$w$的L2正則化項。
令 $1-y_{i}(w \cdot x_{i} + b) = \xi_{i}， \quad \xi_{i} \geq 0$，且 $ \lambda = \frac{1}{2C}$，則：
$$ \min \quad \sum_{i=1}^{N} [1-y_{i}(w \cdot x_{i} +b)]_{+} + \lambda||w||^2 $$
$$ \Downarrow $$
$$ \min \quad \sum_{i=1}^{N} \xi_{i} + \lambda||w||^2 $$
$$ \Downarrow $$
$$ \min \quad \frac{1}{C}(\frac{1}{2}||w||^2 + C\sum_{i=1}^{N}\xi_{i} ) $$
可以看到其最小化的函式和3.4節中的一致。
3.6 非線性支援向量機與核函式：
非線性分類問題，是指通過利用非線性模型才能很好地進行分類的問題。求解非線性問題的一個方法是進行非線性變換，將非線性問題變換為線性問題，通過求解變換後的線性問題的方法求解原來的非線性問題。（類似於將二元的多項式迴歸轉換成五元的線性迴歸）。總結一下，使用線性分類方法求解非線性分類問題分為兩步：
第一步：使用一個變換將原空間的資料對映到新空間
第二步：在新空間裡用線性分類學習方法從訓練資料中學習分類模型
按照上述理論，我們定義一個從輸入空間 $\mathcal{X}$ 到特徵空間 $\mathcal{H}$ 的對映：$\phi(x): \mathcal{X} \rightarrow \mathcal{H}$，將這一對映應用到線性支援向量機的優化函式中，得到：
$$ \begin{align*}
&\min_{\alpha} \quad \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}\phi(x_{i})\cdot \phi(x_{j}) - \sum_{i=1}^{N}\alpha_{i} \\
s.t. \quad & \sum_{i=1}^{N}\alpha_{i}y_{i}=0, \\
& 0 \leq \alpha_{i} \leq C, \quad i=1,2,...,N
\end{align*} $$
但這樣做還存在一個問題，那就是例項的內積 $\phi(x_{i}) \cdot \phi(x_{j})$ 是在對映後的空間$\mathcal{H}$中進行的，這是一個高維空間，計算內積並不容易且費時，於是，我們使用了核技巧，定義了核函式：
$$ K(x,z) = \phi(x) \cdot \phi(z) $$
核技巧的想法是，在學習與預測中只定義核函式$K(x,z)$，而不顯示地定義對映函式$\phi$。通常，計算$K(x,z)$比較容易，因為這是在低維空間上進行的。於是，應用了核函式之後的支援向量機的優化函式為：
$$ \begin{align*}
&\min_{\alpha} \quad \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}K(x_{i},x_{j})- \sum_{i=1}^{N}\alpha_{i} \\
s.t. \quad & \sum_{i=1}^{N}\alpha_{i}y_{i}=0, \\
& 0 \leq \alpha_{i} \leq C, \quad i=1,2,...,N
\end{align*} $$
定義一個核函式，需要各種滿足複雜的數學條件，一般我們使用的常見核函式包括：線性核函式(線性支援向量機中所使用的)、高斯核函式(也稱作徑向基核函式)，其他還包括多項式核函式、Sigmoid核函式：
線性核函式：$ K(x,z) = x \cdot z $
高斯核函式：$ K(x,z) = \exp( -\gamma \| x-z \|^2 ) $
多項式核函式：$ K(x,z) = (\gamma x \cdot z + r)^d $
Sigmoid核函式：$ K(x,z) = \tanh( \gamma x \cdot z + r ) $
以上核函式中的引數 $\gamma$, $r$, $d$ 都需要根據實際情況，調參選取最優值。
4. 實戰：
sklearn中，其實有兩個模組可以應用支援向量機，分別是 sklean.SVM，其中包括了線性 LinearSVC 以及 帶核函式的 SVC；另一個是 sklearn.linear_model.SGDClassifier，只需要將loss設定為 ‘hinge’ 就等同於線性支援向量機。
以下程式碼中，我們分別使用了 LinearSVC()，SGDClassifier，以及 SVC(kernel='rbf') 進行模型訓練，供參考：
1 import pandas as pd 2 import numpy as np 3 import matplotlib.pyplot as plt 4 from sklearn.preprocessing import MinMaxScaler, StandardScaler, OneHotEncoder, OrdinalEncoder 5 from sklearn.impute import SimpleImputer 6 from sklearn.model_selection import StratifiedKFold, GridSearchCV 7 from sklearn.pipeline import Pipeline, FeatureUnion 8 from sklearn.linear_model import SGDClassifier 9 from sklearn.svm import LinearSVC, SVC 10 from sklearn.metrics import accuracy_score, precision_score, recall_score 11 from sklearn.base import BaseEstimator, TransformerMixin 12 13 14 class DataFrameSelector(BaseEstimator, TransformerMixin): 15 def __init__(self, attribute_name): 16 self.attribute_name = attribute_name 17 18 def fit(self, x, y=None): 19 return self 20 21 def transform(self, x): 22 return x[self.attribute_name].values 23 24 25 # Load data 26 data_train = pd.read_csv('train.csv') 27 28 train_x = data_train.drop('Survived', axis=1) 29 train_y = data_train['Survived'] 30 31 # Data cleaning 32 cat_attribs = ['Pclass', 'Sex', 'Embarked'] 33 dis_attribs = ['SibSp', 'Parch'] 34 con_attribs = ['Age', 'Fare'] 35 36 # encoder: OneHotEncoder()、OrdinalEncoder() 37 cat_pipeline = Pipeline([ 38 ('selector', DataFrameSelector(cat_attribs)), 39 ('imputer', SimpleImputer(strategy='most_frequent')), 40 ('encoder', OneHotEncoder()), 41 ]) 42 43 dis_pipeline = Pipeline([ 44 ('selector', DataFrameSelector(dis_attribs)), 45 ('scaler', StandardScaler()), 46 ('imputer', SimpleImputer(strategy='most_frequent')), 47 ]) 48 49 con_pipeline = Pipeline([ 50 ('selector', DataFrameSelector(con_attribs)), 51 ('scaler', StandardScaler()), 52 ('imputer', SimpleImputer(strategy='mean')), 53 ]) 54 55 full_pipeline = FeatureUnion( 56 transformer_list=[ 57 ('con_pipeline', con_pipeline), 58 ('dis_pipeline', dis_pipeline), 59 ('cat_pipeline', cat_pipeline), 60 ] 61 ) 62 63 train_x_cleaned = full_pipeline.fit_transform(train_x) 64 65 cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=2) 66 67 # test1：using Linear SVC 68 clf1 = LinearSVC(tol=1e-4, max_iter=10000) 69 param_grid = [{ 70 'penalty': ['l1'], 71 'loss':['squared_hinge'], 72 'dual':[False], 73 'C':[1e-4, 1e-3, 1e-2, 1e-1, 1, 10], 74 'class_weight':['balanced', None] 75 }, 76 { 77 'penalty': ['l2'], 78 'loss':['hinge', 'squared_hinge'], 79 'dual':[True], 80 'C':[1e-4, 1e-3, 1e-2, 1e-1, 1, 10], 81 'class_weight':['balanced', None] 82 }] 83 84 grid_search = GridSearchCV(clf1, param_grid=param_grid, cv=cv, scoring='accuracy', n_jobs=-1, return_train_score=True) 85 86 grid_search.fit(train_x_cleaned, train_y) 87 predicted_y = grid_search.predict(train_x_cleaned) 88 89 df_cv_results = pd.DataFrame(grid_search.cv_results_) 90 print('-------Result of LinearSVC-------') 91 print(grid_search.best_params_) 92 print(accuracy_score(train_y, predicted_y)) 93 print(precision_score(train_y, predicted_y)) 94 print(recall_score(train_y, predicted_y)) 95 96 # test2: Using SGDClassifier with hinge loss / squared hinge loss 97 clf2 = SGDClassifier(tol=1e-4, max_iter=10000) 98 param_grid2 = [{ 99 'loss': ['hinge', 'squared_hinge'], 100 'penalty': ['l2', 'l1'], 101 'alpha': [1e-4, 1e-3, 1e-2, 1e-1, 1, 10, 100, 1000], 102 'class_weight':[None, 'balanced'] 103 }] 104 105 grid_search2 = GridSearchCV(clf2, param_grid=param_grid2, cv=cv, scoring='accuracy', n_jobs=-1, return_train_score=True) 106 107 grid_search2.fit(train_x_cleaned, train_y) 108 predicted_y2 = grid_search2.predict(train_x_cleaned) 109 110 df_cv_results2 = pd.DataFrame(grid_search2.cv_results_) 111 print('-------Result of SGD-------') 112 print(grid_search2.best_params_) 113 print(accuracy_score(train_y, predicted_y2)) 114 print(precision_score(train_y, predicted_y2)) 115 print(recall_score(train_y, predicted_y2)) 116 117 # test3: Using SVC 118 clf3 = SVC(kernel='rbf', tol=1e-4) 119 param_grid3 = [{ 120 'C': [0.1, 0.5, 0.75, 1, 2.5, 5, 10], 121 'gamma': ['scale', 'auto', 2e-4, 2e-3, 2e-2] 122 }] 123 124 grid_search3 = GridSearchCV(clf3, param_grid=param_grid3, cv=cv, scoring='accuracy', n_jobs=-1, return_train_score=True) 125 126 grid_search3.fit(train_x_cleaned, train_y) 127 predicted_y3 = grid_search3.predict(train_x_cleaned) 128 129 df_cv_results3 = pd.DataFrame(grid_search3.cv_results_) 130 print('-------Result of SVC with rbf-------') 131 print(grid_search3.best_params_) 132 print(accuracy_score(train_y, predicted_y3)) 133 print(precision_score(train_y, predicted_y3)) 134 print(recall_score(train_y, predicted_y3)) 135 136 # 匯入預測資料，預測結果，並生成csv檔案 137 data_test = pd.read_csv('test.csv') 138 submission = pd.DataFrame(columns=['PassengerId', 'Survived']) 139 submission['PassengerId'] = data_test['PassengerId'] 140 141 test_x_cleaned = full_pipeline.fit_transform(data_test) 142 143 submission_LinearSVC = pd.DataFrame(submission, copy=True) 144 submission_LinearSVC['Survived'] = pd.Series(grid_search.predict(test_x_cleaned)) 145 146 submission_SVCrbf = pd.DataFrame(submission, copy=True) 147 submission_SVCrbf['Survived'] = pd.Series(grid_search3.predict(test_x_cleaned)) 148 149 submission_LinearSVC.to_csv('submission_LinearSVC.csv', index=False) 150 submission_SVCrbf.to_csv('submission_SVCrbf.csv', index=False)
4.1 結果分析：
和前幾篇一樣，分別將 LinearSVC() 和 SVC(kernel='rbf')的最優引數用於預測集，並將結果上傳kaggle，結果如下：
可以看到線性SVM和其他線性方法的結果差不多；而使用了核函式後，結果確實有了一定的提升。
|
訓練集 accuracy |
訓練集 precision |
訓練集 recall |
預測集 accuracy（需上傳kaggle獲取結果） |
|
| 樸素貝葉斯最優解 | 0.790 | 0.731 | 0.716 | 0.756 |
| 感知機 | 0.771 | 0.694 | 0.722 | 0.722 |
| 邏輯迴歸 | 0.807 | 0.781 | 0.690 | 0.768 |
| 線性SVM | 0.801 | 0.772 | 0.684 | 0.773 |
| rbf核SVM | 0.834 | 0.817 | 0.731 | 0.785 |