Python · SVM(三)· 核方法

Python · SVM(三)· 核方法

(这里是本章会用到的 Jupyter Notebook 地址)

(考试周写专栏真有种忙里偷闲的感觉 _(:з」∠)_)

关于核方法的比较严谨的叙述,个人建议观众老爷们可以看看这个问题下面的前几个回答;在这里的话,我们就还是只注重直观和应用层面

什么是核方法?

往简单里说,核方法是将一个低维的线性不可分的数据映射到一个高维的空间、并期望映射后的数据在高维空间里是线性可分的。我们以异或数据集为例:在二维空间中、异或数据集是线性不可分的;但是通过将其映射到三维空间、我们可以非常简单地让其在三维空间中变得线性可分。比如定义映射:

\phi(x,y)=\left\{
\begin{aligned}
&(x,y,1),\ \ xy>0 \\
&(x,y,0),\ \ xy\le0
\end{aligned}
\right.

该映射的效果如下图所示:

可以看到,虽然左图的数据集线性不可分、但显然右图的数据集是线性可分的,这就是核工作原理的一个不太严谨但仍然合理的解释

从直观上来说,确实容易想象、同一份数据在越高维的空间中越有可能线性可分,但从理论上是否确实如此呢?1965 年提出的 Cover 定理从理论上解决了这个问题,我们会在文末附上相应的公式,这里暂时按下不表

至此,似乎问题就转化为了如何寻找合适的映射\phi、使得数据集在被它映射到高维空间后变得线性可分。不过可以想象的是,现实任务中的数据集要比上文我们拿来举例的异或数据集要复杂得多、直接构造一个恰当的\phi的难度甚至可能高于解决问题本身。而核方法的巧妙之处就在于,它能将构造映射这个过程再次进行转化、从而使得问题变得简易:它通过核函数来避免显式定义映射\phi。往简单里说,核方法会通过用能够表示成K(x_i,x_j)=\phi(x_i)\cdot\phi(x_j)的核函数K(x_i,x_j)替换各算式中出现的内积x_i\cdot x_j来完成将数据从低维映射到高维的过程。换句话说、核方法的思想如下:


  • 将算法表述成样本点内积的组合(这经常能通过算法的对偶形式实现)
  • 设法找到核函数K(x_i,x_j),它能返回样本点x_ix_j\phi作用后的内积
  • K(x_i,x_j)替换x_i\cdot x_j、完成低维到高维的映射(同时也完成了从线性算法到非线性算法的转换)

当然了,不难想象的是,并不是所有的函数K都能够对应一个映射(亦即不是所有的K(x_i,x_j)都能拆成\phi(x_i)\cdot\phi(x_j);比如说,显然K(x_i,x_j)至少需要是一个对称函数)。幸运的是,1909 年提出的 Mercer 定理解决了这个问题,它的具体叙述会在文末给出

Mercer 定理为寻找核函数带来了极大的便利。可以证明如下两族函数都是核函数:


  • 多项式核
    K(x_i,x_j)=(x_i\cdot x_j+1)^p
  • 径向基(Radial Basis Function,常简称为 RBF)核:
    K(x_i,x_j)=\exp\left(-\gamma\|x_i-x_j\|^2\right)

那么核方法的应用场景有哪些呢?在 2002 年由 Scholkopf 和 Smola 证明的表示定理告诉我们它的应用场景非常广泛。定理的具体内容同样会附在文末,这里就暂时按下不表

核模型的表现

还是用 GIF 来说明问题最为形象。当我们对感知机应用核方法后,它就能对非线性数据集(比如螺旋线数据集)进行分类了,训练过程将如下:

怎么应用核方法?

简单来说,就是把算法中涉及到样本(x_i)的地方都通过某种变换、弄成样本的内积形式(x_i\cdot x_j)。以感知机为例,感知机的原始损失函数为L(D) = \sum_{i=1}^N\left[ -y_i(w\cdot x_i+b)\right]_+

为了让损失函数中的样本都变成内积形式,考虑令

w = \sum_{i=1}^N\alpha_ix_i(也有令w = \sum_{i=1}^N\alpha_iy_ix_i的)

\begin{align}
L(D) &= \sum_{i=1}^N\left[ -y_i\left[\left(\sum_{j=1}^N\alpha_jx_j\right)\cdot x_i+b\right]\right]_+ \\
&= \sum_{i=1}^N\left[ -y_i\left(\sum_{j=1}^N\alpha_j(x_i\cdot x_j)+b\right)\right]_+
\end{align}

在此之上应用核方法是平凡的:设核函数为K,只需把所有的x_i\cdot x_j换成K(x_i,x_j)即可:

L(D) = \sum_{i=1}^N\left[ -y_i\left(\sum_{j=1}^N\alpha_jK(x_i,x_j)+b\right)\right]_+

于是优化问题变为

\min_{\alpha}\sum_{i=1}^N\left[ -y_i\left(\sum_{j=1}^N\alpha_jK(x_i,x_j)+b\right)\right]_+

预测步骤则变为

y_{\text{pred}}=w\cdot x+b=\sum_{i=1}^N\alpha_iK(x_i, x)+b

对于 LinearSVM 而言,用同样的手法不难得出其核形式:

L(D)=\frac12\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jK(x_i,x_j)+C\sum_{i=1}^N\left[ 1-y_i\left(\sum_{j=1}^N\alpha_jK(x_i,x_j)+b\right)\right]_+

预测步骤则仍然是

y_{\text{pred}}=w\cdot x+b=\sum_{i=1}^N\alpha_iK(x_i, x)+b

(有没有发现核形式和对偶形式很像?( σ'ω')σ)

如何训练核模型?

【注意:为简洁,从此往后的推导和实现均以核感知机为例,核 SVM 的相关讨论会放在下一章介绍 SMO 算法时进行】


简洁起见,我们还是用梯度下降法来进行训练,为此我们需要进行求导工作。假设当前模型参数为\alpha=(\alpha_1,\alpha_2,...,\alpha_N)^Tx_i在参数\alpha下的预测值为\hat y_i,则:

\frac{\partial L}{\partial\alpha_i}=-\sum_{y_j\hat y_j<0}y_jK(x_j, x_i)
\frac{\partial L}{\partial b}=-\sum_{y_j\hat y_j<0}y_j

为了加速训练,我们需要将该算式向量化,为此我们需要定义核矩阵。假设现在我们有两组样本:(x^{(1)}_1,x^{(2)}_2,...,x^{(1)}_M)^T(x^{(2)}_1, x^{(2)}_2,...,x^{(2)}_N)^T,那么它们的核矩阵即为

\bold K = \left[\begin{matrix}
K(x^{(1)}_1,x^{(2)}_1) & \ldots & K(x^{(1)}_1,x^{(2)}_N) \\
\vdots & \ddots & \vdots \\
K(x^{(1)}_M,x^{(2)}_1) & \ldots & K(x^{(1)}_M,x^{(2)}_N)
\end{matrix}\right]_{N\times N}

对于训练过程而言,我们关心的是训练样本之间的核矩阵

\bold K = \left[\begin{matrix}
K(x_1,x_1) & \ldots & K(x_1,x_N) \\
\vdots & \ddots & \vdots \\
K(x_N,x_1) & \ldots & K(x_N,x_N)
\end{matrix}\right]_{N\times N}

利用它,不难写出相应的向量化代码:

# 假设 k_mat 存储着原样本之间的核矩阵
# 1、计算损失
err = -y * (k_mat.dot(alpha) + b)
# 2、找出使得损失不小于 0 的样本
mask = err >= 0
# 3、进行相应梯度下降,lr 是学习速率
delta = lr * y[mask]
alpha += np.sum(delta[..., None] * k_mat[mask], axis=0)
b += np.sum(delta)

对于预测过程,我们关心的是原样本和新样本之间的核矩阵。假设新样本为(\tilde x_1,...,\tilde x_n)^T,则

\bold K = \left[\begin{matrix}
K(x_1,\tilde x_1) & \ldots & K(x_1,\tilde x_n) \\
\vdots & \ddots & \vdots \\
K(x_N,\tilde x_1) & \ldots & K(x_N,\tilde x_n)
\end{matrix}\right]_{N\times n}

那么预测过程即为

y_{\text{pred}}=\sum_{i=1}^N\alpha_iK(x_i,x)+b=\alpha^T\bold K+b
于是关键就在于如何定义计算核矩阵的核函数了。对于多项式核来说,核函数的实现是直观的:
@staticmethod
def _poly(x, y, p):
    return (x.dot(y.T) + 1) ** p

但对于 RBF 来说就没那么直观了,用到了 Numpy 的高级实用技巧之一——升维:

@staticmethod
def _rbf(x, y, gamma):
    return np.exp(-gamma * np.sum((x[..., None, :] - y) ** 2, axis=2))

当然直接用 for 来实现也是可以的,不过那将会非常非常慢……

核模型的实现

如果思路能够整理清楚,那么核模型相比原模型来说只有如下两点改变:

  • 需要定义核函数并计算出核矩阵
  • 计算预测值时不是w\cdot x+b=w^Tx+b,而是\alpha^T\bold K+b,其中
    • 在训练时,\bold K为原样本之间的核矩阵
    • 在测试时,\bold K为原样本和新样本的核矩阵

所以实现起来的话会有许多重复代码,这里就只展现其中最核心的部分(仍以核感知机为例):

# 训练代码
def fit(...):
    ...
    self._alpha = np.zeros(len(x))
    self._b = 0.
    self._x = x
    # self._kernel 即为核函数,能够计算两组样本的核矩阵
    k_mat = self._kernel(x, x)
    for _ in range(epoch):
        err = -y * (self._alpha.dot(k_mat) + self._b)
        if np.max(err) < 0:
            continue
        mask = err >= 0
        delta = lr * y[mask]
        self._alpha += np.sum(delta[..., None] * k_mat[mask], axis=0)
        self._b += np.sum(delta)

# 预测代码
def predict(self, x, raw=False):
    x = np.atleast_2d(x).astype(np.float32)
    # 计算原样本与新样本的核矩阵并根据它来计算预测值
    k_mat = self._kernel(self._x, x)
    y_pred = self._alpha.dot(k_mat) + self._b
    if raw:
        return y_pred
    return np.sign(y_pred).astype(np.float32)

相关数学理论

1)Cover 定理

若设 d 维空间中 N 个点线性可分的概率为p(d,N),那么就有:

p(d,N)=\frac{2\sum_{i=0}^mC_{N-1}^i}{2^N}=\left\{
\begin{aligned}
&\frac{\sum_{i=1}^dC^i_{N-1}}{2^{N-1}},\ \ &N>d+1 
\\
&1,\ \ &N\le d+1
\end{aligned}
\right.

其中

m=\min(d,N-1)

证明从略(也就是说我不会)(喂),但是不难从中看出,它证明了当空间的维数 d 越大时、其中的 N 个点线性可分的概率就越大,这构成了核方法的理论基础之一

2)Mercer 定理

K(x_i,x_j)是对称函数(亦即K(x_i,x_j)=K(x_j,x_i))的话,那么它具有 Hilbert 空间中内积形式的充要条件有以下两个:

  • 对任何平方可积的函数g、满足
    \int{K(x_i,x_j)g(x_i)g(x_j)dx_idx_j}\ge0
  • 对含任意 N 个样本的数据集D={x_1,x_2,...,x_N},核矩阵:
    \bold K = \left[\begin{matrix}
K(x_1,x_1) & \ldots & K(x_1,x_N) \\
\vdots & \ddots & \vdots \\
K(x_N,x_1) & \ldots & K(x_N,x_N)
\end{matrix}\right]_{N\times N}
    是半正定矩阵

【注意:通常我们会称满足这两个充要条件之一的函数为 Mercer 核函数而把核函数定义得更宽泛。不过如果不打算在理论上深入太多的话,将 Mercer 核函数简称为核函数是可以的。此外,虽说 Mercer 核函数确实具有 Hilbert 空间中的内积形式、但此时的 Hilbert 空间并不一定具有“维度”这么好的概念(或说、可以认为此时 Hilbert 空间的维度为无穷大;比如说 RBF 核,它映射后的空间就是无穷维的)】

3)表示定理


\mathcal{H}为核函数K对应的映射后的空间(RKHS),\|h\|_\mathcal{H}表示\mathcal{H}h的范数,那么对于任意单调递增的函数C和任意非负损失函数L、优化问题

\min_{h\in\mathcal{H}}L\left(h(x_1),...,h(x_N)\right)+C(\|h\|_{\mathcal{H}})

的解总可以表述为核函数K的线性组合

h^*(x)=\sum_{i=1}^N\alpha_iK(x,x_i)

这意味着对于任意一个损失函数和一个单调递增的正则化项组成的优化问题、我们都能够对其应用核方法


下一篇文章我们则会抛开梯度下降这个有些“偷懒”的做法,并介绍一种叫序列最小最优化(SMO)的算法

希望观众老爷们能够喜欢~

(猛戳我进入下一章!( σ'ω')σ )

文章被以下专栏收录