传统深度模型的uncertainty----Monte Carlo dropout

前言

贝叶斯的一系列方法让我们在对数据建模的时候可以评估模型的认知不确定性，即：我们可以获得一种模型对自己的预测值有多大把握的评价指标，这种评价在许多方法中通常以方差的形式出现。本文介绍2016年在模型uncertainty方面比较有影响力的一篇文章，文章的附录中有大量的推导，考虑到其中有些公式的变换涉及到额外的文献，对一些推导细节不做深入讨论，只大概了解其过程和思想，重点看如何搭建模型，以及如何体现uncertainty。
文章：Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning
附录：Dropout as a Bayesian Approximation: Appendix

背景

Dropout

设有一个只有一层隐藏层的神经网络，隐藏层的神经元数量为$K$，由第一层连接隐藏层的网络权重为$W_1$，大小为$Q\times K$，由隐藏层连接输出的权重为$W_2$，大小为$K\times D$，隐藏层的偏置向量为$b$，大小为$K$。

如果在隐藏层和输出层之前添加dropout应该如何用数学表示？
设${\mathbf{z}}_1,{\mathbf{z}}_2$分别为$Q$和$K$维的二进制向量，向量中的每个元素只能取$0$或$1$，dropout就相当于每个元素服从伯努利分布，即：
$$\begin{gathered} z_{1,q}\sim \text{Bernouli}(p_1)forq=1,...,Q \\ {z}_{2,k}\sim \text{Bernouli}(p_2)fork=1,...,K \end{gathered}$$
其中$p_1,p_2$就是每层神经元激活的概率。将$\mathbf{z}$表示成对角矩阵$\text{diag}(z)$，网络的输出可以表示为:
$$\hat{y}=\sigma (x(z_1{W}_1)+b)(z_2{W}_2)$$
其中$\sigma$为激活函数。网络的损失函数表示为（只考虑回归问题）：
$${\mathcal{L}}_{dropout}=\frac{1}{2N}\sum _{n=1}^N||y_n-{\hat{y}}_n|{|}_2^2+{\lambda }_1||W_1|{|}_2^2+{\lambda }_2||W_2|{|}_2^2+{\lambda }_3||b|{|}_2^2$$
记住dropout所对应的损失函数的表示形式，因为文章的中心意思就是在说变分推断+高斯过程也可以得到相同表示形式的损失函数。

高斯过程

详见高斯过程回归----MC dropout前导篇

变分推断

详见贝叶斯神经网络----从贝叶斯准则到变分推断

如何使用dropout近似贝叶斯

高斯过程

高斯过程是贝叶斯的一种表现形式，现定义一个高斯过程的协方差函数为：
$$K(x,y)=\int p(w)p(b)\sigma (w^{T}x+b)\sigma (w^{T}y+b)\text{d}w\text{d}b$$
使用Monte Carlo近似上式：
$$\hat{K}(x,y)=\frac{1}{K}\sum _{i=1}^K\sigma (w_k^{T}x+b_k)\sigma (w_k^{T}y+b_k)$$
式中$p(\mathbf{w})$和$p(\mathbf{b})$均服从标准的多元高斯分布，即每一元服从均值为$0$，方差为$1$的高斯分布。Monte Carlo采样的次数$K$与隐藏神经元的数量相同，所以此时的$W_1$是一个$Q\times K$的矩阵。

从高斯过程来看，输入到输出的关系如何表示？
$$\begin{gathered} F|X,W_1,b\sim N(0,\hat{K}(X,X)) \\ Y|F\sim N(F,{\tau }^{-1}{I}_N) \end{gathered}$$
这样看来好像只有$W_1$，网络隐藏层去哪了？先接着向下计算，先将$F$视为一个隐变量，通过对上式积分可以得到由输入预测输出的条件概率：
$$p(Y|X)=\int p(Y|F)p(F|W_1,b,X)p(W_1)p(b)\text{d}{W}_1\text{d}b\text{d}F$$
令：
$$\begin{gathered} \varphi (x,W_1,b)=\sqrt{\frac{1}{K}}\sigma (W_1^{T}x+b) \\ \Phi =[\varphi (x_n,W_1,b){]}_{n=1}^N\hat{K}(X,X)=\Phi {\Phi }^T \end{gathered}$$
此时条件概率$p(Y|X)$改写为（高斯分布的性质）：
$$p(Y|X)=\int N(Y;0,\Phi {\Phi }^T+{\tau }^{-1}{I}_N)p(W_1)p(b)\text{d}{W}_1\text{d}b\text{d}F$$
上式的分号表示条件概率，即已知后面的分布，计算$Y$对应的概率值。因为$Y$的每个元素$y_d$代表一个样本的输出，是一个$D$维的向量，所以每个$y_d$是一个多元高斯分布，整个$Y$是一个高斯过程。现引入辅助随机变量$w_d\sim N(0,I_K)$，则有(不要管为什么可以这样改写，作者引入了其他的参考文献)：
$$N(y_d;0,\Phi {\Phi }^T+{\tau }_{-1}{I}_N)=\int N(y_d;\Phi w_d,{\tau }^{-1}{I}_N)N(w_d;0,I_K)\text{d}{w}_d$$
看起来好像是把其中的一个$\Phi$换成了另一种高斯分布的参数，这也是前面的Monte Carlo要采样$K$次的原因：为了与隐藏层的神经元数量保持一致，接下来，令$W_2=[w_d{]}_{d=1}^D$，为$K\times D$的矩阵，条件概率可以改写为：
$$p(Y|X)=\int p(Y|X,W_1,W_2,b)p(W_1)p(W_2)p(b)\text{d}{W}_1\text{d}{W}_2\text{d}b$$
这个公式可以认为是神经网络预测的表达式，但它实际上是通过高斯过程实现的，直接计算高斯分布会牵涉到对于$N\times N$矩阵的逆运算，复杂度很高，所以文章使用变分推断的方式来近似求解。

变分推断

定义一个变分分布：
$$q(W_1,W_2,b):=q(W_1)q(W_2)q(b)$$
$:=$表示用右边的表达式来表达左边。定义$q(W_1)$和$q(W_2)$均为高斯混合分布，$q(b)$为一元高斯分布，则有：
$$\begin{gathered} q(W_1)=\prod _{q=1}^{Q}q(w_q) \\ q(w_q)=p_{1}N(m_q,{\sigma }^2{I}_K)+(1-p_1)N(0,{\sigma }^2{I}_K) \\ q(W_2)=\prod _{k=1}^{K}q(w_k) \\ q(w_k)=p_{2}N(m_k,{\sigma }^2{I}_D)+(1-p_1)N(0,{\sigma }^2{I}_D) \\ q(b)=N(m,{\sigma }^2{I}_K) \end{gathered}$$
令$M_1=[m_q{]}_{q=1}^Q,M_2=[m_k{]}_{k=1}^K$，则原来关于$W_1,W_2,b$的优化变成了对$M_1,M_2,m$的求解，这里的$q(W_1,W_2,b)$是一种条件分布，它是对后验分布$p(W_1,W_2,b|X,Y)$的近似。
我们看$q(w_q)$和$q(w_k)$的表示是由两个高斯分布组成，这很像伯努利的表示形式，即有$p$的概率选$m_q$，有$1-p$的概率为0，不同的地方在于$w_q$从分布中采样，而不是一个确定值。作者用这种分布来近似后验分布也是为了说明dropout的意义，需要注意的是，我们可以选择任意的分布形式来近似待求解的后验分布，这种混合高斯的形式只是一种可用形式而已，作者选它正好解释dropout为什么可以在传统的NN模型中引入不确定性，是因为用这种混合高斯来近似后验分布得到的损失函数与在传统的NN模型中添加dropout得到的损失函数一致。

回归问题的ELBO计算

在贝叶斯推断中，通过优化Evidence Lower Bound（ELBO）来求解后验分布：
$${\mathcal{L}}_{GP-VI}:=\int p(W_1,W_2,b)\log p(Y|X,W_1,W_2,b)-KL(q(W_1,W_2,b)||p(W_1,W_2,b|X,Y))\text{d}{W}_1\text{d}{W}_2\text{d}b$$
根据高斯分布的表达式$\log p(Y|X,W_1,W_2,b)$可以写为：
$$\begin{gathered} \log p(Y|X,W_1,W_2,b)=\sum _{d=1}^D\log N(y_d;\Phi w_d,{\tau }^{-1}{I}_N) \\ =-\frac{ND}{2}\log (2\pi )+\frac{ND}{2}\log (\tau )-\sum _{d=1}^D\frac{\tau }{2}||y_d-\Phi w_d|{|}_2^2 \end{gathered}$$
需要注意的一点，这里的$y_d$是相互独立的，此时$y_d$服从一个高斯分布，因为$Y$是一个$N\times D$的矩阵，如果按照单个样本来看，即$y_n$之间不条件独立，但是从另一个方向来看，即预测值的不同元之间是条件独立的，也就是说$y_d$是条件独立的。令$\hat{Y}=\Phi W_2$，有
$$\sum _{d=1}^D\frac{\tau }{2}||y_d-{\hat{y}}_d|{|}_2^2=\sum _{n=1}^N\frac{\tau }{2}||y_n-{\hat{y}}_n|{|}_2^2$$
此时${\hat{y}}_n=\varphi (x_n,W_1,b)W_2=\sqrt{\frac{1}{K}}\sigma (x_n{W}_1+b)W_2$，有
$$\begin{gathered} \log p(Y|X,W_1,W_2,b)=\sum _{n=1}^N\log N(y_n;\varphi (x_n,W_1,b)W_2,{\tau }^{-1}{I}_D) \\ =\sum _{n=1}^N\log p(y_n|x_n,W_1,W_2,b) \end{gathered}$$
此时的ELBO可以写为：
$$\sum _{n=1}^N\int p(W_1,W_2,b)\log p(y_n|x_n,W_1,W_2,b)-KL(q(W_1,W_2,b)||p(W_1,W_2,b|X,Y))\text{d}{W}_1\text{d}{W}_2\text{d}b$$
令
$$\begin{gathered} q({ϵ}_1)=N(0,I_{Q\times K}),q(z_{1,q})=\text{Bernouli}(p_1)forq=1,...,Q \\ q({ϵ}_2)=N(0,I_{K\times D}),q(z_{1,k})=\text{Bernouli}(p_2)fork=1,...,K \\ q(ϵ)=N(0,I_K) \end{gathered}$$
这里的$ϵ$是为了与$\sigma$相乘时能够得到输出的矩阵大小符合$W$的大小。根据原来变分近似的写法，参数可以分别重新写为：
$$\begin{gathered} W_1=z_1(M_1+\sigma {ϵ}_1)+(1-z_1)\sigma {ϵ}_1 \\ {W}_2=z_2(M_2+\sigma {ϵ}_2)+(1-z_2)\sigma {ϵ}_2 \\ b=m+\sigma ϵ \end{gathered}$$
则有：
$$\begin{gathered} \sum _{n=1}^N\int p(W_1,W_2,b)\log p(y_n|x_n,W_1,W_2,b)\text{d}{W}_1\text{d}{W}_2\text{d}b \\ =\sum _{n=1}^N\int p(z_1,{ϵ}_1,z_2,{ϵ}_2,ϵ)\log p(y_n|x_n,W_1(z_1,{ϵ}_1),W_2(z_2,{ϵ}_2),b(ϵ)) \end{gathered}$$
使用Monte Carlo对上式近似可得：
$$\begin{gathered} \sum _{n=1}^N\int p(W_1,W_2,b)\log p(y_n|x_n,W_1,W_2,b)\text{d}{W}_1\text{d}{W}_2\text{d}b \\ \approx \sum _{n=1}^N\log p(y_n|x_n,{\hat{W}}_1^n,{\hat{W}}_2^n,{\hat{b}}^n) \end{gathered}$$
又有对KL散度的近似（具有足够大的维数且均值随机分布的混合高斯函数的熵趋向于高斯函数体积之和，原文是这样讲的，我也不知道为什么）：
$$\begin{gathered} KL(q(W_1)||p(W_1))\approx QK({\sigma }^2-\log ({\sigma }^2)-1)+\frac{p_1}{2}\sum _{q=1}^Q{m}_q^T{m}_q+C \\ KL(q(W_2)||p(W_2))\approx QK({\sigma }^2-\log ({\sigma }^2)-1)+\frac{p_2}{2}\sum _{k=1}^K{m}_k^T{m}_k+C \\ KL(q(b)||p(b))=\frac{1}{2}(m^{T}m+K({\sigma }^2-\log ({\sigma }^2)-1))+C \end{gathered}$$
我们将原来的ELBO拆分成了近似的形式，忽略常数项，可以得到待优化的目标函数为：
$${\mathcal{L}}_{GP-MC}\propto -\frac{1}{2N}\sum _{n=1}^N||y_n-{\hat{y}}_n|{|}_2^2-\frac{p_1}{2\tau N}||M_1|{|}_2^2-\frac{p_2}{2\tau N}||M_2|{|}_2^2-\frac{1}{2\tau N}||m|{|}_2^2$$
写了这么多，推导过程其实还有很多细节，并且还有分类问题没有写，但是思想已经概括了：高斯过程+变分推断可以得到与NN+dropout相同的损失函数，而因为高斯过程本身就实现了对模型不确定性的评估，所以传统NN在添加dropout也可以引入不确定性。直观来看，在传统的神经网络中添加dropout之后，由于每个输入和每次前向传播dropout选择激活哪个神经元都是随机的，即使是相同的输入样本，得到的输出也有可能不一样，多次输入同一个样本，输出围绕均值上下浮动的程度就可以看成是模型的一种不确定性。

编程实现

import numpy as np
import matplotlib.pyplot as plt
import tensorflow._api.v2.compat.v1 as tf
tf.disable_v2_behavior()

class VariationalDense:
    def __init__(self, n_in, n_out, model_prob, lam):
        self.model_prob = model_prob
        self.lam = lam
        self.bern = tf.compat.v1.distributions.Bernoulli(probs=self.model_prob, dtype=tf.float32)
        self.model_M = tf.Variable(tf.random.truncated_normal([n_in, n_out], stddev=0.01))
        self.model_m = tf.Variable(tf.zeros([n_out]))
        self.model_W = tf.matmul(
            tf.diag(self.bern.sample((n_in, ))), self.model_M
        )

    def __call__(self, X, activation=tf.identity):
        output = activation(tf.matmul(X, self.model_W) + self.model_m)
        if self.model_M.shape[1] == 1:
            output = tf.squeeze(output)

        return output

    @property
    def regularization(self):
        return self.lam * (
                self.model_prob * tf.reduce_sum(tf.square(self.model_W)) + tf.reduce_sum(tf.square(self.model_m)))


n_samples = 100
X = np.random.normal(size=(n_samples, 1))
y = np.random.normal(np.cos(5. * X) / (np.abs(X) + 1.), 0.1).ravel() # np.ravel()函数将多维数组转化为一维数组
X_pred = np.atleast_2d(np.linspace(-3., 3, num=100)).T  # 将数组转化为至少两维
X = np.hstack((X, X**2, X**3)) # 将两个数组按水平方向合起来
X_pred = np.hstack((X_pred, X_pred**2, X_pred**3))

# create the tensorflow model
n_feats = X.shape[1]
n_hidden = 100
model_prob = 0.9
model_lam = 1e-2
model_X = tf.placeholder(tf.float32, [None, n_feats])
model_y = tf.placeholder(tf.float32, [None])
model_L_1 = VariationalDense(n_feats, n_hidden, model_prob, model_lam)
model_L_2 = VariationalDense(n_hidden, n_hidden, model_prob, model_lam)
model_L_3 = VariationalDense(n_hidden, 1, model_prob, model_lam)
model_out_1 = model_L_1(model_X, tf.nn.relu)
model_out_2 = model_L_2(model_out_1, tf.nn.relu)
model_pred = model_L_3(model_out_2)
model_sse = tf.reduce_sum(tf.square(model_pred - model_y))
model_mse = model_sse / n_samples
model_loss = (model_sse + model_L_1.regularization + model_L_2.regularization + model_L_3.regularization)

train_step = tf.train.AdamOptimizer(1e-3).minimize(model_loss)
with tf.Session() as sess:
    sess.run(tf.global_variables_initializer())
    for i in range(1000):
        sess.run(train_step, {model_X: X, model_y:y})
        if i % 100 == 0:
            mse = sess.run(model_mse, {model_X: X, model_y:y})
            print("Iteration {}. Mean squared error: {:.4f}.".format(i, mse))

    # sample from the posterior
    n_post = 1000
    Y_post = np.zeros((n_post, X_pred.shape[0]))
    for i in range(1000):
        Y_post[i] = sess.run(model_pred, {model_X: X_pred})


if __name__ == '__main__':
    plt.figure(figsize=(9, 7))
    for i in range(n_post):
        plt.plot(X_pred[:, 1], Y_post[i], "b-", alpha=1./200)
    plt.plot(X[:, 1], y, "r.")
    plt.grid()
    plt.show()

从实验结果图来看，对于已经训练的样本点，蓝色的实线部分都很密集，说明这部分方差小，模型对这些数据的预测比较有把握，而在预测两边这些未知数据时，蓝色实线明显呈现发散趋势，方差大，反应模型对这部分数据预测的不确定性。
模型的不确定性计算在风险评估，医疗诊断方面是很有意义的，毕竟并不是什么事情都不需要考虑代价，如果没有多大把握且代价很大，就会造成很严重的后果，这个时候评估模型对于预测的不确定程度就变得很重要。