不确定性量化的高精度数值方法和理论-HKBUMATH

英文引用格式:TaoTang,TaoZhou.Recentdevelopmentsinhighordernumericalmethodsforuncertaintyquantification.SciSinMath,2015中国科学:数学评审中稿件不确定性量化的高精度数值方法和理论汤涛¬,周涛­¬香港浸会大学数学系,香港­中科院数学与系统科学研究院,北京100190E-mail:[e-mail]，[e-mail]本文为《中国科学》庆贺林群院士80华诞专辑撰写摘要不确定性量化（UncertaintyQuantification,UQ）是近年来国际上热门的研究课题，其应用领域包括水文学、流体力学、数据同化、天气预测等等。由于UQ问题中的大量随机参数引起的超大计算量，如何设计高效的高精度数值方法变得非常重要，与其相关的计算技术和数学理论也引起人们的高度重视。本文将综述不确定性量化研究中的高精度数值方法和最新进展。我们主要讨论基于正交多项式的逼近方法，其中包括正交多项式Galerkin投影方法和随机配置方法。本文将侧重基于样本（数据）信息的随机配置方法，包括随机抽样，确定性抽样和结构随机样本。我们将重点介绍离散投影算法和压缩感知算法，并介绍相关数值分析进展，即如何确定样本的使用数量M与逼近空间基函数的自由度N的对应关系，以保证算法的稳定性和最优收敛性质。我们还将介绍高维空间中基于任意数量和任意位置节点的插值算法，以及一个相关的研究课题，即正倒向随机微分方程数值方法。最后，我们尝试探讨不确定性量化研究面临的挑战和亟待解决的研究问题。关键词不确定性量化，多项式逼近，随机配置方法，离散L2投影，压缩感知，正倒向随机微分方程MSC(2010)主题分类41A10,60H35,65C30,65C501引言不确定性量化（UQ）研究的是什么？首先，现实世界中许多实际问题的数学模型背后存在很大的不确定性，这些不确定性可能来自于问题中的参数、实验测量值、几何区域的复杂性等等。那么，如何通过量化这些不确定因素，以减少不确定性带来的风险，即是不确定性量化研究的主要目的。近几年，对带有不确定性的实际问题的建模与数值模拟受到了空前的关注，也就是人们常说的不确定性量化（UncertaintyQuantification，UQ）。UQ正是为了处理模型背后的不确定性。UQ方法可以分为概率框架和非概率框架，本文将主要探讨在概率框架下的建模及计算方法。概率论框架为处理不确定性因素的提供了良好的工具。通过使用概率论的工具，可以将大部分的不确定性建模成随机变量（或随机函数），于是，问题归结为如何求解所得到的随机问题。UQ研究的一个核心问题是如何处理由于随机输入带来的高维逼近。对于实际应用问题，比如油藏模拟，天气预测，数据同化，流体动力学，我们通常需要处理30-50维，甚至上百维的随机参数。因此，设计快速有效的不确定性量化方法变得十分迫切。UQ研究在欧美得到了很大的重视和发展，吸引了大量的应用数学家的关注，UQ数值方法也得到了很大的发展。至今，UQ已经在美国成为最重要的应用数学研究方向之一。美国能源部，空军，汤涛等:不确定性量化的高精度数值方法和理论国家实验室都设立专项经费，支持UQ方法的研究。2012年，美国工业与应用数学学会（SIAM）开始组织UQ年会（SIAMConferenceonUQ），两年举办一次，前两届会议每次都有500余人参会。2013年SIAM和美国数理统计协会（ASA）创立联合期刊（SIAM/ASAJournalonUQ），专门发表UQ领域的前沿研究成果。随着UQ研究的深入，对于随机数学模型的计算方法研究有了跨越性的发展，我们这里简单介绍其中的一部分：蒙特卡洛方法（MonteCarloMethods）。蒙特卡洛方法及其改进方法[1]是简单常用的基于样本的方法。在蒙特卡洛方法中，我们利用概率分布随机的产生一些样本，对每一个样本而言，所要解决的问题变成了一个确定的问题。通过求解这些确定问题，可以得到精确解的一些统计量信息，比如均值或者方差。蒙特卡洛方法实施简单，可以利用现成的程序，利于并行计算。但众所周知，蒙特卡洛方法的收敛率非常低：均值收敛率是1/√K（K是所使用的样本的数量），这意味着我们需要使用大量的样本才能得到较精确的数值结果。如果所对应的确定问题求解困难，这将是一个很大的挑战。于是很多改进的蒙特卡洛方法也就应运而生，比如Latin抽样方法[2,3],拟蒙特卡洛方法[4,5]等，但这些方法仍然有一定的应用局限性。摄动方法（PerturbationMethods）。这是一个比较流行的不基于样本的方法。该方法将一个随机函数在其均值附近展开成泰勒级数，然后取一个合理的截断。通常情况下，最多我们可以进行二阶展开的截断，因为对于更高阶的情形，得到的求解系统将会变得非常复杂。此方法被广泛的用于各种工程领域的应用问题[6,7]。然而，该方法还有一个固有缺陷，即对于不确定性的放大，因此一般只应用于小尺度的的随机输入问题，比如小于10%的随机扰动。矩方程方法（MomentEquationMethods）。在矩方程方法中，我们试图直接求解随机解的各阶矩所满足的方程。这些关于矩的方程需要从原始的随机问题出发推导得出。对于一些简单的问题，比如线形问题，该方法比较有效。但通常情况下，当我们推导某阶矩方程的时候，需要使用更高阶矩的信息。这就产生了所谓的“封闭（Closure）”问题，因此，我们通常需要对高阶矩做一些必要的假设。关于此方法的更细节的介绍可以参考文献[8]。多项式逼近方法（GeneralizedPolynomialChaos）。多项式逼近方法是近年来非常流行的计算方法，其基本思想就是将精确解在随机参数空间进行多项式展开。这种想法最早由维纳（Wiener）于1939年引入[9]，维纳使用Hermite多项式来处理带有高斯型随机参数的问题，该方法也被称为“HermiteChaos”，虽然Chaos的含义为“混沌”，我们这里仍将它翻译成多项式逼近，从而减少与动力系统中混沌定义的冲突。Ghanem和Spanos随后将这种逼近方法和空间上有限元方法结合，应用于带有不确定性输入的的动力学模型中[10]，取得了成功。近年来，这种方法得到了迅速发展。Xiu和Karniadakis基于Askey系统将多项式逼近方法推广[11]，从而处理带有任意类型随机参数输入的问题。比如使用Laguerre多项式处理Gamma随机参数输入，使用Jacobi多项式处理Beta分布的随机参数输入等等。并且，Xiu和Karniadakis也给出了相应的离散型随机参数及其对应的离散正交多项式。最初的多项式逼近方法采用Galerkin投影，也就是说，在随机参数空间对精确解做一个有限阶多项式展开，然后将此展开带入原问题，在展开多项式空间中做Galerkin投影，进而得到了一个关于展开系数的联立的方程组，通过求解方程组，我们获得精确解的全部信息。如果精确解对于随机参数有良好的正则性，此方法可以达到指数收敛。然而，求解联立的方程组并不容易，因为方程组的规模远大于原始随机问题的规模，并且我们需要更新现有的程序。而且，对于复杂的数学模型，所形成的联立方程组也将十分复杂，这对于多项式Galerkin投影方法是一个挑战。随机配置法（StochasticCollocation）。随机配点法是当今最流行的UQ计算方法之一，它结2中国科学:数学评审中稿件合了蒙特卡洛方法和传统多项式Galerkin投影方法的优点。随机配置方法试图通过计算一些特殊的样本信息，来构造高精度的多项式逼近（我们这里关注高阶随机配置方法，特别的，基于多项式逼近的随机配置方法）。这种方法避免了求解联立的方程组，只需要求解一定数量、与原随机问题同等规模的确定问题，便于实现，利于并行。随机配置方法最早的想法是基于张量积高斯点和张量多项式插值[12]。这就面临严重的维数祸根问题（Curseofdimensionality）：所需要的张量积节点个数随着维数的增加呈指数快速增长。考虑到对于每一个样本，需要求解与随机问题同等规模的确定问题，这在实际应用中是不能接受的。随后，高维求积中流行的稀疏节点[13,14]被引入到UQ的计算当中，以减弱维数祸根问题。这种方法仍然是一种张量方法，但是稀疏节点的设计试图减轻张量积节点对于维数的严重依赖。基于稀疏节点的插值方法在UQ计算中取得了很大的成功[15-18]。然而，随着UQ需求的不断深入，上述传统的张量型插值方法已经无法满足实际需求，随机配置方法也已经产生了多种变形，并将多个研究领域有机结合，产生了许多重要的数值分析结果，我们将在后面的章节中对此做详细介绍。在本文后面的章节中，我们将主要讨论基于多项式逼近的方法，包括Galerkin投影方法和随机配置方法。这其实可以把UQ问题转化一个随机参数空间的逼近论问题，即我们寻求数学模型（见下节描述）精确解u(x,Z)的多项式逼近：uN(x,Z)=N∑n=1bcn(x)ϕn(Z).其中x是空间变量，Z∈Γ⊂Rd是d维随机参数，ϕn(Z)是选取的多项式基函数。UQ的实际问题通常需要处理很高维数的随机参数，即d≫1。这使得展开项数N=N(d)也相当大。本文中，我们将主要探讨随机配置方法，并分成三种框架：离散最小二乘投影方法，压缩感知方法和基于任意分布高维节点的插值方法。基本思想就是通过一些样本数据构建以上的多项式逼近，假设我们已经获得（通过求解相关的数学模型，详细见后面的讨论）样本信息{u(x,zm)}Mm=1，其中{zm}Mm=1是在随机参数空间抽取的样本。通过样本信息，我们可以构建：N∑n=1ϕn(zm)bcn(x)≈u(x;zm),m=1,...,M.在离散投影方法的框架下，我们要求MN(d)。那么，一个自然并重要的问题是：最少需要多少样本才能保证算法的稳定性？因此，我们将着重讨论如何确定样本数量M与展开项数N(d)的依赖关系，以保证算法具有稳定性和最佳收敛性。特别的，我们将介绍相关的最新进展，包括随机抽样方法，确定性抽样结果，无穷区域问题Γ≡Rd等。在压缩感知的框架下，我们试图寻求精确解的稀疏逼近，这使得我们可以使用较少的样本，即MN(d)，我们同样关心算法的稳定性和最佳逼近性质。通过讨论，我们将呈现压缩感知方法和离散投影方法的联系和对应关系。最后我们讨论高维插值问题，即M=N(d)。为了适应实际需求，我们将强调“无结构”节点重要性，即我们希望使用任意位置，任意数量的节点。本文接下来的组织结构如下：第二节中我们介绍所讨论的模型问题：随机椭圆方程。第三节介绍Galerkin投影方法和随机配置方法。第四节至第六节是本文的重点，我们将系统的分别介绍离散投影方法，压缩感知方法和插值方法，并给出最新的算法和相关的数值分析结果，同时我们将给出一些数值算例。在第七节中，我们将介绍一个相关的研究课题：正倒向随机微分方程的数值方法。最后，我们将在第八节中探讨UQ研究的发展趋势和未来可能的研究方向。3汤涛等:不确定性量化的高精度数值方法和理论2模型问题这一节中我们讨论基于概率框架的UQ数学模型问题，我们将考虑UQ研究中的一个常用数学模型：随机椭圆问题。选取这个模型问题的原因是多方面的。首先，随机椭圆问题结构相对简单，并且有很多重要应用，比如随机振动，油藏管理，复合材料[19-24]；再次，随机椭圆问题的解具有许多特有的数学性质，因此产生了许多漂亮的数值分析结果。我们现在来给出随机椭圆问题−∇x·(κ(x,ω)∇xu)=f(x,ω),u|∂D=0,(2.1)其中D∈Rτ（τ=1,2,3）是一个物理区域。ω是定义在一个完备概率空间(Ω,F,P)上的事件，其中Ω是事件空间，F是σ-代数，P是概率测度。扩散系数κ(x,ω)和f(x,ω)是给定的随机函数，进而，方程的解u(x,ω)也是一个随机函数。为了保证问题的适定性，我们假定∃κmin0,s.t.P{ω∈Ω:κ(x,ω)κm