白给的神经过程Neural Processes笔记 (一)

前言

神经过程neural processes, NP 是2018年ICML会议上提出的一种新的神经网络范式，总的来说是将高斯过程Gaussian Processes, GP 网络化，从而降低GP对计算量的要求。GP的好处是天然地能够对预测的结果给出一个置信度的评分，即uncertainty or confidence， 而且GP模型是非参数的（non parametric ，即模型参数的数量取决于数据集的大小），会根据输入样例自动抽取规律，不断更新自己的预测信息。相比之下，一般的神经网络在训练结束后，模型参数将被固定，无法一直学习。下文首先介绍GP的基础知识，然后过渡到NP模型，中间会插入一些toy example，即小的例子，来更深入地理解GP和NP模型的特点。讲解的过程中会有很多个人的观点且需要有一点点概率论和信号估计理论的基础知识，如有错误欢迎指正和讨论。

1. Gaussian Process (GP)

Gaussian Process的定义：如果对于任意的子集 { x 1 ,   x 2 , ⋯   , x n } ⊂ X \{x_1,~x_2,\cdots,x_n\}\subset\mathcal{X} {x1​, x2​,⋯,xn​}⊂X, 采样一个$\mathcal{X}\to \mathbb{R}$的映射函数$f\sim p(f)$, 记$\mathbf{f}=\left(f(x_1),f(x_2),\cdots ,f(x_n)\right)\in {\mathbb{R}}^n$, 都有$\mathbf{f}\sim \mathcal{N}(\cdot ,\cdot )$, 即 $p(f)$服从多元高斯分布(multivariate Gaussiandistribution)。那么称$p(f)$是高斯过程Gaussian Process

上面的定义应该怎么理解呢？最后一句表明， 高斯过程是一个概率分布$p(f)$，而这个概率分布对应的随机变量是一个函数$f$。这是高斯过程相对于其他机器学习模型最大的不同点，我们所熟知的机器学习模型都是先假设模型能够完全拟合输入输出数据，然后训练或迭代模型参数，这种学习范式对应的函数$f$是确定的，通过样本数据学习模型的参数，优化目标是由先验知识确定的模型的参数。一种另辟蹊径的方法是将拟合输入输出数据的函数作为优化目标，这个目标是函数本身也就是模型本身，是数学上的泛函问题。然而这个函数的候选者数量是无限的，函数形式是未知的。高斯过程模型则是将这种一般化的思路，加入了前提和假设，使得构造的问题能够被求解出来。如下面的例子:
$$\begin{array}{rl}y=& f(x)+\epsilon \\ y=& {\omega }^{T}x+\epsilon \cdots \cdots OLSandridgeregression\\ y=& {\omega }^T\varphi (x)+\epsilon \cdots \cdots kernelregression\end{array}$$高斯过程以概率为框架，将输入$x$输出$y$的回归问题转化为寻找这个未知的函数$f$，其中$y=f(x)$，的概率分布问题，然后以假定观测值$y$经过$f$映射后组成的随机序列满足多元高斯分布。基于此利用贝叶斯定理和多元高斯函数的边缘分布marginalizatiion distribution和条件分布conditioning distribution求解问题。下面我们将对整个求解过程逐一展开。(PS: 从名字上理解，高斯过程和随机过程，这里的过程应该是抽样这个动作的意思，随机过程的定义是在一个变量集合里面随机采样出一个变量的过程; 所以对应地高斯过程可以理解为在一组由无限个可能的均值 $\mu$ 和无限个可能的方差 $\Sigma$ 组合成的多元高斯函数集合中采样出一个符合给定输入样例的函数的过程，而且这个函数并不是一个确定函数，而是一个以给定训练数据集 $\mathcal{D}$ 和给定输入 $x$ 为条件的输出的采样函数 $f=p(y|\mathcal{D},x)\sim \mathcal{G}\mathcal{P}=\mathcal{N}({\mu }_{y|\mathcal{D}},{\Sigma }_{y|\mathcal{D}})$。

1.0 预备知识

预备知识可以参考这条有趣的链接：https://jgoertler.com/visual-exploration-gaussian-processes/

• 多元高斯函数（multi-variate Gaussian function）的数学描述：
  $$\mathbf{x}={\left[\begin{array}{cccc}{x}_1& x_2& \cdots & x_n\end{array}\right]}^T\sim {\mathcal{N}}_{\mathbf{x}}(\mu ,\Sigma )$$多元高斯函数也叫多元正态分布函数，其中
  $${\mathcal{N}}_{\mathbf{x}}(\mu ,\Sigma )=\frac{1}{(2\pi {)}^{n/2}|\Sigma {|}^{1/2}}exp\left(-\frac{1}{2}(\mathbf{x}-\mu {)}^T{\Sigma }^{-1}(\mathbf{x}-\mu )\right)$$ $|\Sigma |$ 表示计算 $\Sigma$ 矩阵的行列式（determination），$\mu =\mathbb{E}[\mathbf{x}]$ 是随机变量 $\mathbf{x}$ 的期望（expectation），在高斯函数中也叫均值。这里我们可以类比联想两个随机变量的协方差矩阵Covariance matrix，计算的是两个随机变量各个维度两两之间的相关程度：
  $$Cov(\mathbf{x},\mathbf{y})=\mathbb{E}\left[(\mathbf{x}-\mathbb{E}[\mathbf{x}])(\mathbf{y}-\mathbb{E}[\mathbf{y}])\right]=\mathbb{E}[\mathbf{xy}]-\mathbb{E}[\mathbf{x}]\mathbb{E}[\mathbf{y}]$$ 而 $\Sigma =Cov(\mathbf{x},\mathbf{x})=\mathbb{E}\left[(\mathbf{x}-\mu )(\mathbf{x}-\mu {)}^T\right]=\mathbb{E}[\mathbf{x}{\mathbf{x}}^T]-\mu {\mu }^T$ 为协方差，计算的是一个自相关矩阵，是一个对称的正定矩阵，对于任意的一个随机向量的协方差矩阵，它都是半正定矩阵。
• 多元高斯函数的边缘分布：
  假设随机变量拆分为两部分 $\mathbf{x}$ 和 $\mathbf{y}$，那么多元高斯函数的数学描述形式可以改写如下：
  $$\left[\begin{array}{c}\mathbf{x}\\ \mathbf{y}\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}{\mu }_{\mathbf{x}}\\ {\mu }_{\mathbf{y}}\end{array}\right],\left[\begin{array}{cc}{\Sigma }_{\mathbf{xx}}& {\Sigma }_{\mathbf{xy}}\\ {\Sigma }_{\mathbf{yx}}& {\Sigma }_{\mathbf{yy}}\end{array}\right]\right)$$记$x$ 和 $y$ 分别为随机变量 $\mathbf{x}$ 和 $\mathbf{y}$ 采样后的具体值，realization sample。假设随机变量序列 $\mathbf{x}$ 和 $\mathbf{y}$（或说样本）之间独立同分布（identity independent distribution, iid.），则有$p_{\mathbf{xy}}(x,y)=p_{\mathbf{y}|\mathbf{x}}(y|x)p_{\mathbf{x}}(x)=p_{\mathbf{x}|\mathbf{y}}(x|y)p_{\mathbf{y}}(y)$：
  $$p_{\mathbf{x}}(x)={\int }_y{p}_{\mathbf{xy}}(x,y)dy={\int }_y{p}_{\mathbf{x}|\mathbf{y}}(x|y)p_{\mathbf{y}}(y)dy$$其中，$p_{\mathbf{xy}}(x,y)$ 表示两个随机变量的联合概率分布， $p_{\mathbf{x}|\mathbf{y}}(x|y)$ 表示给定$y$ 的情况下随机变量 $\mathbf{x}$ 的条件概率分布 （conditioned on $\mathbf{y}$），而 $p_{\mathbf{x}}(x)$ 和 $p_{\mathbf{y}}(y)$ 分别表示随机变量 $\mathbf{x}$ 和 $\mathbf{y}$ 相对于联合概率分布的边缘概率分布。
• 多元高斯函数的条件分布：
  利用贝叶斯条件概率公式和高斯函数运算的一些基本法则：
  $$\begin{array}{c}p(x|y)=\frac{p(x,y)}{p(y)}=\frac{p(y|x)p(x)}{p(y)}\end{array}$$经过整理可以得到条件概率的计算结果如下：
  $$\begin{array}{rl}\mathbf{x}|\mathbf{y}& \sim \mathcal{N}\left({\mu }_{\mathbf{x}}+{\Sigma }_{\mathbf{xy}}{\Sigma }_{\mathbf{yy}}^{-1}(\mathbf{y}-{\mu }_{\mathbf{y}}),{\Sigma }_{\mathbf{xx}}-{\Sigma }_{\mathbf{xy}}{\Sigma }_{\mathbf{yy}}^{-1}{\Sigma }_{\mathbf{yx}}\right)\\ \mathbf{y}|\mathbf{x}& \sim \mathcal{N}\left({\mu }_{\mathbf{y}}+{\Sigma }_{\mathbf{yx}}{\Sigma }_{\mathbf{xx}}^{-1}(\mathbf{x}-{\mu }_{\mathbf{x}}),{\Sigma }_{\mathbf{yy}}-{\Sigma }_{\mathbf{yx}}{\Sigma }_{\mathbf{xx}}^{-1}{\Sigma }_{\mathbf{xy}}\right)\end{array}$$上面的整理结果将会在高斯过程回归（Gaussian process regression, GPR）中经常使用。
• 其他小知识点：
  1）. 如果 $a\sim {\mu }_a,{\Sigma }_a$, $b=Aa+\eta$, $\eta \sim \mathcal{N}(0,{\Sigma }_{\eta })$, 那么 $b\sim \mathcal{N}(A{\mu }_a,A{\Sigma }_a{A}^T+{\Sigma }_{\eta })$。
  2）. 如果 $$\left[\begin{array}{c}x\\ y\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}a\\ b\end{array}\right],\left[\begin{array}{cc}A& C\\ C^T& B\end{array}\right]\right)$$ 那么 (推导过程网上一抓一大把) $$x|y\sim \mathcal{N}(a+C{B}^{-1}(y-b),A-C{B}^{-1}{C}^T)$$
  3）. 贝叶斯公式Bayes formula
  $$posterior=\frac{likelihood\times prior}{marginallikelihood},p(x|y)=\frac{p(x,y)}{p(y)}=\frac{p(y|x)p(x)}{p(y)}$$

1.1 问题描述

为了更好地理解高斯过程，我们首先从实际用例出发现构建一个回归问题：

给定一组输入输出样例数据$(x_s\in {\mathbb{R}}^n,y_s\in {\mathbb{R}}^n)=\{(x_1,y_1),(x_2,y_2),\cdots ,(x_n,y_n)\}$，其中$x_i$为输入样本，$y_i$为输出样本，现给定新的输入$x_t\in {\mathbb{R}}^m$，预测$y_t\in {\mathbb{R}}^m$。

这是典型的回归问题，问题的核心在于通过数据驱动的方式拟合出输入输出的映射关系 $f$，使得 $y^{\prime }=f(x^{\prime })$ 。记 $y_s=(y_1,y_2,\cdots ,y_n)$， $x_s=(x_1,x_2,\cdots ,x_n)$; $\mathcal{G}\mathcal{P}(f|\mu (\mathcal{D}),\Sigma (\mathcal{D}))$ 为多元高斯函数集合。高斯过程首先假设输出序列 $y$ 服从多元高斯分布（现实世界中往往不完全符合，但是为了求解问题，当我们决定使用GP这种方法的时候就已经默认 $y$ 近似服从多元高斯分布）。高斯过程最终的预测输出结果是对映射函数的后验概率进行采样 $y^{\prime }\sim p(f|x,y,x^{\prime })$。由于整个过程都是对多元高斯函数的计算，因此后验概率也是多元高斯分布的形式，高斯过程的高明之处在于映射函数可以用均值和方差参数化表示。

1.2 求解方法

上述问题可以用数学公式描述为：$y=f(x)$。目标是预测映射函数 $f$ 的后验分布 $y_t\sim p(f|\mathcal{D},x_t)$ ，目标是根据输入数据预测模型的均值和方差，最后再根据贝叶斯公式，计算给定样本条件下的输出的条件概率做预测。根据贝叶斯概率框架的基本思想是先根据先验知识prior knowledge猜测变量的概率分布，再根据观测的数据likelihood调整概率分布posterior，这种方法有很多个名字但表达都是同样的概念：最大后验方法 MAP maximum a posterior 或者是岭回归ridge regression.$$p(f|\mathcal{D})=\frac{p(\mathcal{D}|f)p(f)}{p(\mathcal{D})}$$ 解题步骤如下：

• Step 1: 在上述问题中，$f$ 是对 $\mathcal{G}\mathcal{P}$ 的采样，所以 $y$ 是满足高斯分布的，根据上面提到的方法，我们假定输出$y$的方差与输入$x$的方差相等，这个时候，$y$ 的分布如下。 这样做没有毛病是因为我们有贝叶斯公式调整最后的结果，盲猜一下问题不大。
  $$\begin{array}{c}\left[\begin{array}{c}y\\ y_t\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}\mu (y)\\ {\mu }_{(}{y}_t)\end{array}\right],\left[\begin{array}{cc}\mathcal{K}(x,x^{\prime })& \mathcal{K}(x,x_t^{\prime })\\ \mathcal{K}(x_t,x^{\prime })& \mathcal{K}(x_t,x_t^{\prime })\end{array}\right]\right)\end{array}$$其中，$\mathcal{K}$ 表示计算两个向量之间的距离度量 （covariance, kernel），$\mu$ 表示均值（mean）函数。
• Step 2: 当观测数据$y$来了之后，我们先把上式的$x$项去掉，根据上面的多元高斯函数条件概率公式 (1) 可以得到条件概率分布，对步骤1的猜测进行修正：
  $$\begin{array}{c}\begin{array}{rl}{y}_t|\mathcal{D},x_t\sim & \mathcal{N}({\mu }_{y_t|\mathcal{D},x_t},{\Sigma }_{y_t|\mathcal{D},x_t})\\ {\mu }_{y_t|\mathcal{D},x_t}=& \mu (y_t)+\mathcal{K}(x_t,x^{\prime })\mathcal{K}(x,x^{\prime }{)}^{-1}(y-\mu (y))\\ {\Sigma }_{y_t|\mathcal{D},x_t}=& \mathcal{K}(x_t,x_t^{\prime })-\mathcal{K}(x_t,x^{\prime })\mathcal{K}(x,x^{\prime }{)}^{-1}\mathcal{K}(x,x_t^{\prime })\end{array}\end{array}$$

1.3 Question & Answer

以上就是Gaussian Process Regression (GPR) 的整个求解过程。下面我们将会通过Socrates的方式，通过提问和解答的方式来进一步理解公式(3)所告诉我们的内容：

• 1. 为什么高斯过程回归预测输出的同时可以给出预测结果的不确定性值？
     因为公式(3)的输出包含了方差和均值，正常来说，预测值应该通过对（3）结果中均值和方差构成的多元高斯函数进行采样获得。实际上可以直接取均值，而不确定性则是由对应的方差项得出。
• 2. 公式（2）中，为什么不用y直接统计出方差？
     从公式推导来看，是因为如果$y$的分布全都使用$y$的信息，$y$将完全独立，而不是以$x$为条件，$x$的信息也将没有办法融合到整个推导框架中。
• 3. 高斯过程回归能够预测输出的原理是什么？
     当我们仔细观察的公式（3）不难看出，我们可以改写均值的公式为：${\mu }_{y|\mathcal{D},x_t}=\mu (y_t)+\Phi (x,x_t)(y-\mu (y_t))$， 假设我们令$\mu (y)=0$，则有${\mu }_{y|\mathcal{D},x_t}=\Phi (x,x_t)y$。即，预测输出$y_t={\mu }_{y|\mathcal{D},x_t}$是训练集中的输出 $y$ 的线性组合，而整个预测的不确定性 ${\Sigma }_{y_t|\mathcal{D},x_t}$ 完全由 $x$ 和 $x_t$ 决定。再细看 $\mathcal{K}$的计算，更准确地说是由 $x_t$ 处于 $x$ 中的位置决定的。因此，总的来说，高斯过程回归实际上是提供了更加系统的数据插值方法（预测输出实际上就是在做数据的插值interpolation），它根据待预测的输入数据 $x_t$ 是否足够靠近训练集中的样本 $x$ 来确定预测的置信度，也根据自己与训练样本的输入来分配对应输出值的权重分配$\Phi (x,x_t)$，最后输出预测值。思想非常的简单，同时数学上的表述也非常的优美。
• 4. 为什么度量距离的核函数 $\mathcal{K}$ 这么重要？
     核函数之所以这么重要是因为对于训练样本中不存的输入，最终输出的权重分配是由核函数确定的，因此对于没有见过的样本，和距离训练集样本较远的待预测样本，其差值结果完全取决于核函数的选取。这也是为什么教程中屡屡强调核函数的重要性。

1.4 举个栗子写段代码

参照这个链接的用例，魔改了一下gaussian process函数，最终结果如下：

import numpy as np
import scipy
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.gridspec as gridspec
import seaborn as sns

# Set matplotlib and seaborn plotting style
sns.set_style('darkgrid')
np.random.seed(42)

def exponentiated_quadratic(x, y):
    sq_norm = -0.5*scipy.spatial.distance.cdist(x, y, 'sqeuclidean') # 计算两两样本之间的距离
    return np.exp(sq_norm)

class myGP:
    def __init__(self):
        pass
    
    def run_demo(self, domain = (-6, 6)):
        f_sin = lambda x: (np.sin(x)).flatten()
        
        n1 = 8  # number of points to condition on (training point)
        n2 = 75 # number of points in posterior (test points)
        ny = 5  # number of functions that will be sampled from the posterier
        
        X1 = np.random.uniform(domain[0]+2, domain[1]-2, size=(n1, 1))
        y1 = f_sin(X1)
        X2 = np.linspace(domain[0], domain[1], n2).reshape(-1, 1)
        mu2, Sigma2 = self.gaussian_process(X1, y1, X2, exponentiated_quadratic)
        # Compute posterior mean and covariance
        sig2 = np.sqrt(np.diag(Sigma2))

        # Draw some samples of the posterior
        y2 = np.random.multivariate_normal(mean=mu2, cov=Sigma2, size=ny)
        
        # Plot the postior distribution and some samples
        fig, (ax1, ax2) = plt.subplots(
            nrows=2, ncols=1, figsize=(6, 6))
        # Plot the distribution of the function (mean, covariance)
        ax1.plot(X2, f_sin(X2), 'b--', label='$sin(x)$')
        ax1.fill_between(X2.flat, mu2-2*sig2, mu2+2*sig2, color='red', 
                        alpha=0.15, label='$2 \sigma_{2|1}$')
        ax1.plot(X2, mu2, 'r-', lw=2, label='$\mu_{2|1}$')
        ax1.plot(X1, y1, 'ko', linewidth=2, label='$(x_1, y_1)$')
        ax1.set_xlabel('$x$', fontsize=13)
        ax1.set_ylabel('$y$', fontsize=13)
        ax1.set_title('Distribution of posterior and prior data.')
        ax1.axis([domain[0], domain[1], -3, 3])
        ax1.legend()
        # Plot some samples from this function
        ax2.plot(X2, y2.T, '-')
        ax2.set_xlabel('$x$', fontsize=13)
        ax2.set_ylabel('$y$', fontsize=13)
        ax2.set_title('5 different function realizations from posterior')
        ax1.axis([domain[0], domain[1], -3, 3])
        ax2.set_xlim([-6, 6])
        plt.tight_layout()
        plt.show()
    
    def gaussian_process(self, X1, y1, X2, kernel_func):
        sigma11 = kernel_func(X1, X1)
        sigma12 = kernel_func(X1, X2)        
        tmp1112 = np.dot(sigma12.transpose(), np.linalg.pinv(sigma11)) # Sigma21=sigma12.transpose()
        mu2     = np.dot(tmp1112, y1)
        Sigma22 = kernel_func(X2, X2)
        Sigma2  = Sigma22 - np.dot(tmp1112, sigma12)
        return mu2, Sigma2
        
 if __name__ == '__main__':
 	gp = myGP()
    gp.run_demo()

这里我们用到了scipy中的scipy.spatial.distance.cdist(x, y, ‘sqeuclidean’) 来计算两两用例之间的距离关系，核函数类型可以通过最后一个参数给定。

以上，
2023年4月23日
Dianye