北京大学R语言教程(李东风)第44章：支持向量机

分隔超平面

对因变量为两分类的情形，
设自变量xi,i=1,2,…,n在ℝp空间中，
如果存在曲面将ℝp分隔成两部分，
使得两类的样本可以分开，
就有了一种简单的判别准则。
超平面是最简单的分隔曲面。

在平面空间ℝ2={(x1,x2):x1,x2∈ℝ}中，
超平面就是直线，其方程为

β0+β1x1+β2x2=0

其中(β1,β2)T≠0。
在平面空间ℝ2中，直线将空间分为两个部分，
分别满足β0+β1x1+β2x2>0和β0+β1x1+β2x2<0。

可以用其方程的法向β=(β1,β2)T的指向区分两个部分，
设x̃ =(x̃ 1,x̃ 2)T在直线上，
而点x=(x1,x2)T满足β0+β1x1+β2x2>0，
则利用

β0+β1x̃ 1+β2x̃ 2=β0+β1x1+β2x2>00

两式相减得

(β1,β2)(x1−x̃ 1x2−x̃ 2)>0

用⟨⋅,⋅⟩表示ℝ2中的内积，
则x=(x1,x2)T在满足β0+β1x1+β2x2>0的半空间的充分必要条件是

⟨β,x−x̃ ⟩>0

即从平面上的点x̃ 出发到x的向量x−x̃ 与分隔线的法线方向β=(β1,β2)T成锐角；
x̃ =(x̃ 1,x̃ 2)T在满足β0+β1x1+β2x2<0的半空间的充分必要条件是

⟨β,x−x̃ ⟩<0

即即从平面上的点x̃ 出发到x的向量x−x̃ 与分隔线的法线方向β=(β1,β2)T成钝角。

例如，考虑分隔线

1+2x1+3x2=0

法线方向为β=(2,3)T,
上面的一个点x̃ =(1,−1)。
参见图44.1。

对于自变量x∈ℝp的情形，
可以定义分隔超平面：

β0+β1x1+⋯+βpxp=0

将ℝp分隔为满足
β0+β1x1+⋯+βpxp>0的半空间与满足
β0+β1x1+⋯+βpxp<0的半空间，
两个空间的点可以用分隔超平面上的一个点x̃ 到x的矢量
x−x̃ 与超平面的法向量
β=(β1,…,βp)T的夹角是锐角还是钝角来区分，
即内积⟨β,x−x̃ ⟩的正负号来区分。

设两类判别问题的自变量数据保存为矩阵X=(xij)n×p，
其中的每一行为Rn的一个点；
对应于每个点有一个yi∈{−1,1}作为类别标签，
y=(y1,…,yn)T。
设法找到一个超平面将两类的点分开在两边。

例如，
在图44.2中，
要区分蓝色和粉红色的点。
在左图中有三条直线都可以分开两种点。

找到的分隔线必须满足

β0+β1xi1+β2xi2>0⟺yi>0β0+β1xi1+β2xi2<0⟺yi<0

可以统一地写成

yi(β0+β1xi1+β2xi2)>0, i=1,2,…,n.

这种要求可以推广到xi∈ℝp的情形。
即找到超平面使得

yi(β0+β1xi1+⋯+βpxip)>0, i=1,2,…,n.

当分隔超平面存在时，
就可以用β0+β1xi1+⋯+βpxip的正负号来选择y=1还是−1。
对一个待判别的观测x∗，
计算

f(x∗)=β0+β1x∗i1+⋯+βpx∗ip

就可以判别y的类别，
f(x∗)>0时判为1，
f(x∗)<0时判为−1。
f(x∗)的绝对值大小也有意义，
绝对值越大时，
x∗距离分隔边界越远，
判决越可靠。
如果f(x∗)的绝对值接近于零，
则点x∗靠近分隔边界，
对其判决不够可信。

最大分隔边界超平面

如图44.2的左图所示，
只要分隔超平面能分隔开训练集的样本点，
就不止有一个分隔超平面。
需要找到最优的分隔超平面。
想法是让两类的点都尽可能远离分隔面。
对某个候选的分隔超平面，
计算各个观测样本点到超平面的距离并求出最近距离，
称为margin（分隔边界）。
从所有的分隔超平面中求margin最大的一个
（见图44.3），
用这个最大边界的超平面作为判决准则，
这种方法称为最大分隔边界判别法(maximal margin classifier)。

在图44.3中，
有三个观测点到最优超平面的距离相等，都等于最大分隔边界（margin），
这样的点称为支持向量(support vector)。
这些点改变位置会使得分隔线改变位置，
而其它点轻微改变位置则不会使得分隔线改变位置，
除非其它点到分隔线的距离小于分隔边界了。
这是最大分隔边界判别法以及后面讲到的支持向量机方法的重要性质。

最大分隔边界超平面的构造

设法求最大分隔边界超平面。
为了标准化，
令法向量β长度等于1。
求β0,β1,…,βp为如下优化问题的解：

maxβ0,β1,…,βp分隔边界 s.t.yi(β0+β1xi1+⋯+βpxip)≥分隔边界, i=1,2,…,n.β21+⋯+β2p=1

其中的第二个约束条件对可行的分隔超平面是不需要的约束，
而β长度为1的约束也不是必须的，
找到长度不为1的解后除以其长度即可。
不过，
当∥β∥=1成立时，
yi(β0+β1xi1+β2xi2)正是xi到分隔超平面的垂直距离。
有高效的算法可以快速求解上述优化问题。

不能分隔的情形

如果两类的点混杂在一起，
有可能任何的超平面都不能将两个类完全分开。
见图44.4。
这时，
只能退而求其次，
试图找到超平面使得两类尽可能被分开。

问题

要求分隔超平面将两类完全分开有时不能做到，
即使能做到，
因为最大分隔边界超平面的选取过于依赖于边界处的支持向量，
几个点对结果的影响过大，
结果不够稳健，
也容易造成过度拟合。
见图44.5，
右图中增加一个点后，分隔线偏离了很多，
分隔边界变得很窄而且对原有点的判别效果变差了。
各个点距离边界越远，判别效果越好。

解决这个问题的办法是，
不严格要求能够区分所有点，
而是将大多数点区分好，
结果更具有稳健性，
有少数点可以落在边界区域，
甚至于落入分隔超平面的错误一面。
这种方法称为支持向量判别法，
或者软分隔边界法。
见图44.6，
左图中的1号和8号点进入了分隔区域，
右图中的1号和8号点进入了分隔区域，
11和12号点进入了错误的一面。

支持向量判别法

支持向量判别法仍然求出一个分隔超平面，
对待判的观测点，
按其在超平面的那一面判决其类属。
超平面的求法使用软分隔边界，
表述为如下的优化问题：

maxβ0,β1,…,βp,ϵ1,…,ϵn分隔边界 s.t.yi(β0+β1xi1+⋯+βpxip)≥分隔边界×(1−ϵi), i=1,2,…,n.∑i=1pβ2i=1,ϵi≥0,i=1,…,n; ∑i=1nϵi≤C

其中C是一个非负的调节参数，
C越大，
容许进入边界和错判的点的比例越大。
ϵ1,…,ϵn是松弛变量，
使得对点在分割边界以外的要求略微放宽。

求得上述优化问题的解以后，
对新的观测x∗=(x∗1,…,x∗p)T,
只要计算函数
f(x∗)=β0+β1x∗1+…x∗p，
以f(x∗)的正负号决定将因变量判入1或者−1。

优化问题中的ϵi由第i个点与分隔区域的关系决定。
如果ϵi=0，
则软分隔规则对第i个点不起作用，
第i个点仍然被分隔区域分开，
例如图44.6右图的3, 4, 5, 6，10和2, 7, 9号点。
当0<ϵi<1时，
第i个点可以进入分隔边界规定的分隔区域范围内，
如第1，8号点。
当ϵi>1时，
第i个点可以在超平面的错误一侧，
如第11, 12号点。

调节参数C设置了将严格要求所有点在正确一侧且不能进入分隔区域的要求放松到多大程度。
当C=0时对应于最大分隔边界判别法，
不能有任何点在错误一侧，
也不能进入分隔边界区域。
对C>0，
至多有floor(C)个点可以位于分隔超平面的错误一侧。
增大C的值，
解出的边界宽度也会增大。
见图44.7，
各图形对应的调节参数由大到小，
过大的C使得错判点较多，
结果方差较小，偏差较大；
过小的C使得结果不稳健，
预测方差较大，偏差较小。
所以，
这里的调节参数与其它有监督学习方法中的调节参数类似，
应该在偏差与方差之间折衷，
可以用交叉验证方法求最优调节参数值。
C越小，模型复杂度越高。

上面的优化问题的解有一个重要性质：
最终的解仅仅依赖边界点（即落在分隔区域的边界上的点）以及进入边界区域或者进入错误一面的点，
这些点称为支持向量；
那些被分隔区域完全正确地区分在两边的点不起作用。

从这个性质来看条件参数C的影响，
当C大的时候，
分隔边界很宽，
许多点进入了边界区域或者错误一面，
有许多个支持向量，
最终解由许多个点决定，
这时解的方差较低，
但对应于较简单的模型，所以可能会有较大偏差。
当C小的时候，
只有少数几个支持向量，
模型稳定性差，
有过度拟合危险，
方差可能较大。

支持向量判别法只依赖于少数支持向量这个性质与距离判别法很不同，
结果仅由靠近边界的点决定，
那些离边界很远的点则不起作用；
距离判别法则是由所有点决定。
逻辑回归判别法与支持向量判别法有类似的性质，
也是主要依赖那些边界附近的点。

非线性边界方法

分隔超平面是线性边界，
对于图44.8这样的分隔边界非线性的问题无法解决，
其它的线性判别法也不能解决这样的非线性边界判别问题。

解决非线性问题的常用方法是增加非线性项如二次项、三次项。
比如，将自变量由x1,…,xp增加到

x1,…,xp;x21,…,x2p

则支持向量判别法的优化问题变成了

maxβ0,β11,…,βp1,β12,…,βp2,ϵ1,…,ϵn分隔边界 s.t.yi(β0+∑j=1pβj1xij+∑j=1pβj2x2ij)≥分隔边界×(1−ϵi), i=1,2,…,n.∑k=12∑j=1pβ2jk=1ϵi≥0,∑i=1nϵi≤C.

这样得到的边界在(x1,…,xp,x21,…,x2p)所在的ℝ2p空间内是线性的超平面，
但是在原来(x1,…,xp)所在的ℝp空间内则是由
q(x)=0决定的曲面，
q(x)是二次多项式函数。

还可以增加高次项和交叉项，
或者考虑其它的非线性变换。
增加非线性项的方法有无数多种，
一一测试是不现实的，
支持向量机方法则给出了一种增加非线性项的一般方法，
对应的边界可以快速求解。

支持向量机

支持向量机利用了Hilbert空间的方法将线性问题扩展为非线性问题。
线性的支持向量判别法可以通过ℝp的内积转化为判别函数如下的等价表示：

f(x)=β0+∑i=1nαi⟨x,xi⟩.

其中β0,α1,…,αn是待定参数。
记

w=∑i=1nαixi,

则

f(x)=β0+⟨x,∑i=1nαixi⟩=β0+⟨x,w⟩.

为了估计参数，
不需要用到各xi的具体值，
而只需要其两两的内积值，
而且在判别函数中只有支持向量对应的αi才非零，
记为支持向量点集，
则线性判别函数为

f(x)=β0+∑i∈αi⟨x,xi⟩.

支持向量机方法将ℝp中的内积⟨x,x′⟩推广为如下的核函数值：

K(x,x′).

存在从ℝp到某个希尔伯特空间H的映射ϕ(x)，
使得

K(x,x′)=⟨ϕ(x),ϕ(x′)⟩H, ∀x,x′∈ℝp.

核函数K(x,x′)是度量两个观测点x,x′的相似程度的函数。
比如，
取

K(x,x′)=∑j=1pxjx′j,

就又回到了线性的支持向量判别法。

利用核代替内积后，
判别法的判别函数变成

f(x)====β0+∑i∈αiK(x,xi)β0+∑i∈αi⟨ϕ(x),ϕ(xi)⟩β0+⟨ϕ(x),∑i∈αiϕ(xi)⟩β0+⟨ϕ(x),wH⟩.

这在原始观测xi的取值空间中是非线性函数，
但在特征空间H中是线性函数。

核有多种取法。
例如，
取

K(x,x′)={1+∑j=1pxjx′j}d

其中d>1为正整数，
称为多项式核，
则结果是多项式边界的判别法，
本质上是对线性的支持向量方法添加了高次项和交叉项。
前一小节就是d=2时的多项式核。
可以看出，
增加了二次项以后，
将ℝp空间的点x=(x1,…,xp)映射到了ℝ2p空间ϕ(x)=(x1,…,xp,x21,…,x2p)，
考虑在高维的ℝ2p空间的分隔超平面，
但在x所在的ℝp这是分隔曲面。
判决函数可以表示为ℝ2p中ϕ(x)=(x1,…,xp,x21,…,x2p)的线性组合，
这又可以表示为内积形式：

f(x)=β0+∑i∈αi⟨ϕ(x),ϕ(xi)⟩=β0+∑i∈αiK(x,xi).

理论研究表明，给定核函数K(⋅,⋅)后，
必存在一个Hilbert空间H以及ϕ:ℝp→H，使得

K(x,x′)=⟨ϕ(x),ϕ(x′)⟩H, ∀x,x∈ℝp.

于是，
给了一个核函数，
就可以将原来的ℝp空间的判别问题，
通过映射ϕ映射为了H空间的线性判别问题，
这样得到的判别函数在H空间是线性的，
但在原始的ℝp中是非线性的。

训练这样的SVM时只要计算观测两两的(n2)个核函数值。
为什么不直接增加非线性项？
原因是计算这些核函数值计算量是确定的，
而增加许多非线性项，
则可能有很大的计算量，
映射ϕ(x)的值域H可能是高维甚至于无穷维空间，
计算核函数值而不是计算映射ϕ(x)效率更高。

支持向量机的理论基于再生核希尔伯特空间(RKHS)，
可参见44.7，
以及(Trevor Hastie 2009)节5.8和节12.3.3。

图44.9的左图是用多项式核解决非线性边界判别问题的演示。
两条粗实线是分隔曲面，
虚线是边界区域的边缘。

另一种常用的核是径向核(radial kernel)，
定义为

K(x,x′)=exp{−γ∑j=1p(xj−x′j)2}

γ为正常数。
当x和x′分别落在以原点为中心的两个超球面上时，
其核函数值不变。
44.9的右图是用径向核支持向量机进行判别的演示。

使用径向核时，
判别函数为

f(x)=β0+∑i∈αiexp{−γ∑j=1p(xj−xij)2}

对一个待判别的观测x∗，
如果x∗距离训练观测点xi较远，
则K(x∗,xi)的值很小，
xi对x∗的判别基本不起作用。
这样的性质使得径向核方法具有很强的局部性，
只有离x∗很近的点才对其判别起作用。

library(rsample)

Heart <- read_csv(
  "data/Heart.csv",
  show_col_types = FALSE) |>
  dplyr::select(-1) |>
  mutate(
    AHD = factor(AHD, levels=c("Yes", "No"))
  )

Heart <- na.omit(Heart)

set.seed(101)
heart_split <- initial_split(
  Heart, prop = 0.50)
heart_train <- training(heart_split)
heart_test <- testing(heart_split)
test.y <- heart_test$AHD

classifier.error <- function(truth, pred){
  tab1 <- table(truth, pred)
  err <- 1 - sum(diag(tab1))/sum(c(tab1))
  err
}

线性的SVM

支持向量判别法就是SVM取多项式核，
阶数d=1的情形。
需要一个调节参数C，
C越大，
分隔边界越窄，
过度拟合危险越大。

先随便取调节参数C=1试验支持向量判别法：

## 
## Attaching package: 'kernlab'

## The following object is masked from 'package:purrr':
## 
##     cross

## The following object is masked from 'package:ggplot2':
## 
##     alpha

res.svc <- ksvm(
  AHD ~ ., 
  data = heart_train, 
  kernel="vanilladot", 
  C = 1, 
  scaled = TRUE)

##  Setting default kernel parameters

fit.svc <- predict(res.svc)
summary(res.svc)

## Length  Class   Mode 
##      1   ksvm     S4

计算拟合结果并计算错判率：

tab1 <- table(
  truth = heart_train$AHD, 
  fitted = fit.svc); tab1

##      fitted
## truth Yes No
##   Yes  69 10
##   No   12 57

cat("SVC错判率：", 
  round(classifier.error(heart_train$AHD, fit.svc), 2), "\n")

## SVC错判率： 0.15

超参数调优方法参见47.3。

多项式核SVM

res.svm1 <- ksvm(
  AHD ~ ., 
  data = heart_train, 
  kernel="polydot", 
  kpar = list(degree = 2),
  C = 0.1, 
  scale=TRUE)
fit.svm1 <- predict(res.svm1)
summary(res.svm1)

## Length  Class   Mode 
##      1   ksvm     S4

tab1 <- table(truth = heart_train$AHD, fitted = fit.svm1); tab1

##      fitted
## truth Yes No
##   Yes  77  2
##   No    2 67

cat("2阶多项式核SVM错判率：", 
  round(classifier.error(heart_train$AHD, fit.svm1), 2), "\n")

## 2阶多项式核SVM错判率： 0.03

注意这是拟合情况。
在测试集上的预测错误率：

pred.svm1t <- predict(res.svm1, heart_test)
tab1 <- table(truth = heart_test$AHD, pred = pred.svm1t); tab1

##      pred
## truth Yes No
##   Yes  46 12
##   No   26 65

cat("测试集上2阶多项式核SVM错判率：", 
  round(classifier.error(heart_test$AHD, pred.svm1t), 2), "\n")

## 测试集上2阶多项式核SVM错判率： 0.26

径向核SVM

径向核需要的参数为σ值。
取参数sigma = 0.1。

res.svm3 <- ksvm(
  AHD ~ ., 
  data = heart_train, 
  kernel="rbfdot", 
  kpar = list(sigma = 0.1),
  C = 0.1, 
  scaled = TRUE)
fit.svm3 <- predict(res.svm3)
summary(res.svm3)

## Length  Class   Mode 
##      1   ksvm     S4

tab1 <- table(truth = heart_train$AHD, 
  fitted = fit.svm3); tab1

##      fitted
## truth Yes No
##   Yes  74  5
##   No   32 37

cat("径向核（sigma=0.1, cost=0.1）SVM拟合错判率：", 
  round(classifier.error(heart_train$AHD, fit.svm3), 2), "\n")

## 径向核（sigma=0.1, cost=0.1）SVM拟合错判率： 0.25

两两判别

设因变量有K个类，为每两个类都分别用SVM构造判别函数，
得到m=(K2)=K(K−1)/2个判别函数。
对于待判样品x∗，
用这些判别函数每一个都判一遍，
得到m个判别结果，
然后这m个结果投票，
那个类别的票最高就判入哪一类。

一对多判别

设因变量有K个类，
对其中的每一个类，
都以此类以及此类之外的其它类作为一个两类判别问题，
建立K个判别函数。
对于待判样品x∗，
计算K个判别函数值，
以判别函数值最大的类作为最后的判别结果。

多分类支持向量机鸢尾花数据计算实例

考虑著名的鸢尾花数据集的判别问题。
所有数据用作训练集合。

取径向基，用指定的调整参数：

res.msvm1 <- ksvm(
  Species ~ ., 
  data = iris, 
  kernel="rbfdot",
  kpar = list(sigma = 0.01),
  C = 1, 
  scaled=TRUE)
summary(res.msvm1)

## Length  Class   Mode 
##      1   ksvm     S4

fit.msvm1 <- predict(res.msvm1)
tab1 <- table(truth=iris[,"Species"], fitted=fit.msvm1); tab1

##             fitted
## truth        setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         46         4
##   virginica       0         11        39

cat("径向核（gamma=0.01, cost=1）多类SVM错判率：", 
    round(1 - (sum(diag(tab1)))/ sum(c(tab1)), 2), "\n")

## 径向核（gamma=0.01, cost=1）多类SVM错判率： 0.1

正定核