传统的不确定性量化数值方法

输入参数化

当系统的随机输入是系统参数时，参数化过程将很简单，因为输入已经是参数形式。我们关心的问题是这些系统参数之间是否独立， 具体而言：
令 $Y=(Y_1,\cdots ,Y_n),n>1$ 是系统参数, 其对应的分布函数为 $F_Y(y)=P(Y\le y)$, 我们希望寻找一组相互独立的随机变量 $Z=(Z_1,\cdots ,Z_n)\in {\mathbb{R}}^d,1\le d\le n$, 使得$Y=T(Z)$, 其中 $T$ 是一个合适的变换.

我们用一个简单的例子来解释上述事情。 考虑一个带两个随机参数的常微分方程：
$$\frac{du}{dt}(t,\omega )=-\alpha (\omega )u,u(0,\omega )=\beta (\omega ),$$
其中 $\alpha ,\beta$ 是随机的， 也就是说： 输入变量 $Y(\omega )=(\alpha (\omega ),\beta (\omega ))\in {\mathbb{R}}^2$.

如果 $\alpha ,\beta$ 是相互独立的， 则我们令 $Z(\omega )=Y(\omega )$. 则解 $u(t,\omega ):[0,T]\times \Omega \to \mathbb{R}$ 可以被表述为 $$u(T,Z):[0,T]\times {\mathbb{R}}^2\to \mathbb{R},$$ 即方程有一个时间维度和两个随机维度。

如果 $\alpha ,\beta$ 不是相互独立的， 则存在函数 $f$ 使得 $f(\alpha ,\beta )=0$. 则我们可以找到一个随机变量 $Z(\omega )$ 使得
$$\alpha (\omega )=a(Z(\omega )),\beta (\omega )=b(Z(\omega )),f(a,b)=0.$$
或者， 等价地， $\alpha ,\beta$ 相关可以推出存在函数 $g$ 使得 $\beta =g(\alpha ).$ 接着我们令 $Z(\omega )=\alpha (\omega ),\beta (\omega )=g(Z(\omega )).$ 此时， 问题的解便变为
$$u(t,Z):[0,T]\times \mathbb{R}\to \mathbb{R},$$
即只有一个随机维度。
实际问题中， 系统中有许多随机参数， 找到这些随机参数依赖的精确函数形式(即$f,g$) 是很困难的。我们能知道的信息仅仅是他们对应的(联合)概率分布函数， 于是我们的目标便是使用分布函数来找到变换 $T$.

高斯参数

对于高斯参数， 这件事情是比较容易的。 假设 $Y=(Y_1,\cdots ,Y_n)\sim \mathcal{N}(0,\mathbf{C}),$ 其中 $\mathbf{C}\in {\mathbb{R}}^{n\times n}$ 是协方差矩阵。 令 $Z\sim \mathcal{N}(0,\mathbf{I}),\mathbf{I}\in {\mathbb{R}}^{n\times n}$. $\mathbf{A}$ 是一个 $n\times n$ 的矩阵， 由高斯变换的性质我们知 $\mathbf{A}Z\sim \mathcal{N}(0,\mathbf{A}{\mathbf{A}}^{\mathbf{T}})$. 因此， 存在矩阵 $\mathbf{A}$ 使得 $\mathbf{A}{\mathbf{A}}^{\mathbf{T}}=\mathbf{C},$ 则 $Y=\mathbf{A}Z\sim \mathcal{N}(0,\mathbf{C}).$ 实际计算时， 我们可以采用 $Cholesky$ 分解来计算矩阵 $\mathbf{A}$. 具体有
$$\begin{array}{ll}{a}_{i1}=c_{i1}/\sqrt{c_{11}},& 1\le i\le n,\\ a_{ii}=\sqrt{c_{ii}-{\sum }_{k=1}^{i-1}{a}_{ik}^2},& 1<i\le n,\\ a_{ij}=(c_{ij}-{\sum }_{k=1}^{j-1}{a}_{ik}{a}_{jk})/a_{jj},& 1<j<i\le n,\\ a_{ij}=0,& i<j<n.\end{array}$$

非高斯参数

当系统参数是非高斯分布的时候， 直接参数化变变得很困难了。 Rosenblatt 给出了如下结果来理论上解决这一问题：
令 $Y=(Y_1,\cdots ,Y_n)$ 是一个随机向量, 其分布函数为 $F_Y(y)=P(Y\le y),$ 令 $Z=(z_1,\cdots ,z_n)=Ty=T(y_1,\cdots ,y_n)$ 为
$$\begin{array}{l}{z}_1=P(Y_1\le y_1)=F_1(y_1),\\ z_2=P(Y_2\le y_2\mid Y_1=y_1)=F_2(y_2\mid y_1),\\ \cdots ,\\ z_n=P(Y_n\le y_n\mid Y_{n-1}=y_{n-1},\cdots ,Y_1=y_1)=F_n(y_n\mid ,y_{n-1},\cdots ,y_1).\end{array}$$
则简单计算有
P ( Z i ≤ z i ; i = 1 , 2 , ⋯   , n ) = ∫ { Z ∣ Z i ≤ z i } ⋯ ∫ d y n F n ( y n ∣ y n − 1 , ⋯   , y 1 ) ⋯ d y 1 F 1 ( y 1 ) = ∫ 0 z n ⋯ ∫ 0 z 1 d z 1 ⋯ d z n = ∏ i = 1 n z i , \begin{array}{rl} P(Z_i&\leq z_i;i=1,2,\cdots,n)\\ &=\int_{\{Z\mid Z_i\leq z_i\}}\cdots \int d{y_n} F_n(y_n\mid y_{n-1},\cdots,y_1)\cdots d_{y_1} F_1(y_1)\\ &= \int _0^{z_n} \cdots \int_{0}^{z_1}dz_1\cdots dz_n =\prod_{i=1}^{n}z_i, \end{array} P(Zi​​≤zi​;i=1,2,⋯,n)=∫{Z∣Zi​≤zi​}​⋯∫dyn​Fn​(yn​∣yn−1​,⋯,y1​)⋯dy1​​F1​(y1​)=∫0zn​​⋯∫0z1​​dz1​⋯dzn​=∏i=1n​zi​,​
即 $Z=(Z_1,\cdots ,Z_n)\sim [0,1{]}^n$ 独立同分布.

尽管数学上此定理简单有力， 但在实际中基本上不可以使用，原因在于它依赖随机变量的条件概率分布，实际中这样的信息很少完全知道。况且，即使我们知道一些条件概率分布，也没有精确地表示形式。 故实际中我们会采用 Rosenblatt 的一些数值形式， 但这不属于我们讨论的范围。

随机过程与降维

在许多情况下， 随机输入往往是随机过程。 例如， 输入可以是随时间变化的随机强迫项， 此时它为时间上的随机过程； 或者是不确定的材料性质（导电率）， 此时它为空间上的随机过程。 参数化的问题可以被描述为下述形式：

令 $(Y_t,t\in T)$ 是一个随机过程, 求变换 $R$ 使得 $Y_t=R(Z)$, 其中 $Z=(Z_1,\cdots ,Z_d),d\ge 1$ 是相互独立的。
注意到， 指标集 $T$ 可以是时间区域或者空间区域，因而其通常是一个无限维的对象。 因为我们要求 $d$ 是有限的， 故上述变换不是精确的， 即 $Y_t\approx R(Z)$.

一种直接的方式便是考虑 $Y_t$ 的有限维版本， 这便要求将 $T$ 离散为有限维的指标进而考虑
$$(Y_{t_1},Y_{t_2},\cdots ,Y_{t_n}),t_1,t_2,\cdots ,t_n\in T.$$
$Y_t$ 的离散是一个逼近问题，显然，离散越密， 逼近效果越好。但是， 这会导致我们问题的维数变大，计算代价变大。 此处我们讨论三种方式：

Karhunen-Loeve Expansion

KL 展开是表达随机过程中最常用模型降维方法。具体而言， 对于任意一个随机函数 $\kappa (x,\omega )$, 我们定义其均值函数， 即数学期望：
$$\bar{\kappa }(x)=\mathbb{E}[\kappa (x,\cdot )]={\int }_{\Omega }\kappa (x,\omega )dP,x\in D.$$
类似地， 我们可以定义其 $\tau $阶矩($\tau\in \mathbb
N$)
$${\bar{\kappa }}_{\tau }(x)=\mathbb{E}[{\kappa }^{\tau }(x,\cdot )]={\int }_{\Omega }{\kappa }^{\tau }(x,\omega )dP,x\in D.$$
随机函数 $\kappa (x,\omega )$ 的协方差函数定义为
$$Co{v}_{\kappa }(x,y)=\mathbb{E}[(\kappa (x,\cdot )-\bar{\kappa }(x))(\kappa (y,\cdot )-\bar{\kappa }(y))],x,y\in D.$$
根据定义, $Co{v}_{\kappa }(x,y):\overline{D\times D}\to \mathbb{R}$ 是对称正定的协方差函数。 我们引入自共轭紧算子 $T_k:L^2(D)\to L^2(D):$
$$T_{k}v(\cdot )={\int }_{D}Co{v}_{\kappa }(x,\cdot )v(x)dx,\forall v\in L^2(D).$$
令 { λ i } i = 1 ∞ \{\lambda_i\}_{i=1}^{\infty} {λi​}i=1∞​ 和 { ϕ i } i = 1 ∞ \{\phi_i\}_{i=1}^\infty {ϕi​}i=1∞​ 分别为算子 $T_{\kappa }$ 的非负特征值序列和归一化正交特征函数， 即
$$T_{\kappa }{\varphi }_i={\lambda }_i{\varphi }_i,{\lambda }_1\ge {\lambda }_2\ge \cdots >0,⟨{\varphi }_i,{\varphi }_j⟩={\delta }_{ij},i,j\in {\mathbb{N}}^{+}.$$
进一步， 我们定义互不相关、且均值为0方差为1(i.i.d)的随机变量
$$\begin{array}{l}{Z}_i(\omega ):=\frac{1}{\sqrt{{\lambda }_i}}{\int }_D(\kappa (x,\omega )-\bar{\kappa }(x)){\varphi }_i(x)dx,i=1,2,\cdots ,\\ \mathbb{E}[Z_i]=0,\mathbb{E}[Z_i{Z}_j]={\delta }_{ij},i,j\in {\mathbb{N}}^{+}.\end{array}$$
那么，成立如下 Karhunen-Loeve 展开：
$$\kappa (x,\omega )=\bar{\kappa }(x)+\sum _{i=1}^{\infty }\sqrt{{\lambda }_i}{Z}_i(\omega ){\varphi }_i(x).$$
对上述的Karhunen-Loeve 展开作有限项截断， 我们便可以得到关于随机函数 $\kappa (x,\omega )$ 的有限维逼近
$$\kappa (x,\omega )\approx {\kappa }_d(x,\omega )=\bar{\kappa }(x)+\sum _{i=1}^d\sqrt{{\lambda }_i}{Z}_i(\omega ){\varphi }_i(x).$$
并且， 由 Mercer 定理， 我们有
$$\underset{d\to \infty }{\lim }\underset{x\in D}{\sup }\mathbb{E}[(\kappa (x,\cdot )-{\kappa }_d(x,\cdot ){)}^2]=\underset{d\to \infty }{\lim }\underset{x\in D}{\sup }\sum _{i=d+1}^{\infty }{\lambda }_i{\varphi }_i^2(x)=0,$$
并且，在均方意义下， 有限阶截断的 Karhunen-Loeve 展开具有最优误差估计。如果你想较为严格全面的了解KL展开的具体细节， 请参考文献 A brief note on Karhunen-Loeve expansion。

高斯过程与非高斯过程

利用 KL 展开， 我们可以用有限个随机变量来逼近一个随机场， 并且 { Z i ( ω ) } \{Z_i(\omega)\} {Zi​(ω)} 是不相关的。 对于高斯过程来说， 我们可以证明 $Z_i(\omega )$ 也服从高斯分布， 进而我们便知 { Z i ( ω ) } \{Z_i(\omega)\} {Zi​(ω)} 是相互独立的。 这是非常友好的。 对于非高斯过程而言， { Z i ( ω ) } \{Z_i(\omega)\} {Zi​(ω)} 不是相互独立的， 因此 KL 展开不能提供独立参数化。 但是实际中， 我们仍然采用 KL 展开， 并且将 { Z i ( ω ) } \{Z_i(\omega)\} {Zi​(ω)} 看作是相互独立的， 这当然时是不严格的。 此问题到目前仍然是个公开问题，有待进一步的讨论解决。

随机系统的参数化
考虑如下 PDE 系统：
{ u t ( x , t , ω ) = L ( u ) , D × ( 0 , T ] × Ω , B ( u ) = 0 , ∂ D × [ 0 , T ] × Ω , u = u 0 , D × { t = 0 } × Ω , \left\{\begin{array}{ll} u_t(x,t,\omega)=\mathcal{L}(u),& D\times (0,T]\times \Omega,\\ \mathcal{B}(u)=0,& \partial D\times [0,T]\times \Omega,\\ u=u_0,& D\times \{t=0\}\times \Omega, \end{array}\right. ⎩⎨⎧​ut​(x,t,ω)=L(u),B(u)=0,u=u0​,​D×(0,T]×Ω,∂D×[0,T]×Ω,D×{t=0}×Ω,​
其中 $\mathcal{L}$ 是(非线性)微分算子， $\mathcal{B}$ 是边界算子， $u_0$ 是初始条件， $\omega \in \Omega$ 定义了系统的随机输入。 因此解便是一个随机量
$$u(x,t,\omega ):\bar{D}\times [0,T]\times \Omega \to {\mathbb{R}}^{n_u},$$
$n_u\ge 1$ 是 $u$ 的维数。 由之前的知识我们知， 参数化后的系统为
{ u t ( x , t , Z ) = L ( u ) , D × ( 0 , T ] × R d , B ( u ) = 0 , ∂ D × [ 0 , T ] × R d , u = u 0 , D × { t = 0 } × R d , \left\{\begin{array}{ll} u_t(x,t,Z)=\mathcal{L}(u),& D\times (0,T]\times \mathbb{R}^d,\\ \mathcal{B}(u)=0,& \partial D\times [0,T]\times \mathbb{R}^d,\\ u=u_0,& D\times \{t=0\}\times \mathbb{R}^d, \end{array}\right. ⎩⎨⎧​ut​(x,t,Z)=L(u),B(u)=0,u=u0​,​D×(0,T]×Rd,∂D×[0,T]×Rd,D×{t=0}×Rd,​
其中 $Z=(Z_1,\cdots ,Z_d)$ 是一组独立的随机变量，用来参数化随机输入。我们举一个具体的一维随机椭圆问题：
$$\{\begin{array}{l}\nabla (\kappa (x,\omega )\nabla u)=f(x,\omega ),x\in (-1,1),\\ u(-1,\omega )=u_l(\omega ),u(1,\omega )=u_r(\omega ).\end{array}$$
这里 $\kappa ,f$ 是随机场, $u_l,u_r$ 是随机变量。 假定随机场 $\kappa,f $ 有 KL 展开：
$$\begin{array}{l}\kappa (x,\omega )\approx \tilde{\kappa }(x,Z^{\kappa })=\bar{\kappa }(x)+{\sum }_{i=1}^{d_{\kappa }}\sqrt{{\lambda }_i^{\kappa }}{Z}_i^{\kappa }(\omega ){\varphi }_i^{\kappa }(x),\\ f(x,\omega )\approx \tilde{f}(x,Z^f)=\bar{f}(x)+{\sum }_{i=1}^{d_f}\sqrt{{\lambda }_i^f}{Z}_i^f(\omega ){\varphi }_i^f(x).\end{array}$$
我们假定 { Z f } , { Z κ } , u l , u r \{Z^f\},\{Z^\kappa\}, u_l,u_r {Zf},{Zκ},ul​,ur​ 相互独立， 令 $Z=(Z_1,\cdots ,Z_d)=Z_1^{\kappa },\cdots ,Z_{d_{\kappa }}^{\kappa },Z_1^f,\cdots ,Z_{d_f}^f,u_l,u_r),$ 其中 $d=d^{\kappa }+d^f+2,$ 则之前的问题可以表述为
$$\{\begin{array}{l}\nabla (\tilde{\kappa }(x,Z)\nabla u)=\tilde{f}(x,Z),x\in (-1,1),\\ u(-1,Z)=Z_{d-1},u(1,Z)=Z_d.\end{array}$$
方程的解为 $u(x,Z):[-1,1]\times {\mathbb{R}}^d\to \mathbb{R}.$

几种传统的数值方法

我们以之前的常微分方程为例，即
$$\frac{du}{dt}(t,\omega )=-\alpha (\omega )u,u(0,\omega )=\beta (\omega ).$$
上述方程的精确解为 $$u(t,\omega )=\beta (\omega )e^{-\alpha (\omega )t}.$$
如果 $\alpha (\omega ),\beta (\omega )$ 是相互独立的， 则有 $\mathbb{E}[u(t,\omega )]=\mathbb{E}[\beta ]\mathbb{E}[e^{-\alpha t}].$

Monte Carlo Sampling

Monte Carlo 方法. Monte Carlo 方法及其改进方法是简单常用的基于样本的方法. 在 Monte Carlo 方法中, 我们利用概率分布随机地产生一些样本, 对每一个样本而言, 所要解决的问题变成了一 个确定的问题. 通过求解这些确定问题, 可以得到精确解的一些统计量信息, 如均值或者方差. Monte Carlo 方法实施简单, 可以利用现成的程序, 利于并行计算. 但众所周知, Monte Carlo 方法的收敛率非 常低: 均值收敛率是 $1/\sqrt{K}$ ($K$ 是所使用样本的数量), 这意味着我们需要使用大量的样本才能得到较精确的数值结果. 如果所对应的确定问题求解困难, 这将是一个很大的挑战. 具体而言：

• 产生独立同分布的随机实现 { Z i = ( α i , β i ) } ; \{Z_i=(\alpha_i,\beta_i)\}; {Zi​=(αi​,βi​)};
• 求解方程 $u^i(t)=u(t,Z^i);$
• 计算统计量信息 $\bar{u}(t)=\frac{1}{M}{\sum }_{i=1}^{M}u(t,Z^i)\approx \mathbb{E}[u].$

Moment Equation Approach

令 $\mu (t)=\mathbb{E}[u],$ 在方程两边对 $\omega$ 积分有，
$$\frac{d\mu }{dt}(t)=-\mathbb{E}[\alpha u],\mu (0)=\mathbb{E}[\beta ].$$
上述方程需要知道 $\mathbb{E}[\alpha u]$ 的知识， 故我们两边乘 $\alpha$ 再次积分有
$$\frac{d}{dt}\mathbb{E}[\alpha u](t)=-\mathbb{E}[{\alpha }^{2}u],\mathbb{E}[\alpha u](0)=\mathbb{E}[\alpha \beta ],$$
这便出现了一个新的量 $\mathbb{E}[{\alpha }^{2}u]$, 我们对原始系统乘 ${\alpha }^2$ 再积分有
$$\frac{d}{dt}\mathbb{E}[{\alpha }^{2}u]=-\mathbb{E}[{\alpha }^{3}u],\mathbb{E}[{\alpha }^{2}u](0)=\mathbb{E}[{\alpha }^2\beta ].$$
我们定义 ${\tilde{\mu }}_k(t)=\mathbb{E}[{\alpha }^{k}u],$ 接着我们有
$$\frac{d}{dt}{\tilde{\mu }}_k(t)=-{\tilde{\mu }}_{k+1}(t),{\tilde{\mu }}_k(0)=\mathbb{E}[{\alpha }^k\beta ].$$
即原始的系统是不封闭的，北京大学的李若老师在这方面做了很多杰出的工作。最简单的做法是固定一个 $k$, 假定高阶矩可以被低阶矩表示，即
$${\tilde{\mu }}_{k+1}=g({\tilde{\mu }}_0,\cdots ,{\tilde{\mu }}_k).$$

Perturbation method

当 $ϵ=O(\alpha (\omega ))\sim {\sigma }_{\alpha }<<1,$ 我们可以考虑解的Taylor展开：
$$u(t,\omega )=u_0(t)+\alpha (\omega )u_1(t)+{\alpha }^2(\omega )u_2(t)+\cdots ,$$
代入原方程， 我们有
$$\frac{d{u}_0}{dt}+\alpha \frac{d{u}_1}{dt}+{\alpha }^2\frac{d{u}_2}{dt}+\cdots =-\alpha (u_0+\alpha u_1+{\alpha }^2{u}_2+\cdots ).$$
因为 $O({\alpha }^k)={ϵ}^k$， 于是我们有
$$\begin{array}{ll}O(1):& \frac{d{u}_0}{dt}=0,\\ O(ϵ):& \frac{d{u}_1}{dt}=-u_0,\\ O({ϵ}^2):& \frac{d{u}_2}{dt}=-u_1,\\ \cdots & \cdots \\ O({ϵ}^k):& \frac{d{u}_k}{dt}=-u_{k-1},\end{array}$$
且 $$u_0(0)=\beta ,u_k(0)=0,k>1.$$
简单计算有
$$u_0(t)=\beta ,u_1(t)=-\beta t,\cdots ,u_k=\beta (-1{)}^k\frac{t^k}{k!}.$$
故 $$u(t,\omega )=\sum _{k=0}^{\infty }\beta \frac{-t^k}{k!}{\alpha }^k(\omega )=\beta exp(-\alpha t).$$
通常情形下, 最多我们可以进行二 阶展开的截断, 因为对于更高阶的情形, 得到的求解系统将会变得非常复杂. 此方法被广泛地用于各 种工程领域的应用问题. 然而, 该方法还有一个固有缺陷, 即对于不确定性的放大, 因 此一般只应用于小尺度的的随机输入问题, 如小于 10% 的随机扰动.