贝叶斯公式及其在反问题中的应用

贝叶斯公式及其在反问题中的应用1.1 反问题背景有这样一个“盲人听鼓”的问题：蒙上一个人的双眼，让他听鼓的敲击声音来判断这个鼓的形状大小，可能吗？生活经验告诉我们，这也许是可能的。如果一个鼓的形状大小确定了之后，那么它的声音也就随之确定了；如果已知一个鼓的声音，那么能不能反过来确定这个鼓的形状和大小呢？这便是反问题所要研究的范畴。以上这个问题最早是由荷兰物理学家Lorentz以射线理论为背景在1910年提出来的。我们知道，一个鼓的音色可以由它的固有频率来确定，各种鼓的音色综合起来就构成了一串频率谱。“盲人听鼓”这个问题就是想要通过鼓发出的声音的频率来反推鼓的形状和大小等具体情况。经过数学家们一个多世纪的研究发现：根据鼓声，人们确实能得到一些关于鼓的形状的信息并给出了相应的计算公式。例如，鼓的面积可以通过小于的谱数来确定：.但是，这个问题是直到1992年才得到真正解决的。科学家们构造出了两个音色相同，但是形状不同的鼓，从而证明了人们不能仅由鼓的音色就准确判断出鼓的形状和大小，即“盲人听鼓”这个反问题是没有唯一解的。这个经典的问题反映出反问题研究中一个基本的困难，即反问题的不适定性。目前，由于计算机技术的迅猛发展，反问题的研究也突飞猛进，它已成为包含物理学、生物化学、经济学等一系列学科的多学科交叉领域。但是，反问题的研究仍然面临着许多难点，比如上面提到的不适定性。对于反问题的求解，确定性正则化方法已经趋于完善，贝叶斯正则化方法则正处于起步阶段，所以，本文主要讨论了反问题及其贝叶斯求解方法。1.2 反问题的定义 下面我们从数学的角度来理解反问题的定义。定义1.2.1（Banach空间）如果赋范线性空间的度量空间是完备的，即任何柯西列在其中都收敛，则称这个赋范线性空间为Banach空间。记和为两个Banach空间，分别称为“输入空间”，为“输出空间”，假定有一个算子：将“输入空间”映射到“输出空间”，即,则由给定的输出来确定输入或者算子的问题就构成一个反问题。1.3 求解反问题的两种方法求解反问题的最大困难在于大多数反问题都是不适定的，即不满足以下三条中的至少一条：1.存在性：即给定，是否存在相应的解。2.唯一性：即如果相应的解存在，那么是否唯一。3.稳定性：即算子是否连续，即对于的微小变动，相应的解是否也变动较小。 有许多适定性的数学物理正问题，从它们中衍生出来的反问题却是不适定的。例如双曲方程Cauchy问题，这本身是数学物理方程中一个经典的适定问题，但是从它衍生出来的反问题却是不适定的。又比如法国数学家Hadamard给出的椭圆方程Cauchy问题。6除了不适定性之外，求解反问题的困难还有非线性性。有的时候，即使正问题是线性的，它所对应的反问题也有可能是非线性的，这就为求解问题带来了更多困扰。为了求解非线性的反问题，我们通常需要先将问题线性化，然后反复进行迭代，在高维的时候，这种计算量是十分巨大的。从而，如何研究探索出一个能够高效求解反问题的算法也是当今反问题的研究热点之一。2综上我们知道，反问题同时具有不适定性和非线性性的特点。本文主要对不适定的反问题进行讨论。不适定问题通常是不稳定的，即所观测的数据的微小变化，会造成其解的巨大改变。在实际问题之中，观测的数据常常包含误差，因而造成不适定问题的求解变得十分困难。从数学的角度看，假设有算子将映为，即满足.并假设的逆算子不连续，且等式右端受到扰动变为，当满足时，由于的不连续性，我们不能简单地定义对应的解为.从而，我们应该修改逆算子的定义，即引入正则化方法。正则化的理论最初是由前苏联数学家Tikhonov等学者首先提出来的，之后大量的数学家们对这个理论进行了改进和完善。正则化方法可以分为两大类：确定性正则化方法和随机方法。如今，确定性正则化方法的理论体系已经十分完善，但是随机方法仍处于起步阶段，从而吸引着许多数学家的关注。所以本文主要讨论不适定问题的贝叶斯正则化方法。7所谓正则化方法，就是用一组与原不适定问题相“近似”的适定问题的解来逼近原不适定问题的解。沿用上面定义的各个变量，我们用数学的语言来描述正则化方法想要求解的问题：对于有界线性算子，当测量数据满足条件时，我们希望通过正则化方法来求解如下方程：. 因为这个问题是不适定的问题，所以我们不能保证解的存在性，唯一性和稳定性。因此要想利用算子在一般意义上的逆算子来求解方程是不实际的，这通常不能得到令我们满意的解。求解这类问题的一般方法是构造一簇有界算子去逼近的逆算子，即如下定义的正则化方法。定义1.3.1（正则化方法）若依赖于参数的算子，满足：,且满足(1) 对任意以及，都存在正数，使得算子有意义。(2) 对任意，都存在和依赖于的参数，使得算子将的邻域映到的邻域内，并记做。则称为算子的正则逆算子；为正则解；当时，正则解就可以逼近我们要求的精确解；在这里，称为正则化参数；这样的求解方法就称为正则化方法。 1.3.1 确定性正则化方法从上面的定义可知，一个正则化方法是由正则化算子和正则化参数的选取共同决定的。对与正则化算子，一方面要求它有界，从而使得定义的逼近解能够连续依赖于噪音数据，另一方面还要求当噪音逐渐减小时，逼近解能够收敛到精确解。 常用的处理不适定问题的确定性方法是Tikhonov正则化方法。它的基本思想是结合光滑性，有界性等一些先验信息来构造不适定问题的稳定近似解。1 1.3.2 贝叶斯正则化方法 从统计的意义上来看，确定性正则化方法得到的结果是一个点估计，其求解过程忽略了误差的随机性以及模型的不确定性。与此形成对比的是，贝叶斯正则化方法的最突出特点就是将模型中所有的量都看成是随机变量，将解的先验信息看成是概率分布，将它以概率分布函数的形式表现出来，再通过贝叶斯公式，最终得到关于整个空间的条件后验概率分布，而不仅仅是确定性正则化方法所得到的单点估计。通过贝叶斯正则化方法得到的解的后验概率分布之后，利用后验概率分布函数，我们就能得到有关解的各种统计量，比如均值，方法，置信区间等。 利用贝叶斯方法求解反问题的关键点在于如何选择先验概率分布，从而使得该分布能够尽可能完整详细地描述未知参数的全部信息。先验概率分布能够起到正则化参数的效果。在贝叶斯统计中，均匀分布是一种可以选择的先验分布。但是均匀分布通常不能提供足够的信息，达不到很好的正则化效果。另外常见的选取方式还有马尔科夫随机场以及高斯过程等。在许多情况下，这些先验分布往往是主观选择的。 贝叶斯方法的一个难点在于它的巨大的计算量。为了探求后验概率分布函数，我们常常利用马尔科夫链蒙特卡罗抽样方法(Markov Chain Monte Carlo Method)进行大约次的随机抽样。对于每一次随机抽样，都需要重复进行一次似然函数的求解过程。对于偏微分方程的反问题，每一次似然函数的求解就需要计算一次相应的正问题，这样就需要消耗大量的计算时间。为了克服这个难点，设计相应的快速方法在贝叶斯方法处理数学物理反问题中，特别是高维非线性反问题中，起到了至关重要的作用。21.4 本文主要工作本文主要介绍贝叶斯正则化方法的理论和研究思路，并探索相应的贝叶斯随机化算法，以及算法的加速和改进。本文第一章简要介绍了反问题的背景，定义以及求解反问题两种常用的方法。 本文第二章主要介绍了贝叶斯公式，贝叶斯推断，并给出了贝叶斯推断的三种基本模型，最后以一个实际问题的模型选择方法结尾。本文第三章介绍了马尔科夫链蒙特卡罗模拟方法，并给出了其相应的算法和几种常用的取样方法，最后以一个实际的例子做概括。第二章 贝叶斯正则化方法在本章中，我们给出贝叶斯推断的基本思想和计算方法。贝叶斯公式是贝叶斯推断最基础的部分，从而我们首先给出贝叶斯公式的定义。2.1 贝叶斯公式定义2.1.1（贝叶斯公式）两个随机变量和，在随机变量给定时，随机变量的条件概率由以下公式给出： . （1） 在实际问题的贝叶斯推断中，条件密度函数和都是已知的，即在给定时，的分布和的先验分布都已经预先给定。继续沿用之前的假设，我们考虑如下的反问题： . （2）这里代表未知量，代表观测数据。从贝叶斯的角度来看，我们将未知量和观测数据都看做随机变量，并基于经验得出在给定时的条件概率分布，和的先验概率分布。例如，如果我们假设：在给定时，服从一个均值为，方差为的一个高斯分布，那么我们有.的先验概率代表了我们对未知数已有的经验。例如，如果我们已经了解到解十分稀疏，即它有很多零项，那么我们将视为拉普拉斯分布就是十分恰当的，即：.最终求得的后验分布概率密度函数将会包含反问题的解的全部信息。我们只需通过分析这个后验分布函数就能得到我们想知道的解的全部相关信息。2.2 贝叶斯推断 我们把记做，即：,并且我们记：. 因为观测数据的概率密度为常数，从而这个记号表示和是成比例的。 对于基于反问题的模型，我们有服从于条件. 在这里，是一个度量映射，是状态和参量的一个偏微分方程约束。我们直接给出对于这个约束的模型误差和测量误差的联合分布公式如下：. 在这里我们假设误差是服从高斯分布的，从而常数度量了模型的不确定性，常数度量了数据的不确定性。3 除了基于参数的先验分布，我们还会选择一个基于状态的先验分布，从而我们能得到如下的后验分布概率密度函数：这里对应的贝叶斯公式（1）中的分母的定义为： 或者. （3）在这里，叫做的边际密度函数。在上面这个例子中，边际密度函数是在模型类型给定的情况下数据出现的概率，但是并没有假定任何特定的模型参数。代表着为了解释数据而给定的模型的合理性，它通过贝叶斯因数构成了模型选择的基础。这说明，人们可以在依据参数的后验分布得出来的一簇可能的模型中，通过先验模型在公式（3）中对应的势来给先验模型定阶。例如，有如下两个备选模型：.一般来说，我们会通过在所有先验模型上使模型证据最大化的方法，来选择最优模型。2.3 贝叶斯推断的模型 似然函数和先验分布是贝叶斯模型中的两个主要而又最基本的内容。似然函数包含了观察数据的全部信息，也可以说是包含了带有噪音干扰的观察数据中噪音的全部信息；先验分布的概率密度函数则使我们在推断之前就得到了与解相关的一些先验信息。 下面，我们就这两部分分别进行讨论，并给出响应的模型。 2.3.1 噪音模型 在所有实际的反问题中，观察的数据中带有噪音是最突出而困难的一个问题。似然函数完全由噪音的统计量决定。除此之外，还有以下几种原因可能导致模型或者数据中产生误差： (a)如果我们用符号代表我们得到的模型，代表实际模型，代表模型误差，那么显然我们能够得到：. 这个模型误差产生的原因是我们用简化了的模型来替代原本复杂的模型。比如，在物理应用上，我们用连续模型替代完全电极模型。(b)如果我们用符号代表离散误差，那么我们能够得到：. 在实际的反问题算法中，这个模型必须用有限维的模型来近似。离散误差通常可以根据的大小来进行粗略的估计。(c)如果我们用符号代表人们通过测量得到的数据，用代表实际产生的数据，代表数据误差，那么我们能够得到：. 在这里，产生数据误差的一个原因是测量数据时的固有误差，比如测量器械的灵敏度等一些问题；另外一个可能产生数据误差的原因是测量的不准确；除此之外，数据误差还有可能是在数据传播过程中产生的，例如，在嘈杂的频道中传播的数据常有可能带有十分严重的误差问题。在实际问题中，因为各种各样原因产生的误差都会集中在我们最终得到的数据中，因而如何从带有误差的观察数据中提取我们需要的有用信息就变得至关重要。现在最常用且简便的噪音模型是加性高斯模型：.这里是维空间中的一个随机变量，并且它的所有分量都是独立同分布的（i.i.d.），均值为零，方差为的高斯分布。此外，如果噪音与真实的观测数据相互独立，那么它的似然函数可以由下式给出：.在下面的公式中，我们将把简记为，叫做反方差。 2.3.2 先验模型从概率统计的角度来看，先验概率密度函数包含了解的全部先验信息。一般来说，先验信息是通过历史资料，或者已有的统计知识等得到的。由于先验信息匮乏和噪音等问题，我们之前也已经陈述过，反问题通常是不适定的，从而如何尽最大程度结合所有可利用的先验信息就变得十分重要。在贝叶斯模型中也是如此，并且先验在贝叶斯模型中扮演着正则化的作用。因此我们可以这样说，先验模型是贝叶斯模型的核心与关键。马尔科夫随机场是一个非常常用的先验模型，用公式表示如下：.这里，是指示随机场中各部分之间能量相互作用的势函数，并且可以是任何正则化函数；这里是一个标量参数，指示出局部作用的强度，它起到正则化参数的作用，从而的选取十分重要。 2.3.3 分层模型似然函数和先验都可以包含附加的未知参数，例如：和.这里代表精确度，代表标量参数。我们一般将这些参数叫做超参数。在贝叶斯估计中，马尔科夫随机场的先验概率密度函数中的会影响后验分布的密度函数，因此，它会影响后验分布密度函数的点估计和置信区间等性质。分层贝叶斯模型提供了一个自动选择这些参数的一个简便有效的方法：我们将标量参数和反方差看做为随机变量，并通过观测数据来决定它们。利用共轭分布是选择这些先验的一个简便的方法，这样会使得后验分布和先验分布在同一簇中。对于和而言，共轭分布由Gamma分布给出：,.这里的参数对和决定着参数和的先验信息的广度。对于参数来说，一个不给出任何信息的先验分布也是可以接受的，这相当于在上式中令。对于反方差来说，通过重复试验是能够得到一个相当合理的估计的，因此人们可以使用比较精确的先验分布。一旦我们选好了先验参数对和，我们就能用（1）中的贝叶斯公式来计算带有这些参数的后验分布的概率密度函数：. 2.3.4 贝叶斯推断中的创新之处后验分布概率密度函数代表了（2）中反问题的完整的贝叶斯解，并且包含了这个反问题的解的全部信息。与确定性正则化方法相比较，贝叶斯方法有如下的特点：首先，贝叶斯方法的解是一个概率分布，从而这个解包含了所有与观察数据一致的合理的解，例如我们可以考虑均值和方差：.特别地，我们可以将均值看作具有代表性的反问题的解，方差让我们能够量化一个特定解的不确定性。例如，的第个分量的的置信区间由以均值为中心的区间给出，即：.从而我们能够得到：.这里表示。对于超参数和来说，我们可以将它们固定为一些选定的值，或者直接将它们排除在外。相比之下，决定的技巧是从整体中孤立出一个解。例如，Tikhonov正则化方法只着眼于寻找使得后验分布最大化的点，而完全忽略了解在统计上的浮动性，这就很容易使人们误入歧途。其次，恰当的统计模型对于提取出观测数据中的全部信息来说是十分必要的。例如，对于脉冲噪音来说，拉普拉斯分布或者分布都可以作为合适的似然函数，因为它们对于有大量异常点的观测数据而言仍然是稳健的。最后，贝叶斯推断提供了一个灵活的正则化方法，因为分层模型能解决选择一个合适的正则化参数的非平凡问题。因为有这些特点，贝叶斯推断在近年来得到了数学工作者们的许多关注。然而，就像我们之前谈到的那样，贝叶斯推断需要大量的计算，这对数学工作者们来说同样是一个不小的挑战。2.4 实际问题的模型选择通常来说，对于相同一组观测数据，我们会有许多可以用来解释它们的模型。我们建立一个模型，这里的这个可以指代不同的先验模型，或者指代依赖于一个参数的一簇先验模型，比如。现在就有一个重要的问题摆在我们面前：我们应当怎样选择恰当的模型来解释观测数据呢？依赖于模型选择的边缘概率包含了给定模型的很多有价值的信息。特别地，它衡量了模型解释观测数据的能力的大小。为了特别指出边缘概率密度对模型的依赖，我们将它记做。在给定了噪音数据的情况下，模型选择的贝叶斯方法是选择使后验概率最大的那个模型。由贝叶斯公式（1）得知，模型的后验概率由下式给出：. 这里是观测数据的边际似然函数，是模型的先验概率函数。 以下我们介绍两个简单例子的模型：在这里我们分别用和来表示，并用贝叶斯公式在和之间选择较优的模型。 首先，我们给出贝叶斯因数的定义如下： 这里的数值可以用来衡量模型与模型对观测数据的拟合优劣的比较。如果因数，那么这意味着模型对观测数据的拟合较优；如果，那么这意味着模型对观测数据的拟合较优；如果，那么这意味着两个模型对观测数据的拟合情况相同。 除此之外，贝叶斯因数的大小是指示模型与模型的拟合情况相比较的程度。根据Kass and Raftery，当贝叶斯系数超过320时，我们就能据此说：模型对观测数据的拟合远优于模型。4,5 现在我们来考虑多元模型的情况。如果所有的备选模型都是相似的，那么在观测数据给定时，要使模型的后验概率最大化就等价于要使以下的边际似然函数最大化：. 这里为模型的参数向量，为的后验概率分布密度函数。 为了得到更加简洁有效的法则，我们用拉普拉斯方法来估计边际似然。我们首先指出： 接下来我们在后验周围展开，这就是MAP估计；在这里后验概率分布密度函数如果存在，就能达到它严格的最大值，并且后验的对数值如果存在，则也会达到它严格的最大值。然后我们就可以在周围用对数的泰勒公式将展开：.这里是Hessian矩阵 黑塞矩阵（Hessian Matrix），又译作海森矩阵、海瑟矩阵、海塞矩阵等，是一个多元函数的二阶偏导数构成的方阵，描述了函数的局部曲率。黑塞矩阵最早于19世纪由德国数学家Ludwig Otto Hesse提出，并以其名字命名。黑塞矩阵常用于牛顿法解决优化问题。通过的最优性我们知道，一定为负值，并且。我们记，从而用下式来估计：这个过程叫做拉普拉斯方法。如果数据的观测量足够大的话，那么这个后验分布就会在周围集中，从而我们的这个近似就是合理的。因此，我们能够得到：.为了更好地理解这个近似的含义，我们考虑的情况。这样的话，后验最大概率与最大似然估计(MLE)就会完全相同。不仅如此，Hessian中的每一项都能用如下式子来表达：.这说明Heissian矩阵是所观测得到的Fisher信息矩阵。不仅如此，如果我们观测到的数据是独立同分布的，且观测到的数据量足够大，那么弱大数定律告诉我们：.因此，观察到的Fisher信息矩阵可以化简为： 所以我们有，这里是参数的值，与单个数据点的Fisher信息矩阵一致。 总结来说，现在我们有： 对于足够大的， 我们可以忽略与无关的项，而只保留下与有关的项，从而可以得到：.第三章 马尔科夫链蒙特卡罗方法MCMC 我们之前介绍过先验概率密度函数，下面我们将简记为，这个分布本身并不能给我们提供直接的信息。因此，计算摘要统计量是十分必要的，例如，均值和方差：.下面，我们将均值记做，将方差记做。这些都是高维的积分，并且因为维数的原因，标准的正交规则对此并不适用。例如，对于一个一百维空间中的概率密度函数，在每个方向上有两个正交点的张量乘积类型的正交规则需要估计这么多的点，这个数值实在是太大了。当前避免这个计算问题的常用方法是蒙特卡罗方法，特别是马尔科夫链蒙特卡罗(MCMC)方法。下面我们就MCMC方法具体展开讨论。3.1 蒙特卡罗模拟 蒙特卡罗模拟的基本步骤是从定义在高维空间的目标分布中抽取一系列独立同分布的样本。然后再计算样本均值： 然后通过样本均值来估计任意函数的期望。 我们通过下面的式子来定义蒙特卡罗积分误差：.从而偏差为，均方误差根为。当充分大时，根据中心极限定理，我们能由下式得到蒙特卡罗积分误差为：.这里为服从标准正态分布的随机变量，为以为自变量的函数的方差。下面我们再给出一个更加精确的算式： 因此，蒙特卡罗积分中的误差是阶数为的常数，并且仅依赖于被积函数的方差。不仅如此，误差的统计分布粗略来看就是正态分布。并且我们指出这个估计与空间的维度无关。 对于我们感兴趣的非线性反问题和非高斯模型，从高维联合分布中生成一系列独立同分布的样本这种做法是非常典型的，并且是高度非平凡的。因此，对高级蒙特卡罗方法的研究多种多样。对我们有用的想法是重要性样本。假设是一个接近于后验分布概率密度函数的对于样本来说比较简单的概率密度分布函数。那么我们就可以据此估计函数的期望： . （4）这里的独立同分布的样本是从我们定义的辅助分布中抽取的，并且它的重要性权重下式给出：.在（4）式的估计中，我们假定密度为正态分布的，否则正则化常数就应该由下式计算出来：.3.2 MCMC算法 MCMC算法是探索后验分布密度函数最常用的方法。在这部分中，我们介绍MCMC算法的基本知识。 MCMC抽样的想法最早是由Metropolis提出来的，当时是作为能有效模拟原子的能量水平的一种方法，后来经由Hastings改造并一般化为解决统计问题的方法。9,10这个方法的基本思想十分简单，即给定一个复杂的目标分布函数，在状态空间上，我们构造一个非周期的，不可约的马尔科夫链，使得它的稳定分布为。通过运行足够长的马尔科夫链，从这个马尔科夫链中模拟的数值就能被看做是从目标分布中抽取的相关样本，从而能被用来计算摘要统计量。 Metropolis-Hastings算法是最基本的MCMC方法。Metropolis-Hastings算法如下：Metropolis-Hastings 算法1：赋初值，设定N；2：for i=0：N do3： 样本；4： 样本；5： if ， then6： ；7： else8： ；9： end if10：end for 在这里面，是从标准正态分布中产生的随机变量，是目标分布，是样本的一个建议分布。从分布中产生了一个新的状态之后，我们将这个新生成的点看做马尔科夫链的新状态，这个新状态的概率由下式给出：. 然而，如果我们拒绝状态，那么这个马尔科夫链保持当前的状态。我们指出目标分布只有通过比率才通过进入算法，所以我们只要知道分布中一个关于乘积的常数，对实现程序就足够了。不仅如此，在建议分布是对称的情况下，即，接受概率函数可以减小到： . （5）这个Metropolis-Hastings算法保证了马尔科夫链对任何合理的建议分布都收敛到