从零单排PRML 第6章 Kernel Methods

从零单排PRML 第6章 Kernel Methods

作者: Houye

单位: 北京邮电大学博士生

知乎主页: https://www.zhihu.com/people/Jhy1993

公众号 图与推荐

PRML一般指代Pattern Recognition and Machine Learning(模式识别与机器学习)一书。该书出版于2006年，是贝叶斯机器学习领域的经典之作。

更重要的是,算法工程面试的真题很多都来自这本书. 作者在面试阿里巴巴研究型实习生的时候就被问到了这本书的很多知识.

马春鹏同学的中译版本给了我们很好的参考,但是其并没有对PRMl的内容进行解读.

从零单排PRML系列旨在用浅显易懂的白话对机器学习相关知识进行梳理(个人理解,可能有不精确的地方,欢迎批评指教).

同时,作者对一些公式和模型罗列了一些简单的Python代码示例以促进理解.代码示例主要来自https://github.com/ctgk/PRML

引言

机器学习里面对待训练数据有2种,抛弃或者保留

• 有的是训练完得到参数后就可以抛弃了，比如神经网络；
• 有的是还需要原来的训练数据比如KNN，SVM也需要保留一部分数据–支持向量。

许多线性参数模型可以转化为对偶形式,转化后依然是线性模型.

只是线性组合的对象变成了一组非线性基函数.

这里的非线性基函数就是在基于各个数据点计算的核函数.

对偶表示

最小均方误差形式
$$J(\mathbf{w})=\frac{1}{2}\sum _{n=1}^N{\left\{{\mathbf{w}}^T\mathbf{\varphi }\left({\mathbf{x}}_n\right)-t_n\right\}}^2+\frac{\lambda }{2}{\mathbf{w}}^T\mathbf{w}$$
对偶形式:核函数形式
$$J(\mathbf{a})=\frac{1}{2}{\mathbf{a}}^T\mathbf{K}\mathbf{K}\mathbf{a}-{\mathbf{a}}^T\mathbf{K}\mathbf{t}+\frac{1}{2}{\mathbf{t}}^T\mathbf{t}+\frac{\lambda }{2}{\mathbf{a}}^T\mathbf{K}\mathbf{a}$$
这里K是著名的Gram矩阵$\mathbf{K}=\mathbf{\Phi }{\mathbf{\Phi }}^T$.上式中的$\boldsymbol { a } $又有闭式解
$$\mathbf{a}={\left(\mathbf{K}+\lambda {\mathbf{I}}_N\right)}^{-1}\mathbf{t}$$
这样,问题的解$\boldsymbol { a } $就只与核函数$\boldsymbol { K } $有关了.

最小均方转为核函数的组合$k\left(x,x^{\prime }\right)$后,我们就只关心$x , x ^ { \prime } $基于核函数的结果,对偶表示可以避免了显式的使用特征向量$\phi(x)$.这样做不关心中间过程$\boldsymbol { \phi} ( x ) ^ { \mathrm { T } } \boldsymbol { \phi } \left( x ^ { \prime } \right) $ ,只关心结果$k\left(x,x^{\prime }\right)$.

核函数形式$k \left( x , x ^ { \prime } \right) =\boldsymbol { \phi} ( x ) ^ { \mathrm { T } } \boldsymbol { \phi } \left( x ^ { \prime } \right) $ ,也就是映射后高维特征空间的内积可以通过原来低维的特征得到, 这样做的好处是可以使用隐式的高维特征,所以kernel method得到了广泛应用. 以SVM为例,低维不可分的数据,经过kernel将特征映射到高维后,就可以分开了.

以下式为例, 这里 $k ( \boldsymbol { x } , \boldsymbol { z } ) $ 的形式很简洁,但是如果想找到$\varphi (x)$是比较难的.这里的$\varphi (x)$只是中间需求,最终我们需要只是$k ( \boldsymbol { x } , \boldsymbol { z } ) $.
$$\begin{gathered} k(\mathbf{x},\mathbf{z})={\left({\mathbf{x}}^T\mathbf{z}\right)}^2,形式简洁,好算 \\ \varphi (\mathbf{x})={\left(x_1^2,\sqrt{2}{x}_1{x}_2,x_2^2\right)}^{\mathrm{T}},很难找到,计算量大 \end{gathered}$$

构建kernel的方法

特征映射法

选择一个特征映射函数$\boldsymbol { \phi} ( x ) $ ,并通过这个映射来寻找kernel.
$$k\left(x,x^{\prime }\right)=\mathbf{\varphi }(x{)}^{\mathrm{T}}\mathbf{\varphi }\left(x^{\prime }\right)=\sum _{i=1}^M{\varphi }_i(x){\varphi }_i\left(x^{\prime }\right)$$

直接构造

直接写出kernel的形式,
$$k(\mathbf{x},\mathbf{z})={\left({\mathbf{x}}^{\mathrm{T}}\mathbf{z}\right)}^2$$
这里kernel的形式是写出来了,但是这个形式合法吗?

合法性检验有2个方法:

1. 找到对应的特征映射函数

2. Gram矩阵在所有可能的{xn}的选择下都是半正定的.

下面是通过方法1来对kernel进行合法性检查.上式可以认为是二维输入$\mathbf{x}=\left(x_1,x_2\right)$, 经过特征映射法的结果,原因是可以将上式展开了为2个经过特征映射后的$\boldsymbol { \phi} ( x ) $的内积.
$$\begin{array}{rl}k(\mathbf{x},\mathbf{z})& ={\left({\mathbf{x}}^{\mathrm{T}}\mathbf{z}\right)}^2={\left(x_1{z}_1+x_2{z}_2\right)}^2\\ & =x_1^2{z}_1^2+2x_1{z}_1{x}_2{z}_2+x_2^2{z}_2^2\\ & =\left(x_1^2,\sqrt{2}{x}_1{x}_2,x_2^2\right){\left(z_1^2,\sqrt{2}{z}_1{z}_2,z_2^2\right)}^{\mathrm{T}}\\ & =\varphi (\mathbf{x}{)}^{\mathrm{T}}\varphi (\mathbf{z})\end{array}$$
这里的特征映射函数形式为$\varphi (\mathbf{x})={\left(x_1^2,\sqrt{2}{x}_1{x}_2,x_2^2\right)}^{\mathrm{T}}$ ,可以验证$k(\mathbf{x},\mathbf{z})={\left({\mathbf{x}}^{\mathrm{T}}\mathbf{z}\right)}^2$ 是一个合法的kernel.

基于现有的合法kernel进行构建

假设 $k _ { 1 } \left( \mathbf { x } , \mathbf { x } ^ { \prime } \right) $和$ k _ { 2 } \left( \mathbf { x } , \mathbf { x } ^ { \prime } \right)$都是合法kernel,那么通过以下形式组合构建的新的kernel也是合法的.
$$\begin{gathered} k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)=c{k}_1\left(\mathbf{x},{\mathbf{x}}^{\prime }\right),c>0 \\ k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)=k_1\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)+k_2\left(\mathbf{x},{\mathbf{x}}^{\prime }\right) \\ k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)=\exp \left(k_1\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)\right) \\ k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)=f(\mathbf{x})k_1\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)f\left({\mathbf{x}}^{\prime }\right),fisanyfunction \\ ... \end{gathered}$$
这里以经典的 高斯核 为例, 通过其构造过程来验证其合法性.

首先高斯核函数是一种径向基核
$$k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)=\exp \left(-{\Vert \mathbf{x}-{\mathbf{x}}^{\prime }\Vert }^2/2{\sigma }^2\right)$$
把exp中的平方项展开得到
$$\begin{gathered} {\Vert \mathbf{x}-{\mathbf{x}}^{\prime }\Vert }^2={\mathbf{x}}^{\mathrm{T}}\mathbf{x}+{\left({\mathbf{x}}^{\prime }\right)}^{\mathrm{T}}{\mathbf{x}}^{\prime }-2{\mathbf{x}}^{\mathrm{T}}{\mathbf{x}}^{\prime } \\ k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)=\exp \left(-{\mathbf{x}}^{\mathrm{T}}\mathbf{x}/2{\sigma }^2\right)\exp \left({\mathbf{x}}^{\mathrm{T}}{\mathbf{x}}^{\prime }/{\sigma }^2\right)\exp \left(-{\left({\mathbf{x}}^{\prime }\right)}^{\mathrm{T}}{\mathbf{x}}^{\prime }/2{\sigma }^2\right) \end{gathered}$$
这里线性kernel $k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)={\mathbf{x}}^{\mathrm{T}}{\mathbf{x}}^{\prime }$ 是合法的,然后基于上述的构造性质,可以看出高斯核是合法的.同时,由于$exp$的泰勒展开有无穷多项,所以高斯核将特征向量映射到无穷多维.

通过概率生成模型来构建

给定一个生成模型$p(x)$ , 我们可以定义kernel
$$k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)=p(\mathbf{x})p\left({\mathbf{x}}^{\prime }\right)$$
这里的$p(x)$可以认为是$\varphi (x)$.

高斯过程

什么是高斯过程?

首先高斯分布定义为函数$y(\mathbf{x})$上的一个高斯分布,并且这个函数在一系列点 ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_N$ 的响应值 $y\left({\mathbf{x}}_1\right),\dots ,y\left({\mathbf{x}}_N\right)$ 联合概率分布也是高斯的.

高斯过程=高斯+过程

• 高斯: 对于$y(\mathbf{x})$,为什么它是高斯的呢?我们令
  $$\begin{gathered} y(\mathbf{x})={\mathbf{w}}^{\mathrm{T}}\mathbf{\varphi }(\mathbf{x}) \\ p(\mathbf{w})=\mathcal{N}(\mathbf{w}|0,{\alpha }^{-1}\mathbf{I}) \end{gathered}$$
  这里的$y$是一些$\mathbf{w}$的混合,而$\mathbf{w}$ 我们给了一个高斯分布的先验.然后利用性质:高斯分布的线性组合依然是高斯分布.所以$y(\mathbf{x})$也是高斯分布. 这里我们只要确定$y(\mathbf{x})$的均值和方差就可以确定这个高斯分布.

• 过程:时间或者空间上一系列点${\mathbf{x}}_1,\dots ,{\mathbf{x}}_N$ 的响应值$y\left({\mathbf{x}}_1\right),\dots ,y\left({\mathbf{x}}_N\right)$ .

考虑一系列的点,我们可以得到
$$\mathbf{y}={\mathbf{\Phi }}_{\mathbf{W}}$$
下面就是这个多维高斯分布的均值和协方差.
$$\begin{array}{rl}\mathbb{E}[\mathbf{y}]& =\mathbf{\Phi }\mathbb{E}[\mathbf{w}]=0\\ \mathrm{cov}[\mathbf{y}]& =\mathbb{E}\left[{\mathbf{y}\mathbf{y}}^{\mathrm{T}}\right]=\mathbf{\Phi }\mathbb{E}\left[\mathbf{w}{\mathbf{w}}^{\mathrm{T}}\right]{\mathbf{\Phi }}^{\mathrm{T}}=\frac{1}{\alpha }\mathbf{\Phi }{\mathbf{\Phi }}^{\mathrm{T}}=\mathbf{K}\end{array}$$
这里我们进一步做了简化,将$y(\mathbf{x})$的均值设为0,那么一个高斯过程就仅仅受协方差矩阵的控制.更重要的是,这个协方差矩阵我们用一个kernel来确定,kernel形式为特征映射法的形式. 这也是为什么高斯过程放在kernel method这一章! 综上,我们得到确定一个高斯过程的步骤:

1. 确定核函数的形式
2. 求解核形式的协方差矩阵
3. 确定高斯过程

有了前几章的基础,书中这里举了2个核函数例子(高斯核,指数核).给定了$\mathbf { y } $的分布,我们可以进行采样,得到不同的$\mathbf { y } $.下图就是基于不同的核采样得到的$\mathbf { y } $.

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-eoEuGGeW-1580383052689)(https://ws2.sinaimg.cn/large/006tNbRwly1fy38weej2tj30nt0813zk.jpg)]

到这里什么是高斯过程及其如何确定一个高斯过程就清楚了.接下来是一些高斯过程应用(回归,分类).

基于高斯过程的回归

考虑观测值$y_n$上的一个高斯噪声${ϵ}_n$,
$$\begin{gathered} t_n=y_n+{ϵ}_n \\ {y}_n=y({\mathbf{x}}_n) \end{gathered}$$
这样$\mathbf{t}$可以认为是从以$\mathbf{y}$为均值的高斯分布中采样得到.似然函数如下
$$p(\mathbf{t}|\mathbf{y})=\mathcal{N}(\mathbf{t}|\mathbf{y},{\beta }^{-1}{\mathbf{I}}_N)$$
基于高斯过程的定义, 下式为高斯过程先验
$$p(\mathbf{y})=\mathcal{N}(\mathbf{y}|0,\mathbf{K})$$
可以看出:核函数$K$控制先验.

为了得到$\mathbf{t}$的分布,我们需要对其联合概率分布求积分,即边缘化
$$p(\mathbf{t})=\int p(\mathbf{t}|\mathbf{y})p(\mathbf{y})\mathrm{d}\mathbf{y}=\mathcal{N}(\mathbf{t}|0,\mathbf{C})$$
这里的$p ( \mathbf { t } | \mathbf { y } ) 和 p ( \mathbf { y } ) $都是高斯分布,其中
$$\begin{gathered} C\left({\mathbf{x}}_n,{\mathbf{x}}_m\right)=k\left({\mathbf{x}}_n,{\mathbf{x}}_m\right)+{\beta }^{-1}{\delta }_{nm} \\ {\delta }_{nm}=1,ifm=n \end{gathered}$$
这里意思是只有协方差矩阵的对角线叠加了噪声.这是因为每次的噪声都是相互独立,所以非对角线的位置相关性都为0.这里对核函数$K$ 的唯一限制就是协方差矩阵是正定的.假设${\lambda }_i$ 是核函数$K$ 的特征值,那么协方差矩阵$\mathbf{C}$的特征值就是${\lambda }_i+{\beta }^{-1}$.核函数需要满足半正定,即${\lambda }_i\ge 0$,加上一个正数$\beta$就可以${\lambda }_i+{\beta }^{-1}>0$, 所以$\mathbf{C}$是正定.

回归问题可以认为是基于之前的输入${\mathbf{x}}_1,\dots ,{\mathbf{x}}_N$和输出${\mathbf{t}}_N={\left(t_1,\dots ,t_N\right)}^{\mathrm{T}}$ ,对新的输入${\mathbf{x}}_{N+1}$ 进行预测(求$t_{N+1}$).写成条件概率的形式就是
$$\begin{gathered} p\left(t_{N+1}|{\mathbf{t}}_N,{\mathbf{x}}_1,\dots ,{\mathbf{x}}_N,{\mathbf{x}}_{N+1}\right),(省略输入) \\ =p\left(t_{N+1}|{\mathbf{t}}_N\right) \end{gathered}$$
前面已经得到了$p(\mathbf{t})=\mathcal{N}(\mathbf{t}|0,\mathbf{C})$ (高斯分布),那么
$$p\left({\mathbf{t}}_{N+1}\right)=\mathcal{N}\left({\mathbf{t}}_{N+1}|0,{\mathbf{C}}_{N+1}\right)$$
这里的$\mathbf { C } _ { N + 1 } $和前面的$\mathbf { C }_N$ 的关系为,为什么要写成这种形式呢?
$${\mathbf{C}}_{N+1}=\left(\begin{array}{cc}{\mathbf{C}}_N& \mathbf{k}\\ {\mathbf{k}}^{\mathrm{T}}& c\end{array}\right)$$
写成上述形式,可以利用公式
$$\begin{array}{rl}{\mathbf{\mu }}_{a|b}& ={\mathbf{\mu }}_a+{\mathbf{\Sigma }}_{ab}{\mathbf{\Sigma }}_{bb}^{-1}\left({\mathbf{x}}_b-{\mathbf{\mu }}_b\right)\\ {\mathbf{\Sigma }}_{a|b}& ={\mathbf{\Sigma }}_{aa}-{\mathbf{\Sigma }}_{ab}{\mathbf{\Sigma }}_{bb}^{-1}{\mathbf{\Sigma }}_{ba}\end{array}$$
这里的$p\left({\mathbf{x}}_a|{\mathbf{x}}_b\right)=p\left(t_{N+1}|{\mathbf{t}}_N\right)$ , b对应数据${\mathbf{t}}_N={\left(t_1,\dots ,t_N\right)}^{\mathrm{T}}$ , a对应$t_{N+1}$.代入得到
$$\begin{gathered} p\left(t_{N+1}|\mathbf{t}\right)=\mathcal{N}\left({\mathbf{t}}_{N+1}|m\left({\mathbf{x}}_{N+1}\right),{\sigma }^2\left({\mathbf{x}}_{N+1}\right)\right) \\ \begin{array}{rl}m\left({\mathbf{x}}_{N+1}\right)& ={\mathbf{k}}^{\mathrm{T}}{\mathbf{C}}_N^{-1}\mathbf{t}\\ {\sigma }^2\left({\mathbf{x}}_{N+1}\right)& =c-{\mathbf{k}}^{\mathrm{T}}{\mathbf{C}}_N^{-1}\mathbf{k}\end{array} \end{gathered}$$

这里有了均值和方差就定了:具有任意核函数$k\left({\mathbf{x}}_n,{\mathbf{x}}_m\right)$ 的高斯过程回归的预测分布.

这里高斯过程回归的预测是一堆基函数的表示,也就是在函数空间

对应之前的线性回归的预测是一堆参数的表示,也就是在参数空间

学习超参数

第一次听到GP是用于神经网络超参数的搜索,这里终于看到了真身.

GP依赖于协方差函数的选择 ,而实际中用一组带有参数的函数来代替协方差函数,这组函数可以描述协方差矩阵的参数.举例如下:
$$k\left({\mathbf{x}}_n,{\mathbf{x}}_m\right)={\theta }_0\exp \left\{-\frac{{\theta }_1}{2}{\Vert {\mathbf{x}}_n-{\mathbf{x}}_m\Vert }^2\right\}+{\theta }_2+{\theta }_3{\mathbf{x}}_n^{\mathrm{T}}{\mathbf{x}}_m$$
这里的参数${\theta }_0,{\theta }_1...$控制了协方差矩阵.不同参数的效果如下

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-NFwk75Pa-1580383052689)(https://ws3.sinaimg.cn/large/006tNbRwly1fy8nz32snmj30ny0hr40v.jpg)]

这里利用最大似然来估计参数.
$$\ln p(\mathbf{t}|\mathbf{\theta })=-\frac{1}{2}\ln \left|{\mathbf{C}}_N\right|-\frac{1}{2}{\mathbf{t}}^{\mathrm{T}}{\mathbf{C}}_N^{-1}\mathbf{t}-\frac{N}{2}\ln (2\pi )$$
求导,并利用矩阵求导公式
$$\begin{array}{rl}\frac{\partial }{\partial \mathbf{x}}\left(A^{-1}\right)& =-A^{-1}\frac{\partial A}{\partial \mathbf{x}}{A}^{-1}\\ \frac{\partial }{\partial \mathbf{x}}\ln |A|& =\mathrm{Tr}\left(A^{-1}\frac{\partial A}{\partial \mathbf{x}}\right)\end{array}$$
可以得到
$$\frac{\partial }{\partial {\theta }_i}\ln p(\mathbf{t}|\mathbf{\theta })=-\frac{1}{2}\mathrm{Tr}\left({\mathbf{C}}_N^{-1}\frac{\partial {\mathbf{C}}_N}{\partial {\theta }_i}\right)+\frac{1}{2}{\mathbf{t}}^{\mathrm{T}}{\mathbf{C}}_N^{-1}\frac{\partial {\mathbf{C}}_N}{\partial {\theta }_i}{\mathbf{C}}_N^{-1}\mathbf{t}$$

自动相关性确定

⾼斯过程中的⾃动相关性确定: 通过最⼤似然⽅法进⾏的参数最优化，能 够将不同输⼊的相对重要性从数据中推断出来。

假定一个二维高斯过程$\mathbf{x}=\left(x_1,x_2\right)$,其核函数如下
$$k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)={\theta }_0\exp \left\{-\frac{1}{2}\sum _{i=1}^2{\eta }_i{\left(x_i-x_i^{\prime }\right)}^2\right\}$$
不同的参数${\eta }_i$的影响如下: 当${\eta }_i$变小的时候,核函数变的对相应的输入$x_i$ 不敏感了.

我们首先看看,根据参数${\eta }_i$来输入$x_i$的影响.下右图: 参数${\eta }_2$变小了,那么在其相应的输入$x_2$下,输出变得较为稳定了.直观理解就是上式中$\eta _ { i } \left( x _ { i } - x _ { i } ^ { \prime } \right) ^ { 2 } $变小了,影响自然也下降.

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-vtvhbWwB-1580383052690)(https://ws4.sinaimg.cn/large/006tNbRwly1fy8ohlkmgpj31010d1gnq.jpg)]

另一个角度,根据输入$x_i$来看看参数${\eta }_i$.

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-RExVYdfr-1580383052690)(https://ws3.sinaimg.cn/large/006tNbRwly1fy8ovn5sunj30ss0g4ju7.jpg)]

高斯过程分类

GP分类的预测在整个实数轴上,如果分类需要进行放缩（sigmoid来放缩).

这里通过$\sigma$函数将高斯过程$a(\mathbf{x})$ 变成了非高斯的随机过程$y$.

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-UvTPYdFo-1580383052690)(https://ws2.sinaimg.cn/large/006tNbRwly1fy8p1so66rj311d0emq4n.jpg)]

与前面的回归类似,分类问题可以认为是基于之前的输入${\mathbf{x}}_1,\dots ,{\mathbf{x}}_N$和输出${\mathbf{t}}_N={\left(t_1,\dots ,t_N\right)}^{\mathrm{T}}$ ,对新的输入${\mathbf{x}}_{N+1}$ 进行预测(求$t_{N+1}$).写成条件概率的形式就是
$$p\left(t_{N+1}|{\mathbf{t}}_N\right)$$
只不过这里的高斯过程换成了$a$,而预测的$t_{N+1}$ 是非高斯过程.
$$p\left({\mathbf{a}}_{N+1}\right)=\mathcal{N}\left({\mathbf{a}}_{N+1}|0,{\mathbf{C}}_{N+1}\right)$$
同时,在分类时,我们假定所有label都是正确的,所以协方差矩阵没有像分类一样叠加了一个噪声项${\beta }^{-1}{\delta }_{nm}$. 这里的$\nu$看起来像噪声,实际是为了保证协方差正定引入的.

对于2分类,将$t_{N+1}$预测为1的概率
$$\begin{array}{rl}p\left(t_{N+1}=1|{\mathbf{t}}_N\right)& =\int p\left(t_{N+1}=1|a_{N+1}\right)p\left(a_{N+1}|{\mathbf{t}}_N\right)\mathrm{d}{a}_{N+1}\\ p\left(t_{N+1}=1|a_{N+1}\right)& =\sigma \left(a_{N+1}\right)\end{array}$$
那么问题来了:上面这个积分怎么解决? 目前有3种方法:

• variational inference
• expectation propagation
• Laplace approximation (下一节重点)

拉普拉斯近似

关于拉普拉斯近似的详细介绍见 PRML 4.4 . 简而言之,laplace approximation就是使用Gaussian去近似一个连续变量的概率密度函数,这个Gaussian是以后验分布的概率的众数为中心的。

注: 本节的${\mathbf{a}}_N$就是前面的 ${\mathbf{t}}_N$.

上节说到, $p\left(t_{N+1}=1|{\mathbf{t}}_N\right)$ 不好求解.
$$\begin{array}{rl}p\left(a_{N+1}|{\mathbf{t}}_N\right)& =\int p\left(a_{N+1},{\mathbf{a}}_N|{\mathbf{t}}_N\right)\mathrm{d}{\mathbf{a}}_N\\ & =\frac{1}{p\left({\mathbf{t}}_N\right)}\int p\left(a_{N+1},{\mathbf{a}}_N\right)p\left({\mathbf{t}}_N|a_{N+1},{\mathbf{a}}_N\right)\mathrm{d}{\mathbf{a}}_N\\ & =\frac{1}{p\left({\mathbf{t}}_N\right)}\int p\left(a_{N+1}|{\mathbf{a}}_N\right)p\left({\mathbf{a}}_N\right)p\left({\mathbf{t}}_N|{\mathbf{a}}_N\right)\mathrm{d}{\mathbf{a}}_N\\ & =\int p\left(a_{N+1}|{\mathbf{a}}_N\right)p\left({\mathbf{a}}_N|{\mathbf{t}}_N\right)\mathrm{d}{\mathbf{a}}_N\end{array}$$
其中$p\left({\mathbf{t}}_N|a_{N+1},{\mathbf{a}}_N\right)=p\left({\mathbf{t}}_N|{\mathbf{a}}_N\right)$ .上述积分中

• 第一项是高斯过程,$p\left(a_{N+1}|{\mathbf{a}}_N\right)=\mathcal{N}\left(a_{N+1}|{\mathbf{k}}^{\mathrm{T}}{\mathbf{C}}_N^{-1}{\mathbf{a}}_N,c-{\mathbf{k}}^{\mathrm{T}}{\mathbf{C}}_N^{-1}\mathbf{k}\right)$
• 第二项是后验分布,如果我们能通过拉普拉斯近似找到一个高斯分布

那么,$p\left(t_{N+1}=1|{\mathbf{t}}_N\right)$就变成对2个高斯分布的积分,可以求解.

现在核心问题变成了第二项:后验怎么求?

后验=先验*似然(省略归一化)

• 先验:前面已经得到,

$$\ln p(\mathbf{t}|\mathbf{\theta })=-\frac{1}{2}\ln \left|{\mathbf{C}}_N\right|-\frac{1}{2}{\mathbf{t}}^{\mathrm{T}}{\mathbf{C}}_N^{-1}\mathbf{t}-\frac{N}{2}\ln (2\pi )$$

​ 把上式中的$t$换成$a$ ,就得到 $p({\mathbf{a}}_N)$

• 似然: 由现有的数据及其预测得到,
  $$\begin{gathered} p\left({\mathbf{t}}_N|{\mathbf{a}}_N\right)=\prod _{n=1}^N\sigma {\left(a_n\right)}^{t_n}{\left(1-\sigma \left(a_n\right)\right)}^{1-t_n}=\prod _{n=1}^N{e}^{a_n{t}_n}\sigma \left(-a_n\right) \\ \sigma {\left(a_n\right)}^{t_n}{\left(1-\sigma \left(a_n\right)\right)}^{1-t_n}=(\frac{e^{a_n}}{1+e^{a_n}}{)}^{t_n}(\frac{1}{e^{a_n}+1}{)}^{1-t_n}=e^{a_n{t}_n}\frac{1}{1+e^{a_n}} \end{gathered}$$

• 结合先验和似然,得到后验的表达形式

$$\begin{array}{rl}\Psi \left({\mathbf{a}}_N\right)=& \ln p\left({\mathbf{a}}_N\right)+\ln p\left({\mathbf{t}}_N|{\mathbf{a}}_N\right)\\ =& -\frac{1}{2}{\mathbf{a}}_N^T{\mathbf{C}}_N^{-1}{\mathbf{a}}_N-\frac{N}{2}\ln (2\pi )-\frac{1}{2}\ln \left|{\mathbf{C}}_N\right|+{\mathbf{t}}_N^T{\mathbf{a}}_N\\ & -\sum _{n=1}^N\ln \left(1+e^{a_n}\right)\end{array}$$

目前我们已经有了后验的表达形式,如何来拉普拉斯近似来解呢?

1. 找到后验分布的众数

   后验求导,并代入

$$[\ln (1+e^x){]}^{\prime }=\frac{e^x}{1+e^x}=\frac{1}{1+e^{-x}}=\sigma (x)$$

​ 得到
$$\nabla \Psi \left({\mathbf{a}}_N\right)={\mathbf{t}}_N-{\mathbf{\sigma }}_N-{\mathbf{C}}_N^{-1}{\mathbf{a}}_N$$
​ 其中,${\mathbf{\sigma }}_N$是元素为$\sigma (a_n)$的向量. 这里由于${\mathbf{\sigma }}_N$和${\mathbf{a}}_N$的非线性关系, 所以令一阶导为零没有闭式解.利用二阶导和IRLS算法我们可以找到众数 ${\mathbf{a}}_N^{\star }$.

2. 计算以众数 ${\mathbf{a}}_N^{\star }$为中心的高斯分布$q \left( \mathbf { a } _ { N } \right) $ .
   $$q\left({\mathbf{a}}_N\right)=\mathcal{N}\left({\mathbf{a}}_N|{\mathbf{a}}_N^{\star },{\mathbf{H}}^{-1}\right)$$

到这里, $p\left(a_{N+1}|{\mathbf{t}}_N\right)$已经可以解了,分别将2个高斯分布代入得到,最终$$p\left(a_{N+1}|{\mathbf{t}}_N\right)$$ 是一个高斯分布,均值方差如下
$$\begin{array}{rl}\mathbb{E}\left[a_{N+1}|{\mathbf{t}}_N\right]& ={\mathbf{k}}^{\mathrm{T}}\left({\mathbf{t}}_N-{\mathbf{\sigma }}_N\right)\\ \mathrm{var}\left[a_{N+1}|{\mathbf{t}}_N\right]& =c-{\mathbf{k}}^{\mathrm{T}}{\left({\mathbf{W}}_N^{-1}+{\mathbf{C}}_N\right)}^{-1}\mathbf{k}\end{array}$$
接下来,我们需要确定 协方差函数的参数$\theta$. 最大似然即可

参考

1. https://www.zhihu.com/question/46631426
2. PRML
3. http://blog.videolectures.net/100-most-popular-machine-learning-talks-at-videolectures-net/
4. MLAPP
5. http://www.52nlp.cn/prml读书会第六章-kernel-methods