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)$ 的符号来确定类别。

$$\begin{gathered} f(x)=\mathbf{wx}+b \\ y=sign(f(x)) \end{gathered}$$

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

$$distance=\frac{|\mathbf{wx}+b|}{||w|{|}_2}$$

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

$$y(\mathbf{wx}+b)=1$$

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

$$\begin{array}{rl}maximize& \frac{1}{||\mathbf{w}||}=minimize\frac{1}{2}||\mathbf{w}|{|}^2\\ s.t.& y_i(\mathbf{w}{x}_i+b)\ge 1fori=1,2,...n\end{array}$$

求解 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=\left[\begin{array}{cc}{I}_d& 0\\ 0& 0\end{array}\right]$$ $$q=[0]\ast (d+1)$$ $$\begin{gathered} G_{nd}=-y\odot \left[\begin{array}{c}\mathbf{x}|1\end{array}\right] \\ h=[-1]\ast n \end{gathered}$$ $$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{array}{rl}\text{minimize}& f_0(\mathbf{x})\\ \text{subject to}& f_i(\mathbf{x})\le 0\forall i=1,2,\dots ,m\\ & h_i(\mathbf{x})=0\forall i=1,2,\dots ,p\end{array}$$

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

$$\begin{array}{rl}\text{maximize}& g(\lambda ,\nu )\\ \text{subject to}& \lambda \ge 0\end{array}$$

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

$$g(\lambda ,\nu )=\underset{\mathbf{x}}{\inf }\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^{\ast }$ ，原问题的最优解记为 $p^{\ast }$ , 对偶问题的最优解构成了原问题的下解（证明略），即

$$d^{\ast }\le p^{\ast }$$

我们称为弱对偶性，$p^{\ast }-d^{\ast }$ 称为对偶间隙。即使原问题不是凸问题，弱对偶性仍然成立。若

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条件包括以下几个方面： 1. 对于最小化问题，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$$ 2. 不等式约束的松弛条件（complementary slackness）： $$\lambda_i^* g_i(x^*) = 0, \quad i = 1, 2, ..., m$$ 其中，$\lambda_i^*$是不等式约束对应的拉格朗日乘子，$g_i(x^*)$是不等式约束函数在最优解$x^*$处的值。这个条件表示，在最优解中，不等式约束和对应的拉格朗日乘子之间有一个松弛条件。 3. 非负约束条件： $$\lambda_i^* \geq 0, \quad i = 1, 2, ..., m$$ 这个条件要求不等式约束的拉格朗日乘子非负。 4. 等式约束条件： $$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 条件得到： 1. 目标函数的梯度与约束的梯度线性相关 $$\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$ 求偏导数的操作。 2. 互补松弛 $$\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$ 3. 非负条件 $$\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实现： ```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](https://github.com/cjl99/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 都在 [这里](https://github.com/cjl99/svm/)。 我使用 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 1. [中科大凸优化课程](https://www.youtube.com/watch?v=yGmAnZ1iNtA) 2. [OpenAI ChatGPT](https://chat.openai.com/) 3. 沈春华，陈昊 ZJU机器学习课程PPT 4. 李航《统计学习方法》 5. Boyd《Convex Optimization》