SVM 学习

问题引入

假设现在有一个二分类任务，我们有一些标签数据 $ (X_i, y_i) $ ，其中 $i=1,2,…n$，$X_i$代表了数据特征，设其维度为 $d$, 即 $X_i \in \mathbf{R}^d$。 $y_i$是标签，由于是二分类任务，其取值为 $1$ 或 $-1$。二分类任务的目标是训练一个模型，可以模型准确的预测与训练数据满足同一分布的新数据。

Hard-margin SVM

由于数据都是 d-维空间上的点，Hard-margin SVM 首先假设这些点可以被一个超平面 hyperplane分割，将这样的超平面作为模型完成分类任务。超平面将空间分为两侧，每一侧分别作为二分类问题中的一类。

模型可以简单表述为，给定数据 $x \in \mathbf{R}^d$，计算 $f(x)$，查看 $f(x)$ 的符号来确定类别。

\[f(x) = \mathbf{w} \mathbf{x} + b \\ y = sign(f(x))\]

如下图所示，直线H1、H2、H3都是超平面，可以将空间分隔成两部分。但我们应该选取哪条直线作为模型呢，即如何选取最优分隔平面呢？根据直觉（也是SVM的思想），我们认为 训练样本中的点离超平面的距离的最小值越大越好。得到的超平面离边界点距离大，可以让模型遇到轻微扰动的点时候能够正确分类，模型抗干扰性更强。点到直线的公式:

\[distance = \frac{| \mathbf{w} \mathbf{x} + b|}{||w||_2}\]

在此基础上，我们加上一些归一化操作，认为离超平面最近的点满足

\[y (\mathbf{w} \mathbf{x} + b) = 1\]

综上，SVM问题变为了（称之为 primal hard margin SVM problem）

\[\begin{align} maximize& \frac{1}{||\mathbf{w}||} = minimize \frac{1}{2} ||\mathbf{w}||^2 \\ s.t.& \ y_i(\mathbf{w} x_i + b) \geq 1 \quad for \ i=1,2,...n \end{align}\]

求解 Hard-margin SVM 问题

前面我们给出了 Hard-margin SVM 问题描述和定义，我们现在想办法求解改问题。这个问题是凸二次规划问题，凸二次规划问题在可行域上具有全局最优解。有许多算法可以求解凸二次规划问题，如内点法（Interior Point Method）、椭圆法等等。标准库如LibSVM中一般采取 SMO（Sequential Minimal Optimization）算法来解。我了解到 CVXOPT工具包中的solver.qp就使用了内点法求解二次规划问题，下面展示了它的简单用法 from cvxopt docs。我们需要做一些简单的变化让前面提到的Hard-Margin Primal SVM 问题适应这个API的调用：对于偏执项 $b$，我们将其附加到 $\mathbf{w}$ 最后一维，即 $ \hat{\mathbf{w}} = (\mathbf{w} | b)$

调用的参数为：

\[P = \begin{bmatrix} I_{d} & 0\\ 0 & 0 \end{bmatrix}\] \[q = [0] * (d+1)\] \[G_{nd} = -y \odot \begin{bmatrix} \mathbf{x} | \mathbf{1} \end{bmatrix} \\ h = [-1] * n\] \[A = b = None \\\]

python 实现(使用 numpy 库):

  
class Primal_Hard_SVM():
    """Primal Hard-margin SVM problem solved by cvxopt.solvers.qp."""
    def __init__(self) -> None:
        # bias is the last term of w.
        self.w = None  # shape (d+1, ) 
    
    def fit(self, X, y):
        n, d = X.shape[0], X.shape[1]

        P =  np.zeros((d+1, d+1)) # shape (d+1, d+1)
        P[:d, :d] = np.eye(d)

        q = np.zeros((d+1, ))  # shape (d+1, )

        G = np.zeros((n, d + 1))  # shape (n, d+1)
        G[:, :d] = -y.reshape(-1, 1) * X # broadcast & pointwise product
        G[:, d] = -y 
        h = np.array([-1.0] * n)  # (n,)
        
        sol = solvers.qp(matrix(P), matrix(q), matrix(G), matrix(h))
        self.w = np.array(sol['x']).flatten()

        return self
    
    def predict(self, X):
        """Predict given data."""
        n = X.shape[0]
        y_values = np.concatenate((X, np.ones((n, 1))), axis=1) @ self.w  # (n, d+1) @ (d+1, ) = (n, )
        y_labels = np.sign(y_values)  # (n, )
        
        return y_values, y_labels

对偶 dual Hard-Margin SVM 问题

没有凸优化背景可能不了解对偶问题，先补充一些对偶问题的基本知识。

对偶问题的来源和定义

对偶问题最早的想法来源于数学家拉格朗日（Joseph-Louis Lagrange）在18世纪的研究工作。拉格朗日提出了一种在优化问题中引入拉格朗日乘子的方法，以解决约束条件下的极值问题。考虑一个 General 的凸问题的 primal 形式：

\[\begin{align} \text{minimize} & \quad f_0(\mathbf{x}) \\ \text{subject to} & \quad f_i(\mathbf{x}) \leq 0 \quad \forall i = 1, 2, \ldots, m \\ & \quad h_i(\mathbf{x}) = 0 \quad \forall i = 1, 2, \ldots, p \\ \end{align}\]

其中，$\mathbf{x}$是优化变量，$f_0(\mathbf{x})$是凸目标函数，$f_i(\mathbf{x})$是凸不等式约束函数，$h_i(\mathbf{x})$是等式约束函数。这是一个求解使目标函数最小化的优化问题。对偶问题是通过对原始问题的约束条件引入拉格朗日乘子（也称为对偶变量）来构建的。通过引入拉格朗日乘子$\mathbf{\lambda}$和$\mathbf{\nu}$，对偶问题可以表示为：

\[\begin{align*} \text{maximize} & \quad g(\mathbf{\lambda}, \mathbf{\nu}) \\ \text{subject to} & \quad \mathbf{\lambda} \geq 0 \end{align*}\]

其中，$g(\mathbf{\lambda}, \mathbf{\nu})$是拉格朗日对偶函数。定义为：

\[g(\mathbf{\lambda}, \mathbf{\nu}) = \inf_{\mathbf{x}} \left( f_0(\mathbf{x}) + \sum_{i=1}^{m} \lambda_i f_i(\mathbf{x}) + \sum_{i=1}^{p} \nu_i h_i(\mathbf{x}) \right)\]

可以证明（这里略去），拉格朗日对偶函数是原问题的下界。而 maximize 拉格朗日对偶函数的 idea 就是去找到一个最优下界。

强/弱对偶与 Slater 条件

我们把对偶问题的最优解记为 $d^*$ ，原问题的最优解记为 $p^$ , 对偶问题的最优解构成了原问题的下解（证明略），即 $d^* \leq p^*$ 我们称为弱对偶性，$p^ - d^*$ 称为对偶间隙。即使原问题不是凸问题，弱对偶性仍然成立。若 $d^* = p^*$ 我们称为强对偶性成立，即对偶间隙为0。显然，强对偶性质很好，满足强对偶性意味着对偶问题和原问题的最优解是相同的。但对一般情况（大多数问题不是凸问题）下，强对偶条件是不成立的。那满足什么条件可以让强对偶性成立，帮助我们构造对偶问题来求解原问题呢？

Slater 条件给出了强对偶条件的一个成立的充分条件（证明略），即Slater条件成立，强对偶一定成立，但强对偶成立，Slater条件不一定成立。Slater条件要求存在一个向量$x’$，满足

不等式约束：$f_i(x’) < 0, \quad i = 1, 2, …, m$
等式约束：$h_j(x’) = 0, \quad j = 1, 2, …, p$

换句话说，存在一个可行解$x’$，使得所有的不等式约束严格满足，而等式约束恰好满足。幸运的是，对于大部分（不是所有）凸问题，Slater 条件都成立，即强对偶性成立。 还有一个弱化版本的Slater 条件更有利于我们解决简单问题：

若不等式约束为仿射函数，仿射不等式不需要严格成立。只要可行域非空，强对偶性成立。

如何求解对偶问题 KKT 条件

我们已经知道了对偶问题的解对应于原问题的最优解。那如何求解对偶问题呢？我们引入 KKT 条件。 KKT条件（Karush-Kuhn-Tucker条件）是优化问题中判断最优解的一组必要条件，即如果某个解是最优解，那它一定满足KKT条件。但满足KKT条件的不一定是最优解。KKT条件包括以下几个方面：

对于最小化问题，KKT条件的第一部分是目标函数的梯度与约束函数的梯度线性相关：

\[\nabla f(x^*) + \sum_{i=1}^m \lambda_i^* \nabla g_i(x^*) + \sum_{j=1}^p \mu_j^* \nabla h_j(x^*) = 0\]

不等式约束的松弛条件（complementary slackness）：

\[\lambda_i^* g_i(x^*) = 0, \quad i = 1, 2, ..., m\]

其中，$\lambda_i^$是不等式约束对应的拉格朗日乘子，$g_i(x^)$是不等式约束函数在最优解$x^*$处的值。这个条件表示，在最优解中，不等式约束和对应的拉格朗日乘子之间有一个松弛条件。

非负约束条件：

\[\lambda_i^* \geq 0, \quad i = 1, 2, ..., m\]

这个条件要求不等式约束的拉格朗日乘子非负。

等式约束条件：

\[h_j(x^*) = 0, \quad j = 1, 2, ..., p\]

上面提到的优化问题不单指凸优化问题，如果某些问题满足一些额外的条件（如Slater条件），KKT条件也是最优解的充分条件，即充要条件(证明略)。 这意味着，我们求解一个满足Slater条件的问题的对偶问题时，只需要让KKT条件中的所有等式和不等式都成立。

求解对偶问题的优势

有了这么多背景知识，应该会有疑问——我们为什么要废这么大劲去解对偶问题，而不是直接接原问题呢？

对于强对偶性成立的问题，对偶问题往往复杂度会更低、更容易求解。时间 & 空间复杂度
求解对偶问题的优化算法收敛更快
将原问题的约束（可能是非线性的）转化成了对偶问题的线性约束可能更容易求解。
可以引入 Kernel 函数（对SVM来说）
…

Dual hard-margin SVM problem

根据上面的背景知识，我们知道 SVM 问题原问题是一个凸问题且满足weak Slater条件，所以我们只需要求解对偶问题的KKT条件，就能获得和原问题一样的最优解。我们首先写出 Hard-margin SVM 的拉格朗日函数。

拉格朗日函数的形式为：

\[\mathcal{L}(w, b, \alpha) = \frac{1}{2} |w|^2 - \sum_{i=1}^{n} \alpha_i [y_i (w^T x_i + b) - 1]\]

其中，$w$ 是超平面的法向量，$b$ 是超平面的截距，$\alpha_i$ 是拉格朗日乘子。

使用 KKT 条件得到：

目标函数的梯度与约束的梯度线性相关 $\nabla_w \mathcal{L}(w, b, \alpha) = 0 \\ \nabla_b \mathcal{L}(w, b, \alpha) = 0$ 这里的 $\mathcal{L}(w, b, \alpha)$ 是拉格朗日函数，$\nabla_w$ 和 $\nabla_b$ 分别表示对 $w$ 和 $b$ 求偏导数的操作。
互补松弛 $\alpha_i [y_i (w^T x_i + b) - 1] = 0, \quad i = 1, 2, ..., n$对于 $\alpha_i > 0$ 的训练样本，需要满足 $y_i (w^T x_i + b) = 1$
非负条件 $\alpha_i \geq 0, \quad i = 1, 2, ..., n$

化简后可以得到对应的对偶问题：

\[\begin{align} \text{minimize} & \quad -\sum_{i=1}^n \alpha_i + \frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n \alpha_i \alpha_j y_i y_j \mathbf{x}_i^T \mathbf{x}_j \\ \text{subject to} & \quad \alpha_i \geq 0 \quad i = 1,2,...,n\\ & \quad \sum_{i=1}^n \alpha_i y_i = 0 \end{align}\]

对偶问题同样是一个二次规划问题，只不过所求的参数从 $\mathbf{w}, b$变为了拉格朗日乘子 $\alpha_i$，通过使用工具包 cvxopt中的solvers.qp 方法，我们可以求解出拉格朗日乘子 $\alpha_i$。和 Primal 问题类似的，我们只需要构造方程适配API调用。求得乘子 $\alpha_i$ 后，我们可以根据公式

\[\begin{align} w^* &= \sum_{i=1}^{n} \alpha_i* y_i x_i \\ b^* &= y_j - w^T x_j \quad \exists \alpha_j > 0 \end{align}\]

恢复法向量 $w^*$，以及截距 $b^*$，得到 hyperplane，实现了原问题和对偶问题最优解的对应关系。

具体的python实现：

  
class Dual_Hard_SVM(SVMModel):
    """Primal Soft-margin SVM problem solved by cvxopt.solvers.qp."""
    def __init__(self) -> None:
        super().__init__()
        self.alpha = None  # shape (n, )
        self.w = None
    
    def fit(self, X, y):
        n, d = X.shape[0], X.shape[1]

        P = np.zeros((n, n)) # shape (n, n)
        for i in range(n):
            for j in range(n):
                P[i, j] = y[i] * y[j] * X[i] @ X[j].reshape(-1, 1)  # (1, d) @ (d, 1) = (1, 1)
        
        q = -np.ones((n, ))  # shape (n, )
        G = -np.eye(n)  # shape (n, n)
        h = np.zeros((n, ))  # shape (n, )
        A = y.reshape(1, -1)  # shape (1, n)
        b = 0.

        sol = solvers.qp(matrix(P), matrix(q), matrix(G), matrix(h), matrix(A), matrix(b))
        self.alpha = np.array(sol['x']).flatten()

        self.w = np.zeros((d, ))
        for i in range(n):
            self.w += self.alpha[i] * y[i] * X[i] # (1, ) * (1, ) * (1, d) = (d, ) 

        for i in range(n):
            if self.alpha[i] > 1e-8:
                b = y[i] - X[i] @ self.w
        self.w = np.concatenate((self.w, np.array([b])), axis=0)  # (d, ) -> (d+1, )

        return self

Soft-margin SVM

现在已经成功解决了 Hard-margin SVM 问题。我们知道，Hard-margin SVM 是基于假设：数据线性可分。但现实中数据往往不是线性可分的，会有 Outliers 的存在，如果直接使用Hard-margin SVM算法会失败。所以我们稍微改进一下 Hard-margin，增加一些可调整的small variables $\xi$。让 SVM 在任何情况下都可以求解，更加robust。下面简单罗列了 primal/hard soft-margin svm 问题的数学公式，具体基于ovxopt.solvers.qp实现见svm。

Primal Soft-margin Proble

\[\begin{align} \text{minimize} & \quad \frac{1}{2} \|\mathbf{w}\|^2 + C \sum_{i=1}^n \xi_i \\ \text{subject to} & \quad y_i(\mathbf{w}^T\mathbf{x}_i + b) \geq 1 - \xi_i \quad \forall i \\ & \quad \xi_i \geq 0 \quad \forall i \end{align}\]

$C$是正则系数。若C设置很大（如1e6），$\xi$就会接近0，和 hard-margin 问题相似了。

对比hard-margin 原问题soft-margin 问题多了 n 个约束，需要求解的向量变为了 $\hat{\mathbf{w}} = \begin{bmatrix}\mathbf{w} & b & \xi \end{bmatrix}$ 维度为 n + d + 1。

Dual Soft-margin Problem

\[\begin{align*} \text{minimize} & \quad -\sum_{i=1}^n \alpha_i + \frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n \alpha_i \alpha_j y_i y_j \mathbf{x}_i^T \mathbf{x}_j \\ \text{subject to} & \quad 0 \leq \alpha_i \leq C \quad \forall i \\ & \quad \sum_{i=1}^n \alpha_i y_i = 0 \end{align*}\]

复杂度分析

可以看到，对于 Soft Margin SVM 中的primal 问题求解QP时，P矩阵的大小为 (n, d+1)， G 矩阵（代表约束）的大小为 (2n, d+1) 而 dual 问题求解QP问题时候 P矩阵大小为 (n, n)，G矩阵的大小为(2n, n)。两者可能有不同的应用场景，比如维度较多时使用 dual 问题更有优势，数据量较多时使用 primal 问题直接求解比较有优势。

支持向量

最后再回顾一下，为什么这样的模型被叫做 Support Vector Machine 呢。出处主要来自于 SVM 的对偶问题，由 KKT 条件中的互补松弛条件，对于 $\alpha_i > 0$ 的情况下，必有 $y (\mathbf{w} \mathbf{x} + b) = 1$。“支持向量“就指的是满足这个条件的样本，所有margin上的向量被称为支持向量。可以发现，最终的模型的复杂度其实与样本数量无关，而与最后的“支持”向量有关。

结果比较与分析

Addition version 只需要修改一下 C，把 C都变成 C/n，然后再关注一下 bias 的符号，基本和 Soft margin SVM 模型没有区别。Code 都在这里。

我使用 7:3 顺序划分了训练集和验证集。下面展示了一些对于不同的C结果对比，主要对比了方法的正确率和 weights 之间的L2 distance。

====== results for C=0.01 ======
distance between primal weights and dual weights: 0.009500775235839143 
distance between primal weights and libsvm weights: 0.005876487666099964 
distance between dual weights and libsvm weights: 0.015372744251807083
primal soft-margin svm acc: 0.9181286549707602
dual soft-margin svm acc: 0.9181286549707602
libsvm acc: 0.9181286549707602
====== results for C=1.0 ======
distance between primal weights and dual weights: 0.10224759687652629 
distance between primal weights and libsvm weights: 0.19837006598001697 
distance between dual weights and libsvm weights: 0.22549480613486989
primal soft-margin svm acc: 0.9473684210526315
dual soft-margin svm acc: 0.935672514619883
libsvm acc: 0.9415204678362573
====== results for C=5 ======
distance between primal weights and dual weights: 0.1863953725032128 
distance between primal weights and libsvm weights: 8.738408498121661 
distance between dual weights and libsvm weights: 8.555750491132391
primal soft-margin svm acc: 0.9473684210526315
dual soft-margin svm acc: 0.935672514619883
libsvm acc: 0.9590643274853801
====== results for C=100000.0 ======
distance between primal weights and dual weights: 0.04013768203724823 
distance between primal weights and libsvm weights: 48297.19872750908 
distance between dual weights and libsvm weights: 48297.18696223501
primal soft-margin svm acc: 0.9415204678362573
dual soft-margin svm acc: 0.9415204678362573
libsvm acc: 0.9239766081871345

可能我的实现与 libsvm 有些不同，选取的 C 越大，我的实现和libsvm的实现的 weights & bias 差距就越大。但是我自己实现的 primal 和 dual 方法得到的 weights and bias 基本一致。可能是 SMO 算法的问题和我使用的 QP solver 结果的精度不同（不确定）。可以看到，选取不同的 C 对结果还是有较大的影响。

TODO

SMO 算法
Kernel Function

Refenences

中科大凸优化课程
OpenAI ChatGPT
沈春华，陈昊 ZJU机器学习课程PPT
李航《统计学习方法》
Boyd《Convex Optimization》