Chapter 6：Kernel Methods

一：导论

1：在之前我们介绍的线性参数回归和分类模型或者是神经网络模型中，我们利用训练集(training set)给出参数$\vec w$的点推测或者是后验分布，当我们做预测时，训练集将不再被使用到，我们纯粹是基于参数$\vec w$的点推测或后验分布来做出预测的。同时还有另外一类技术，训练集不仅被用来训练模型参数，同时当我们对新的输入矢量$\vec x$做预测时，训练集也会被使用到，例如使用kernel function的Gaussian process。

2：如果模型是基于fixed nonlinear feature space mapping $\vec \phi{(\vec x)}$，则kernel function $k(\vec x,\vec x')=\vec \phi{(\vec x)}^{T}\vec \phi{(\vec x)}$，$k(\vec x,\vec x')=\vec x^{T}\vec x'$称为linear kernel。注意kernel function是特征向量的内积，kernel trick或者是kernel substitution恰好是利用了这一性质。同时我们也注意到kernel function 是一个对称函数$k(\vec x,\vec x')=k(\vec x',\vec x)$。有着各种各样形式的内核函数，比如stationary kernel $k(\vec x,\vec x')=k(\vec x-\vec x')$，又比如homogeneour kernel 或者是radial basis functions $k(\vec x,\vec x')=k(||\vec x-\vec x'||)$。

3：正如上面提到的，kernel trick或者是kernel substitution利用了内核函数是特征矢量的内积这一性质。kernel trick其基本思路是如果一个算法中的输入矢量$\vec x$仅仅是以标量积的形式存在，则我们可以用其它的kernel function去替换这一标量积(scalar product)。

二：Constructing Kernels

1：valid kernel的评判标准：内核函数能否可以写成特征矢量标量积的形式$\{\vec \phi{(\vec x)}^{T}\vec \phi{(\vec x)}\}$；或者是判断Gram矩阵K$\{K_{nm}=k(\vec x_{n},\vec x_{m})\}$对于所有可能的矢量$\vec x_{n}$是否是正半定矩阵。

2：构造kernel function的三种方法：

（1）：选择一个feature space mapping $\vec \phi{(\vec x)}$，然后令$k(\vec x,\vec x')=\vec \phi{(\vec x)}^{T}\vec \phi{(\vec x)}$；

（2）：直接构造一个kernel function，比如$k(\vec x,\vec z)=(\vec x^{T}\vec z)^{2}$，但必须要确保该kernel function是一个有效内核；

（3）：以简单内核为基础，通过下列技术来构造新的内核：假设我们有两个有效内核，$k_{1}(\vec x,\vec x')，k_{2}(\vec x,\vec x')$，则

$$k(\vec x,\vec x')=ck_{1}(\vec x,\vec x') \quad c为大于0的常数\\ k(\vec x,\vec x')=f(\vec x)k_{1}(\vec x,\vec x')f(\vec x') \quad f()为任意函数 \\k(\vec x,\vec x')=q(k_{1}(\vec x,\vec x')) \quad q()多项式，其coefficient为非负值 \\k(\vec x,\vec x')=exp(k_{1}(\vec x,\vec x')) \\k(\vec x,\vec x')=k_{1}(\vec x,\vec x')+k_{2}(\vec x,\vec x')\\k(\vec x,\vec x')=k_{1}(\vec x,\vec x')k_{2}(\vec x,\vec x')\\k(\vec x,\vec x')=k_{3}(\vec \phi{(\vec x),\vec \phi{(\vec x'))}} \quad \vec \phi{(\vec x)}为把\vec x映射到R^{M}的函数，k_{3}()为在R^{M}空间中的有效kernel\\ k(\vec x,\vec x')=\vec x^{T}A\vec x' \quad A是对称正半定矩阵\\k(\vec x,\vec x')=k_{a}(\vec x_{a},\vec x_{a}^{'})+k_{b}(\vec x_{b},\vec x_{b}^{'})\quad \vec x=(\vec x_{a},\vec x_{b})，k_{a}和k_{b}是有效kernel \\ k(\vec x,\vec x')=k_{a}(\vec x_{a},\vec x_{a}^{'})k_{b}(\vec x_{b},\vec x_{b}^{'})$$
我们希望构造出来的kernel $k(\vec x,\vec x')$对称并且是正半定，同时它应该表达着 $\vec x,\vec x'$的相似性，当然这相似性要视具体的应用而定。

3：经常会用到的一些kernel function：

1）：$k(\vec x,\vec x')=(\vec x^{T}\vec x'+c)^{M}$；

2）：Gaussian kernel： $k(\vec x,\vec x')=exp(-||\vec x-\vec x'||^{2}/(2\sigma^{2}))$

因为$||\vec x-\vec x'||^{2}=\vec x^{T}\vec x+(\vec x')^{T}\vec x'-2\vec x^{T}\vec x'$，则如果我们运用kernel substitution技术，用非线性kernel $k(\vec x,\vec x')$替换$\vec x^{T}\vec x'$，则我们可以获得一个新的kernel $k(\vec x,\vec x')=exp\{-\frac{1}{2\sigma^{2}}(k(\vec x,\vec x)+k(\vec x',\vec x')-2k(\vec x,\vec x'))\}$

3）：利用generative model $p(\vec x)$来定义kernel

方法1：我们可以构造简单kernel $k(\vec x,\vec x')=p(\vec x,\vec x')$，这个kernel表示如果一对$\vec x,\vec x'$对应着较高的概率，则它们是相似的。我们能够在此比较简单kernel的基础上扩展一个新的内核$k(\vec x,\vec x')=\sum_{i}p(\vec x|i)p(\vec x'|i)p(i)$，其中p(i)为正的权重因子，如果采取无限和的极限，则$k(\vec x,\vec x')=\int p(\vec x|\vec z)p(\vec x'|\vec z)p(\vec z)d\vec z$，其中$\vec z$称为latent variable。

方法二：我们可以利用gengerative model $p(\vec x)$构造Fisher kernel。假设我们有一个parametric generative model $p(\vec x|\vec \theta)$，则$k(\vec x,\vec x')=\vec g(\vec \theta,\vec x)^{T}F^{-1}\vec g(\vec \theta,\vec x)$，其中$\vec g(\vec \theta,\vec x)=\bigtriangledown_{\vec \theta} \ln{p(\vec x|\vec \theta)}$，F为Fisher information matrix，表达式为$E_{\vec x}[\vec g(\vec \theta,\vec x)\vec g(\vec \theta,\vec x)^{T}]=\int \vec g(\vec \theta,\vec x)\vec g(\vec \theta,\vec x)^{T} p(\vec x|\vec \theta)d\vec x$

但是通常Fisher information matrix是不可解的，因此我们可以采取样品平均近似，表达式为$F \approx \frac{1}{N}\sum_{n=1}^{N}\vec g(\vec \theta,\vec x_{n})\vec g(\vec \theta,\vec x_{n})^{T}$；或者我们可以略去Fisher information matrix，获得Fisher kernel $k(\vec x,\vec x')=\vec g(\vec \theta,\vec x)^{T}\vec g(\vec \theta,\vec x)$

4）：$k(\vec x,\vec x')=tanh(a\vec x^{T}\vec x'+b)$，虽然它的Gram matrix不是正半定的，但是此kernel经常被使用。

三：Gaussian Processes

1：求解回归或者是分类问题的两种思路：

第一种思路：basis function Model：我们令$p(t|\vec x,\vec w)=N(t|y(\vec x,\vec w),\beta^{-1})=N(t|\vec w^{T}\vec \phi(\vec x),\beta^{-1})$或者是$p(C_{k}|\vec x,\vec w)=\sigma(a)=\sigma(\vec w^{T}\vec \phi(\vec x))$，通过训练集(training set)训练模型，我们给出$\vec w$的后验分布，然后再给出预测分布；

第二种思路：Gaussian Processes：我们不定义basis function，我们直接给函数$y(\vec x)$或者是$a(\vec x)$定义一个概率分布，使得一些列数据点$\vec x_{1},...,\vec x_{N}$对应的$y(\vec x)$或者是$a(\vec x)$值的联合分布$p(y(\vec x_{1}),...,p(\vec x_{N}))$或者是$p(a(\vec x_{1}),...,a(\vec x_{N}))$为高斯分布，其均值为0，方差为Gram矩阵K($K_{nm}=k(\vec x_{n},\vec x_{m})$)，$k(\vec x_{n},\vec x_{m})$为kernel function，所以在Gaussian Processes这一方法中，合适的kernel function的选取显得尤为重要。

2：Gaussian processes for regression：

1）：Calculating the prediction distribution：

假设$p(t_{n}|y_{n})=N(t_{n}|y_{n},\beta^{-1})$，$\beta$为超参数并且表示的是具有高斯分布噪音的precision，因为每个观察数据点的噪音是相互独立的，则$p(\vec t|\vec y)=N(\vec t|\vec y,\beta^{-1}I_{N})$，其中$\vec t=(t_{1},...,t_{N})^{T}，\vec y=(y_{1},...,y_{N})，I_{N}为N*N的单位矩阵$。从Gaussian processes的定义我们知道$p(\vec y)=N(\vec y|0,K)$，决定Gram矩阵K的kernel function应该要表达这样一条性质：对于相似点的$\vec x_{n}和\vec x_{m}$，它们对应的值$y_{n}和y_{m}$应该要比非相似点对更加强烈的关联，当然相似的定义视具体实际情况而定。

则边际分布:

$$p(\vec t)=\int p(\vec t|\vec y)p(\vec y)d\vec y=N(\vec t|0,C) \\ C(\vec x_{n},\vec x_{m})=k(\vec x_{n},\vec x_{m})+\beta^{-1}\delta_{nm}$$ 对于regression问题，一个广泛使用的kernel function是 $k(\vec x_{n},\vec x_{m})=\theta_{0}exp\{-\frac{\theta_{1}}{2}||\vec x_{n}-\vec x_{m}||^{2}\}+\theta_{2}+\theta_{3}\vec x_{n}^{T}\vec x_{m}$。

我们最终想要得到的预测分布:

$$p(t_{N+1}|\vec x_{N+1},\vec t_{N})=N(t_{N+1}|m(\vec x_{N+1}),\sigma^{2}(\vec x_{N+1})) \\ m(\vec x_{N+1})=\vec k^{T}C_{N}^{-1}\vec t \\ \sigma^{2}(\vec x_{N+1})=c-\vec k^{T}C_{N}^{-1}\vec k$$
$\vec k$由一些列元素 $k(\vec x_{n},\vec x_{N+1})$组成，并且 $m(\vec x_{N+1})$ $可以表达成$ $\sum_{n=1}^{N}a_{n}k(\vec x_{n},\vec x_{N+1})$，在这里 $a_{n}$是 $C_{N}^{-1}\vec t$的第n个元素。

Gaussian processes的计算包括了计算一个N*N矩阵的逆，所需的时间为$O(N^{3})$，以及矢量与矩阵的相乘，所需时间为$O(N^{2})$；而在basis function的模型中，所需的时间分别为$O(M^{3})和O(M^{2})$，一般来说训练集中数据点的个数要多于basis function的个数，使用Gaussian processes求解效率更低一点，但是Gaussian processes的一个优点是 we can consider kernel function that only can be expressed in terms of an infinite number of basis functions。

2）：Learning the hyperparameters：

正如我们上面看到的，Gaussian processes的预测分布部分地依赖于kernel function的选择，我们一般会选择带有参数的kernel function，然后从训练集中推出这些参数的值或者是后验分布，我们称这些参数为hyperparameters，推断出参数相关信息有极大似然法和贝叶斯方法这两种方法:

极大似然法：我们用基于梯度的最优化算法，比如conjugate gradients最大化似然函数$p(\vec t|\vec \theta)$的对数，其中$\theta$代表着Gaussian processes的hyperparameters，其中:

$$\ln p(\vec t|\vec \theta)=-\frac{1}{2}\ln|C_{N}|-\frac{1}{2}\vec t^{T}C_{N}^{-1}\vec t-\frac{N}{2}\ln(2\pi) \\ \frac{\partial \ln p(\vec t|\vec \theta)}{\partial \theta_{i}}=-\frac{1}{2}Tr(C_{N}^{-1}\frac{\partial C_{N}}{\partial \theta_{i}})+\frac{1}{2}\vec t^{T}C_{N}^{-1}\frac{\partial C_{N}}{\partial \theta_{i}}C_{N}^{-1}\vec t$$
但必须要说的是 $\ln{p(\vec t|\vec \theta)}$一般来说是nonconvex函数，因此它可能有着多个最大值；

贝叶斯方法：我们能够定义一个关于$\theta$的先验分布，然后用基于梯度的算法最大化其后验分布的log，获得一组$\vec \theta_{MAP}$；在全贝叶斯处理中，我们需要对$\vec \theta$进行边际化处理，但这一般没有解析解，因此我们需要做近似处理。

3）：Automatic relevance determination（ARD）：

我们之前给出的kernel function为

$$k(\vec x_{n},\vec x_{m})=\theta_{0}exp\{-\frac{\theta_{1}}{2}||\vec x_{n}-\vec x_{m}||^{2}\}+\theta_{2}+\theta_{3}\vec x_{n}^{T}\vec x_{m}$$
在这里参数 $\theta_{1}$对应着所有的输入变量；

现在运用ARD方法，我们能够对每个输入变量分别定义一个参数，则上述的kernel function显示如下：

$$k(\vec x_{n},\vec x_{m})=\theta_{0}exp\{-\frac{1}{2}\sum_{i=1}^{D}\eta_{i}(x_{ni}-x_{mi})^{2}\}+\theta_{2}+\theta_{3}\sum_{i=1}^{D}x_{ni}x_{mi}$$
此方法的好处是the optimization of these parameters by maximum likelihood allows the relative importance of different inputs to be inferred from the data。

3：Gaussian processes for classification

1）：Calculating the predictive distribution：

在之前的分类问题中，我们使用公式$p(t=1|\vec x)=\sigma(a)$，当我们用Gaussian processes处理分类问题时，我们能够对$a(\vec x)$使用Gaussian processes，然后使用logistic sigmoid把$a(\vec x)$转变成y{$\sigma(a)$}，这样我们就获得了一个关于函数$y(\vec x)$的non-Gaussian stochastic process。对于两类问题，我们知道$p(t|a)=\sigma(a)^{t}(1-\sigma(a))^{1-t}$。

假设我们有一个训练集，其输入为$\vec x_{1},...,\vec x_{N}$，对应的目标变量$\vec t_{N}=(t_{1},...,t_{N})^{T}$，并且我们也有一个test point $\vec x_{N+1}$，对应的目标值为$t_{N+1}$。由于我们对$a(\vec x)$定义了一个gaussian processes，因此$p(\vec a_{N+1})=N(\vec a_{N+1}|0,C_{N+1})$，其中$C(\vec x_{n},\vec x_{m})=k(\vec x_{n},\vec x_{m})+\upsilon\delta_{nm}$，$\upsilon$的值会预先确定好，在这里要说的是这里的kernel function由一组参数$\theta$控制。

预测分布为$p(t_{N+1}=1|\vec t_{N})=\int p(t_{N+1}=1|a_{N+1})p(a_{N+1}|\vec t_{N})da_{N+1}$，这个没有解析解，因此我们可以对$p(a_{N+1}|\vec t_{N})$采取拉普拉斯近似，近似的结果为$p(a_{N+1}|\vec t_{N})=N(a_{N+1}|E[a_{N+1}|\vec t_{N}],var[a_{N+1}|\vec t_{N}])$，这样我们能够使用结果(4.153)计算预测分布$p(t_{N+1}=1|\vec t_{N})$。拉普拉斯近似的具体计算过程见课本。

2）：Learning the hyperparameters：

正如我们上面说的，kernel function是由一组参数控制的，因此我们能够通过最大化似然函数$p(\vec t_{N}|\vec \theta)$来获得参数值，$p(\vec t_{N}|\vec \theta)=\int p(\vec t_{N}|\vec a_{N})p(\vec a_{N}|\vec \theta)d\vec a_{N}$