支持向量机对偶问题求解：从数学原理到Python实现

1次阅读

没有评论

共计 4545 个字符，预计需要花费 12 分钟才能阅读完成。

支持向量机 (SVM) 最初的形式是一个凸二次规划问题，其原始优化目标为：

$$
\min_{w,b} \frac{1}{2}||w||^2 \quad \text{s.t.} \quad y_i(w^Tx_i + b) \geq 1, \forall i
$$

为了将这个带约束的优化问题转化为无约束问题，我们引入拉格朗日乘子法。构建拉格朗日函数：

$$
L(w,b,\alpha) = \frac{1}{2}||w||^2 – \sum_{i=1}^n \alpha_i[y_i(w^Tx_i + b) – 1]
$$

其中 $\alpha_i \geq 0$ 是拉格朗日乘子。根据 KKT 条件，最优解需要满足：

梯度为零：$\frac{\partial L}{\partial w} = 0$ 和 $\frac{\partial L}{\partial b} = 0$
原始可行性：$y_i(w^Tx_i + b) \geq 1$
对偶可行性：$\alpha_i \geq 0$
互补松弛性：$\alpha_i[y_i(w^Tx_i + b) – 1] = 0$

通过求解梯度为零的条件，我们得到：

$$
w = \sum_{i=1}^n \alpha_i y_i x_i
$$

$$
\sum_{i=1}^n \alpha_i y_i = 0
$$

将这两个关系式代回拉格朗日函数，就得到了对偶问题：

$$
\max_{\alpha} \sum_{i=1}^n \alpha_i – \frac{1}{2} \sum_{i,j=1}^n \alpha_i \alpha_j y_i y_j x_i^T x_j
$$

$$
\text{s.t.} \quad \sum_{i=1}^n \alpha_i y_i = 0, \alpha_i \geq 0
$$

原始问题：需要优化 $w \in \mathbb{R}^d$ 和 $b \in \mathbb{R}$，变量数量与特征维度 $d$ 相关
对偶问题：需要优化 $\alpha \in \mathbb{R}^n$，变量数量与样本数量 $n$ 相关

当特征维度 $d$ 远大于样本数量 $n$ 时，对偶问题更高效。此外，对偶形式天然支持核技巧，可以通过核函数 $K(x_i,x_j)$ 隐式映射到高维空间。

import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_classification
import matplotlib.pyplot as plt

# 生成数据
X, y = make_classification(n_samples=100, n_features=2, n_redundant=0, 
                          n_clusters_per_class=1, random_state=42)
y[y==0] = -1  # 将标签转换为 - 1 和 1

# 数据标准化
scaler = StandardScaler()
X = scaler.fit_transform(X)

# 核函数实现
def linear_kernel(x1, x2):
    return np.dot(x1, x2)

def rbf_kernel(x1, x2, gamma=0.1):
    return np.exp(-gamma * np.linalg.norm(x1 - x2)**2)

# SVM 对偶问题求解
class SVM:
    def __init__(self, kernel=linear_kernel, C=1.0):
        self.kernel = kernel
        self.C = C

    def fit(self, X, y, max_iter=1000, tol=1e-3):
        n_samples, n_features = X.shape
        self.alpha = np.zeros(n_samples)
        self.b = 0

        # 计算核矩阵
        K = np.zeros((n_samples, n_samples))
        for i in range(n_samples):
            for j in range(n_samples):
                K[i,j] = self.kernel(X[i], X[j])

        # 使用 SMO 算法求解
        for _ in range(max_iter):
            alpha_prev = np.copy(self.alpha)

            for i in range(n_samples):
                # 计算预测误差
                E_i = np.sum(self.alpha * y * K[:,i]) + self.b - y[i]

                # 选择第二个变量
                j = np.random.choice(list(range(i)) + list(range(i+1, n_samples)))
                E_j = np.sum(self.alpha * y * K[:,j]) + self.b - y[j]

                # 计算边界
                if y[i] != y[j]:
                    L = max(0, self.alpha[j] - self.alpha[i])
                    H = min(self.C, self.C + self.alpha[j] - self.alpha[i])
                else:
                    L = max(0, self.alpha[i] + self.alpha[j] - self.C)
                    H = min(self.C, self.alpha[i] + self.alpha[j])

                if L == H:
                    continue

                # 计算 eta
                eta = 2 * K[i,j] - K[i,i] - K[j,j]
                if eta >= 0:
                    continue

                # 更新 alpha_j
                alpha_j = self.alpha[j] - y[j] * (E_i - E_j) / eta
                alpha_j = max(L, min(H, alpha_j))

                # 更新 alpha_i
                alpha_i = self.alpha[i] + y[i]*y[j]*(self.alpha[j] - alpha_j)

                # 更新 b
                b1 = self.b - E_i - y[i]*(alpha_i - self.alpha[i])*K[i,i] - y[j]*(alpha_j - self.alpha[j])*K[i,j]
                b2 = self.b - E_j - y[i]*(alpha_i - self.alpha[i])*K[i,j] - y[j]*(alpha_j - self.alpha[j])*K[j,j]

                if 0 < alpha_i < self.C:
                    self.b = b1
                elif 0 < alpha_j < self.C:
                    self.b = b2
                else:
                    self.b = (b1 + b2)/2

                self.alpha[i] = alpha_i
                self.alpha[j] = alpha_j

            # 检查收敛
            if np.linalg.norm(self.alpha - alpha_prev) < tol:
                break

        # 获取支持向量
        idx = self.alpha > 1e-5
        self.support_vectors = X[idx]
        self.support_vector_labels = y[idx]
        self.support_vector_alphas = self.alpha[idx]

        # 计算 w (仅对线性核有效)
        if self.kernel == linear_kernel:
            self.w = np.sum((self.alpha * y).reshape(-1,1) * X, axis=0)
        else:
            self.w = None

    def predict(self, X):
        if self.w is not None:
            return np.sign(np.dot(X, self.w) + self.b)
        else:
            y_pred = np.zeros(len(X))
            for i in range(len(X)):
                s = 0
                for alpha, sv_y, sv in zip(self.support_vector_alphas, 
                                          self.support_vector_labels, 
                                          self.support_vectors):
                    s += alpha * sv_y * self.kernel(X[i], sv)
                y_pred[i] = s
            return np.sign(y_pred + self.b)

# 训练模型
svm = SVM(kernel=linear_kernel, C=1.0)
svm.fit(X, y)

# 可视化
def plot_decision_boundary(model, X, y):
    h = 0.02
    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h),
                         np.arange(y_min, y_max, h))

    Z = model.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)

    plt.contourf(xx, yy, Z, alpha=0.8)
    plt.scatter(X[:, 0], X[:, 1], c=y, edgecolors='k')
    plt.scatter(model.support_vectors[:, 0], model.support_vectors[:, 1], 
                s=100, facecolors='none', edgecolors='k')
    plt.title('SVM Decision Boundary')
    plt.show()

plot_decision_boundary(svm, X, y)

from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV

# 使用线性核
svm_linear = SVC(kernel='linear', C=1.0)
svm_linear.fit(X, y)

# 使用 RBF 核
svm_rbf = SVC(kernel='rbf', gamma=0.1, C=1.0)
svm_rbf.fit(X, y)

# 超参数调优
param_grid = {'C': [0.1, 1, 10, 100], 'gamma': [0.01, 0.1, 1, 10]}
grid = GridSearchCV(SVC(kernel='rbf'), param_grid, cv=5)
grid.fit(X, y)

print("Best parameters:", grid.best_params_)
print("Best score:", grid.best_score_)