第五章 支持向量机(SVM)

支持向量机(SVM)详解

文章目录

from sklearn.datasets import load_iris
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC  # 线性SVM模型
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import font_manager

my_font = font_manager.FontProperties(fname="C:/Windows/Fonts/simhei.ttf")

支持向量机

支持向量机(简称SVM)是一个功能强大并且全面的机器学习模型,它能够执行线性或非线性分类、回归,甚至是异常值检测任务。它是机器学习领域最受欢迎的模型之一,任何对机器学习感兴趣的人都应该在工具箱中配备一个。SVM特别适用于中小型复杂数据集的分类。

本章将会介绍不同SVM的核心概念,怎么使用它们以及它们的工作原理。

线性SVM分类

大间隔分类

SVM的基本思想可以用一些图来说明。下图所示的数据集来自引用的鸢尾花数据集的一部分。两个类别可以轻松地被一条直线(它们是线性可分离的)分开。左图显示了三种可能的线性分类器的决策边界。其中虚线所代表的模型表现非常糟糕,甚至都无法正确实现分类。其余两个模型在这个训练集上表现堪称完美,但是它们的决策边界与实例过于接近,导致在面对新实例时,表现可能不会太好。相比之下,右图中的实线代表SVM分类器的决策边界,这条线不仅分离了两个类别,并且尽可能远离了最近的训练实例。你可以将SVM分类器视为在类别之间拟合可能的最宽的街道(平行的虚线所示)。因此这也叫作大间隔分类(large margin classification)

from sklearn.svm import SVC
from sklearn import  datasets

iris = datasets.load_iris()
X = iris["data"][:, (2, 3)]  # petal length,petal width
y = iris["target"]

setosa_or_versicolor = (y == 0) | (y == 1)
X = X[setosa_or_versicolor]
y = y[setosa_or_versicolor]

# SVM分类模型
svm_clf = SVC(kernel="linear", C=float("inf"))
svm_clf.fit(X, y)
SVC(C=inf, cache_size=200, class_weight=None, coef0=0.0,
  decision_function_shape='ovr', degree=3, gamma='auto', kernel='linear',
  max_iter=-1, probability=False, random_state=None, shrinking=True,
  tol=0.001, verbose=False)
# 不好的模型
import numpy as np

x0 = np.linspace(0, 5.5, 200)
pred_1 = 5*x0 - 20
pred_2 = x0 - 1.8
pred_3 = 0.1*x0 + 0.5

def plot_svc_decision_boundary(svm_clf, xmin, xmax):
    w = svm_clf.coef_[0]
    b = svm_clf.intercept_[0]
    
    # w0*x0 + w1*x1 + b = 0
    # => x1 = -w0/w1 *x0 -b/w1
    x0 = np.linspace(xmin, xmax, 200)
    decision_boundary = -w[0]/w[1] * x0 - b/w[1]
    
    margin = 1/w[1]
    gutter_up = decision_boundary + margin
    gutter_down = decision_boundary - margin
    
    svs = svm_clf.support_vectors_
    plt.scatter(svs[:, 0], svs[:, 1], s=180, facecolors="#FFAAAA")
    plt.plot(x0, decision_boundary, "k-", linewidth=2)
    plt.plot(x0, gutter_up, "k--", linewidth=2)
    plt.plot(x0, gutter_down, "k--", linewidth=2)

plt.figure(figsize=(12, 2.7))
plt.suptitle("图1:大间隔分类", fontproperties=my_font, fontsize=14)
plt.subplot(121)
plt.plot(x0, pred_1, "g--", linewidth=2)
plt.plot(x0, pred_2, "m-", linewidth=2)
plt.plot(x0, pred_3, "r-", linewidth=2)

plt.plot(X[:, 0][y==1], X[:, 1][y==1], "bs", label="Iris-Versicolor")
plt.plot(X[:, 0][y==0], X[:, 1][y==0], "yo", label="Iris-Setosa")
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.legend(loc="upper left", fontsize=14)
plt.axis([0, 5.5, 0, 2])

plt.subplot(122)
plot_svc_decision_boundary(svm_clf, 0, 5.5)
plt.plot(X[:, 0][y==1], X[:, 1][y==1], "bs", label="Iris-Versicolor")
plt.plot(X[:, 0][y==0], X[:, 1][y==0], "yo", label="Iris-Setosa")
plt.xlabel("Petal length", fontsize=14)
plt.axis([0, 5.5, 0, 2])


[0, 5.5, 0, 2]

第五章 支持向量机(SVM)

注意:在街道以外的地方增加更多训练实例,不会对决策边界产生影响:也就是说它完全由位于街道边缘的实例所决定(或者称之为“支持”)。这些实例被称为支持向量(在上图中已圈出)。

特征缩放的敏感度

SVM对特征的缩放非常敏感,如下图所示,在左图中,垂直刻度比水平刻度大得多,因此可能的最宽的街道接近于水平。在特征缩放(例如使用Scikit-Learn的StandardScaler)后,决策边界看起来好很多(见右图)。

from sklearn.svm import SVC
Xs = np.array([[1, 50], [5, 20], [3, 80], [5, 60]]).astype(np.float64)
ys = np.array([0, 0, 1, 1])
svm_clf = SVC(kernel="linear", C=100)
svm_clf.fit(Xs, ys)
SVC(C=100, cache_size=200, class_weight=None, coef0=0.0,
  decision_function_shape='ovr', degree=3, gamma='auto', kernel='linear',
  max_iter=-1, probability=False, random_state=None, shrinking=True,
  tol=0.001, verbose=False)
from sklearn.svm import SVC
Xs = np.array([[1, 50], [5, 20], [3, 80], [5, 60]]).astype(np.float64)
ys = np.array([0, 0, 1, 1])
svm_clf = SVC(kernel="linear", C=100)
svm_clf.fit(Xs, ys)

plt.figure(figsize=(12, 3.2))
plt.suptitle("图2;对特征缩放的敏感度", fontproperties=my_font, fontsize=14)

plt.subplot(121)
plt.plot(Xs[:, 0][ys == 1], Xs[:, 1][ys == 1], "bo")
plt.plot(Xs[:, 0][ys == 0], Xs[:, 1][ys == 0], "ms")
plot_svc_decision_boundary(svm_clf, 0, 6)
plt.xlabel("$x_0$", fontsize=20,rotation=0)
plt.ylabel("$x_1$", fontsize=20,rotation=0)
plt.title("Unscaled", fontsize=12)
plt.axis([0, 6, 0, 90])

from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_scaled = scaler.fit_transform(Xs)
svm_clf.fit(X_scaled, ys)

plt.subplot(122)
plt.plot(X_scaled[:, 0][ys==1], X_scaled[:, 1][ys==1], "bo")
plt.plot(X_scaled[:, 0][ys==0], X_scaled[:, 1][ys==0], "ms")
plot_svc_decision_boundary(svm_clf, -2, 2)
plt.xlabel("$x_0$", fontsize=20,rotation=0)
plt.ylabel("$x_1$", fontsize=20,rotation=0)
plt.title("scaled", fontsize=12)
plt.axis([-2, 2, -2, 2])
plt.savefig("./plot_data/SVM/sensitivity_to_feature_scales_plot.png")

第五章 支持向量机(SVM)

软间隔分类

硬间隔对异常值的敏感度

如果我们严格地让所有实例都不在街道上,并且位于正确的一边,这就是硬间隔分类

硬间隔分类有两个主要问题,首先,它只在数据是线性可分离的时候才有效;其次,它对异常值非常敏感。图3显示了有一个额外异常值的鸢尾花数据:左图的数据根本找不出硬间隔,而右图最终显示的决策边界与我们在图1中所看到的无异常值时的决策边界也大不相同,可能无法很好地泛化。

X_outliers = np.array([[3.4, 1.3], [3.2, 0.8]])
y_outliers = np.array([0, 0])
Xo1 = np.concatenate((X, X_outliers[:1]), axis=0)  # array添加元素,都必须为list, 注意:输入必须为tuple
yo1 = np.concatenate((y, y_outliers[:1]), axis=0)
Xo2 = np.concatenate((X, X_outliers[1:]), axis=0)
yo2 = np.concatenate((y, y_outliers[1:]), axis=0)

svm_clf2 = SVC(kernel="linear", C=10**9)
svm_clf2.fit(Xo2, yo2)

plt.figure(figsize=(12, 2.7))
plt.suptitle("图3:硬间隔对异常值的敏感度", fontproperties=my_font, fontsize=14)

plt.subplot(121)
plt.plot(Xo1[:, 0][yo1==1], Xo1[:, 1][yo1==1], "bs")
plt.plot(Xo1[:, 0][yo1==0], Xo1[:, 1][yo1==0], "yo")
plt.text(0.3, 1.0, "Impossible!", fontsize=24, color="red")
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.annotate("Outlier",
             xy=(X_outliers[0][0], X_outliers[0][1]),
             xytext=(2.5, 1.7),
             ha="center",
             arrowprops=dict(facecolor="black", shrink=0.1),
             fontsize=16,
            )
plt.axis([0, 5.5, 0, 2])

plt.subplot(122)
plt.plot(Xo2[:, 0][yo2==1], Xo2[:, 1][yo2==1], "bs")
plt.plot(Xo2[:, 0][yo2==0], Xo2[:, 1][yo2==0], "yo")
plot_svc_decision_boundary(svm_clf2, 0, 5.5)
plt.xlabel("Petal length", fontsize=14)
plt.annotate("Outlier",
             xy=(X_outliers[1][0], X_outliers[1][1]),
             xytext=(3.2, 0.08),
             ha="center",
             arrowprops=dict(facecolor="black", shrink=0.1),
             fontsize=16,
            )
plt.axis([0, 5.5, 0, 2])

plt.savefig("./plot_data/SVM/sensiticity_to_outliers_plot.png")
plt.show()

第五章 支持向量机(SVM)

街道宽阔 vs 限制间隔违例
(之间找到良好的平衡,即:软间隔分类)

要避免这些问题,最好使用更灵活的模型。目标是尽可能在保持街道宽阔和限制间隔违例(即位于街道之上,甚至在错误的一边的实例)之间找到良好的平衡,这就是软间隔分类。在Scikit-Learn的SVM类中,可以通过超参数C来控制这个平衡:C值越小,则街道越宽,但是间隔违例也会越多。图4显示了在一个非线性可分离数据集上,两个软间隔SVM分类器各自的决策边界和间隔。左边使用了高C值,分类器的间隔违例较少,但是间隔也较小。右边使用了低C值,间隔大了很多,但是位于街道上的实例也更多。看起来第二个分类器的泛化效果更好,因为大多数间隔违例实际上都位于决策边界正确的一边,所以即便是在该训练集上,它做出的错误预测也会更少。

如果你的SVM模型过度拟合,可以试试通过降低C来进行正则化。

加载鸢尾花数据集,缩放特征,然后训练一个线性SVM模型(使用LinearSVC类,C=0.1,用即将介绍的hinge损失函数)用来检测Virginica鸢尾花。得到的模型如图4右图所示。

import numpy as np
from sklearn import datasets
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC

iris = datasets.load_iris()
X = iris["data"][:, (2,3)]  # petal length, petal width
y = (iris["target"]==2).astype(np.float64)  # Iris-Virginica

svm_clf = Pipeline([("scaler", StandardScaler()),
                    ("linear_svc", LinearSVC(C=0.1, loss="hinge", random_state=42)),
                   ])

svm_clf.fit(X,y)
Pipeline(memory=None,
     steps=[('scaler', StandardScaler(copy=True, with_mean=True, with_std=True)), ('linear_svc', LinearSVC(C=0.1, class_weight=None, dual=True, fit_intercept=True,
     intercept_scaling=1, loss='hinge', max_iter=1000, multi_class='ovr',
     penalty='l2', random_state=42, tol=0.0001, verbose=0))])
# 模型做出预测
svm_clf.predict([[5.5, 1.7]])
array([1.])

与Logistic回归分类器不同的是,SVM分类器不会输出每个类别的概率。

或者,你还可以选择SVC类,使用SVC(kernel=“linear”,C=1),但是这要慢得多,特别是对于大型训练集而言,因此不推荐使用。另一个选择是SGDClassifier类,使用SGDClassifier(loss=“hinge”,alpha=1/(m*C))。这适用于常规随机梯度下降(参见第4章)来训练线性SVM分类器。它不会像LinearSVC类那样快速收敛,但是对于内存处理不了的大型数据集(核外训练)或是在线分类任务,它非常有效。

现在让我们生成比较不同正则化设置的图表:

LinearSVC类会对偏置项进行正则化,所以你需要先减去平均
值,使训练集集中。如果使用StandardScaler会自动进行这一步。此
外,请确保超参数loss设置为"hinge",因为它不是默认值。最后,为
了获得更好的性能,还应该将超参数dual设置为False,除非特征数量
比训练实例还多。

scaler = StandardScaler()

svm_clf1 = LinearSVC(C=1, loss="hinge", random_state=42)
svm_clf2 = LinearSVC(C=100, loss="hinge", random_state=42)

scaled_svm_clf1 = Pipeline([("scaler", scaler),
                            ("linear_svc", svm_clf1),
                           ])

scaled_svm_clf2 = Pipeline([("scaler", scaler),
                           ("linear_svc", svm_clf2),
                           ])

scaled_svm_clf1.fit(X, y)
scaled_svm_clf2.fit(X, y)
Pipeline(memory=None,
     steps=[('scaler', StandardScaler(copy=True, with_mean=True, with_std=True)), ('linear_svc', LinearSVC(C=100, class_weight=None, dual=True, fit_intercept=True,
     intercept_scaling=1, loss='hinge', max_iter=1000, multi_class='ovr',
     penalty='l2', random_state=42, tol=0.0001, verbose=0))])
# 转换为未缩放的参数
# svm_clf1.coef_ 权重系数
# svm_clf1.intercept_ 截距
# scaler.scale_ 缩放比例

b1 = svm_clf1.decision_function([-scaler.mean_/scaler.scale_])
b2 = svm_clf2.decision_function([-scaler.mean_/scaler.scale_])
w1 = svm_clf1.coef_[0] / scaler.scale_
w2 = svm_clf2.coef_[0] / scaler.scale_
svm_clf1.intercept_ = np.array([b1])
svm_clf2.intercept_ = np.array([b2])
svm_clf1.coef_ = np.array([w1])
svm_clf2.coef_ = np.array([w2])

# 查找支持向量(LinearSVC不会自动执行此操作)
t = y * 2 - 1
support_vectors_idx1 = (t * (X.dot(w1) + b1) < 1).ravel()
support_vectors_idx2 = (t * (X.dot(w2) + b2) < 1).ravel()
svm_clf1.support_vectors_ = X[support_vectors_idx1]
svm_clf2.support_vectors_ = X[support_vectors_idx2]
plt.figure(figsize=(12, 3.2))
plt.suptitle("图4:较少间隔违例核大间隔对比", fontproperties=my_font, fontsize=14)
plt.subplot(121)
plt.plot(X[:, 0][y==1], X[:, 1][y==1], "g^",label="Iris-Virginica")
plt.plot(X[:, 0][y==0], X[:, 1][y==0], "bs",label="Iris-Versicolor")
plot_svc_decision_boundary(svm_clf1, 4, 6)
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.legend(loc="upper left", fontsize=14)
plt.title("$C = {}$".format(svm_clf1.C), fontsize=16)
plt.axis([4, 6, 0.8, 2.8])

plt.subplot(122)
plt.plot(X[:, 0][y==1], X[:, 1][y==1], "g^",label="Iris-Virginica")
plt.plot(X[:, 0][y==0], X[:, 1][y==0], "bs",label="Iris-Versicolor")
plot_svc_decision_boundary(svm_clf2, 4, 6)
plt.xlabel("Petal length", fontsize=14)
plt.title("$C = {}$".format(svm_clf2.C), fontsize=16)
plt.axis([4, 6, 0.8, 2.8])

plt.savefig("./plot_data/SVM/regularization_plot.png")
plt.show()

第五章 支持向量机(SVM)

非线性SVM分类

虽然在许多情况下,线性SVM分类器是有效的,并且通常出人意料的好,但是,有很多数据集远不是线性可分离的。处理非线性数据集的方法之一是添加更多特征,比如多项式特征,某些情况下,这可能导致数据集变得线性可分离。参见图5的左图:这是一个简单的数据集,只有一个特征x1x_1x1​,可以看出,数据集线性不可分。但是如果添加第二个特征x2=(x1)2x_2=(x_1)^2x2​=(x1​)2,生成的2D数据集则完全线性可分离。

X1D = np.linspace(-4, 4, 9).reshape(-1, 1)
X2D = np.c_[X1D, X1D**2]
y = np.array([0, 0, 1, 1, 1, 1, 1, 0, 0])

plt.figure(figsize=(11, 4))
plt.suptitle("图5:通过添加特征使数据集线性可分离", fontproperties=my_font, fontsize=14)

plt.subplot(121)
plt.grid(True, which="both")
plt.axhline(y=0,color="k")
plt.plot(X1D[:, 0][y==0], np.zeros(4), "bs")
plt.plot(X1D[:, 0][y==1], np.zeros(5), "g^")
plt.gca().get_yaxis().set_ticks([])
plt.axis([-4.5, 4.5, -0.2, 0.2])

plt.subplot(122)
plt.grid(True, which="both")
plt.axhline(y=0,color="k")
plt.axvline(x=0, color="k")
plt.plot(X2D[:, 0][y==0], X2D[:, 1][y==0], "bs")
plt.plot(X2D[:, 0][y==1], X2D[:, 1][y==1], "g^")
plt.xlabel("$x_1$", fontsize=20, rotation=0)
plt.ylabel("$x_2$", fontsize=20)
plt.gca().get_yaxis().set_ticks([0, 4, 8, 12, 16])
plt.plot([-4.5, 4.5], [6.5, 6.5], "r-.", linewidth=3)
plt.axis([-4.5, 4.5, -1, 17])

plt.show()

第五章 支持向量机(SVM)

非线性SVM分类(卫星数据集)

from sklearn.datasets import make_moons
X, y = make_moons(n_samples=100, noise=0.15, random_state=42)

def plot_dataset(X, y, ases):
    plt.plot(X[:, 0][y==0], X[:, 1][y==0], "bs")
    plt.plot(X[:, 0][y==1], X[:, 1][y==1], "g^")
    plt.axis(axes)
    plt.grid(True, which="both")
    plt.xlabel("$x_1$", fontsize=20, rotation=0)
    plt.ylabel("$x_2$", fontsize=20)

axes = [-1.5, 2.5, -1, 1.5]
plot_dataset(X, y, axes)
plt.show()

第五章 支持向量机(SVM)

from sklearn.datasets import make_moons
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.svm import LinearSVC
polynomial_svm_clf = Pipeline([
    ("poly_fetures", PolynomialFeatures(degree=3)), # 生成2次多项式为(1,a,b,a^2,ab, b^2)
    ("scaler", StandardScaler()),
    ("svm_clf", LinearSVC(C=10, loss = "hinge", random_state=42))
])
polynomial_svm_clf.fit(X,y)
Pipeline(memory=None,
     steps=[('poly_fetures', PolynomialFeatures(degree=3, include_bias=True, interaction_only=False)), ('scaler', StandardScaler(copy=True, with_mean=True, with_std=True)), ('svm_clf', LinearSVC(C=10, class_weight=None, dual=True, fit_intercept=True,
     intercept_scaling=1, loss='hinge', max_iter=1000, multi_class='ovr',
     penalty='l2', random_state=42, tol=0.0001, verbose=0))])
def plot_predictions(clf, axes):
    x0s = np.linspace(axes[0], axes[1], 100)
    x1s = np.linspace(axes[2], axes[3], 100)
    x0, x1 = np.meshgrid(x0s, x1s)  # 从坐标向量中返回坐标矩阵
    X = np.c_[x0.ravel(), x1.ravel()]
    y_pred = clf.predict(X).reshape(x0.shape)
    y_decision = clf.decision_function(X).reshape(x0.shape)
    plt.contourf(x0, x1, y_pred, cmap=plt.cm.brg, alpha=0.2)  #  等高线图
    plt.contourf(x0, x1, y_decision, cmap=plt.cm.brg, alpha=0.1)

plt.suptitle("图6:使用多项式特征的线性LVM分类器", fontproperties=my_font, fontsize=14)
plot_predictions(polynomial_svm_clf, [-1.5, 2.5, -1, 1.5])
plot_dataset(X, y, [-1.5, 2.5, -1, 1.5])

# save_fig("moons_polynomial_svc_plot")
plt.show()

第五章 支持向量机(SVM)

多项式核

添加多项式特征实现起来非常简单,并且对所有的机器学习算法(不只是SVM)都非常有效。但问题是,如果多项式太低阶,处理不了非常复杂的数据集,而高阶则会创造出大量的特征,导致模型变得太慢。

幸运的是,使用SVM时,有一个魔术般的数学技巧可以应用,这就是核技巧(稍后解释)。它产生的结果就跟添加了许多多项式特征,甚至是非常高阶的多项式特征一样,但实际上并不需要真的添加。因为实际没有添加任何特征,所以也就不存在数量爆炸的组合特征了。这个技巧由SVC类来实现,我们看看在卫星数据集上的测试:

from sklearn.svm import SVC

poly_kernel_svm_clf = Pipeline([("scaler", StandardScaler()),
                               ("svm_clf", SVC(kernel="poly", degree=3, coef0=1, C=5)),
                               ])
poly_kernel_svm_clf.fit(X, y)
Pipeline(memory=None,
     steps=[('scaler', StandardScaler(copy=True, with_mean=True, with_std=True)), ('svm_clf', SVC(C=5, cache_size=200, class_weight=None, coef0=1,
  decision_function_shape='ovr', degree=3, gamma='auto', kernel='poly',
  max_iter=-1, probability=False, random_state=None, shrinking=True,
  tol=0.001, verbose=False))])
from sklearn.svm import SVC

poly100_kernel_svm_clf = Pipeline([("scaler", StandardScaler()),
                               ("svm_clf", SVC(kernel="poly", degree=10, coef0=100, C=5)),
                               ])
poly100_kernel_svm_clf.fit(X, y)
Pipeline(memory=None,
     steps=[('scaler', StandardScaler(copy=True, with_mean=True, with_std=True)), ('svm_clf', SVC(C=5, cache_size=200, class_weight=None, coef0=100,
  decision_function_shape='ovr', degree=10, gamma='auto', kernel='poly',
  max_iter=-1, probability=False, random_state=None, shrinking=True,
  tol=0.001, verbose=False))])
plt.figure(figsize=(11, 4))
plt.suptitle("图7:多项式核的SVM分类器", fontproperties=my_font, fontsize=14)

plt.subplot(121)
plot_predictions(poly_kernel_svm_clf, [-1.5, 2.5, -1, 1.5])
plot_dataset(X, y, [-1.5, 2.5, -1, 1.5])
plt.title("$d=3, r=1, C=5$", fontsize=18)

plt.subplot(122)
plot_predictions(poly100_kernel_svm_clf, [-1.5, 2.5, -1, 1.5])
plot_dataset(X, y, [-1.5, 2.5, -1, 1.5])
plt.title("$d=10, r=100, C=5$", fontsize=18)

plt.savefig("./plot_data/SVM/moons_kernelized_polynormial_svc_plot.png")
plt.show()

第五章 支持向量机(SVM)

如图7的左图使用了一个3阶多项式内核训练SVM分类器。而右图是另一个使用了10阶多项式核的SVM分类器。显然,如果模型过度拟合,你应该降低多项式阶数;反过来,如果拟合不足,则可以尝试使之提升。超参数coef0控制的是模型受高阶多项式还是低阶多项式影响的程度。

寻找正确的超参数:

寻找正确的超参数值的常用方法是网格搜索。先进行一次粗略的网格搜索,然后在最好的值附近展开一轮更精细的网格搜索,这样通常会快一些。多了解每个超参数实际上是用来做什么的,有助于你在超参数空间层正确搜索。

添加相似特征

解决非线性问题的另一种技术是添加相似特征。这些特征经过相似函数计算得出,相似函数可以测量每个实例与一个特定地标(landmark)之间的相似度。以前面提到过的一维数据集为例,在x1=-2和x1=1处添加两个地标(见下图中的左图)。接下来,我们采用高斯径向基函数(RBF)作为相似函数,y=0.3(见下等式)。

公式1:高斯RBF:
ϕγ(X,)=exp(γx2) \phi\gamma(X, \ell) = exp(-\gamma\parallel x-\ell \parallel^{2}) ϕγ(X,ℓ)=exp(−γ∥x−ℓ∥2)

这是一个从0(离地标差得非常远)到1(跟地标一样)变化的钟形函数。现在我们准备计算新特征。例如,我们看实例x1=1x_1=-1x1​=−1:它与第一个地标的距离为1,与第二个地标的距离为2。因此它的新特征为x2=eps(0.3×12)0.74x3=eps(0.3×22)0.30x_2=eps(-0.3×1^2)≈0.74,x_3=eps(-0.3×2^2)≈0.30x2​=eps(−0.3×12)≈0.74,x3​=eps(−0.3×22)≈0.30。下图的右图显示了转换后的数据集(去除了原始特征),现在你可以看出,数据呈线性可分离的了。

def gaussian_rbf(x, landmark, gamma):
    return np.exp(-gamma * np.linalg.norm(x - landmark, axis=1)**2)

gamma = 0.3

x1s = np.linspace(-4.5, 4.5, 200).reshape(-1, 1)
x2s = gaussian_rbf(x1s, -2, gamma)
x3s = gaussian_rbf(x1s, 1, gamma)

Xk = np.c_[gaussian_rbf(X1D, -2, gamma), gaussian_rbf(X1D, 1, gamma)]
yk = np.array([0, 0, 1, 1, 1, 1, 1, 0, 0])

plt.figure(figsize=(11, 4))
plt.suptitle("图8:使用高斯RBF的相似特征", fontproperties=my_font, fontsize=14)

plt.subplot(121)
plt.grid(True, which="both")
plt.axhline(y=0, color="k")
plt.scatter(x=[-2, 1], y=[0, 0], s=150, alpha=0.5, c='red')
plt.plot(X1D[:, 0][yk==0], np.zeros(4), "bs")
plt.plot(X1D[:, 0][yk==1], np.zeros(5), "g^")
plt.plot(x1s, x2s, "g-.")
plt.plot(x1s, x3s, "b:")
plt.gca().get_yaxis().set_ticks([0, 0.25, 0.5, 0.75, 1])
plt.xlabel("$x_1$",fontsize=20)
plt.xlabel("$Similarity$",fontsize=14)
plt.annotate("$\mathbf{x}$",
             xy=(X1D[3, 0], 0),
             xytext=(-0.5, 0.2),
             ha="center",
             arrowprops=dict(facecolor="black", shrink=0.1),
             fontsize=18,
            )
plt.text(-2, 0.9, "$x_2$", ha="center", fontsize=20)
plt.text(1, 0.9, "$x_3$", ha="center", fontsize=20)
plt.axis([-4.5, 4.5, -0.1, 1.1])

plt.subplot(122)
plt.grid(True, which="both")
plt.axhline(y=0, color="k")
plt.axvline(x=0, color="k")
plt.plot(Xk[:, 0][yk==0], Xk[:, 1][yk==0], "bs")
plt.plot(Xk[:, 0][yk==1], Xk[:, 1][yk==1], "g^")
plt.xlabel("$x_2$", fontsize=20, rotation=0)
plt.ylabel("$x_2$", fontsize=20)
plt.annotate(r"$\phi\left(\mathbf{x}\right)$",
             xy=(Xk[3, 0], Xk[3, 1]),
             xytext=(0.65, 0.50),
             ha="center",
             arrowprops=dict(facecolor="black", shrink=0.1),
             fontsize=18,
            )
plt.plot([-0.1, 1.1], [0.57, -0.1], "r-.", linewidth=3)
plt.axis([-0.1, 1.1, -0.1, 1.1])

plt.subplots_adjust(right=1)
plt.savefig("./plot_data/SVM/kernel_method_plot.png")
plt.show()

第五章 支持向量机(SVM)

你可能想问怎么选择地标呢?最简单的方法是在数据集里每一个实例的位置上创建一个地标。这会创造出许多维度,因而也增加了转换后的训练集线性可分离的机会。缺点是,一个有m个实例n个特征的训练集会被转换成一个m个实例m个特征的训练集(假设抛弃了原始特征)。如果训练集非常大,那就会得到同样大数量的特征。

# 计算新的特征值
x1_example = X1D[3, 0]
for landmark in (-2,1):
    k = gaussian_rbf(np.array([[x1_example]]), np.array([[landmark]]), gamma)
    print("Pho({}, {}) = {}".format(x1_example, landmark, k))
Pho(-1.0, -2) = [0.74081822]
Pho(-1.0, 1) = [0.30119421]

高斯RBF核函数

与多项式特征方法一样,相似特征法也可以用任意机器学习算法,但是要计算出所有附加特征,其计算代价可能非常昂贵,尤其是对大型训练集来说。然而,核技巧再一次施展了它的SVM魔术:它能够产生的结果就跟添加了许多相似特征一样,但实际上也并不需要添加。我们来使用SVC类试试高斯RBF核:

rbf_kernel_svm_clf = Pipeline([("scaler", StandardScaler()),
                              ("svm_clf", SVC(kernel="rbf", gamma=5, C=0.01)),
                              ])
rbf_kernel_svm_clf.fit(X, y)
Pipeline(memory=None,
     steps=[('scaler', StandardScaler(copy=True, with_mean=True, with_std=True)), ('svm_clf', SVC(C=0.01, cache_size=200, class_weight=None, coef0=0.0,
  decision_function_shape='ovr', degree=3, gamma=5, kernel='rbf',
  max_iter=-1, probability=False, random_state=None, shrinking=True,
  tol=0.001, verbose=False))])

下图的左下方显示了这个模型。其他图显示了超参数gamma()和C使用不同值时的模型。增加gamma值会使钟形曲线变得更窄(上图8的左图),因此每个实例的影响范围随之变小:决策边界变得更不规则,开始围着单个实例绕弯。反过来,减小gamma值使钟形曲线变得更宽,因而每个实例的影响范围增大,决策边界变得更平坦。所以就像是一个正则化的超参数:模型过度拟合,就降低它的值,如果拟合不足则提升它的值(类似超参数C)。

还有一些其他较少用到的核函数,例如专门针对特定数据结构的核函数。字符串核常用于文本文档或是DNA序列(如使用字符串子序列核或是基于莱文斯坦距离的核函数)的分类。

from sklearn.svm import SVC

gamma1, gamma2 = 0.1, 5
C1, C2 = 0.001, 1000
hyperparams = (gamma1, C1), (gamma1, C2), (gamma2, C1), (gamma2, C2)

svm_clfs = []
for gamma, C in hyperparams:
    rbf_kernel_svm_clf = Pipeline([("scaler", StandardScaler()),
                                  ("svm_clf", SVC(kernel="rbf", gamma=gamma, C=C)),
                                  ])
    rbf_kernel_svm_clf.fit(X, y)
    svm_clfs.append(rbf_kernel_svm_clf)

plt.figure(figsize=(20, 11))
plt.suptitle("图9:使用RBF核的SVM分类器", fontproperties=my_font, fontsize=20)

for i,svm_clf in enumerate(svm_clfs):
    plt.subplot(221 + i)
    plot_predictions(svm_clf, [-1.5, 2.5, -1, 1.5])
    plot_dataset(X, y, [-1.5, 2.5, -1, 1.5])
    gamma, C = hyperparams[i]
    plt.title("$\gamma={}, C ={}$".format(gamma, C), fontsize=16)
    
plt.savefig("./plot_data/SVM/moons_rbf_svc_plot.png")
plt.show()

第五章 支持向量机(SVM)

有这么多的核函数,该如何决定使用哪一个呢?有一个经验法则是,永远先从线性核函数开始尝试(要记住,LinearSVC比SVC(kernel=“linear”)快得多),特别是训练集非常大或特征非常多的时候。如果训练集不太大,你可以试试高斯RBF核,大多数情况下它都非常好用。如果你还有多余的时间和计算能力,你可以使用交叉验证和网格搜索来尝试一些其他的核函数,特别是那些专门针对你的数据集数据结构的核函数。

计算复杂度

liblinear库为线性SVM实现了一个优化算法,LinearSVC正是基于该库的。这个算法不支持核技巧,不过它与训练实例的数量和特征数量几乎呈线性相关:其训练时间复杂度大致为O(m×n)。

如果你想要非常高的精度,算法需要的时间更长。它由容差超参数(在Scikit-Learn中为tol)来控制。大多数分类任务中,默认的容差就够了。

SVC则是基于libsvm库的,这个库的算法支持核技巧。训练时间复杂度通常在O(m2×n)和O(m3×n)之间。很不幸,这意味着如果训练实例的数量变大(例如上十万个实例),它将会慢得可怕,所以这个算法完美适用于复杂但是中小型的训练集。但是,它还是可以良好适应地特征数量的增加,特别是应对稀疏特征(即,每个实例仅有少量的非零特征)。在这种情况下,算法复杂度大致与实例的平均非零特征数成比例。下表比较了Scikit-Learn的SVM分类器类别。

时间复杂度 是否制程核外 是否需要缩放 核技巧
LinearSVC O(m*n)
SDGClassifier O(m*n)
SVC O(m2n)toO(m3n)O(m^{2}*n)to O(m^{3}*n)O(m2∗n)toO(m3∗n)

SVM回归

正如前面提到的,SVM算法非常全面:它不仅支持线性和非线性分类,而且还支持线性和非线性回归。诀窍在于将目标反转一下:不再是尝试拟合两个类别之间可能的最宽的街道的同时限制间隔违例,SVM回归要做的是让尽可能多的实例位于街道上,同时限制间隔违例(也就是不在街道上的实例)。街道的宽度由超参数ε控制。图5-10显示了用随机线性数据训练的两个线性SVM回归模型,一个间隔较大(ε=1.5),另一个间隔较小(ε=0.5)。

np.random.seed(42)  # 设置相同的seed,每次生成的随机数也相同,否则每次生成的随机数都会不一样
m = 50
X = 2 * np.random.rand(m, 1) # np.random.rand 随机生成[0,1)之间的数
y = (4 + 3 * X + np.random.randn(m, 1)).ravel()  # np.random.randn 从标准正态分布中返回一个或多个样本值
from sklearn.svm import LinearSVR

svm_reg = LinearSVR(epsilon=1.5, random_state=42)  # random_state 是一个随机种子,是在任意带有随机性的类或函数里作为参数来控制随机模式,当其取某一值,也就确定了一种规则
svm_reg.fit(X, y)
LinearSVR(C=1.0, dual=True, epsilon=1.5, fit_intercept=True,
     intercept_scaling=1.0, loss='epsilon_insensitive', max_iter=1000,
     random_state=42, tol=0.0001, verbose=0)
svm_reg1 = LinearSVR(epsilon=1.5, random_state=42)
svm_reg2 = LinearSVR(epsilon=0.5, random_state=42)
svm_reg1.fit(X, y)
svm_reg2.fit(X, y)

def find_support_vectors(svm_reg, X, y):
    y_pred = svm_reg.predict(X)
    off_margin = (np.abs(y - y_pred) >= svm_reg.epsilon)
    return np.where(off_margin)

svm_reg1.support_ = find_support_vectors(svm_reg1, X, y)
svm_reg2.support_ = find_support_vectors(svm_reg2, X, y)

eps_x1 = 1
eps_y_pred = svm_reg1.predict([[eps_x1]])
def plot_svm_regression(svm_reg, X, y, axes):
    x1s = np.linspace(axes[0], axes[1], 100).reshape(100,-1)
    y_pred = svm_reg.predict(x1s)
    plt.plot(x1s, y_pred, "k-", linewidth=2, label="$\hat{y}$")
    plt.plot(x1s, y_pred + svm_reg.epsilon, "k--")
    plt.plot(x1s, y_pred - svm_reg.epsilon, "k--")
    plt.scatter(X[svm_reg.support_], y[svm_reg.support_], s=180, facecolors="#FFAAAA")
    plt.plot(X, y, "bo")
    plt.xlabel("$x_1$", fontsize=18)
    plt.legend(loc="upper left", fontsize=18)
    plt.axis(axes)

plt.figure(figsize=(20, 8))
plt.suptitle("图10:SVM回归", fontproperties=my_font, fontsize=18)

plt.subplot(121)
plot_svm_regression(svm_reg1, X, y, [0, 2, 3, 11])
plt.title("$\epsilon={}$".format(svm_reg1.epsilon), fontsize=18)
plt.ylabel("$y$", fontsize=18)
plt.annotate("",
             xy=(eps_x1, eps_y_pred),
             xycoords="data",
             xytext=(eps_x1, eps_y_pred-svm_reg1.epsilon),
             arrowprops={"arrowstyle":"<->", "linewidth":1.5}
            )
plt.text(0.91, 5.6, "$\epsilon$", fontsize=20)

plt.subplot(122)
plot_svm_regression(svm_reg2, X, y, [0, 2, 3, 11])
plt.title("$\epsilon={}$".format(svm_reg2.epsilon), fontsize=18)

plt.savefig("./plot_data/SVM/svm_regression_plot.png")
plt.show()

第五章 支持向量机(SVM)

在间隔内添加更多的实例不会影响模型的预测,所以这个模型被称为ε不敏感

我们可以使用Scikit-Learn的LinearSVR类来执行线性SVM回归。以下代码生成如图10左图所示的模型(训练数据需要先缩放并集中):

np.random.seed(42)
m = 100
X = 2 * np.random.rand(m ,1) - 1
y = (0.2 + 0.1 * X + 0.5 * X **2 + np.random.randn(m, 1)/10).ravel()

警告:在版本0.22中,gamma的默认值将从“auto”更改为“scale”,以更好地考虑未缩放的功能。 为了保留与书中相同的结果,我们明确地将其设置为“auto”,但您可能只需在自己的代码中使用默认值。

要解决非线性回归任务,可以使用核化的SVM模型。例如,图11显示了在一个随机二次训练集上,使用二阶多项式核的SVM回归。左图几乎没有正则化(C值很大),右图则过度正则化(C值很小)。

使用Scikit-Learn的SVR类(支持核技巧)生成如图11左图所示的模型。SVR类是SVC类的回归等价物,LinearSVR类也是LinearSVC类的回归等价物。LinearSVR与训练集的大小线性相关(跟LinearSVC一样),而SVR则在训练集变大时,变得很慢(SVC也是一样)。

from sklearn.svm import SVR

svm_poly_reg = SVR(kernel="poly", degree=2, C=100, epsilon=0.1, gamma="auto")
svm_poly_reg.fit(X, y)
SVR(C=100, cache_size=200, coef0=0.0, degree=2, epsilon=0.1, gamma='auto',
  kernel='poly', max_iter=-1, shrinking=True, tol=0.001, verbose=False)
from sklearn.svm import SVR

svm_poly_reg1 = SVR(kernel="poly", degree=2, C=100, epsilon=0.1, gamma="auto")
svm_poly_reg2 = SVR(kernel="poly", degree=2, C=0.01, epsilon=0.1, gamma="auto")
svm_poly_reg1.fit(X, y)
svm_poly_reg2.fit(X, y)

plt.figure(figsize=(20, 8))
plt.suptitle("图11:使用二阶多项式核的SVM回归", fontproperties=my_font, fontsize=18)

plt.subplot(121)
plot_svm_regression(svm_poly_reg1, X, y, [-1, 1, 0, 1])
plt.title("$degree={}, C={}, \epsilon={}$".format(svm_poly_reg1.degree, svm_poly_reg1.C, svm_poly_reg1.epsilon))
plt.ylabel("$y$", fontsize=18)

plt.subplot(122)
plot_svm_regression(svm_poly_reg2, X, y, [-1, 1, 0, 1])
plt.title("$degree={}, C={}, \epsilon={}$".format(svm_poly_reg2.degree, svm_poly_reg2.C, svm_poly_reg2.epsilon))
plt.ylabel("$y$", fontsize=18)
plt.show()

第五章 支持向量机(SVM)

SVM也可用于异常值检测:详细信息请参考Scikit-Learn文档。

工作原理

在处理SVM时它更为方便(也更常见):偏置项表示为b,特征权重向量表示为w,同时输入特征向量中不添加偏置特征。

决策函数和预测

线性SVM分类器通过简单地计算决策函数WTX+b=w1x1++wnxn+bW^{T}·X+b=w_{1}x_{1}+…+w_{n}x_{n}+bWT⋅X+b=w1​x1​+…+wn​xn​+b来预测新实例x的分类。如果结果为正,则预测类别是正类(1),不然则预测其为负类(0),见x下面公式。
公式2:线性SVM分类器预测

y^={0ifWTX+b&lt;01ifWTX+b0 \hat{y}= \left\{ \begin{aligned} 0 &amp;&amp; if\quad\quad W^{T}·X + b &lt; 0\\ 1 &amp;&amp; if\quad\quad W^{T}·X + b \geq 0 \end{aligned} \right. y^​={01​​ifWT⋅X+b<0ifWT⋅X+b≥0​

图12显示了图4右侧的模型所对应的决策函数:数据集包含两个特征(花瓣宽度和长度),所以是一个二维平面。决策边界是决策函数等于0的点的集合:它是两个平面的交集,也就是一条直线(加粗实线所示)

iris = datasets.load_iris()
X = iris["data"][:, (2, 3)]  # petal length, petal width
y = (iris["target"]==2).astype(np.float64)  # Iris-Virginica
%matplotlib inline
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def plot_3D_decision_function(ax, w, b, x1_lim=[4, 6], x2_lim=[0.8, 2.8]):
    x1_in_bounds = (X[:, 0] > x1_lim[0]) & (X[:, 0] < x1_lim[1])
    X_crop = X[x1_in_bounds]
    y_crop = y[x1_in_bounds]
    x1s = np.linspace(x1_lim[0], x1_lim[1], 20)
    x2s = np.linspace(x2_lim[0], x2_lim[1], 20)
    x1, x2 = np.meshgrid(x1s, x2s)
    xs = np.c_[x1.ravel(), x2.ravel()]
    df = (xs.dot(w) + b).reshape(x1.shape)
    m = 1 / np.linalg.norm(w)
    
    boundary_x2s = -x1s*(w[0]/w[1])-b/w[1]
    margin_x2s_1 = -x1s*(w[0]/w[1])-(b-1)/w[1]
    margin_x2s_2 = -x1s*(w[0]/w[1])-(b+1)/w[1]
    ax.plot_surface(x1s, x2, np.zeros_like(x1), color="b", alpha=0.2, cstride=100, rstride=100)
    ax.plot(x1s, boundary_x2s, 0, "k-", linewidth=2, label=r"$h=0$")
    ax.plot(x1s, margin_x2s_1, 0, "k--", linewidth=2, label=r"$h=\pm 1$")
    ax.plot(x1s, margin_x2s_2, 0, "k--", linewidth=2)
    ax.plot(X_crop[:, 0][y_crop==1], X_crop[:, 1][y_crop==1], 0, "g^")
    ax.plot_wireframe(x1, x2, df, alpha=0.3, color="k")  # 绘制网格的图形
    ax.plot(X_crop[:, 0][y_crop==0], X_crop[:, 1][y_crop==0], 0, "bs")
    ax.axis(x1_lim + x2_lim)
    ax.text(4.5, 2.5, 3.8, "Decision function $h$", fontsize=15)
    ax.set_xlabel(r"Petal length", fontsize=15)
    ax.set_ylabel(r"Petal width", fontsize=15)
    ax.set_zlabel(r"$h = \mathbf{w}^T \mathbf{x} + b$", fontsize=18)
    ax.legend(loc="upper left", fontsize=16)

fig = plt.figure(figsize=(20, 8))
ax1 = fig.add_subplot(111, projection="3d")
plot_3D_decision_function(ax1, w=svm_clf2.coef_[0], b=svm_clf2.intercept_[0])

plt.title("图12:鸢尾花数据集的决策函数", fontproperties=my_font, fontsize=18)
plt.show()

第五章 支持向量机(SVM)

虚线表示决策函数等于1或-1的点:它们互相平行,并且与决策边界的距离相等,从而形成了一个间隔。训练线性SVM分类器即意味着找到w和b的值,从而使这个间隔尽可能宽的同时,避免(硬间隔)或是限制(软间隔)间隔违例。

更概括地说,当有n个特征时,决策函数是一个n维的超平面,决策边界是一个(n-1)维的超平面。

训练目标

决策函数的斜率:它等于权重向量的范数,即||w||。如果我们将斜率除以2,那么决策函数等于±1的点也将变得离决策函数两倍远。也就是说,将斜率除以2,将会使间隔乘以2。也许2D图更容易将其可视化,见图13。权重向量w越小,间隔越大。

def plot_2D_decision_function(w, b, ylabel=True, x1_lim=[-3, 3]):
    x1 = np.linspace(x1_lim[0], x1_lim[1], 200)
    y = w * x1 + b
    m = 1/w
    
    plt.plot(x1,y)
    plt.plot(x1_lim, [1, 1], "k:")
    plt.plot(x1_lim, [-1, -1], "k:")
    plt.axhline(y=0, color="k")
    plt.axvline(x=0, color="k")
    plt.plot([m, m], [0, 1], "k--")
    plt.plot([-m, -m], [0, -1], "k--")
    plt.plot([-m, m], [0, 0], "k-o", linewidth=3)
    plt.axis(x1_lim + [-2, 2])
    plt.xlabel("$x_1$", fontsize=16)
    if ylabel:
        plt.ylabel("$w_1 x_1$", fontsize=16)
    plt.title("$w_1={}$".format(w), fontsize=16)
    
plt.figure(figsize=(20, 8))
plt.suptitle("图13:权重向量大小,间隔越大", fontproperties=my_font, fontsize=18)

plt.subplot(121)
plot_2D_decision_function(1, 0)
plt.subplot(122)
plot_2D_decision_function(0.5, 0, ylabel=False)
plt.savefig("./plot_data/SVM/small_w_large_margin_plot.png")
plt.show()

第五章 支持向量机(SVM)

from sklearn.svm import SVC
from sklearn import datasets

iris = datasets.load_iris()
X = iris["data"][:, (2,3)]  # petal length, petal width
y = (iris["target"]==2).astype(np.float64)  # Iris-Virginia

svm_clf = SVC(kernel="linear", C=1)
svm_clf.fit(X, y)
svm_clf.predict([[5.3, 1.3]])
array([1.])

所以我们要最小化||w||来得到尽可能大的间隔。但是,如果我们想避免任何间隔违例(硬间隔),那么就要使所有正类训练集的决策函数大于1,负类训练集的决策函数小于-1。如果我们定义,实例为负类(如果y(i)=0)时,t(i)=-1;实例为正类(如果y(i)=1)时,t(i)=1。那么我们就可以将这个约束条件表示为:对所有实例来说,t(i)(wT·x(i)+b)≥1。

因此,我们可以将硬间隔线性SVM分类器的目标,看作一个约束优化问题,如公式3所示。
公式3:硬间隔线性SVM分类器的目标
12WTW使t(i)(WTX(i)+b)1(i=1,2,...,m) 最小化\frac{1}{2}W^{T}·W\\ 使得t^{(i)}(W^{T}·X^{(i)}+b) \geq 1 (i=1, 2,..., m) 最小化21​WT⋅W使得t(i)(WT⋅X(i)+b)≥1(i=1,2,...,m)

我们最小化的是12WTW\frac{1}{2}W^{T}·W21​WT⋅W,它等于12W2\frac{1}{2}||W||^221​∣∣W∣∣2,而不是最小化W||W||∣∣W∣∣。
这是因为,二者虽然会得到同样的结果(因为让某个值最小的w和b,同样也使其平方的一半最小),但是12W2\frac{1}{2}||W||^221​∣∣W∣∣2有一个简单好用的导数(就是W),而W||W||∣∣W∣∣在W=0W=0W=0时,是不可微的。优化算法在可微函数上的工作效果要好得多。

要达到软间隔的目标,我们需要为每个实例引入一个松弛变量ζ(i)0ζ^{(i)}≥0ζ(i)≥0,ζ(i)ζ^{(i)}ζ(i)衡量的是第i个实例多大程度上允许间隔违例。那么现在我们有了两个互相冲突的目标:使松弛变量越小越好从而减少间隔违例,同时还要使wT·w/2最小化以增大间隔。这正是超参数C的用武之地:允许我们在两个目标之间权衡。公式5-4给出了这个约束优化问题。

公式4:软间隔线性SVM分类器目标
12WTW+Ci=1mζ(i)使t(i)(WTX(i)+b)1ζ(i)ζ(i)0(i=1,2,...,m) 最小化\frac{1}{2}W^{T}·W + C\sum_{i=1}^{m}\zeta^{(i)}\\ 使得t^{(i)}(W^{T}·X^{(i)}+b) \geq 1 - \zeta^{(i)}且 \zeta^{(i)} \geq 0 (i=1, 2,..., m) 最小化21​WT⋅W+Ci=1∑m​ζ(i)使得t(i)(WT⋅X(i)+b)≥1−ζ(i)且ζ(i)≥0(i=1,2,...,m)

二次规划

硬间隔和软间隔问题都属于线性约束的凸二次优化问题。这类问题被称为二次规划(QP)问题。要解决二次规划问题有很多现成的求解器,使用到的技术各不相同,这些不在本书的讨论范围之内。公式5给出的是问题的一般形式。

公式5-5:二次规划问题

12PTHP+fTP使APb{Pnp(np)HnpnpfnpAncnp(nc)bnp 最小化\frac{1}{2}P^{T}·H·P + f^{T}·P\\ 使得 A·P \leq b\\ 其中\left\{ \begin{aligned} &amp;P是一个n_{p}维向量(n_{p}为参数数量)\\ &amp;H是一个n_{p}*n_{p}矩阵\\ &amp;f是一个n_{p}维向量\\ &amp;A是一个n_{c}*n_{p}矩阵(n_{c}为约束数量)\\ &amp;b是一个n_{p}维向量 \end{aligned} \right. 最小化21​PT⋅H⋅P+fT⋅P使得A⋅P≤b其中⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧​​P是一个np​维向量(np​为参数数量)H是一个np​∗np​矩阵f是一个np​维向量A是一个nc​∗np​矩阵(nc​为约束数量)b是一个np​维向量​

注意表达式APbA·P \leq bA⋅P≤b实际上定义了ncn_cnc​个约束:对于i=12ncPTa(i)b(i)i=1,2,…,n_c,P^T·a^{(i)}≤b^{(i)}i=1,2,…,nc​,PT⋅a(i)≤b(i),其中a(i)a^{(i)}a(i)是包含A的第i行元素的向量,而b^{(i)}是b的第i个元素。

可以轻松验证一下,如果你把二次规划参数按以下方式设置,是否能够实现硬间隔线性SVM分类器的目标:

·np=n+1n_p=n+1np​=n+1,其中n为特征数量(+1是偏置项)。

·nc=mn_c=mnc​=m,其中m是训练实例的数量。

·H是np×npn_p×n_pnp​×np​的单位矩阵,但是顶左单元格为零(为了忽略偏置项)。

·f=0,一个全是0的n_p维向量。

·b=1,一个全是1的n_c维向量。

·a(i)=t(i)x^(i)a^{(i)}=-t^{(i)}\hat{x}^{(i)}a(i)=−t(i)x^(i),其中x^(i)\hat{x}^{(i)}x^(i)等于x(i)x^{(i)}x(i),除了一个额外的偏置特征。

所以,要训练硬间隔线性SVM分类器,有一种办法是直接将上面的参数用在一个现成的二次规划求解器上。得到的向量p将会包括偏置项b=p0b=p_0b=p0​,以及特征权重wi=pii=12mw_i=p_i,i=1,2,…,mwi​=pi​,i=1,2,…,m。类似地,你也可以用二次规划求解器来解决软间隔问题。

但是,为了运用核技巧,接下来我们将要看一个不同的约束优化问题。

对偶问题

针对一个给定的约束优化问题,称之为原始问题,我们常常可以用另一个不同的,但是与之密切相关的问题来表达,这个问题我们称之为对偶问题。通常来说,对偶问题的解只能算是原始问题的解的下限,但是在某些情况下,它也可能跟原始问题的解完全相同。幸运的是,SVM问题刚好就满足这些条件,所以你可以选择是解决原始问题还是对偶问题,二者解相同。公式6给出了线性SVM目标的对偶形式(如果你对如何从原始问题导出对偶问题感兴趣,请参阅附录C)。

公式6:线性SVM目标的对偶形式

12i=1mj=1mα(i)α(j)t(i)t(j)x(i)Tx(j)=i=1mα(i)使α(i)0(i=1,2,...,m) 最小化\frac{1}{2}\sum_{i=1}^{m}\sum_{j=1}^{m}\alpha^{(i)}\alpha^{(j)}t^{(i)}t^{(j)}x^{(i)T}x^{(j)} = \sum_{i=1}^{m}\alpha^{(i)}\\ 使得\alpha^{(i)}\geq 0 (i=1, 2, ..., m) 最小化21​i=1∑m​j=1∑m​α(i)α(j)t(i)t(j)x(i)Tx(j)=i=1∑m​α(i)使得α(i)≥0(i=1,2,...,m)

一旦得到使得该等式最小化(使用二次规划求解器)的向量,就可以使用公式7来计算使原始问题最小化的和。

公式7:从对偶问题到原始问题

w^=i=1mα(i)t(i)x(i)b^=1nsi=1,α^(i)&gt;0m(1t(i)(w^Tx(i))) \hat{w} = \sum_{i=1}^{m}\alpha^{(i)}t^{(i)}x^{(i)}\\ \hat{b} = \frac{1}{n_s}\sum_{i=1, \hat{\alpha}^{(i)&gt;0}}^{m}(1- t^{(i)}(\hat{w}^{T}·x^{(i)})) w^=i=1∑m​α(i)t(i)x(i)b^=ns​1​i=1,α^(i)>0∑m​(1−t(i)(w^T⋅x(i)))

当训练实例的数量小于特征数量时,解决对偶问题比原始问题更快速。更重要的是,它能够实现核技巧,而原始问题不可能实现。

核化SVM

假设你想要将一个二阶多项式转换为一个二维训练集(例如卫星训练集),然后在转换训练集上训练线性SVM分类器。这个二阶多项式的映射函数如公式8所示。

公式8:二阶多项式映射

ϕ(x)=ϕ((x1x2))=(x122x1x2x22).\phi(x) = \phi\left(\left( \begin{aligned} x_1\\ x_2 \end{aligned} \right)\right) = \left(\begin{aligned} x_{1}^{2}\\ \sqrt{2}x_{1}x_{2}\\ x_{2}^{2} \end{aligned} \right). ϕ(x)=ϕ((x1​x2​​))=⎝⎜⎜⎛​x12​2​x1​x2​x22​​⎠⎟⎟⎞​.
注意转换后的向量是三维的而不是二维的。现在我们来看看,如果我们应用这个二阶多项式映射,两个二维向量a和b会发生什么变化,然后计算转换后两个向量的点积(参见公式9)。

公式9:二阶多项式映射的核技巧
ϕ(a)Tϕb=(a122a1a2a22)T(b122b1b2b22)=a12b12+2a1b1a2b2+a22b22=(a1b1+a2b2)2=((a1a2)T(b1b2))2=(aTb)2.\phi(a)^{T}·\phi{b} = \left(\begin{aligned} a_{1}^{2}\\ \sqrt{2}a_{1}a_{2}\\ a_{2}^{2} \end{aligned} \right)^{T}·\left(\begin{aligned} b_{1}^{2}\\ \sqrt{2}b_{1}b_{2}\\ b_{2}^{2} \end{aligned} \right) = a_{1}^{2}b_{1}^{2}+2a_{1}b_{1}a_{2}b_{2}+a_{2}^{2}b_{2}^{2}\\ =(a_{1}b_{1}+a_{2}b_{2})^{2} = \left(\left( \begin{aligned} a_1\\ a_2 \end{aligned} \right)^{T}·\left( \begin{aligned} b_1\\ b_2 \end{aligned} \right)\right)^{2} =(a^{T}·b)^{2}. ϕ(a)T⋅ϕb=⎝⎜⎜⎛​a12​2​a1​a2​a22​​⎠⎟⎟⎞​T⋅⎝⎜⎜⎛​b12​2​b1​b2​b22​​⎠⎟⎟⎞​=a12​b12​+2a1​b1​a2​b2​+a22​b22​=(a1​b1​+a2​b2​)2=⎝⎛​(a1​a2​​)T⋅(b1​b2​​)⎠⎞​2=(aT⋅b)2.

怎么样?转换后向量的点积等于原始向量的点积的平方:ϕ(a)Tϕ(b)=(aTb)2\phi(a)^{T}·\phi(b) = (a^{T}·b)^{2}ϕ(a)T⋅ϕ(b)=(aT⋅b)2.

关键点来了:如果将转换映射应用于所有训练实例,那么对偶问题(公式6)将包含点积ϕ(x(i))Tϕ(x(j))\phi(x^{(i)})^T·\phi(x^{(j)})ϕ(x(i))T⋅ϕ(x(j))的计算。如果ϕ\phiϕ是公式8所定义的二阶多项式转换,那么可以直接用KaTeX parse error: Double superscript at position 9: (x^{(i)}^̲{T}·x^{(j)})^2来替代这个转换向量的点积。所以你根本不需要转换训练实例,只需将公式6里的点积换成点积的平方即可。如果你不嫌麻烦,可以动手将训练集进行转换,然后拟合线性SVM算法,你会发现,结果一模一样。但是这个技巧大大提高了整个过程的计算效率。这就是核技巧的本质。

函数K(ab)=(aTb)2K(a,b)=(a^T·b)^2K(a,b)=(aT⋅b)2被称为二阶多项式核。在机器学习里,核是能够仅基于原始向量aaa和bbb来计算点积\phi(a)^T·\phi(b)$的函数,它不需要计算(甚至不需要知道)转换函数\phi。公式10列出了一些最常用的核函数。

公式10:常用核函数
线K(a,b)=aTbK(a,b)=(γaTb+r)dK(a,b)=exp(γab2)SigmoidK(a,b)=tanh(γaTb+r) 线性核函数:K(a,b) = a^T· b\\ 多项式核函数:K(a,b) = (\gamma a^T·b + r)^d\\ 高斯核函数:K(a,b) = exp(-\gamma||a-b||^2)\\ Sigmoid核函数:K(a,b) = tanh(\gamma a^T·b + r) 线性核函数:K(a,b)=aT⋅b多项式核函数:K(a,b)=(γaT⋅b+r)d高斯核函数:K(a,b)=exp(−γ∣∣a−b∣∣2)Sigmoid核函数:K(a,b)=tanh(γaT⋅b+r)


Mercer定理

根据Mercer定理,如果函数K(a,b)K(a,b)K(a,b)符合几个数学条件——也就是MercerMercerMercer条件(K必须是连续的,并且在其参数上对称,所以K(a,b)=K(b,a),等等),则存在函数将a和b映射到另一空(可能是更高维度的空间),使得K(ab)ϕ(a)Tϕ(b)K(a,b)= \phi(a)^T·\phi(b)K(a,b)=ϕ(a)T⋅ϕ(b)。所以你可以将K用作核函数,因为你知道是存在的,即使你不知道它是什么。对于高斯RBF核函数,可以看出, 实际上将每个训练实例映射到了一个无限维空间,幸好不用执行这个映射。

注意,也有一些常用的核函数(如Sigmoid核函数)不符合Mercer条件的所有条件,但是它们在实践中通常也表现不错。


还有一个未了结的问题我们需要说明。公式7显示了用线性SVM分类器如何从对偶解走到原始解,但是如果你应用了核技巧,最终得到的是包含ϕ(x(i))\phi(x^{(i)})ϕ(x(i))的方程。而w^\hat{w}w^的维度数量必须与ϕ(x(i))\phi(x^{(i)})ϕ(x(i))相同,后者很有可能是巨大甚至是无穷大的,所以你根本没法计算。可是不知道w^\hat{w}w^该如何做出预测呢?你可以将公式7中w^\hat{w}w^的公式插入新实例x(n)x^{(n)}x(n)的决策函数中,这样就得到了一个只包含输入向量之间点积的公式。这时你就可以再次运用核技巧了(见公式11)。

公式11:使用核化SVM做出预测

hW^,b^(ϕ(X(n)))=W^Tϕ(X(n))+b^=(i=1mα^(i)t(i)ϕ(X(i)))Tϕ(X(n))+b^=i=1mα^(i)t(i)(ϕ(X(i))Tϕ(X(n)))+b^=i=1,α^(i)&gt;0mα^(i)t(i)K(X(i),X(n))+b^ h_{\hat{W},\hat{b}}(\phi(X^{(n)})) = \hat{W}^{T}·\phi(X^{(n)}) + \hat{b} = (\sum_{i=1}^{m}\hat{\alpha}^{(i)}t^{(i)}\phi(X^{(i)}))^{T}·\phi(X^{(n)}) + \hat{b}\\ =\sum_{i=1}^{m}\hat{\alpha}^{(i)}t^{(i)}(\phi(X^{(i)})^{T}·\phi(X^{(n)}))+ \hat{b}\\ =\sum_{i=1,\hat{\alpha}^{(i)}&gt;0}^{m}\hat{\alpha}^{(i)}t^{(i)}K(X^{(i)}, X^{(n)})+ \hat{b} hW^,b^​(ϕ(X(n)))=W^T⋅ϕ(X(n))+b^=(i=1∑m​α^(i)t(i)ϕ(X(i)))T⋅ϕ(X(n))+b^=i=1∑m​α^(i)t(i)(ϕ(X(i))T⋅ϕ(X(n)))+b^=i=1,α^(i)>0∑m​α^(i)t(i)K(X(i),X(n))+b^

注意,因为仅对于支持向量才有α(i)0\alpha^{(i)}\neq 0α(i)̸​=0,所以预测时,计算新输入向量x(n)x^{(n)}x(n)的点积,使用的仅仅是支持向量而不是全部训练实例。当然,你还需要使用同样的技巧来计算偏置项(见公式12)。

公式12:使用核技巧计算偏置项
b^=1nsi=1,α^(i)m(1t(i)W^Tϕ(X(i)))=1nsi=1,α^(i)m(1t(i)(j=1m(^α)jtjϕ(X(j)))Tϕ(X(i)))=1nsi=1,α^(i)m(1t(i)j=1m(^α)jtjK(X(i),X(j))) \hat{b} = \frac{1}{n_s}\sum_{i=1, \hat{\alpha}^{(i)}}^{m}(1-t^{(i)}\hat{W}^T·\phi(X^{(i)})) = \frac{1}{n_s}\sum_{i=1, \hat{\alpha}^{(i)}}^{m}(1-t^{(i)}(\sum_{j=1}^{m}\hat(\alpha)^{j}t^{j}\phi(X^{(j)}))^T·\phi(X^{(i)}))\\ = \frac{1}{n_s}\sum_{i=1, \hat{\alpha}^{(i)}}^{m}(1-t^{(i)}\sum_{j=1}^{m}\hat(\alpha)^{j}t^{j}K(X^{(i)},X^{(j)})) b^=ns​1​i=1,α^(i)∑m​(1−t(i)W^T⋅ϕ(X(i)))=ns​1​i=1,α^(i)∑m​(1−t(i)(j=1∑m​(^​α)jtjϕ(X(j)))T⋅ϕ(X(i)))=ns​1​i=1,α^(i)∑m​(1−t(i)j=1∑m​(^​α)jtjK(X(i),X(j)))

如果你现在开始觉得头痛,完全正常:这正是核技巧的副作用。

损失函数

hinge损失函数(max(0, 1-t))

当t≥1时,函数等于0。如果t<1,其导数(斜率)等于-1,如果t>1,则导数(斜率)为0,t=1,时,函数不可导。但是,在t=0处可以使用任意次导数(即-1到0之间的任意值)

t = np.linspace(-2, 4, 200)
h = np.where(1-t<0, 0, 1-t)  # max(0, 1-t)

plt.figure(figsize=(20, 8))
plt.plot(t, h, "b-", linewidth=2, label="$max(0, 1-t)$")
plt.grid(True, which="both")
plt.axhline(y=0, color="k")
plt.axvline(x=0, color="k")
plt.yticks(np.arange(-1, 2.5, 1))
plt.xlabel("$t$", fontsize=16)
plt.axis([-2, 4, -1, 2.5])
plt.legend(loc="upper right", fontsize=16)
plt.savefig("./plot_data/SVM/hinge_plot.png")
plt.show()

第五章 支持向量机(SVM)

还可以使用梯度下降

梯度下降的线性SVM分类器实现

对线性SVM分类器来说,方法之一是使用梯度下降,使从原始问题导出的成本函数(见公式13)最小化。但不幸的是,这种方法的收敛速度比二次规划方法要慢得多。

公式13:线性SVM分类器成本函数
J(W,b)=12WTW+Ci=1mmax(0,1t(i)(WTX(i)+b)) J(W,b) = \frac{1}{2}W^{T}\cdot{W} + C\sum_{i=1}^{m}max(0, 1-t^{(i)}(W^{T}\cdot{X^{(i) + b}})) J(W,b)=21​WT⋅W+Ci=1∑m​max(0,1−t(i)(WT⋅X(i)+b))
成本函数中的第一项会推动模型得到一个较小的权重向量w,从而使间隔更大。第二项则计算全部的间隔违例。如果没有一个示例位于街道之上,并且都在街道正确的一边,那么这个实例的间隔违例为0;如不然,则该实例的违例大小与其到街道正确一边的距离成正比。所以将这个项最小化,能够保证模型使间隔违例尽可能小,也尽可能少。

# Training set
X = iris["data"][:, (2, 3)] # petal length, petal width
y = (iris["target"] == 2).astype(np.float64).reshape(-1, 1) # Iris-Virginica

自己手动编译的梯度下降方法

from sklearn.base import BaseEstimator

class MyLinearSVC(BaseEstimator):
    def __init__(self, C=1, eta0=1, eta_d=10000, n_epochs=1000, random_state=None):
        self.C = C
        self.eta0 = eta0
        self.n_epochs = n_epochs
        self.random_state = random_state
        self.eta_d = eta_d

    def eta(self, epoch):
        return self.eta0 / (epoch + self.eta_d)
        
    def fit(self, X, y):
        # Random initialization
        if self.random_state:
            np.random.seed(self.random_state)
        w = np.random.randn(X.shape[1], 1) # n feature weights
        b = 0

        m = len(X)
        t = y * 2 - 1  # -1 if t==0, +1 if t==1
        X_t = X * t
        self.Js=[]

        # Training
        for epoch in range(self.n_epochs):
            support_vectors_idx = (X_t.dot(w) + t * b < 1).ravel()
            X_t_sv = X_t[support_vectors_idx]
            t_sv = t[support_vectors_idx]

            J = 1/2 * np.sum(w * w) + self.C * (np.sum(1 - X_t_sv.dot(w)) - b * np.sum(t_sv))
            self.Js.append(J)

            w_gradient_vector = w - self.C * np.sum(X_t_sv, axis=0).reshape(-1, 1)
            b_derivative = -C * np.sum(t_sv)
                
            w = w - self.eta(epoch) * w_gradient_vector
            b = b - self.eta(epoch) * b_derivative
            

        self.intercept_ = np.array([b])
        self.coef_ = np.array([w])
        support_vectors_idx = (X_t.dot(w) + t * b < 1).ravel()
        self.support_vectors_ = X[support_vectors_idx]
        return self

    def decision_function(self, X):
        return X.dot(self.coef_[0]) + self.intercept_[0]

    def predict(self, X):
        return (self.decision_function(X) >= 0).astype(np.float64)

C=2
svm_clf = MyLinearSVC(C=C, eta0 = 10, eta_d = 1000, n_epochs=6000, random_state=2)
svm_clf.fit(X, y)
svm_clf.predict(np.array([[5, 2], [4, 1]]))
array([[1.],
       [0.]])
plt.plot(range(svm_clf.n_epochs), svm_clf.Js)
# plt.axis([0, svm_clf.n_epochs, 0, 1000])
[<matplotlib.lines.Line2D at 0x203d5deb550>]

第五章 支持向量机(SVM)

print(svm_clf.intercept_, svm_clf.coef_)
[-22.69442785] [[[3.71213881]
  [2.79939099]]]
svm_clf2 = SVC(kernel="linear", C=C)
svm_clf2.fit(X, y.ravel())
print(svm_clf2.intercept_, svm_clf2.coef_)
[-15.51721253] [[2.27128546 2.71287145]]
yr = y.ravel()
plt.figure(figsize=(12,3.2))
plt.subplot(121)
plt.plot(X[:, 0][yr==1], X[:, 1][yr==1], "g^", label="Iris-Virginica")
plt.plot(X[:, 0][yr==0], X[:, 1][yr==0], "bs", label="Not Iris-Virginica")
plot_svc_decision_boundary(svm_clf, 4, 6)
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.title("MyLinearSVC", fontsize=14)
plt.axis([4, 6, 0.8, 2.8])

plt.subplot(122)
plt.plot(X[:, 0][yr==1], X[:, 1][yr==1], "g^")
plt.plot(X[:, 0][yr==0], X[:, 1][yr==0], "bs")
plot_svc_decision_boundary(svm_clf2, 4, 6)
plt.xlabel("Petal length", fontsize=14)
ax1 = plt.twinx()
plt.title("SVC", fontsize=14)
plt.axis([4, 6, 0.8, 2.8])
[4, 6, 0.8, 2.8]

第五章 支持向量机(SVM)

使用scikit-learn自带的 梯度下降的方法

from sklearn.linear_model import SGDClassifier

sgd_clf = SGDClassifier(loss="hinge", alpha = 0.017, max_iter = 50, tol=-np.infty, random_state=42)
sgd_clf.fit(X, y.ravel())

m = len(X)
t = y * 2 - 1  # -1 if t==0, +1 if t==1
X_b = np.c_[np.ones((m, 1)), X]  # Add bias input x0=1
X_b_t = X_b * t
sgd_theta = np.r_[sgd_clf.intercept_[0], sgd_clf.coef_[0]]
print(sgd_theta)
support_vectors_idx = (X_b_t.dot(sgd_theta) < 1).ravel()
sgd_clf.support_vectors_ = X[support_vectors_idx]
sgd_clf.C = C

plt.figure(figsize=(5.5,3.2))
plt.plot(X[:, 0][yr==1], X[:, 1][yr==1], "g^")
plt.plot(X[:, 0][yr==0], X[:, 1][yr==0], "bs")
plot_svc_decision_boundary(sgd_clf, 4, 6)
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.title("SGDClassifier", fontsize=14)
plt.axis([4, 6, 0.8, 2.8])
[-14.062485     2.24179316   1.79750198]





[4, 6, 0.8, 2.8]

第五章 支持向量机(SVM)

在线SVM也可以实现核技巧,可参考“Incremental and Decremental SVM Learning”(http://goo.gl/JEqVui ),以及“Fast Kernel Classifiers with Online and Active Learning”(https://goo.gl/hsoUHA )。但是这些是在Matlab和C++上实现的。对于大规模非线性问题,你可能需要使用神经网络.

上一篇:Vue系列四:class与style绑定


下一篇:CocosCreator -- 根据字体样式获取Label宽高