Press "Enter" to skip to content

使用Python轻松从头实现多类支持向量机

免费提供SVM的深入概述

在本故事中,我们将以一般软间隔和核化形式实现支持向量机学习算法。我们将首先提供SVM的简要概述以及其训练和推断方程式,然后将其转化为代码以开发SVM模型。随后,我们将扩展我们的实现以处理多类别情况,并通过使用Sci-kit Learn测试我们的模型。

因此,通过本故事的结束:

  • 您将对各种重要的SVM概念有清晰的了解。
  • 您将能够从头开始实现二进制和多类别情况下的SVM模型,并真正理解。
Van Gogh Starry Night with a Line of Symmetry and Stars— Generated by author using DALLE 2

目录

· 简要概述硬间隔支持向量机软间隔支持向量机核软间隔支持向量机· 实现基本导入定义Kernels和SVM超参数定义预测方法定义预测方法测试实现将适应多类别情况泛化将预测泛化到多类别情况测试实现

简要概述

硬间隔支持向量机

SVM的目标是拟合超平面,以获得最大化间隔(两类别中最接近的点之间的距离)。可以证明,并且直观上认为,这样的超平面(A)具有更好的泛化性能,并且比不能最大化间隔的超平面(B)更能抵抗噪音。

Figure by Ennepetaler86 on Wikimedia. CC BY-SA 3.0 Unported.

为了实现这一目标,SVM通过解决以下优化问题来找到超平面的W和b:

使用Python轻松从头实现多类支持向量机 四海 第3张

它试图找到使最接近点的距离最大化且将所有事物正确分类的W,b(如y取±1的约束条件中)。这可以等价地表示为以下优化问题:

使用Python轻松从头实现多类支持向量机 四海 第4张

对于这个问题,可以写出等价的对偶优化问题

使用Python轻松从头实现多类支持向量机 四海 第5张

这个问题的解产生了数据集中每个点的拉格朗日乘子,我们假设数据集大小为m:(α, α₂, …, α_N)。目标函数在α上明显是二次的,约束是线性的,这意味着可以很容易地用二次规划求解。一旦找到解决方案,根据对偶的推导可知:

(xₛ,yₛ) is any point with α>0

请注意,只有α>0的点定义了超平面(对总和有贡献)。这些被称为支持向量。

因此,预测方程在给定新示例x时返回其预测y=±1的结果为:

Involves plugging then doing some algebraic simplification

这种基本形式的SVM称为硬间隔SVM,因为它解决的优化问题(如上所定义)强制要求训练中的所有点都必须被正确分类。在实际情况下,可能会由于一些噪声而阻止或限制存在完全分离数据的超平面,这种情况下优化问题将返回无解或低质量解。

软间隔SVM

Fit by Soft Margin SVM by Mangat et al on Research Gate. CC BY-SA 4.0 International

为了推广硬间隔SVM,软间隔SVM通过引入一个C常数(用户给定的超参数)来调整问题的优化,该常数控制它的“硬度”。具体来说,它修改原始优化问题变为以下形式:

modifications in blue

该问题允许每个点产生一些违规ϵₙ(例如,位于超平面的错误一侧),但仍旨在通过在目标函数中加权它们的总和来减少它们。当C接近无穷大时,它变为硬间隔的等价问题(通常是在接近无穷大之前)。与此同时,较小的C允许更多违规(作为更宽的间隔的代价;即更小的wᵗw)。

相当令人惊讶的是,可以证明等效的对偶问题只需将每个点的α约束为≤C。

使用Python轻松从头实现多类支持向量机 四海 第10张

由于允许违规,支持向量(α>0的点)不再全部位于边缘边缘。可以证明,任何已经违规的支持向量都将具有α=C,并且非支持向量(α=0)不能违规。我们称可能已经违规的支持向量(α=C)为“非边际支持向量”,而其他纯粹的支持向量(没有违规;位于边缘上)称为“边际支持向量”(0<α<C)。

可以证明推断方程不会改变:

使用Python轻松从头实现多类支持向量机 四海 第11张

然而,现在(xₛ,yₛ)现在必须是未违规的支持向量,因为方程假定它位于边缘的边缘(以前,可以使用任何支持向量)。

核函数软间隔支持向量机

软间隔支持向量机扩展了硬间隔支持向量机以处理噪声,但是通常由于噪声之外的因素(如自然非线性),数据无法由超平面分开。在这种情况下可以使用软间隔支持向量机,但是最优解可能涉及比现实中可容忍的错误更多的超平面。

Machine Learner在维基百科上的图像。 CC BY-SA 4.0国际许可协议

核函数软间隔支持向量机将软间隔支持向量机推广到处理数据自然非线性的情况。例如,在左侧显示的示例中,无论设置C为何值,软间隔支持向量机都无法找到线性超平面,以合理地分离数据。

然而,可以通过对数据集中的每个点x进行一些变换函数z = Φ(x),将数据映射到更高维度,使其在新的高维空间中更线性(或完全线性)。这相当于在对偶问题中将x替换为z来得到:

使用Python轻松从头实现多类支持向量机 四海 第13张

实际上,特别是当Φ转换为一个非常高维的空间时,计算z可能需要很长时间。这通过核技巧得以解决,该技巧用一种数学函数(称为核函数)的“等效”计算替换了zᵗz,并且速度更快(例如,zᵗz的代数简化)。例如,这里是一些流行的核函数(每个核函数对应于某个变换Φ到更高维空间):

多项式的次数(Q)和RBF的γ是超参数(由用户决定)

有了这个,对偶优化问题变为:

使用Python轻松从头实现多类支持向量机 四海 第15张

直观地说,推断方程变为(代数推导后):

使用Python轻松从头实现多类支持向量机 四海 第16张

你可以在这里找到基于数学背景的上述方程的完整推导:数学背景

来自Unsplash上Scott Graham的照片

实施

我们将使用以下内容进行实施:

基本导入

首先导入一些基本库:

import numpy as np                  # 用于数组的基本操作
from scipy.spatial import distance  # 用于计算高斯核
import cvxopt                       # 用于解决对偶优化问题
import copy                         # 用于复制numpy数组 

定义核函数和SVM超参数

我们首先使用它们各自的函数定义三个核函数:

多项式的次数(Q)和RBF的γ是超参数(由用户决定)
class SVM:    linear = lambda x, xࠤ , c=0: x @ xࠤ.T    polynomial = lambda x, xࠤ , Q=5: (1 + x @ xࠤ.T)**Q    rbf = lambda x, xࠤ, γ=10: np.exp(-γ*distance.cdist(x, xࠤ,'sqeuclidean'))    kernel_funs = {'linear': linear, 'polynomial': polynomial, 'rbf': rbf}

为了与其他核函数保持一致,线性核函数使用一个额外无用的超参数。很明显,kernel_funs接受一个字符串表示的核函数,并返回相应的核函数。

现在让我们继续定义构造函数:

class SVM:    linear = lambda x, xࠤ , c=0: x @ xࠤ.T    polynomial = lambda x, xࠤ , Q=5: (1 + x @ xࠤ.T)**Q    rbf = lambda x, xࠤ, γ=10: np.exp(-γ*distance.cdist(x, xࠤ,'sqeuclidean'))    kernel_funs = {'linear': linear, 'polynomial': polynomial, 'rbf': rbf}        def __init__(self, kernel='rbf', C=1, k=2):        # 设置超参数        self.kernel_str = kernel        self.kernel = SVM.kernel_funs[kernel]        self.C = C                  # 正则化参数        self.k = k                  # 核函数参数                # 训练数据和支持向量(稍后设置)        self.X, y = None, None        self.αs = None                # 多类分类(稍后设置)        self.multiclass = False        self.clfs = []                          

SVM有三个主要超参数:核函数(我们存储给定的字符串和相应的核函数),正则化参数C和核函数超参数(将传递给核函数);它表示多项式核函数的次数Q和RBF核函数的γ。

定义适应方法

为了在不同的单元格中扩展这个类以定义fitpredict函数,我们将定义以下函数,并将其用作装饰器:

SVMClass = lambda func: setattr(SVM, func.__name__, func) or func

回顾一下,适应SVM对应于通过解决对偶优化问题找到每个点的支持向量α:

使用Python轻松从头实现多类支持向量机 四海 第19张

α 成为一个变量列向量 (α α₂ … α_N)ᵗ,y 成为一个常量列向量对于标签 (y y₂ … y_N)ᵗ,并且 K 成为一个常量矩阵,其中 K[n,m] 计算核函数在 (xₙ, xₘ) 处的值。回忆以下基于索引的点积、外积和二次型的等价关系:

使用Python轻松从头实现多类支持向量机 四海 第20张

以便能够将对偶优化问题写成如下的矩阵形式:

使用Python轻松从头实现多类支持向量机 四海 第21张

根据之前的暗示,我们知道这是一个二次规划问题,我们可以阅读 CVXOPT 文档中关于 Quadratic Programming 的部分:

From CVXOPT Documentation. GNU General Public License

方括号告诉你,你可以只使用 (P,q),或者 (P,q,G,h),或者 (P, q, G, h, A, b) 等等(未给出的部分将使用默认值,如 1)。

为了了解在我们的情况中 (P, q, G, h, A, b) 的值,我们进行以下比较:

使用Python轻松从头实现多类支持向量机 四海 第23张

通过重写第一个式子,使比较更容易了解:

Notice that we changed the max to min by multiplying the function by -1

现在明显了(注意 0≤α 等价于 -α≤0):

使用Python轻松从头实现多类支持向量机 四海 第25张

因此,我们可以编写以下拟合函数:

@SVMClassdef fit(self, X, y, eval_train=False):    # 如果有超过两个唯一的标签,调用多类别版本    if len(np.unique(y)) > 2:        self.multiclass = True        return self.multi_fit(X, y, eval_train)        # 如果标签是 {0,1},将其更改为 {-1,1}    if set(np.unique(y)) == {0, 1}: y[y == 0] = -1    # 确保 y 是一个 Nx1 列向量(CVXOPT 需要)    self.y = y.reshape(-1, 1).astype(np.double) # 必须是列向量    self.X = X    N = X.shape[0]  # 点的数量        # 通过 Numpy 的向量化计算数据中所有可能的 (x, x') 对的核函数    # 将得到矩阵 K    self.K = self.kernel(X, X, self.k)        ### 设置优化参数    # 对于 1/2 x^T P x + q^T x    P = cvxopt.matrix(self.y @ self.y.T * self.K)    q = cvxopt.matrix(-np.ones((N, 1)))        # 对于 Ax = b    A = cvxopt.matrix(self.y.T)    b = cvxopt.matrix(np.zeros(1))    # 对于 Gx <= h    G = cvxopt.matrix(np.vstack((-np.identity(N),                                 np.identity(N))))    h = cvxopt.matrix(np.vstack((np.zeros((N,1)),                                 np.ones((N,1)) * self.C)))    # 求解        cvxopt.solvers.options['show_progress'] = False    sol = cvxopt.solvers.qp(P, q, G, h, A, b)    self.αs = np.array(sol["x"])            # 我们的解决方案            # 一个布尔数组,标记支持向量    self.is_sv = ((self.αs-1e-3 > 0)&(self.αs <= self.C)).squeeze()    # 一些边界支持向量的索引    self.margin_sv = np.argmax((0 < self.αs-1e-3)&(self.αs < self.C-1e-3))        if eval_train:        print(f"训练完成,准确率为{self.evaluate(X, y)}")

我们确保这是一个二元问题,并且二元标签被设定为SVM(±1),并且y是一个具有维度(N,1)的列向量。然后我们解决优化问题,以找到(α α₂ … α_N)ᵗ。

我们使用(α α₂ … α_N)ᵗ来获得一个标志数组,在与支持向量对应的任何索引处为1,以便我们以后只需对支持向量求和并为(xₛ,yₛ)寻找一个边缘支持向量的索引应用预测方程。请注意,在我们的检查中,我们假设非支持向量可能没有完全 α=0,如果 α≤1e-3,则近似为零(我们知道 CVXOPT 的结果可能不是完全精确的)。同样,我们假设非边缘支持向量可能没有 α=C 完全。

定义预测方法

回顾预测方程:

使用Python轻松从头实现多类支持向量机 四海 第16张

@SVMClassdef predict(self, X_t):    if self.multiclass: return self.multi_predict(X_t)    # 计算 (xₛ, yₛ)    xₛ, yₛ = self.X[self.margin_sv, np.newaxis], self.y[self.margin_sv]    # 寻找支持向量    αs, y, X= self.αs[self.is_sv], self.y[self.is_sv], self.X[self.is_sv]    # 计算第二项    b = yₛ - np.sum(αs * y * self.kernel(X, xₛ, self.k), axis=0)    # 计算分数    score = np.sum(αs * y * self.kernel(X, X_t, self.k), axis=0) + b    return np.sign(score).astype(int), score

就是这样。我们还可以实现一个 evaluate 方法来计算准确率(在上面的拟合中使用)。

@SVMClassdef evaluate(self, X,y):      outputs, _ = self.predict(X)    accuracy = np.sum(outputs == y) / len(y)    return round(accuracy, 2)

测试实现

from sklearn.datasets import make_classificationimport numpy as np# 加载数据集np.random.seed(1)X, y = make_classification(n_samples=2500, n_features=5,                            n_redundant=0, n_informative=5,                            n_classes=2,  class_sep=0.3)# 测试已实现的 SVMsvm = SVM(kernel='rbf', k=1)svm.fit(X, y, eval_train=True)y_pred, _ = svm.predict(X)print(f"准确率:{np.sum(y==y_pred)/y.shape[0]}")  #0.9108# 使用 Scikit 进行测试from sklearn.svm import SVCclf = SVC(kernel='rbf', C=1, gamma=1)clf.fit(X, y)y_pred = clf.predict(X)print(f"准确率:{sum(y==y_pred)/y.shape[0]}")    #0.9108

您可以更改数据集和超参数,以进一步确保它们是相同的。理想情况下,通过比较决策函数而不是准确率来进行比较。

将拟合泛化为多类别

@SVMClassdef multi_fit(self, X, y, eval_train=False):    self.k = len(np.unique(y))      # 类别数    # 对于每一对类别    for i in range(self.k):        # 获取该类别的数据        Xs, Ys = X, copy.copy(y)        # 将标签更改为 -1 和 1        Ys[Ys!=i], Ys[Ys==i] = -1, +1        # 拟合分类器        clf = SVM(kernel=self.kernel_str, C=self.C, k=self.k)        clf.fit(Xs, Ys)        # 保存分类器        self.clfs.append(clf)    if eval_train:          print(f"训练完成,准确率为 {self.evaluate(X, y)}")

为了将模型泛化到多类别,我们对每个存在的类别训练一个二元 SVM 分类器,其中我们循环遍历每个类别,并将属于该类别的点重新标记为 +1,将其他所有类别的点标记为 -1。

给定k类别,训练结果是k个分类器,其中第i个分类器是根据第i类别标记为+1,其他类别标记为-1的数据进行训练。

将预测推广到多类别

然后,在新样本上进行预测时,我们选择相应分类器最有信心(得分最高)的类别。

@SVMClassdef multi_predict(self, X):    # 获取所有分类器的预测结果    N = X.shape[0]    preds = np.zeros((N, self.k))    for i, clf in enumerate(self.clfs):        _, preds[:, i] = clf.predict(X)        # 获取最大值及相应的得分    return np.argmax(preds, axis=1), np.max(preds, axis=1)

测试实现

from sklearn.datasets import make_classificationimport numpy as np# 加载数据集np.random.seed(1)X, y = make_classification(n_samples=500, n_features=2,                            n_redundant=0, n_informative=2,                            n_classes=4, n_clusters_per_class=1,                             class_sep=0.3)# 测试 SVMsvm = SVM(kernel='rbf', k=4)svm.fit(X, y, eval_train=True)y_pred = svm.predict(X)print(f"准确率: {np.sum(y==y_pred)/y.shape[0]}") # 0.65# 使用Scikit测试from sklearn.multiclass import OneVsRestClassifierfrom sklearn.svm import SVCclf = OneVsRestClassifier(SVC(kernel='rbf', C=1, gamma=4)).fit(X, y)y_pred = clf.predict(X)print(f"准确率: {sum(y==y_pred)/y.shape[0]}")    # 0.65

绘制每个决策区域的结果如下图所示:

作者的图

请注意,尽管Sci-kit Learn的SVM默认支持OVR(无需像上面显示的明确调用),但这可能还具有特定于SVM的进一步优化。

完整代码

import numpy as np                  # 用于数组的基本操作from scipy.spatial import distance  # 用于计算高斯核import cvxopt                       # 用于求解对偶优化问题import copy                         # 用于复制numpy数组 class SVM:    linear = lambda x, xࠤ , c=0: x @ xࠤ .T    polynomial = lambda x, xࠤ , Q=5: (1 + x @ xࠤ.T)**Q    rbf = lambda x, xࠤ , γ=10: np.exp(-γ * distance.cdist(x, xࠤ,'sqeuclidean'))    kernel_funs = {'linear': linear, 'polynomial': polynomial, 'rbf': rbf}        def __init__(self, kernel='rbf', C=1, k=2):        # 设置超参数        self.kernel_str = kernel        self.kernel = SVM.kernel_funs[kernel]        self.C = C                  # 正则化参数        self.k = k                  # 核参数                # 训练数据和支持向量        self.X, y = None, None        self.αs = None                # 多类别分类        self.multiclass = False        self.clfs = []                                  # 这在这个地方是无用的(仅用于笔记本)SVMClass = lambda func: setattr(SVM, func.__name__, func) or func@SVMClassdef fit(self, X, y, eval_train=False):    if len(np.unique(y)) > 2:        self.multiclass = True        return self.multi_fit(X, y, eval_train)        # 如果需要重新标记    if set(np.unique(y)) == {0, 1}: y[y == 0] = -1    # 确保y的维度为 Nx1    self.y = y.reshape(-1, 1).astype(np.double) # 必须是列向量    self.X = X    N = X.shape[0]        # 计算数据中所有可能的(x, x')对的核    self.K = self.kernel(X, X, self.k)        # 为 1/2 x^T P x + q^T x    P = cvxopt.matrix(self.y @ self.y.T * self.K)    q = cvxopt.matrix(-np.ones((N, 1)))        # 为 Ax = b    A = cvxopt.matrix(self.y.T)    b = cvxopt.matrix(np.zeros(1))    # 为 Gx <= h    G = cvxopt.matrix(np.vstack((-np.identity(N),                                 np.identity(N))))    h = cvxopt.matrix(np.vstack((np.zeros((N,1)),                                 np.ones((N,1)) * self.C)))    # 解决        cvxopt.solvers.options['show_progress'] = False    sol = cvxopt.solvers.qp(P, q, G, h, A, b)    self.αs = np.array(sol["x"])            # 映射到支持向量    self.is_sv = ((self.αs > 1e-3) & (self.αs <= self.C)).squeeze()    self.margin_sv = np.argmax((1e-3 < self.αs) & (self.αs < self.C - 1e-3))        if eval_train:        print(f"训练完成,准确率为 {self.evaluate(X, y)}")@SVMClassdef multi_fit(self, X, y, eval_train=False):    self.k = len(np.unique(y))      # 类别数    # 对于每一对类别    for i in range(self.k):        # 获取该类别对应的数据        Xs, Ys = X, copy.copy(y)        # 将标签变为-1和1        Ys[Ys!=i], Ys[Ys==i] = -1, +1        # 拟合分类器        clf = SVM(kernel=self.kernel_str, C=self.C, k=self.k)        clf.fit(Xs, Ys)        # 保存分类器        self.clfs.append(clf)    if eval_train:        print(f"训练完成,准确率为 {self.evaluate(X, y)}")@SVMClassdef predict(self, X_t):    if self.multiclass: return self.multi_predict(X_t)    xₛ, yₛ = self.X[self.margin_sv, np.newaxis], self.y[self.margin_sv]    αs, y, X= self.αs[self.is_sv], self.y[self.is_sv], self.X[self.is_sv]    b = yₛ - np.sum(αs * y * self.kernel(X, xₛ, self.k), axis=0)    score = np.sum(αs * y * self.kernel(X, X_t, self.k), axis=0) + b    return np.sign(score).astype(int), score@SVMClassdef multi_predict(self, X):    # 获取所有分类器的预测结果    preds = np.zeros((X.shape[0], self.k))    for i, clf in enumerate(self.clfs):        _, preds[:, i] = clf.predict(X)        # 获取最大值及相应的得分    return np.argmax(preds, axis=1)@SVMClassdef evaluate(self, X,y):      outputs, _ = self.predict(X)    accuracy = np.sum(outputs == y) / len(y)    return round(accuracy, 2)from sklearn.datasets import make_classificationimport numpy as np# 加载数据集np.random.seed(1)X, y = make_classification(n_samples=500, n_features=2,           n_redundant=0, n_informative=2, n_classes=4,           n_clusters_per_class=1,  class_sep=0.3)# 测试 SVMsvm = SVM(kernel='rbf', k=4)svm.fit(X, y, eval_train=True)y_pred = svm.predict(X)print(f"准确率: {np.sum(y==y_pred)/y.shape[0]}")# 使用Scikit测试from sklearn.multiclass import OneVsRestClassifierfrom sklearn.svm import SVCclf = OneVsRestClassifier(SVC(kernel='rbf', C=1, gamma=4)).fit(X, y)y_pred = clf.predict(X)print(f"准确率: {sum(y==y_pred)/y.shape[0]}")
Nathan Van Egmond在Unsplash上的照片

总之,我们实现了支持向量机(SVM)学习算法,涵盖了其一般的软间隔和核方法。我们提供了SVM的概述,编写了代码中的模型,扩展了多类别情况下的实现,并使用Sci-kit Learn验证了我们的实现。

希望您能从这个故事中学到一些对您的工作有用的知识。下次见,再见。

资源:

代码大部分是对此处的修改(MIT许可证):Persson, Aladdin. “SVM from Scratch — Machine Learning Python (Support Vector Machine).” YouTube

Leave a Reply

Your email address will not be published. Required fields are marked *