PRML第六章读书笔记——Kernel Methods 静止核、对偶表示、构造核、高斯核函数、Fisher核、径向基函数网络/Nadaraya-Watson模型、高斯过程回归/分类/自动相关确定

（这一章感觉比较抽象，而且和之前学的线性回归和逻辑回归之间有奇妙的联系。最后一页高斯过程用于分类时选模型的公式推导我略去了，没有细看）

许多线性参数模型可以被转化为一个等价的“对偶表示”dual representation，其中，预测是基于训练数据点处的核函数线性组合
对于非线性特征空间映射$\varphi (\mathbf{x})$，核函数为
$$k(\mathbf{x},{\mathbf{x}}^{\prime })=\varphi (\mathbf{x}{)}^T\varphi ({\mathbf{x}}^{\prime })$$
注意核函数是对称的
核函数有许多扩展，例如核技巧。如果有一个算法，只和输入向量的标量积有关，则可以用其它的核函数形式来替换这个标量积。例如把PCA扩展到非线性PCA，把knn扩展到非线性knn，以及带核的Fisher判别分析等。
（核技巧之前写过一篇博客：核技巧的一些知识点——CVMLI Prince读书随笔第7章，这里要融会贯通一下）

• 核技巧能把特征升到高维，从而避免避免维度灾难吗？
  实际上，升到高维后，数据仍然是高维空间中的低维流形。本身维度没有发生变化。

静止核stationary kernel

只和差有关，具有平移不变性
$$k(\mathbf{x},{\mathbf{x}}^{\prime })=k(\mathbf{x}-{\mathbf{x}}^{\prime })$$

径向基函数radial basis function（同质核homogeneous kernel）

只和距离有关，具有平移不变性
$$k(\mathbf{x},{\mathbf{x}}^{\prime })=k(\Vert \mathbf{x}-{\mathbf{x}}^{\prime }\Vert )$$

6.1 Dual Representations

考虑线性回归MAP解，损失函数为
J ( w ) = 1 2 ∑ n = 1 N { w T ϕ ( x n ) − t n } 2 + λ 2 w T w J(\bm w)=\frac{1}{2}\sum_{n=1}^N \{ \bm w^T\bm\phi(\bm x_n)-t_n \}^2+\frac{\lambda}{2}\bm w^T\bm w J(w)=21​n=1∑N​{wTϕ(xn​)−tn​}2+2λ​wTw
回顾第3章，这个解为
$$\mathbf{w}=(\lambda {\mathbf{I}}_M+{\mathbf{\Phi }}^T\mathbf{\Phi }{)}^{-1}{\mathbf{\Phi }}^T\mathbf{t}$$
其中$\mathbf{\Phi }=[\mathbf{\varphi }(x_1),\mathbf{\varphi }(x_2),\cdots ,\mathbf{\varphi }(x_N){]}^T$
用Woodbury恒等式，右侧变为
$$\begin{array}{rl}& [{\lambda }^{-1}{\mathbf{\Phi }}^T-{\lambda }^{-1}{\mathbf{\Phi }}^T({\mathbf{I}}_N+\mathbf{\Phi }{\lambda }^{-1}{\mathbf{\Phi }}^T{)}^{-1}\mathbf{\Phi }{\lambda }^{-1}{\mathbf{\Phi }}^T]\mathbf{t}\\ =& [{\lambda }^{-1}{\mathbf{\Phi }}^T-{\lambda }^{-1}{\mathbf{\Phi }}^T(\lambda {\mathbf{I}}_N+\mathbf{\Phi }{\mathbf{\Phi }}^T{)}^{-1}\mathbf{\Phi }{\mathbf{\Phi }}^T]\mathbf{t}\\ =& [{\lambda }^{-1}{\mathbf{\Phi }}^T(\lambda {\mathbf{I}}_N+\mathbf{\Phi }{\mathbf{\Phi }}^T{)}^{-1}(\lambda {\mathbf{I}}_N+\mathbf{\Phi }{\mathbf{\Phi }}^T-\mathbf{\Phi }{\mathbf{\Phi }}^T)]\mathbf{t}\\ =& {\mathbf{\Phi }}^T(\lambda {\mathbf{I}}_N+\mathbf{\Phi }{\mathbf{\Phi }}^T{)}^{-1}\mathbf{t}\end{array}$$
这波操作有点神奇。。。
（书上用的另一种推法，我没有细看，但是结果一样的，我为了保持思路连续，沿用了第3章的结论）
从而
$$y(\mathbf{x})=\mathbf{\varphi }(\mathbf{x}{)}^T\mathbf{w}=\mathbf{\varphi }(\mathbf{x}{)}^T{\mathbf{\Phi }}^T(\lambda {\mathbf{I}}_N+\mathbf{\Phi }{\mathbf{\Phi }}^T{)}^{-1}\mathbf{t}=\mathbf{k}(\mathbf{x}{)}^T(\mathbf{K}+\lambda {\mathbf{I}}_N{)}^{-1}\mathbf{t}$$

其中$\mathbf{K}=\mathbf{\Phi }{\mathbf{\Phi }}^T$是Gram矩阵，元素为$K_{nm}=\mathbf{\varphi }({\mathbf{x}}_n{)}^T\mathbf{\varphi }({\mathbf{x}}_m)=k({\mathbf{x}}_n,{\mathbf{x}}_m)$，向量$\mathbf{k}$由$k_n(\mathbf{x})=k({\mathbf{x}}_n,\mathbf{x})$组成

• 这里的对偶形式，完全由核函数表达
• 结果是新观测和训练集计算核函数的结果的线性组合，满足$y(\mathbf{x})={\sum }_{n=1}^N{\alpha }_{n}k({\mathbf{x}}_n,\mathbf{x})$
• 实际上$(\lambda {\mathbf{I}}_N+\mathbf{K}{)}^{-1}\mathbf{t}$是向量$\mathbf{\varphi }({\mathbf{x}}_n)$的线性组合（没证……）
• 这里对偶形式求逆的维度为$N$，当数据量大时，比原形式的复杂度要大很多
• 不过这里可以直接定义核函数$k$，而不用考虑特征映射$\mathbf{\varphi }$

感知机也可以写出对偶形式

6.2 Constructing Kernels

最简单的构造核方法是从
$$k(x,x^{\prime })=\varphi (x{)}^T\varphi (x^{\prime })$$
如图所示

但我们可以不需要直接设计$\varphi$，直接构造$k$，但是要保证它对应一个标量积（可以是无穷维）

关于核函数的判定参见最上方贴出的之前的博客，这里再回顾一下：

• 对于任意 { x n } \{\bm x_n\} {xn​}，$\mathbf{K}$需要半正定，其中元素由$k({\mathbf{x}}_n,{\mathbf{x}}_m$)组成

一种构造核的办法是根据已有的核构造新的

P296 高斯核

$$k(\mathbf{x},{\mathbf{x}}^{\prime })=\exp \left(-\Vert \mathbf{x}-{\mathbf{x}}^{\prime }{\Vert }^2/(2{\sigma }^2)\right)$$
可以认为构造方法为
$$\exp (-{\mathbf{x}}^T\mathbf{x}/(2{\sigma }^2)\exp ({\mathbf{x}}^T{\mathbf{x}}^{\prime }/{\sigma }^2)\exp (-{\mathbf{x}}^{\prime T}{\mathbf{x}}^{\prime }/(2{\sigma }^2))$$
使用上述图中(6.14)和(6.16)进行构造
注意特征向量是无穷维的
高斯核函数可以不局限于欧氏距离
$$k(\mathbf{x},{\mathbf{x}}^{\prime })=\exp (-\frac{1}{2{\sigma }^2}(\kappa (\mathbf{x},\mathbf{x})+\kappa ({\mathbf{x}}^{\prime },{\mathbf{x}}^{\prime })-2\kappa (\mathbf{x},{\mathbf{x}}^{\prime }))$$
（注意，高斯核函数是径向基核函数的一种）

P297 其他数据类型的核函数

核函数中的自变量不一定要是实数或向量，可以是其他符号（如图、集合、字符串等等）

P297 从生成式模型构造核函数

我们可以在判别式框架中使用生成式模型，生成式模型可以自然处理缺失数据，以及隐马尔科夫模型下的边长序列。判别式模型则在判别式任务中表现更好。
一种将结合的方式是使用生成式模型定义核，在判别式方法中使用
例如
$$k(\mathbf{x},{\mathbf{x}}^{\prime })=p(\mathbf{x})p({\mathbf{x}}^{\prime })$$
这里实际上把$p$看作是特征向量
进一步，可以构造出
$$k(\mathbf{x},{\mathbf{x}}^{\prime })=\sum _{i}p(\mathbf{x}|i)p({\mathbf{x}}^{\prime }|i)p(i)$$
或
$$k(\mathbf{x},{\mathbf{x}}^{\prime })=\int p(\mathbf{x}|\mathbf{z})p({\mathbf{x}}^{\prime }|\mathbf{z})p(\mathbf{z})d\mathbf{z}$$
在隐马尔科夫模型生成的序列数据中 X = { x 1 , ⋯   , x L } \bm X=\{\bm x_1, \cdots, \bm x_L\} X={x1​,⋯,xL​}，隐状态为 Z = { z 1 , ⋯   , z L } \bm Z=\{\bm z_1, \cdots, \bm z_L\} Z={z1​,⋯,zL​}，从而可以度量两个序列的相似性
$$k(\mathbf{X},{\mathbf{X}}^{\prime })=\sum _{\mathbf{Z}}p(\mathbf{X}|\mathbf{Z})p({\mathbf{X}}^{\prime }|\mathbf{Z})p(\mathbf{Z})$$
这种方法可以扩展到变长序列

P298 Fisher核

也是用于生成式模型。这里不展开讲，仅Mark一下
其想法是找到一个核，度量分布中两个变量$\mathbf{x},{\mathbf{x}}^{\prime }$之间的相似性（分布由参数$\mathbf{\theta }$控制）。Fisher得分Fisher score定义为
$$\mathbf{g}(\mathbf{\theta },\mathbf{x})={\nabla }_{\mathbf{\theta }}\ln p(\mathbf{x}|\mathbf{\theta })$$
Fisher核定义为
$$k(\mathbf{x},{\mathbf{x}}^{\prime })=\mathbf{g}(\mathbf{\theta },\mathbf{x}{)}^T{\mathbf{F}}^{-1}\mathbf{g}(\mathbf{\theta },{\mathbf{x}}^{\prime })$$
其中$\mathbf{F}$为Fisher信息矩阵Fisher information matrix
$$\mathbf{F}={\mathbb{E}}_{\mathbf{x}\sim p_{\mathbf{x}}(\cdot |\mathbf{\theta })}[g(\mathbf{\theta },\mathbf{x})g(\mathbf{\theta },\mathbf{x}{)}^T]$$
这种方法的动机来自信息几何information geometry
对于高斯分布$p(\mathbf{x}|\mu )$，均值是参数，方差固定为$\mathbf{S}$，Fisher核的形式可以算出为
$$k(\mathbf{x},{\mathbf{x}}^{\prime })=(\mathbf{x}-\mu {)}^T{\mathbf{S}}^{-1}({\mathbf{x}}^{\prime }-\mu )$$
这是马氏距离Mahalanobis distance的平方

6.3 Radial Basis Function Networks

P299 输入特征带噪声

${\varphi }_j(\mathbf{x})=h(\Vert \mathbf{x}-{\mu }_j\Vert )$

考虑输入$\mathbf{x}$带有噪声的回归问题，噪声分布为$\nu (\mathbf{\xi })$
则平方和损失变成
E = 1 2 ∑ n = 1 N ∫ { y ( x n + ξ ) − t n } 2 ν ( ξ ) d ξ E=\frac{1}{2}\sum_{n=1}^N\int \{ y(\bm x_n+\bm\xi)-t_n\}^2 \nu(\bm\xi)d\bm\xi E=21​n=1∑N​∫{y(xn​+ξ)−tn​}2ν(ξ)dξ
用变分法，可以得到
$$y(\mathbf{x})=\sum _{n=1}^N{t}_{n}h(\mathbf{x}-{\mathbf{x}}_n)$$
其中
$$h(\mathbf{x}-{\mathbf{x}}_n)=\frac{\nu (\mathbf{x}-{\mathbf{x}}_n)}{{\sum }_{i=1}^N\nu (\mathbf{x}-{\mathbf{x}}_i)}$$
如果$\nu (\xi )$各向同性，即只和$\Vert \xi \Vert$有关，那么基函数就是径向的
径向基函数出现了！
${\sum }_{n}h(\mathbf{x}-{\mathbf{x}}_n)=1$，这种标准化的作用如图所示。这可以避免所有的基函数都很小。

如果每个数据点都关联一个基函数，那么开销很大。可以只让部分数据关联基函数。例如随机选择一个子集，或者采用正交最小平方orthogonal least squares，顺序选择每一个样本点，使平方和损失下降最快。也可以基于K-means算法，得到一组基中心，这是基中心不与数据点重合

P301 Nadaraya-Watson模型

Nadaraya-Watson模型和第3章的等价核模型，可以对比起来看
假定有训练集 { x n , t n } \{\bm x_n, t_n\} {xn​,tn​}，采用第2章的Parzen密度估计来建模联合分布
$$p(\mathbf{x},t)=\frac{1}{N}\sum _{n=1}^{N}f(\mathbf{x}-{\mathbf{x}}_n,t-t_n)$$
其中$f$是分量密度函数，则
$$\begin{array}{rl}y(\mathbf{x})& =\mathbb{E}[t|\mathbf{x}]={\int }_{-\infty }^{\infty }tp(t|\mathbf{x})dt\\ & =\frac{\int tp(\mathbf{x},t)dt}{\int p(\mathbf{x},t)dt}\\ & =\frac{{\sum }_n\int tf(\mathbf{x}-{\mathbf{x}}_n,t-t_n)dt}{{\sum }_m\int f(\mathbf{x}-{\mathbf{x}}_m,t-t_m)dt}\end{array}$$
假定$${\int }_{-\infty }^{\infty }f(\mathbf{x},t)tdt=0$$
则$$y(\mathbf{x})=\sum _{n}k(\mathbf{x},{\mathbf{x}}_n)t_n$$
其中
$$\begin{array}{rl}k(\mathbf{x},{\mathbf{x}}_n)=& \frac{g(\mathbf{x}-{\mathbf{x}}_n)}{{\sum }_{m}g(\mathbf{x}-{\mathbf{x}}_m)}\\ g(\mathbf{x})=& {\int }_{-\infty }^{\infty }f(\mathbf{x},t)dt\end{array}$$
而且易得，这个结果称之为Nadaraya-Watson模型，或者核回归kernel regression
$$\sum _{n}k(\mathbf{x},{\mathbf{x}}_n)=1$$
这和等价核类似
注意这里不仅仅是$t$的期望，$p(t|\mathbf{x})$是可以求出来的
$$p(t|\mathbf{x})=\frac{{\sum }_{n}f(\mathbf{x}-{\mathbf{x}}_n,t-t_n)}{{\sum }_m\int f(\mathbf{x}-{\mathbf{x}}_m,t-t_m)dt}$$
如果把$f$建模成各向同性的高斯，则结果如图所示

此外，可以把$p(t,\mathbf{x})$建模成混合高斯，这样核函数不在训练点处，但是成分少了，推断时计算效率高了。

总之，在（1）输入变量被污染；（2）Nadaraya-Watson模型，这两种假设下，都能得到同样的径向基函数模型

6.4 Gaussian Processes

直接定义函数$y$上的先验概率分布，而不是定义在参数上。（这句话很深刻，也是高斯过程的主要特点。另外，高斯过程是一个随机过程）
很多模型都能看作是高斯过程的形式，包括ARMA、卡尔曼滤波、径向基函数网络等。

高斯过程是定义在函数$y(\mathbf{x})$上的一个概率分布，使得在任意点集${\mathbf{x}}_1,\cdots ,{\mathbf{x}}_N$计算的$y$的集合联合起来服从高斯分布。更一般的情况下可以用更合理的方式为$y({\mathbf{x}}_1),\cdots ,y({\mathbf{x}}_N)$赋予联合概率分布，来确定一个随机过程$y(\mathbf{x})$，高斯过程中，$N$个变量$y_1,\cdots ,y_N$上的联合概率分布完全由二阶统计确定，大部分应用中，直接把$y$的均值设为0，方差则由核函数确定
$$\mathbb{E}[y({\mathbf{x}}_n)y({\mathbf{x}}_m)]=k({\mathbf{x}}_n,{\mathbf{x}}_m)$$
对于线性回归而言
$$\mathbf{y}=\mathbf{\Phi }\mathbf{w}$$
$$\begin{array}{rl}\mathbb{E}[\mathbf{y}]=& \mathbf{\Phi }\mathbb{E}[\mathbf{w}]=0\\ cov[\mathbf{y}]=& \mathbb{E}[\mathbf{y}{\mathbf{y}}^T]=\mathbf{\Phi }\mathbb{E}[\mathbf{w}{\mathbf{w}}^T]{\mathbf{\Phi }}^T=K\end{array}$$
也可以直接定义核函数，而不用选择基函数。如图所示（图中横轴是$x$，纵轴是$y$，每一条线是一个$y(x)$的采样。这里应该可以理解为有无穷个 { x , y } \{x,y\} {x,y}样本），第一个是高斯核，第二个是指数核$k(x,x^{\prime })=\exp (-\theta |x-x^{\prime }|)$，后者对应Ornstein-Uhlenbeck过程，描述布朗运动

P306 高斯过程回归

在回归中，要考虑噪声。$y=y(\mathbf{x})，t=y+ϵ$，也即
$$\begin{array}{rl}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})\end{array}$$
根据线性高斯模型，得到
$$p(\mathbf{t})=\mathcal{N}(\mathbf{t}|0,\mathbf{C})$$
其中
$$C({\mathbf{x}}_n,{\mathbf{x}}_m)=k({\mathbf{x}}_n,{\mathbf{x}}_m)+{\beta }^{-1}{\delta }_{nm}$$

一种常用的核函数是指数-二次核exponential-quadratic kernel
$$k({\mathbf{x}}_n,{\mathbf{x}}_m)={\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^T{\mathbf{x}}_m$$

推断时，考虑一个新样本点$p(t_{N+1}|\mathbf{t})$（这里略去了$\mathbf{x}$为了书写方便）。先考虑$p({\mathbf{t}}_{N+1})$，有
$$p({\mathbf{t}}_{N=1})=\mathcal{N}({\mathbf{t}}_{N+1}|0,{\mathbf{C}}_{N+1})$$
其中${\mathbf{C}}_{N+1}$是$(N+1)\times (N+1)$的协方差矩阵，满足
$${\mathbf{C}}_{N+1}=\left[\begin{array}{cc}{\mathbf{C}}_N& \mathbf{k}\\ {\mathbf{k}}^T& c\end{array}\right]$$
其中$\mathbf{k}$中第$n$个元素为$k({\mathbf{x}}_n,{\mathbf{x}}_{N+1})$，$c=k({\mathbf{x}}_{N+1},{\mathbf{x}}_{N+1})+{\beta }^{-1}$
利用第2章求高斯分布的条件概率，得$p(t_{N+1}|\mathbf{t})$的均值、方差为
$$\begin{array}{rl}m({\mathbf{x}}_{N+1})& ={\mathbf{k}}^T{\mathbf{C}}_N^{-1}\mathbf{t}\\ {\sigma }^2({\mathbf{x}}_{N+1})& =c-{\mathbf{k}}^T{\mathbf{C}}_N^{-1}\mathbf{k}\end{array}$$

另外，如果考虑有限个样本，则这里的结果与贝叶斯线性回归完全一致（方差用woodbury恒等式直接替换就可以验证）
所以这里也有核函数的线性组合
$$m({\mathbf{x}}_{N+1})=\sum _{n=1}^N{a}_{n}k({\mathbf{x}}_n,{\mathbf{x}}_{N+1})$$
其中$a_n$是${\mathbf{C}}_N^{-1}\mathbf{t}$的第$n$个分量

• 对比高斯过程和线性基函数模型，高斯过程需要对维度为$N\times N$的矩阵求逆，而线性基函数只需要考虑特征有关的$M\times M$矩阵求逆（这也和woodbury恒等式的变换有关）。所以效率低。但是高斯过程的一个好处是可以考虑协方差只能被表达为无穷维的情况
• 对于大数据集，直接应用高斯过程很困难，所以有很多种近似方法
• 高斯过程可以扩展到高维目标观测上

P311 学习超参数

高斯过程的超参数$\mathbf{\theta }$在协方差矩阵的核函数中（也就是说核函数带参数了！）
所以目标可以是最大化似然$p(\mathbf{t}|\mathbf{\theta })$
$$\ln p(\mathbf{t}|\mathbf{\theta })=-\frac{1}{2}\ln |{\mathbf{C}}_N|-\frac{1}{2}{\mathbf{t}}^T{\mathbf{C}}_N^{-1}\mathbf{t}-\frac{N}{2}\ln (2\pi )$$
求梯度得到
$$\frac{\partial }{\partial {\theta }_i}\ln p(\mathbf{t}|\mathbf{\theta })=-\frac{1}{2}Tr\left({\mathbf{C}}_N^{-1}\frac{\partial {\mathbf{C}}_N}{\partial {\theta }_i}\right)+\frac{1}{2}{\mathbf{t}}^T{\mathbf{C}}_N^{-1}\frac{\partial {\mathbf{C}}_N}{\partial {\theta }_i}{\mathbf{C}}_N^{-1}\mathbf{t}$$
如果用完全贝叶斯的想法。可以直接在$\mathbf{\theta }$上引入先验，这样能直接积分出$p(\mathbf{t})$，但精确积分困难，这里进行近似（所以从贝叶斯的角度来看，求极大似然是一种近似）

P312 自动相关确定automatic relevance determination(ARD)

对于高维输入，ARD用来选特征的重要性。例如采用二维输入对应的核函数
$$k(\mathbf{x},{\mathbf{x}}^{\prime })={\theta }_0\exp \left\{-\frac{1}{2}\sum _{i=1}^2{\eta }_i(x_i-x_i^{\prime }{)}^2\right\}$$
从先验采样的$y$如图所示

当${\eta }_n$减小时，$y$和$x_i$之间的关系变得不敏感。
第二个例子是从$y=\sin (2\pi x_1)$中采样，$t$在$y$之上加入随机噪声；$x_2$则在$x_1$基础上加入随机噪声；$x_3$则从另一个独立的高斯分布中采样，优化结果如图所示

ARD可以整合到上文的指数-二次核中，其中$D$是输入维度
$$k({\mathbf{x}}_n,{\mathbf{x}}_m)={\theta }_0\exp \left\{-\frac{1}{2}\sum _{i=1}^D{\eta }_i(x_{ni}-x_{mi}{)}^2\right\}+{\theta }_2+{\theta }_3{\mathbf{x}}_n^T{\mathbf{x}}_m$$

P313 高斯过程分类

由于高斯过程中$y$是联合高斯分布，所以为了适应分类任务，要在输出$\mathbf{a}$套上激活函数得到$\mathbf{t}$。例如二分类问题，如图所示

$$p(t|a)=\sigma (a{)}^t(1-\sigma (a){)}^{1-t}$$
仿照回归
$$p({\mathbf{a}}_{N+1})=\mathcal{N}({\mathbf{a}}_{N+1}|0,{\mathbf{C}}_{N+1})$$
其中$C({\mathbf{x}}_n,{\mathbf{x}}_m)=k({\mathbf{x}}_n,{\mathbf{x}}_m)+\nu {\delta }_{nm}$，$\nu$是一个提前固定的参数，保证$\mathbf{C}$正定（注意$\mathbf{K}$半正定）
$$p(t_{N+1}=1|{\mathbf{t}}_N)=\int p(t_{N+1}=1|a_{N+1})p(a_{N+1}|{\mathbf{t}}_N)d{a}_{N+1}$$
其中$p(t_{N+1}=1|a_{N+1})=\sigma (a_{N+1})$
这积分没有解析解，所以采用近似。三种方法：

1. 变分推断variational inference
2. 期望传播expectation propagation
3. 拉普拉斯近似laplace approximation

这里考虑第三种方法，把$p(a_{N+1}|{\mathbf{t}}_N)$近似成高斯分布，这样就可以用第4章的高斯分布和logistic函数卷积进行近似（probit函数代替logistic函数）的方法求解
$$p(a_{N+1}|{\mathbf{t}}_N)=\int p(a_{N+1}|{\mathbf{a}}_N)p({\mathbf{a}}_N|{\mathbf{t}}_N)d{\mathbf{a}}_N\text{\hspace{1em}(1)}$$
这里$a_{N+1}$与${\mathbf{t}}_N$在给定${\mathbf{a}}_N$的条件下独立。其中已知
$$p(a_{N+1}|{\mathbf{a}}_N)=(a_{N+1}|{\mathbf{k}}^T{\mathbf{C}}_N^{-1}{\mathbf{a}}_N,c-{\mathbf{k}}^T{\mathbf{C}}_N^{-1}\mathbf{k})$$
这里再对$p({\mathbf{a}}_N|{\mathbf{t}}_N)$运用拉普拉斯近似
考虑到$p({\mathbf{a}}_N|{\mathbf{t}}_N)\propto p({\mathbf{a}}_N)p({\mathbf{t}}_N|{\mathbf{a}}_N)$，其中
$$\begin{array}{rl}p({\mathbf{a}}_N)& =\mathcal{N}({\mathbf{a}}_N|0,{\mathbf{C}}_N)\\ p({\mathbf{t}}_N|{\mathbf{a}}_N)& =\prod _{n=1}^N\sigma (a_n{)}^{t_n}(1-\sigma (a_n){)}^{1-t_n}=\prod _{n=1}^N{e}^{a_n{t}_n}\sigma (-a_n)\end{array}$$
所以$p({\mathbf{a}}_N|{\mathbf{t}}_N)$的概率对数为
$$\begin{array}{rl}\Psi ({\mathbf{a}}_N|{\mathbf{t}}_N)& =\ln p({\mathbf{a}}_N)+\ln p({\mathbf{t}}_N|{\mathbf{a}}_N)\\ & =-\frac{1}{2}{\mathbf{a}}_N^T{\mathbf{C}}_N^{-1}{\mathbf{a}}_N-\frac{N}{2}\ln (2\pi )-\frac{1}{2}\ln |{\mathbf{C}}_N|+{\mathbf{t}}_N^T{\mathbf{a}}_N-\sum _{n=1}^N\ln (1+e^{a_n})+\text{const}\end{array}$$

$$\nabla \Psi ({\mathbf{a}}_N)={\mathbf{t}}_N-{\mathbf{\sigma }}_N-{\mathbf{C}}_N^{-1}{\mathbf{a}}_N$$
因为${\mathbf{\sigma }}_N$和${\mathbf{a}}_N$有关，导致上式零点不好求，用泰勒展开近似（这第三波近似了）
$${\nabla }^2\Psi ({\mathbf{a}}_N)=-{\mathbf{W}}_N-{\mathbf{C}}_N^{-1}$$
其中${\mathbf{W}}_N$是对角矩阵，元素为$\sigma (a_n)(1-\sigma (a_n))$. 能看出${\nabla }^2\Psi ({\mathbf{a}}_N)$负定，所以泰勒展开近似的$p({\mathbf{a}}_N|{\mathbf{t}}_N)$上凸，有最大值
所以二次迭代公式为
a N n e w = C N ( I + W N C N ) − 1 { t N − σ N + W N a N } \bm a_N^{new}=\bm C_N(\bm I+\bm W_N\bm C_N)^{-1}\{ \bm t_N-\bm \sigma_N+\bm W_N\bm a_N \} aNnew​=CN​(I+WN​CN​)−1{tN​−σN​+WN​aN​}
这样一波算下来，能得到$q({\mathbf{a}}_N)$去近似$p({\mathbf{a}}_N|{\mathbf{t}}_N)$，进而求出
$$\begin{array}{rl}\mathbb{E}[a_{N+1}|{\mathbf{t}}_N]=& {\mathbf{k}}^T({\mathbf{t}}_N-{\mathbf{\sigma }}_N)\\ var[a_{N+1}|{\mathbf{t}}_N]=& c-{\mathbf{k}}^T({\mathbf{W}}_N^{-1}+{\mathbf{C}}_N{)}^{-1}\mathbf{k}\end{array}$$
得到$p(a_{N+1}|{\mathbf{t}}_N)$的近似，从而得到$p(t_{N+1}|{\mathbf{t}}_N)$
如果只关心决策边界，那么只需要看均值就已经结束了，不用再处理式（1）的积分

此外，我们还可以用上文的方法学习超参数$\mathbf{\theta }$，即优化$p({\mathbf{t}}_N|\mathbf{\theta })$，这一页公式太繁琐，我跳过了没看。。。（饶恕了自己。。）

可视化的结果如图

P319 高斯过程和神经网络的关系

当神经网络的隐含单元$M\to \infty$时（把隐含单元看作基函数，也即基函数的数量有无穷多的极限情况下），对于参数$\mathbf{w}$的一大类恰当先验分布，神经网络产生的函数分布会趋于高斯过程。不过，在极限情况下，神经网络的输出变量会相互独立，不再共享隐含单元的“统计优势”

参考文献：
[1] Christopher M. Bishop. Pattern Recognition and Machine Learning. 2006