機器學習：支持向量機

支持向量機（Support Vector Machine）是一種二類分類模型，其基本模型定義為特征空間上的間隔最大的廣義線性分類器，其學習策略便是間隔最大化，最終可轉化為一個凸二次規劃問題的求解。

在這里插入圖片描述

假設兩類數據可以被 ${x:w^Tx + b \ge c}$ 分離，垂直于法向量 $w$ ，移動 $H$ 直到碰到某個訓練點，可以得到兩個超平面 $H_1$ 和 $H_2$ ，兩個平面稱為支撐超平面，題目分別支撐兩類數據。而位于 $H_1$ 和 $H_2$ 正中間的超平面是分離這兩類數據的最好選擇。支持向量就是離分隔超平面最近的那些點。

法向量 $w$ 有很多種選擇，超平面 $H_1$ 和 $H_2$ 之間的距離稱為間隔，這個間隔是 $w$ 的函數，**目的就是尋找這樣的 $w$ 使得間隔達到最大。

在這里插入圖片描述

在求解最優化問題中，拉格朗日乘子法（Lagrange Multiplier）和KKT（Karush Kuhn Tucker）條件是兩種最常用的方法。在有等式約束時使用拉格朗日乘子法，在有不等約束時使用KKT條件。

拉格朗日乘子法

拉格朗日乘子法是一種尋找多元函數在一組約束下的極值的方法。通過引入拉格朗日乘子，可將有 $d$ 個變量與 $k$ 個約束條件的最優化問題轉化為具有 $d + k$ 個變量的無約束優化問題求解。
二次規劃

二次規劃是一類典型的優化問題，包括凸二次優化和非凸二次優化。在此類問題中，目標函數是變量的二次函數，而約束條件是變量的線性不等式。
$\frac{1} {2} x^T Q x + c^T x \\ s.t. \vec{A} \vec{x} \le \vec{b}$

具體公式證明：【整理】深入理解拉格朗日乘子法（Lagrange Multiplier) 和KKT條件 - mo_wang - 博客園 (cnblogs.com)

序列最小優化（Sequential Minimal Optimization，SMO）

序列最小優化是將大優化問題分界成多個小優化問題來求解。

SMO算法工作原理：每次循環中選擇兩個變量進行優化處理。一旦找到一對合適的變量，那么就增大其中一個同時減小另一個。這里的“合適”指的是兩個變量必須要符合一定的條件，條件之一就是這兩個變量必須要在間隔邊界之外，而其第二個條件則是這兩個變量還沒有進行過區間化處理或者不在邊界上。

在這里插入圖片描述

代碼實現

參考《機器學習實戰》，代碼鏈接：https://github.com/golitter/Decoding-ML-Top10/tree/master/SVM

這里采用簡化的SMO代碼，數據集是https://blog.caiyongji.com/assets/mouse_viral_study.csv。

data_processing.py：

import numpy as np
import pandas as pd# https://zhuanlan.zhihu.com/p/350836534
def data_processing():data_csv = pd.read_csv('mouse_viral_study.csv')data_csv = data_csv.dropna()# print(data_csv)X = data_csv.iloc[:-1, 0:2].values# print(X)Y = data_csv.iloc[:-1, 2].map({0: -1, 1: 1}).valuesY = Y.reshape(-1, 1)# print(Y.shape)return X, Y# X, Y = data_processing()
# print(X)

工具模塊，smo_assist.py：

import random
def select_Jrandom(i:int, m:int) -> int:"""隨機選擇一個不等于 i 的整數"""j = iwhile j == i:j = int(random.uniform(0, m))return jdef clip_alpha(alpha_j:float, H:float, L:float) -> float:"""修剪 alpha_j"""if alpha_j > H:alpha_j = Hif alpha_j < L:alpha_j = Lreturn alpha_j

簡化SMO的代碼實現，smoSimple.py：

from smo_assist import (select_Jrandom, clip_alpha)import numpy as np
import pdbdef smoSimple(data_mat_in:np.ndarray, class_labels:np.ndarray, C:float, toler:float, max_iter:int):"""data_mat_in: 數據集class_labels: 類別標簽C: 松弛變量toler: 容錯率max_iter: 最大迭代次數"""b = 0; # 初始化bm, n = np.shape(data_mat_in) # m: 樣本數, n: 特征數alphas = np.zeros((m, 1)) # 初始化alphaiter = 0 # 迭代次數while iter < max_iter:alphaPairsChanged = 0for i in range(m):fXi = float(np.multiply(alphas, class_labels).T @ (data_mat_in @ data_mat_in[i, :].T)) + b"""(1 , m) * (m, n) * (n, 1) = (1, 1) = 標量再 加上 b 就是 f(x) 的值"""Ei = fXi - float(class_labels[i])"""Ei = f(x) - y 預測誤差"""if (# 第一種情況：樣本被誤分類，且權重可以增加((class_labels[i] * Ei < -toler) # 預測誤差與標簽方向相反，且誤差大于容忍度and (alphas[i] < C)) # 當前權重小于正則化參數 C，可以增加權重or # 第二種情況：樣本被誤分類，且權重需要調整((class_labels[i] * Ei > toler) # 預測誤差與標簽方向相同，且誤差大于容忍度and (alphas[i] > 0)) # 當前權重大于 0，需要調整權重):j = select_Jrandom(i, m)fxj = float(np.multiply(alphas, class_labels).T @ (data_mat_in @ data_mat_in[j, :].T)) + bEj = fxj - float(class_labels[j])alpha_j_old = alphas[j].copy(); alpha_i_old = alphas[i].copy()if (class_labels[i] != class_labels[j]):L = max(0, alphas[j] - alphas[i]) # 左邊界H = min(C, C + alphas[j] - alphas[i]) # 右邊界else:L = max(0, alphas[j] + alphas[i] - C)H = min(C, alphas[j] + alphas[i])if L == H: continue # 跳出本次循環eta = 2.0 * data_mat_in[i, :] @ data_mat_in[j, :].T - data_mat_in[i, :] @ data_mat_in[i, :].T - data_mat_in[j, :] @ data_mat_in[j, :].T"""計算 eta = K11 + K22 - 2 * K12 = 2 * x_i * x_j - x_i * x_i - x_j * x_j """     if eta >= 0:continuealphas[j] -= class_labels[j] * (Ei - Ej) / eta # 更新權重alphas[j] = clip_alpha(alphas[j], H, L) # 調整權重if abs(alphas[j] - alpha_j_old) < 0.00001:continue # 跳出本次循環，不更新 ialphas[i] += class_labels[j] * class_labels[i] * (alpha_j_old - alphas[j]) # 更新權重b1 = b - Ei - class_labels[i] * (alphas[i] - alpha_i_old) * data_mat_in[i, :] @ data_mat_in[i, :].T - class_labels[j] *(alphas[j] - alpha_j_old) * data_mat_in[i, :] @ data_mat_in[j, :].Tb2 = b - Ej - class_labels[i] * (alphas[i] - alpha_i_old) * data_mat_in[i, :] @ data_mat_in[j, :].T - class_labels[j] *(alphas[j] - alpha_j_old) * data_mat_in[j, :] @ data_mat_in[j, :].T"""更新 b"""     if 0 < alphas[i] < C:b = b1elif 0 < alphas[j] < C:b = b2else:b = (b1 + b2) / 2.0alphaPairsChanged += 1if alphaPairsChanged == 0:iter += 1else:iter = 0return b, alphasif __name__ == '__main__':print(  smoSimple(np.array([[1, 2], [3, 4]]), np.array([[-1],[1]]), 0.6, 0.001, 40))

test.py：

from data_processing import *
from smoSimple import *
import numpy as np
import matplotlib.pyplot as plt# 數據處理和 SVM 訓練
data_mat_in, class_labels = data_processing()
b, alphas = smoSimple(data_mat_in, class_labels, 0.6, 0.001, 40)# 打印結果
print("Bias (b):", b)
print("Non-zero alphas:", alphas[alphas > 0])# 打印數據形狀
print("Shape of data_mat_in:", np.shape(data_mat_in))
print("Shape of class_labels:", np.shape(class_labels))# 將 Y 轉換為一維數組（如果它是二維的）
Y = class_labels
# 提取不同類別的索引
class_1_indices = np.where(Y == 1)[0]  # 類別為 1 的樣本索引
class_2_indices = np.where(Y == -1)[0]  # 類別為 -1 的樣本索引
X = data_mat_in# 繪制散點圖
plt.figure(figsize=(8, 6))
plt.scatter(X[class_1_indices, 0], X[class_1_indices, 1], c='blue', label='Class 1', alpha=0.5)
plt.scatter(X[class_2_indices, 0], X[class_2_indices, 1], c='red', label='Class -1', alpha=0.5)# 計算權重向量 w
w = np.dot((alphas * Y).T, X).flatten()
# print(f"w: {w}")
print("Shape of X:", X.shape)  # 應該是 (m, n)
print("Shape of Y:", Y.shape)  # 應該是 (m, 1)
print("Shape of alphas:", alphas.shape)  # 應該是 (m, 1)# 繪制超平面
# 超平面方程：w[0] * x1 + w[1] * x2 + b = 0
# 解出 x2: x2 = -(w[0] * x1 + b) / w[1]
x1 = np.linspace(np.min(X[:, 0]), np.max(X[:, 0]), 100)
x2 = -(w[0] * x1 + b) / w[1]
print(f"w_shape: {w.shape}")
# 繪制超平面
plt.plot(x1, x2, label='SVM Hyperplane', color='green', linewidth=2)# 標出支持向量
support_vectors_indices = np.where(alphas > 0)[0]  # 找到所有支持向量的索引
plt.scatter(X[support_vectors_indices, 0], X[support_vectors_indices, 1], facecolors='none', edgecolors='k', s=50, label='Support Vectors')# 添加圖例和標簽
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('Scatter Plot of Data with SVM Hyperplane')
plt.legend()# 顯示圖形
plt.show()