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

1次阅读

没有评论

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

支持向量机（SVM）最初的形式是一个带约束的优化问题，称为原始问题（Primal Problem）。它试图找到一个超平面，使得分类间隔最大化。数学上，原始问题可以表示为：

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

为了解决这个问题，我们引入拉格朗日乘子，将其转换为对偶问题（Dual Problem）。对偶问题的关键优势在于：

将原始的高维参数优化转换为对偶变量（拉格朗日乘子）的优化
自然地引入了核函数，使 SVM 能够处理非线性分类问题
减少了优化变量的数量，特别是在高维特征空间的情况下

对偶问题的最终形式为：

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

其中 $K(x_i,x_j)$ 是核函数，$\alpha_i$ 是拉格朗日乘子，需要满足 KKT 条件：

$\alpha_i \geq 0$
$y_i(w^T x_i + b) – 1 \geq 0$
$\alpha_i[y_i(w^T x_i + b) – 1] = 0$

虽然对偶问题有诸多优势，但在实际求解时仍会面临几个关键挑战：

计算复杂度：当样本量很大时，核矩阵的大小为 $n \times n$，存储和计算都变得非常昂贵
数值稳定性：某些核函数可能导致矩阵条件数很差，影响求解精度
收敛速度：特别是当数据具有特殊结构时，某些算法可能收敛很慢
内存瓶颈：大规模数据集下，存储整个核矩阵可能超出内存容量

最直接的方法是使用通用 QP 求解器（如 CVXOPT）：

优点：实现简单，可靠性高
缺点：难以扩展到大规模问题，核矩阵存储是瓶颈

SMO 是专门为 SVM 设计的算法：

每次只优化一对 $\alpha$ 变量
解析求解这两个变量的最优值
优点：内存效率高，实现相对简单
缺点：收敛速度可能较慢

内点法是一类有效的凸优化算法：

优点：理论收敛性好
缺点：每步迭代计算量大，难以处理大规模问题

以下是简化版 SMO 算法的 Python 实现：

import numpy as np

class SVM:
    def __init__(self, C=1.0, kernel='linear', tol=1e-3, max_iter=1000):
        self.C = C  # 正则化参数
        self.kernel = kernel
        self.tol = tol  # 容忍度
        self.max_iter = max_iter  # 最大迭代次数

    def fit(self, X, y):
        n_samples, n_features = X.shape
        self.alpha = np.zeros(n_samples)
        self.b = 0

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

        # SMO 主循环
        iter = 0
        while iter < self.max_iter:
            num_changed_alphas = 0
            for i in range(n_samples):
                # 计算误差
                E_i = self._decision_function(X[i]) - y[i]

                # 检查是否违反 KKT 条件
                if ((y[i]*E_i < -self.tol) and (self.alpha[i] < self.C)) or \
                   ((y[i]*E_i > self.tol) and (self.alpha[i] > 0)):

                    # 随机选择第二个 alpha
                    j = np.random.choice(list(range(i)) + list(range(i+1, n_samples)))
                    E_j = self._decision_function(X[j]) - y[j]

                    # 保存旧值
                    alpha_i_old = self.alpha[i]
                    alpha_j_old = self.alpha[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 * self.K[i,j] - self.K[i,i] - self.K[j,j]
                    if eta >= 0:
                        continue

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

                    # 裁剪到边界
                    self.alpha[j] = np.clip(self.alpha[j], L, H)

                    # 检查变化是否显著
                    if abs(self.alpha[j] - alpha_j_old) < 1e-5:
                        continue

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

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

                    if 0 < self.alpha[i] < self.C:
                        self.b = b1
                    elif 0 < self.alpha[j] < self.C:
                        self.b = b2
                    else:
                        self.b = (b1 + b2) / 2

                    num_changed_alphas += 1

            if num_changed_alphas == 0:
                iter += 1
            else:
                iter = 0

        # 保存支持向量
        self.support_vectors = X[self.alpha > 0]
        self.support_vector_labels = y[self.alpha > 0]
        self.support_vector_alphas = self.alpha[self.alpha > 0]

    def _kernel(self, x1, x2):
        if self.kernel == 'linear':
            return np.dot(x1, x2)
        elif self.kernel == 'rbf':
            gamma = 1.0 / x1.shape[0]
            return np.exp(-gamma * np.linalg.norm(x1 - x2)**2)
        else:
            raise ValueError("Unknown kernel")

    def _decision_function(self, x):
        return np.dot(self.alpha * y, self.K[:, np.where((X == x).all(axis=1))[0][0]]) + self.b

    def predict(self, X):
        return np.sign([self._decision_function(x) for x in X])