Certified Offline-Free Reduced Basis (COFRB) Methods for Stochastic Differential Equations Driven by Arbitrary Types of Noise

Abstract

In this paper, we propose, analyze, and implement a new reduced basis method (RBM) tailored for the linear (ordinary and partial) differential equations driven by arbitrary (i.e. not necessarily Gaussian) types of noise. There are four main ingredients of our algorithm. First, we propose a new space-time-like treatment of time in the numerical schemes for ODEs and PDEs. The second ingredient is an accurate yet efficient compression technique for the spatial component of the space-time snapshots that the RBM is adopting as bases. The third ingredient is a non-conventional “parameterization” of a non-parametric problem. The last is a RBM that is free of any dedicated offline procedure yet is still efficient online. The numerical experiments verify the effectiveness and robustness of our algorithms for both types of differential equations.

Appendices

Appendix A: Implementation of Algorithm 2

In this appendix, we present some details for efficient implementation of Algorithm 2 whose steps 5, 6, and 10 are key.

1.1 A.1 Step 5

For the forward Euler time discretization, we rewrite the right hand side of (3.8) in the matrix form

$$\begin{aligned} \vec {f}_\alpha&=\left( \begin{array}{c} u_\alpha ^0+\Delta t u_\alpha ^0+\Delta t \sum _{\varepsilon _k \le \alpha }m_k^0u_{\alpha -\varepsilon _k}^0\\ 0\\ \vdots \\ 0\\ \end{array}\right) +\Delta t \sum _{\varepsilon _k \le \alpha }\left( \begin{array}{cccc} 0 &{} &{} &{} \\ &{}m_k^1 &{} &{}\\ &{} &{} \ddots &{}\\ &{} &{} &{} m_k^{n-1} \end{array}\right) \left( \begin{array}{c} 0\\ u_{\alpha -\varepsilon _k}^1\\ \vdots \\ u_{\alpha -\varepsilon _k}^{n-1} \end{array}\right) \nonumber \\&=\vec {f}_{\alpha }^0+\Delta t \sum _{\varepsilon _k \le \alpha } M_k^0 \vec {U}^0 \vec {c}_{\alpha -\varepsilon _k} \end{aligned}$$

(6.1)

where $M_k^0=\text {diag}\{0, m_k^1, \ldots ,m_k^{n-1}\}$ and

$$\begin{aligned} \vec {U}^0&=\left( \begin{array}{cccc} 0 &{}0 &{}\ldots &{}0\\ u_{\alpha _1}^1 &{}u_{\alpha _2}^1 &{}\ldots &{}u_{\alpha _i}^1\\ \vdots &{}\vdots &{} \vdots &{}\vdots \\ u_{\alpha _1}^{n-1} &{}u_{\alpha _2}^{n-1} &{} \ldots &{} u_{\alpha _i}^{n-1} \end{array}\right) =\left( \begin{array}{cccc} \vec {U}_{\alpha _1}^0&\vec {U}_{\alpha _2}^0&\ldots&\vec {U}_{\alpha _i}^0 \end{array}\right) . \end{aligned}$$

(6.2)

For the Crank–Nicolson method, we can also rewrite (3.10) similarly

$$\begin{aligned} \vec {f}_\alpha&=\left( \begin{array}{c} u_\alpha ^0+\frac{1}{2}\Delta t u_\alpha ^0+\frac{1}{2}\Delta t \sum _{\varepsilon _k \le \alpha }m_k^0u_{\alpha -\varepsilon _k}^0\\ 0\\ \vdots \\ 0\\ \end{array}\right) +\frac{1}{2}\Delta t \sum _{\varepsilon _k \le \alpha }\left( \begin{array}{cccc} 0 &{} &{} &{} \\ &{}m_k^1 &{} &{}\\ &{} &{} \ddots &{}\\ &{} &{} &{} m_k^{n-1} \end{array}\right) \left( \begin{array}{c} 0\\ u_{\alpha -\varepsilon _k}^1\\ \vdots \\ u_{\alpha -\varepsilon _k}^{n-1} \end{array}\right) \nonumber \\&\quad +\frac{1}{2}\Delta t \sum _{\varepsilon _k \le \alpha }\left( \begin{array}{cccc} m_k^1 &{} &{} &{} \\ &{}m_k^2 &{} &{}\\ &{} &{} \ddots &{}\\ &{} &{} &{} m_k^{n} \end{array}\right) \left( \begin{array}{c} u_{\alpha -\varepsilon _k}^1\\ u_{\alpha -\varepsilon _k}^2\\ \vdots \\ u_{\alpha -\varepsilon _k}^{n} \end{array}\right) \nonumber \\&=\vec {f}_{\alpha }^0+\frac{1}{2}\Delta t \sum _{\varepsilon _k \le \alpha } \left( M_k^0 \vec {U}^0+ {M_k} \vec {U}\right) \vec {c}_{\alpha -\varepsilon _k} \end{aligned}$$

(6.3)

where $M_k=\text {diag}\{m_k^1, m_k^2, \ldots ,m_k^{n}\}$.

Due to the hierarchical nature of the RB space, the RB stiffness matrix $A_{RB}=\vec {U}^TA^TA\vec {U}$ can be formed by appending a row and a column each time a new reduced basis $\vec {U}_{\alpha _{i+1}}$ is identified. That is, we exploit the following identity.

$$\begin{aligned} \left( \begin{array}{c} \vec {U}^T\\ (\vec {U}_{\alpha _{i+1}})^{T} \end{array}\right) A^T A \left( \vec {U},\vec {U}_{\alpha _{i+1}}\right) =\left( \begin{array}{cc} \vec {U}^TA^TA\vec {U}&{} \vec {U}^TA^TA\vec {U}_{\alpha _{i+1}}\\ (\vec {U}_{\alpha _{i+1}})^{T}A^TA\vec {U}&{}(\vec {U}_{\alpha _{i+1}})^{T}A^TA\vec {U}_{\alpha _{i+1}} \end{array}\right) . \end{aligned}$$

(6.4)

For the RB right hand side $\vec {f}_{RB}=\vec {U}^{T}A^T\vec {f}_\alpha$, recognizing that

$$\begin{aligned} \vec {f}_{RB}&=\vec {U}^{T}A^T\vec {f}_\alpha =\left\{ \begin{array}{cl} \vec {U}^{T}A^T\vec {f}_{\alpha }^0+\Delta t \sum _{\varepsilon _k \le \alpha } \vec {U}^{T}A^TM_k^0 \vec {U}^0 \vec {c}_{\alpha -\varepsilon _k} &{}\text {for FE}\\ \vec {U}^{T}A^T\vec {f}_{\alpha }^0+\frac{\Delta t}{2}\sum _{\varepsilon _k \le \alpha }\vec {U}^{T}A^T(M_k^0 \vec {U}^0+ M_k \vec {U})\vec {c}_{\alpha -\varepsilon _k} &{}\text {for CN } \end{array} \right. \end{aligned}$$

(6.5)

we can also exploit the hierarchical nature to gradually build up $\vec {U}^{T}A^T$, $\vec {U}^{T}A^TM_k^0 \vec {U}^0$ and $\vec {U}^{T}A^TM_k \vec {U}$.

1.2 A.2 Step 6

Efficient and accurate evaluation of the error estimator is critical for the correct identification of the key multi-indices and thus the convergence of the COFRB algorithms. The classical approach of computing the square norm of $A\vec {U} \vec {c}_{\alpha }-\vec {f}_{\alpha }$ and then expanding it to enable an offline-online decomposition leads to numerical instability [11]. In our setting, it will result in this norm being negative. To detail the numerically stable method, we follow [11]. Noting that we can assume $\vec {f}_{\alpha }^0=0$ since $u_{\alpha -\varepsilon _k}^0=0$ when $|\alpha | \ge 2$ due to the initial condition being deterministic and all $\alpha$ with $|\alpha | = 1$ are usually chosen meaning the residual will be zero, we can rewrite the residual $A\vec {U} \vec {c}_{\alpha }-\vec {f}_{\alpha }$ as

$$\begin{aligned} A\vec {U} \vec {c}_{\alpha }-\vec {f}_{\alpha }= {\mathcal {B}}{\tilde{C}}_{\alpha } \end{aligned}$$

(6.6)

where

$$\begin{aligned} {\mathcal {B}}=&\left( A\vec {U}_{\alpha _1}, M_{1}^0\vec {U}_{\alpha _1}^0, M_2^0\vec {U}_{\alpha _1}^0,\ldots ,M_K^0\vec {U}_{\alpha _1}^0,A\vec {U}_{\alpha _2},M_1^0\vec {U}_{\alpha _2}^0, \ldots ,M_K^0\vec {U}_{\alpha _2}^0,\ldots ,A\vec {U}_{\alpha _i},\right. \nonumber \\&\left. \qquad M_1^0\vec {U}_{\alpha _i}^0,\ldots ,M_K^0\vec {U}_{\alpha _i}^0\right) \end{aligned}$$

(6.7)

$$\begin{aligned} {\tilde{C}}_{\alpha }=&\left( c_1({\alpha }),-\Delta t \delta _{\alpha -\varepsilon _1}{c}_1({\alpha -\varepsilon _1}),-\Delta t \delta _{\alpha -\varepsilon _2}{c}_1({\alpha -\varepsilon _2}),\ldots ,-\Delta t \delta _{\alpha -\varepsilon _K}{c}_1({\alpha -\varepsilon _K}), \right. \nonumber \\&\left. {c}_2({\alpha }),-\Delta t \delta _{\alpha -\varepsilon _1}{c}_2({\alpha -\varepsilon _1}),-\Delta t \delta _{\alpha -\varepsilon _2}{c}_2({\alpha -\varepsilon _2}),\ldots ,-\Delta t \delta _{\alpha -\varepsilon _K}{c}_2({\alpha -\varepsilon _K}),\ldots ,\nonumber \right. \\&\left. {c}_i({\alpha }),-\Delta t \delta _{\alpha -\varepsilon _1}{c}_i({\alpha -\varepsilon _1}),-\Delta t \delta _{\alpha -\varepsilon _2}{c}_i({\alpha -\varepsilon _2})\ldots , -\Delta t \delta _{\alpha -\varepsilon _K}{c}_i({\alpha -\varepsilon _K})\right) ^T \end{aligned}$$

(6.8)

and $\delta _{\alpha -\varepsilon _k}$ is defined as follows

$$\begin{aligned} \delta _{\alpha -\varepsilon _k}=\left\{ \begin{aligned}&1, \quad \text { if } \varepsilon _k\le \alpha \text { is true }\\&0, \quad \text { otherwise. } \end{aligned}\right. \end{aligned}$$

(6.9)

We adopt the rank-revealing QR factorization through modified Gram–Schmidt for matrix ${\mathcal {B}}$.

$$\begin{aligned} {\mathcal {B}}={\mathcal {Q}}{\mathcal {R}},\quad {\mathcal {Q}}\in {\mathbb {R}}^{n \times \text {rank}({\mathcal {B}})},\quad {\mathcal {R}} \in {\mathbb {R}}^{\text {rank}({\mathcal {B}}) \times M(K+1)}. \end{aligned}$$

(6.10)

where $\text {rank}({\mathcal {B}})\le M(K+1)$ is the rank of matrix ${\mathcal {B}}$. Then

$$\begin{aligned} \Vert A\vec {U} \vec {c}_{\alpha }-\vec {f}_{\alpha }\Vert ^2={\tilde{C}}_{\alpha }^ T{\mathcal {R}}^T{\mathcal {R}}{\tilde{C}}_{\alpha }. \end{aligned}$$

(6.11)

The cost of computing this term is independent of n if we pre-compute matrix ${\mathcal {R}}$. Notice that the matrices ${\mathcal {Q}}$ and ${\mathcal {R}}$ must also be gradually expanded similar to the way RB stiffness matrix is handled recognizing that QR factorization can be expanded in a hierarchical fashion as the data matrix is expanded.

1.3 A.3 Step 10

When a new multi-index $\alpha _{i+1}$ is deemed a candidate for addition to the RB space i.e. $\Delta _i(\alpha _{i+1}) > \varepsilon _{\text {tol}}$, we take $v_0=\vec {U}_{\alpha _{i+1}}$. Then the modified Gram–Schmidt follows

$$\begin{aligned}&{\tilde{v}}_j=v_{j-1}-<v_{j-1},\vec {U}_{\alpha _i}>\vec {U}_{\alpha _i}, \end{aligned}$$

(6.12)

$$\begin{aligned}&v_j=\frac{{\tilde{v}}_j}{\Vert {\tilde{v}}_j\Vert }, \quad j=1,2,\ldots ,i \end{aligned}$$

(6.13)

where $<\cdot ,\cdot>$ denotes the inner product. If $\Vert {\tilde{v}}_j\Vert =0$ for some $j\le i$, we discard this candidate and set the (RB) coefficients for this multi-index as

$$\begin{aligned} \vec {c}_{\alpha _{i+1}}=\left(<v_0,\vec {U}_{\alpha _1}>,<v_1, \vec {U}_{\alpha _2}>\Vert {\tilde{v}}_1\Vert ,\ldots , <v_{j-1},\vec {U}_{\alpha _j}>\Pi _{l=1}^{j-1}\Vert {\tilde{v}}_l\Vert ,0,\ldots ,0\right) . \end{aligned}$$

Otherwise, we update $\vec {U}_{\alpha _{i+1}}=v_i$ and augment the basis matrix, and set the (RB) coefficients

$$\begin{aligned} \vec {c}_{\alpha _{i+1}}=\left(<v_0,\vec {U}_{\alpha _1}>,<v_1,\vec {U}_{\alpha _2}>\Vert {\tilde{v}}_1\Vert ,\ldots , <v_{i-1},\vec {U}_{\alpha _i}>\Pi _{j=1}^{i-1}\Vert {\tilde{v}}_j\Vert , \Pi _{j=1}^{i}\Vert {\tilde{v}}_j\Vert \right) . \end{aligned}$$

Appendix B: Implementation of Algorithm 3

In this appendix, we present some details for the efficient implementation of Algorithm 3. The RB stiffness matrix and vector $A_{RB}$ and $\vec {f}_{RB}=\vec {W}^T\vec {V}^TA^T\vec {f}_{\alpha }$ are expanded in the same fashion as in Appendix A. For calculating $\Delta _i(\alpha )$, slight changes are necessary. Under the same assumption $\vec {f}^0_{\alpha }=0$, we have

$$\begin{aligned} \Vert A\vec {V}\vec {W} \vec {c}_{\alpha }-\vec {f}_{\alpha }\Vert ^2&= \Vert {\mathcal {B}}{\tilde{C}}_{\alpha }\Vert ^2 \end{aligned}$$

(6.14)

$$\begin{aligned}&={\tilde{C}}_{\alpha }^T{\mathcal {B}}^T{\mathcal {B}}{\tilde{C}}_{\alpha }, \end{aligned}$$

(6.15)

where

$$\begin{aligned} {\mathcal {B}}=&\left( A\vec {V}_{\alpha _1}^m\vec {W}_{\alpha _1}, M_{1}^0\vec {V}_{\alpha _1}^m\vec {W}_{\alpha _1}^0, M_2^0\vec {V}_{\alpha _1}^m\vec {W}_{\alpha _1}^0,\ldots , M_K^0\vec {V}_{\alpha _1}^m\vec {W}_{\alpha _1}^0,A\vec {V}_{\alpha _2}^m \vec {W}_{\alpha _2},M_1^0\vec {V}_{\alpha _2}^m\vec {W}_{\alpha _2}^0, \ldots ,\right. \nonumber \\&\left. M_K^0\vec {V}_{\alpha _2}^m\vec {W}_{\alpha _2}^0, \ldots ,A\vec {V}_{\alpha _i}^m\vec {W}_{\alpha _i},M_1^0 \vec {V}_{\alpha _i}^m\vec {W}_{\alpha _i}^0,\ldots ,M_K^0 \vec {V}_{\alpha _i}^m\vec {W}_{\alpha _i}^0\right) . \end{aligned}$$

(6.16)

We first use hierarchical expansion to form the matrix ${\mathcal {B}}^T{\mathcal {B}}$ which is of size $i(K+1)$-by-$i(K+1)$. Since it is symmetrical and positive definite, there exists an orthogonal decomposition. Hence,

$$\begin{aligned} {\tilde{C}}_{\alpha }^T{\mathcal {B}}^T{\mathcal {B}}{\tilde{C}}_{\alpha }= {\tilde{C}}_{\alpha }^TP^T\Lambda P {\tilde{C}}_{\alpha }=\omega _\alpha ^T\Lambda \omega _{\alpha } \end{aligned}$$

(6.17)

where P is an $i(K+1)$-by-$i(K+1)$ orthogonal matrix, $\Lambda$ is a diagonal matrix whose diagonal elements are square of singular values of matrix ${\mathcal {B}}$ denoted by $s_j$ and $\omega _{\alpha }=P{\tilde{C}}_{\alpha }$. Thus the residue

$$\begin{aligned} \Vert A\vec {U} {\vec {c}}_{\alpha }-\vec {f}_{\alpha }\Vert ^2=\sum _{j=1}^{i(K+1)}s_j(\omega _{\alpha })_j^2. \end{aligned}$$

(6.18)

In this way, we only need to store the $i(K+1)$-by-$i(K+1)$ matrices ${\mathcal {B}}^T{\mathcal {B}}$, P and $i(K+1)$ diagonal elements of $\Lambda$.