L1-ROC and R2-ROC: L1- and R2-based Reduced Over-Collocation methods for parametrized nonlinear partial differential equations

TL;DR

This paper introduces L1-ROC and R2-ROC, innovative reduced over-collocation methods that enhance efficiency and stability in solving parametrized nonlinear PDEs, overcoming limitations of traditional empirical interpolation methods.

Contribution

The paper develops two new reduced over-collocation methods, L1-ROC and R2-ROC, which improve online efficiency and stability for nonlinear pPDEs by extending EIM with strategic collocation and parameter selection.

Findings

High efficiency and accuracy demonstrated in numerical tests

Methods outperform traditional EIM-based approaches

Stable performance across different nonlinear problems

Abstract

The onerous task of repeatedly resolving certain parametrized partial differential equations (pPDEs) in, e.g. the optimization context, makes it imperative to design vastly more efficient numerical solvers without sacrificing any accuracy. The reduced basis method (RBM) presents itself as such an option. With a mathematically rigorous error estimator, RBM seeks a surrogate solution in a carefully-built subspace of the parameter-induced high fidelity solution manifold. It can improve efficiency by several orders of magnitudes leveraging an offline-online decomposition procedure. However, this decomposition, usually through the empirical interpolation method (EIM) when the PDE is nonlinear or its parameter dependence nonaffine, is either challenging to implement, or severely degrading to the online efficiency. In this paper, we augment and extend the EIM approach in the context of solving pPDEs in two different ways, resulting in the Reduced Over-Collocation methods (ROC). These are stable and capable of avoiding the efficiency degradation inherent to a direct application of EIM. There are two ingredients of these methods. First is a strategy to collocate at about twice as many locations as the number of bases for the surrogate space. The second is an efficient approach for the strategic selection of the parameter values to build the reduced solution space for which we study two choices, a recent empirical L1 approach and a new indicator based on the reduced residual. Together, these two ingredients render the schemes, L1-ROC and R2-ROC, online efficient and immune from the efficiency degradation of EIM for nonlinear and nonaffine problems offline and online. Numerical tests on three different families of nonlinear problems demonstrate the high efficiency and accuracy of these new algorithms and their superior stability performance.

Tables

Table 1. Table 1: Notation and terminology used throughout this article.

$𝝁=({\mu }_1,\dots ,{\mu }_p)$	Parameter in $p(=2,\text{in this paper})$-dimensional parameter domain $𝒟\subseteq I{R}^p$
${\Xi }_{\mathrm{train}}$	Parameter training set, a finite subset of $𝒟$
$u\left(𝝁\right)$	Function-valued solution of a parameterized PDE on and $\Omega \subset {ℝ}^d$
$𝒫(u\left(𝝁\right);𝝁)$	A (nonlinear) PDE operator
$K$	Number of finite difference intervals per direction of the physical domain
$𝒩\approx K^d$	Degrees of freedom (DoF) of a high-fidelity PDE discretization, called “truth” solver
$u^{𝒩}\left(𝝁\right)$	Finite-dimensional truth solution
$N$	Number of reduced basis snapshots, $N\ll 𝒩$
${𝝁}^j$	“Snapshot” parameter values, $j=1,\dots ,N$
${\hat{u}}_n\left(𝝁\right)$	Reduced basis solution in the $n$-dimensional RB space spanned by $\{u^{𝒩}\left({𝝁}^1\right),\dots ,u^{𝒩}\left({𝝁}^n\right)\}$
$e_n\left(𝝁\right)$	Reduced basis solution error, equals $u^{𝒩}\left(𝝁\right)-{\hat{u}}_n\left(𝝁\right)$
${\Delta }_N\left(𝝁\right)$	A residual-based error estimate (upper bound) for $\Vert e_N\left(𝝁\right)\Vert$ or an error indicator
$X^{𝒩}$	A size-$𝒩$ (full) collocation grid
$X_r^{N-1}=\{{𝒙}_{\ast \ast }^1,\dots ,{𝒙}_{\ast \ast }^N\}$	A size-$N-1$ reduced collocation grid. It is a subset of $X^{𝒩}$ determined based on residuals
$X_s^N=\{{𝒙}_{\ast }^1,\dots ,{𝒙}_{\ast }^N\}$	An additional size-$N$ reduced collocation grid, a subset of $X^{𝒩}$ determined based on the solutions
$X^M$	A reduced collocation grid of size $M$ that is $X_r^{N-1}\cup X_s^N$
${ϵ}_{\mathrm{tol}}$	Error estimate stopping tolerance in greedy sweep
Offline component	The pre-computation phase, where we produce our surrogate solver with a greedy selection of bases from the solution space
Online component	The process of solving the reduced problem, yielding the surrogate solution

Table 2. Table 2: Online computational times at different partition numbers K 𝐾 K when V = 3.85 , D = 0.15 2 formulae-sequence 𝑉 3.85 𝐷 superscript 0.15 2 V=3.85,D=0.15^{2} , M = 2 ​ N − 1 𝑀 2 𝑁 1 M=2N-1 , and N = 30 𝑁 30 N=30 .

$K$	Residual-based ROC	L1-ROC	R2-ROC	Direct FDM
200	0.000678	0.000688	0.000703	1.439812
400	0.000770	0.000646	0.000705	6.492029
800	0.000728	0.000625	0.000800	33.722112

Table 3. Table 3: Online computational times at different partition numbers K 𝐾 K when N = 40 𝑁 40 N=40 .

$({\mu }_1,{\mu }_2)$	$K$	Residual-based ROC	L1-ROC	R2-ROC	Direct FDM
$(4.55,0.42)$	200	0.003150	0.003159	0.004781	2.310034
	400	0.003067	0.003136	0.003931	11.779558
	800	0.003258	0.003162	0.004185	53.727031
$(1,1.82)$	200	0.001125	0.001060	0.001416	0.662095
	400	0.001141	0.001205	0.001299	3.338956
	800	0.001207	0.001261	0.001732	15.173460

Table 4. Table 4: Online computational times at different partition numbers K 𝐾 K when N = 20 𝑁 20 N=20 . The first two lines are for μ 1 = 32 , μ 2 = 3 formulae-sequence subscript 𝜇 1 32 subscript 𝜇 2 3 \mu_{1}=32,\mu_{2}=3 .

$K$	Residual-based ROC	L1-based ROC	R2-based ROC	Direct FDM
200	0.000422	0.000428	0.000523	0.569732
400	0.000397	0.000410	0.000515	2.838783
800	0.000424	0.000425	0.000512	12.582593

Equations

$$\mathcal{P}(u(\mathbf{x};\mathbf{\mu });\mathbf{\mu })-f(\mathbf{x})=0,\mathbf{x}\in \Omega ,$$

$${\mathcal{P}}_{\mathcal{N}}(u^{\mathcal{N}}(X^{\mathcal{N}},\mathbf{\mu });\mathbf{\mu })-f(X^{\mathcal{N}})=0,$$

$$X^M=\{{\mathbf{x}}_{\ast }^1,\dots ,{\mathbf{x}}_{\ast }^M\},$$

$$W_{n,M}=[u_1(X^M),\dots ,u_n(X^M)]\in {\mathbb{R}}^{M\times n},\text{mbox}forn=1,\dots ,N.$$

$${\hat{u}}_n(\mathbf{\mu })=W_n{\mathbf{c}}_n(\mathbf{\mu }).$$

$${\mathcal{P}}_{\mathcal{N}}(W_{n,M}{\mathbf{c}}_n(\mathbf{\mu });\mathbf{\mu })-f(X^M)=0.$$

$${\nabla }_h{\hat{u}}_n(\mathbf{\mu })$$

$${\Delta }_h{\hat{u}}_n(\mathbf{\mu })$$

$${\mathbf{\mu }}^{n+1}=\underset{\mathbf{\mu }\in {\Xi }_{\mathrm{train}}}{\mathrm{argmax}}{\Delta }_n(\mathbf{\mu }),W_{n+1}=W_n\oplus \{u^{\mathcal{N}}({\mathbf{\mu }}^{n+1})\}.$$

$${\Delta }_n^R(\mathbf{\mu })=\frac{\Vert f-{\mathcal{P}}_{\mathcal{N}}({\hat{u}}_n;\mathbf{\mu }){\Vert }_2}{\sqrt{{\beta }_{LB}(\mathbf{\mu })}},$$

$${\Delta }_n^L(\mathbf{\mu })=||{\mathbf{c}}_n(\mathbf{\mu })|{|}_1.$$

$${\Delta }_n^{RR}(\mathbf{\mu })=\Vert f-{\mathcal{P}}_{\mathcal{N}}({\hat{u}}_n;\mathbf{\mu }){\Vert }_{L^{\infty }(X^M)}.$$

$$r_{n-1}={\mathcal{P}}_{\mathcal{N}}(X^{\mathcal{N}},{\hat{u}}_{n-1}({\mathbf{\mu }}^n);{\mathbf{\mu }}^n)-f(X^{\mathcal{N}}),n\in \{2,\dots ,N\}$$

$$M=2N-1\text{mbox}and{X}^M=X_s^N\cup X_r^{N-1},$$

$$u^{\mathcal{N}}(\mathbf{\mu })=\sum _{j=1}^J{\alpha }_j(\mathbf{\mu }){\zeta }_j.$$

$$({\zeta }_i,{\zeta }_j{)}_{{\mathcal{P}}_{\mathcal{N}}}={\delta }_{ij},$$

$$c_1(\mathbf{\mu })=\underset{c\in \mathbb{R}}{\mathrm{argmin}}\sum _{j=1}^J{\left({\alpha }_j(\mathbf{\mu })-c{\alpha }_j({\mathbf{\mu }}^1)\right)}^2=\frac{{\sum }_{j=1}^J{\alpha }_j(\mathbf{\mu }){\alpha }_j({\mathbf{\mu }}^1)}{{\sum }_{j=1}^J{\left({\alpha }_j({\mathbf{\mu }}^1)\right)}^2}.$$

$${\mathbf{\mu }}^2=\underset{\mathbf{\mu }={\Xi }_{\mathrm{train}}}{\mathrm{argmax}}|c_1(\mathbf{\mu })|=\underset{\mathbf{\mu }={\Xi }_{\mathrm{train}}}{\mathrm{argmax}}\left|\frac{{\sum }_{j=1}^J{\alpha }_j(\mathbf{\mu }){\alpha }_j({\mathbf{\mu }}^1)}{{\sum }_{j=1}^J{\left({\alpha }_j({\mathbf{\mu }}^1)\right)}^2}\right|=\underset{\mathbf{\mu }={\Xi }_{\mathrm{train}}}{\mathrm{argmax}}\left|\sum _{j=1}^J{\alpha }_j(\mathbf{\mu }){\alpha }_j({\mathbf{\mu }}^1)\right|.$$

$${\mathbf{\mu }}^2=\underset{({\mu }_1,{\mu }_2)\in \mathcal{D}}{\mathrm{argmax}}\left|{\mu }_1^1{\mu }_1+{\mu }_2^1{\mu }_2\right|$$

$$D{\nabla }^{2}u=\sinh u+g(\mathbf{x})$$

$$u(x_1=-1,x_2)=0,$$

$$u(x_1=+1,x_2)=V,$$

$${\partial }_{x_2}u(x_1,x_2=\pm 1)=0.$$

$${\Xi }_{\mathrm{train}}=(0.08:0.02:0.4)\times (0:0.25:5),$$

$${\Xi }_{\mathrm{test}}=(0.085:0.01:0.395)\times (0.4:0.5:4.4),$$

$$E(n)=\underset{\mathbf{\mu }\in {\Xi }_{\mathrm{test}}}{\max }\{\Vert u(\mathbf{\mu })-{\hat{u}}_n(\mathbf{\mu }){\Vert }_{\infty }/\Vert u{\Vert }_{L^{\infty }({\Xi }_{\mathrm{test}},L^{\infty }(\Omega ))}\}$$

$$||u|{|}_{L^{\infty }({\Xi }_{\mathrm{test}},L^{\infty }(\Omega ))}=\underset{\mathbf{\mu }\in {\Xi }_{\mathrm{test}}}{\max }\Vert u(\mathbf{\mu }){\Vert }_{\infty }.$$

$${\Xi }_{test2}=(0.09:0.01:0.31)\times (0.7:0.2:1.5),$$

$$-{\mu }_2\Delta u+u{(u-{\mu }_1)}^2=f(\mathbf{x})$$

$${\Xi }_{\mathrm{train}}=(0.2:4h_1:5)\times (0.2:4h_2:2),$$

$${\Xi }_{\mathrm{test}}=((0.2+2h_1):4h_1:(5-2h_1))\times ((0.2+2h_2):4h_2:(2-2h_2)),$$

$${\Xi }_{test2}=((0.2+3h_1):4h_1:(5-3h_1))\times ((0.2+3h_2):2h_2:(2-3h_2)).$$

Full text

In this paper, we augment and extend the EIM approach in the context of solving pPDEs in two different ways, resulting in the Reduced Over-Collocation methods (ROC). These are stable and capable of avoiding the efficiency degradation inherent to a direct application of EIM. There are two ingredients of these methods. First is a strategy to collocate at about twice as many locations as the number of bases for the surrogate space. Half of these points come from bases while the other half are from the residuals when these bases are adopted to solve the pPDE. The second is an efficient approach for the strategic selection of the parameter values to build the reduced solution space for which we propose two choices. In addition to the recently introduced empirical L1 approach which is further analyzed and tested, we propose and test a new indicator that is based on the $L^{\infty}$ norm of the reduced residual (R2), the residual sampled at these reduced collocation points. Together, these two ingredients render the schemes, L1-ROC and R2-ROC, online efficient (i.e. online cost is independent of the number of degrees of freedom of the high-fidelity truth approximation) and immune from the efficiency degradation of EIM for nonlinear and nonaffine problems. Moreover, they are highly efficient offline in that they require very little computation in addition to the bare minimum of acquiring the bases for the surrogate space. Numerical tests on three different families of nonlinear problems demonstrate the high efficiency and accuracy of these new algorithms and its superior stability performance.

1 Introduction

Design of fast numerical algorithms with certifiable accuracies for parametrized systems arising from various engineering and applied science disciplines has continued to attract researchers’ attention. The parameters delineating these systems may include boundary conditions, material properties, geometric settings, source properties etc. The wide variety, the complicated dependence of the system on these parameters, and their potential high dimensionality consist of some of the major challenges. The reduced basis method (RBM) has proved an effective option for this purpose [34, 24].

RBM was first introduced for nonlinear structure problem [1, 32] in 1970s and has proven to be effective for linear evolutionary equations [23], viscous Burgers equation [36], Navier-Stokes equations [16], and harmonic Maxwell’s equation [13, 14], just to name a few. The key to RBM’s success in realizing the efficiency gain per parameter instance, is an offline-online decomposition process where the basis selection is performed offline by a greedy algorithm, see review papers [35, 22] and monographs [34, 24] for details. During the offline process, the necessary preparations for the online reduced solver are performed. The ultimate goal is that the complexity of the reduced solver, called upon in a potentially real-time fashion online, is independent of the degrees of freedom of the high-fidelity approximation of the basis functions, solutions to the system at certain judiciously selected configurations. For the nonaffine and nonlinear equations, the Empirical Interpolation Method (EIM) [3, 21, 10, 33] is usually used to achieve the online independence of the degrees of freedom. In practice, EIM is often not feasible due to severe nonlinearity and/or nonaffinity of the problem. Unfortunately performing EIM, even when feasible, severely degrades this online efficiency when either the parameter dependence or the nonlinearity is complicated such as when it involves geometric parametrization [14, 4]. The reason is that the online complexity is dependent on the number of terms resulting from the EIM decomposition.

In this paper, we design two reduced over-collocation (ROC) methods achieving full online-efficiency. They are stable and much more efficient than the typical RBM adopting directly EIM, thanks to an augmentation of EIM and further leveraging of the collocation philosophy originally explored in [11]. There are two ingredients of the ROC methods.

First is a strategy to fully explore the EIM framework and partially circumvent its efficiency degradation by adopting the collocation approach as opposed to a variational (i.e. Galerkin or Petrov-Galerkin) approach [5, 7, 6] when seeking the reduced solution. This so-called reduced collocation method is proposed and documented to work well in circumventing the EIM degradation for the reduced solver in our previous work [11]. However, its stability is lacking [12]. Our reduced over-collocation methods mitigate this stability defect by collocating at about twice as many locations as the number of bases for the surrogate space. Half of these points come from the bases. They interpolate the reduced solution (a linear combination of these bases) well. The other half are from the judiciously selected residuals when these bases are determined during the offline procedure. They are present to ensure a good interpolation of the residual corresponding to an arbitrary parameter value when the bases are adopted to solve the pPDE.

This ingredient alone is not enough to achieve online and offline efficiency as the efficient calculation of the error estimators, critical for the construction of the reduced solution space, still relies on direct application of EIM. This is now resolved by the second ingredient of our ROC methods, an efficient alternative for guiding the strategic selection of parameter values to build the reduced solution space. We examine two choices toward that end. In addition to the recently introduced empirical L1 approach [15] which is further analyzed and tested (producing L1-ROC), we propose and test a new indicator that is based on the $L^{\infty}$ norm of the reduced residual (R2, producing R2-ROC), the residual sampled at these reduced collocation points.

Together, these two ingredients render the schemes online efficient (i.e. online cost is independent of the number of degrees of freedom of the high-fidelity truth approximation) successfully avoiding the efficiency degradation of a direct EIM for nonlinear and nonaffine problems. Moreover, the ROC methods are highly efficient offline in that they require minimal computation in addition to that for acquiring the basis snapshots of the surrogate space. As a consequence, the “break-even” number of simulations for the pPDE (minimum number of simulations that make the offline preparation stage worthwhile) is significantly smaller than traditional RBM and, in fact, comparable to the dimension of the surrogate space, the minimum possible break-even number. We test the algorithms on the Poisson-Boltzmann equation (PBE) [20, 9, 2, 17] and two additional nonaffine and nonlinear PDEs with severely nonlinear reaction or convection terms. We note that PBE is a boundary layer problem and plays important roles in understanding the electrostatic phenomenon in physical, biological and materials sciences [18, 37, 29, 31] at the nano/micro scale. We show the new L1-ROC and R2-ROC methods improve upon the performance of RBM in our paper [27] which achieves partially order reduction of fully nonlinear PBE. For all test problems, both ROC methods are shown to have accuracy on par with the classical RBM while possessing much better efficiency due to the independence of the number of expansion terms resulting from the EIM decomposition.

The paper is organized as follows. In Section 2, we introduce the L1-ROC and R2-ROC methods. Section 3 is devoted to theoretical and numerical understandings of the reliability of the L1 approach. Numerical results for three test problems to demonstrate the accuracy and efficiency of our L1-ROC and R2-ROC methods are shown in Section 4. Finally, concluding remarks are drawn in Section 5.

2 The Reduced over-collocation (ROC) methods

In this section, we introduce the L1-ROC and R2-ROC methods. Toward that end, we first describe the problem we are solving. The framework of the online algorithm is then presented in Section 2.1. Specification of part of the algorithm is postponed until the introduction of two versions of the reduced over collocation offline algorithm in Section 2.2 which repeatedly calls the online solver to construct a surrogate solution space. The design of the main algorithm, the reduced over-collocation (ROC) approach, is detailed in Section 2.2.2. To facilitate the reading of this and the following sections, we list our notations in Table 1.

We let $\mathcal{D}\subset\mathbb{R}^{p}$ be the domain for a $p$-dimensional parameter $\boldsymbol{\mu}$, and $\Omega\subset\mathbb{R}^{d}$ (for $d=2~{}\text{or}~{}3$) be a bounded physical domain. Given $\bm{\mu}\in\mathcal{D}$, the goal is to compute $u(\bm{\mu}):=u(\boldsymbol{x};\boldsymbol{\mu})\in H^{1}(\Omega)$ satisfying

[TABLE]

with $\mathcal{P}$ encoding a parametric partial differential operator that may include linear and nonlinear functions of $u(\boldsymbol{x};\boldsymbol{\mu})$, $\nabla u(\boldsymbol{x};\boldsymbol{\mu})$, and $\Delta u(\boldsymbol{x};\boldsymbol{\mu})$. We further discretize this equation by a high-fidelity scheme, termed “truth solver”. In this paper, we adopt Finite Difference Methods (FDM) as the truth solver. However, extension to point-wise schemes such as spectral collocation is obvious, and to Finite Element Method is possible. Indeed, we let $X^{\cal N}$ be a set of (roughly ${\cal N}$) collocation points on $\Omega$ at which the equation is enforced on a discrete level. The discretized equation then becomes to find $u^{\cal N}(X^{\cal N},\boldsymbol{\mu})$, a discretization of the $H^{1}$ function $u(\boldsymbol{\mu})$, such that we have

[TABLE]

with $\nabla u(X^{\cal N};\boldsymbol{\mu})$, and $\Delta u(X^{\cal N};\boldsymbol{\mu})$ approximated by their numerical counterparts $\nabla_{h}u(X^{\cal N};\boldsymbol{\mu})$, and $\Delta_{h}u(X^{\cal N};\boldsymbol{\mu})$. With a slight abuse of notation, we are adopting ${\cal N}$ for the degrees of freedom as well, even though the ${\cal N}$ points in $X^{\cal N}$ might include, e.g. the points on the Dirichlet boundary. In fact, for simplicity we will generally drop the superscript ${\cal N}$ for the solution $u^{\cal N}$ in the remainder of the paper as we will not make any reference to the exact solution of the PDE.

2.1 Online algorithm

The online component of the L1-ROC is essentially the same as the previously-introduced reduced collocation method [11] with the critical difference being that the number of collocation points is larger than the number of reduced bases. This over-collocation feature gives the method its name and provides additional stabilization of the online solver as we will observe in the numerical results. To describe the online algorithm, given $N$ selected parameters $\{\boldsymbol{\mu}^{1},\dots,\boldsymbol{\mu}^{N}\}$, the corresponding high fidelity truth approximations $\{u_{n}\equiv u^{\cal N}(\boldsymbol{\mu}^{n}),1\leq n\leq N\}$, and $M$ collocation points

[TABLE]

we are able to perform the online algorithm. Note that, whenever there is no confusion, we are adopting the same notation for a function and its discrete representation in the form of a vector of its values at the grid points. These vectors $\{u_{n},1\leq n\leq N\}$ constitute the basis spaces/matrices $W_{n}\in\mathbb{R}^{{\cal N}\times n}$ for $n\in\{1,\dots,N\}$. Furthermore, we denote the corresponding reduced representation of the basis space on the set $X^{M}$, by a matrix of the following form,

[TABLE]

Reduced approximations of the solution for any given parameter $\boldsymbol{\mu}$ is sought in the form of

[TABLE]

Substituting this into equation (2), we will obtain a system of equations for the unknown coefficients $\boldsymbol{c}_{n}(\boldsymbol{\mu})$ at the reduced collocation nodes $X^{M}$,

[TABLE]

We note that this is a nonlinear system of equations for $\boldsymbol{c}_{n}$ with $\nabla_{h}\widehat{u}_{n}(\boldsymbol{\mu})$ and $\Delta_{h}\widehat{u}_{n}(\boldsymbol{\mu})$ computed on the full grid and then evaluated on the reduced grid $X^{M}$ according to

[TABLE]

whose right hand sides are simply denoted as $\nabla_{h}(W_{n,M})\boldsymbol{c}_{n}$ and $\Delta_{h}(W_{n,M})\boldsymbol{c}_{n}$ respectively. Iterative methods, such as Newton’s method, will be used to solve for the coefficients $\boldsymbol{c}_{n}(\boldsymbol{\mu})$.

Once the offline preparation is under its way and the snapshot locations $\boldsymbol{\mu}^{j}$ are gradually determined, we precompute as many quantities as possible so that minimal update is performed at each iteration of the iterative method. The online procedure of the nonlinear solve for obtaining $\boldsymbol{c}_{n}(\boldsymbol{\mu})$ from equation (3) is independent of the degrees of freedom ${\cal N}$ of the underlying truth solver, and involves:

realizing/updating $W_{n,M}\boldsymbol{c}_{n}$, $\nabla_{h}(W_{n,M})\boldsymbol{c}_{n}$, and $\Delta_{h}(W_{n,M})\boldsymbol{c}_{n}$ at each iteration taking $O(Mn)$ operations;

• 2)

calculating the forcing term $f(X^{M})$ taking $O(M)$ operations; and

• 3)

solving the reduced linear systems at each iteration of the nonlinear solve taking $O(N^{3})$ operations.

2.2 Offline algorithm

In this section, we describe the offline procedure of the reduced over collocation framework resulting in two different approaches depending on how the reduced collocation set $X^{M}$ is determined. We are going to use the L1-approach proposed in [15] and reviewed briefly below or a newly proposed R2-approach for the critical greedy algorithm executed offline to construct $W_{N}$. The remaining ingredients of the offline procedure is identical with the traditional RBM algorithm [35, 22, 34, 24].

2.2.1 L1- and R2-based greedy algorithm

We first briefly describe the procedure for selecting the representative parameters $\boldsymbol{\mu}^{1},\ldots,\boldsymbol{\mu}^{N}$ for constructing the solution space $W_{N}$. RBM utilizes a greedy scheme to iteratively construct $W_{N}$ relying on an efficiently-computable error estimates that quantify the discrepancy between the dimension-$n$ RBM surrogate solution $\widehat{u}_{n}(\boldsymbol{\mu})$ and the truth solution $u^{\mathcal{N}}(\boldsymbol{\mu})$. Denoted $\Delta_{n}$, this error estimate traditionally satisfies $\Delta_{n}(\boldsymbol{\mu})\geq\left\|\widehat{u}_{n}(\boldsymbol{\mu})-u^{\mathcal{N}}(\boldsymbol{\mu})\right\|$. Assuming existence of this error estimate, the greedy procedure for constructing $W_{N}$ then starts by selecting the first parameter $\boldsymbol{\mu}^{1}$ randomly from $\Xi_{\rm train}$ (a discretization of the parameter domain $\mathcal{D}$) and obtaining its corresponding high-fidelity truth approximation $u^{\mathcal{N}}(\boldsymbol{\mu}^{1})$ to form a (one-dimensional) RB space $W_{1}=\{u^{\mathcal{N}}(\boldsymbol{\mu}^{1})\}$. Next, we obtain an RB approximation $\widehat{u}_{n}(\boldsymbol{\mu})$ for each parameter in $\Xi_{\rm train}$ together with an error bound $\Delta_{n}(\boldsymbol{\mu})$. The greedy choice for the $(n+1)$th parameter $(n=1,\cdots,N-1)$ is made and the RB space augmented by

[TABLE]

The design and efficient implementation of the error bound $\Delta_{n}$ is usually accomplished with a residual-based a posteriori error estimate from the truth discretization. Mathematical rigor and implementational efficiency of this estimate are crucial for the accuracy of the reduced basis solution and its efficiency gain over the truth approximation. When $\mathcal{P}(u;\boldsymbol{\mu})$ is a linear operator, the Riesz representation theorem and a variational inequality imply that $\Delta_{n}$ can be taken as

[TABLE]

which is a rigorous bound (with the R-superscript denoting it is based on the full residual). Here $\beta_{LB}(\boldsymbol{\mu})$ is a lower bound for the smallest eigenvalue of ${P}_{\mathcal{N}}(\boldsymbol{\mu})^{T}{P}_{\mathcal{N}}(\boldsymbol{\mu})$ with ${P}_{\mathcal{N}}(\boldsymbol{\mu})$ being the matrix corresponding to the discretized linear operator $\mathcal{P}_{\mathcal{N}}(\cdot;\boldsymbol{\mu})$.

Deriving the counterpart of this estimation for the general nonlinear equation is far from trivial. Moreover, even for the linear equations, the robust evaluation of the residual norm in the numerator is delicate [8, 15]. We would also have to resort to an offline-online decomposition to retain efficiency which usually means application of EIM for nonlinear or nonaffine terms. This complication degrades, sometimes significantly [4], the online efficiency due to the large number of resulting EIM terms. What exacerbates the situation further is that the (parameter-dependent) stability factor $\beta_{LB}(\boldsymbol{\mu})$ must be calculated by a computationally efficient procedure such as the Successive Constraint Method [26, 25]. For these reasons, we are going to adopt the following two empirical alternatives.

• •

[L1-based greedy] A much simpler importance indicator proposed in [15] in place of $\Delta_{n}^{R}$:

[TABLE]

The L-superscript denotes that it is based on the L1-norm making our scheme L1-based reduced over collocation method. We note that this is not an error indicator because $\Delta_{n}^{L}$ does not decrease as we increase $n$ since $\Delta_{n}^{L}(\boldsymbol{\mu}^{i})=1$ for $i\in\{1,\dots,n\}$. Nevertheless, we demonstrate that it is a reliable quantity to monitor when deciding which representative parameters $\boldsymbol{\mu}^{1},\dots,\boldsymbol{\mu}^{N}$ will form the surrogate space.

• •

[R2-based greedy] An equally simple error indicator in place of $\Delta_{n}^{R}$:

[TABLE]

Note that we are not evaluating the residual over the entire collocation grid of the truth approximation, just the reduced set $X^{M}$. It is therefore based on the reduced residual, thus the method called R2-based reduced over collocation method. In addition to being a reliable quantity to monitor when deciding the representative parameters $\boldsymbol{\mu}^{1},\dots,\boldsymbol{\mu}^{N}$ as demonstrated by our numerical results, a further advantage is that $\Delta_{n}^{RR}$ does decrease as $n$ increases. In fact, the numerical results seem to indicate the effectivity index is rather constant and small, an aspect of the algorithm that we are further investigating.

We finish this subsection by pointing out that the calculation of $\Delta_{n}^{L}$ and $\Delta_{n}^{RR}$ is independent of $\mathcal{N}$ while that for the traditional $\Delta_{n}^{R}$ is. This difference leads to the dramatic efficiency gain of the L1-ROC and R2-ROC, as numerically confirmed in Section 4.

2.2.2 Construction of the reduced over-collocation set $X^{M}$

Let us now describe how we determine the reduced collocation set $X^{M}$ to complete the offline algorithm. Toward that end, we first describe the construction of two sets. The first one, denoted by $X^{N}_{s}$, consists of the maximizers of the EIM-orthonormalized basis functions of $W_{N}$. The second one examines the residual of the RB solution at $\boldsymbol{\mu}^{n}$ when only $n-1$ bases are used,

[TABLE]

It takes these $N-1$ residual vectors and performs an EIM orthonormalization. The $N-1$ maximizers from this orthonormalization form the second set $X^{N-1}_{r}$. The reduced collocation approach in [11] amounts to simply taking $M=N$ and $X^{M}=X^{N}_{s}$. The resulting reduced scheme can be unstable particularly when high accuracy (i.e. large $N$) of the reduced solution is wanted. It can be resolved in special cases by the analytical preconditioning approach [12]. The second obvious choice of $X^{M}$ is to append $X^{N-1}_{r}$ with one more point such as the maximizer of the first basis. Numerical tests (not reported in this paper) also reveals instability of the scheme.

The stabilization mechanism and name of the reduced over-collocation methods, outlined in Algorithm 1, come from that we take

[TABLE]

and solve a least squares problem on the reduced level by collocating on about twice as many points as the number of bases in the RB space. The first basis has no accompanying residual vector (5). From the second onward, there are two collocation points selected whenever a new parameter is identified by the greedy algorithm.

We emphasize that Algorithm 1 leads to L1-ROC if we take $\Delta_{n}\equiv\Delta_{n}^{L}$ and it leads to R2-ROC if we take $\Delta_{n}\equiv\Delta_{n}^{RR}$.

3 Analysis of the L1-ROC method

In this section, we first provide some theoretical intuition of the reliability of the L1-based importance indicator $\Delta_{n}^{L}$ that was originally proposed in [15] and now serves as one of the two major components of our L1-ROC approach. Toward that end, we assume that equation (2) is linear and its truth approximation can be expressed in the following form

[TABLE]

Moreover, we assume that $\{\zeta_{j}:j=1,\dots,J\}$ are othornormal with respect to the inner product induced by the (linear operator) $\mathcal{P}_{\cal N}$ (i.e. the “energy” inner product):

[TABLE]

where $\delta_{ij}$ is the Kronecker delta function. Under this setting, the first reduced basis can be expressed as $u^{\cal N}(\boldsymbol{\mu}^{1})=\sum_{j=1}^{J}\alpha_{j}(\boldsymbol{\mu}^{1})\zeta_{j}$ for a given $\boldsymbol{\mu}^{1}$. The first round of greedy sweep (i.e. solving the reduced problems) under the Galerkin projection setting is to find, for each $\boldsymbol{\mu}\in\Xi_{\rm train}$,

[TABLE]

This means that the greedy choice guided by the L1-approach is

[TABLE]

The last equality stands due to the previous denominator being a constant. Therefore, we are in fact solving two (constrained) linear programing problems approximately to identify $\boldsymbol{\mu}^{2}$ if we release the condition $\boldsymbol{\mu}\in\Xi$ to $\boldsymbol{\mu}\in\mathcal{D}$. The solution is a vertex in the $(\alpha_{1},\dots,\alpha_{J})$-hypercube. The greedy algorithm may stall if, for example,

• •

$\boldsymbol{\mu}^{1}$ was also this vertex,

• •

there is a $j_{*}$ such that $\alpha_{j_{*}}(\boldsymbol{\mu}^{1})=0$ and the linear program is solved in such a way that $\alpha_{j_{*}}(\boldsymbol{\mu})=0$.

The set of $\boldsymbol{\mu}^{1}$ is a [math]-measure set in the parameter domain $\mathcal{D}$ and the probability of $\boldsymbol{\mu}^{1}$ being from this set is [math] if it was chosen randomly. We therefore conclude that the random start of the greedy offline algorithm is key to the success of the L1-approach.

Example for the importance of the random start. Taking $J=2$ (thus $\boldsymbol{\mu}=(\mu_{1},\mu_{2})$) and $\alpha_{j}(\boldsymbol{\mu})=\mu_{j}$, we are then solving

[TABLE]

We see that the procedure may stall (producing $u^{\cal N}(\boldsymbol{\mu}^{2})$ that is linearly dependent on $u^{\cal N}(\boldsymbol{\mu}^{1})$) if, for example, $\mu_{1}^{1}=0$ or $\mu_{2}^{1}=0$ or $(\mu_{1}^{1},\mu_{2}^{1})$ happens to be the vertex of this linear program solution. However, these happen with zero probability if $\boldsymbol{\mu}=(\mu_{1}^{1},\mu_{2}^{1})$ is chosen at random.

Numerical comparison with POD and random generation. There are two extremal means of building the reduced basis space. On one end, the proper orthogonal decomposition (POD) [19, 28, 38, 30] based on an exhaustive selection of snapshots produces the best reduced solution space and thus the most accurate, albeit costly, surrogate solution. On the other end, random selection of $N$ parameters as our RB snapshots is a fast but crude method.

To establish numerically the reliability of our L1-approach, we measure its convergence against these two extremal algorithms. To guarantee that the POD solution is the most accurate possible, we include all solutions $u^{\cal N}(\boldsymbol{\mu})$ for $\boldsymbol{\mu}\in\Xi_{\rm train}$. On the other end, we perform $20$ random selections of $N$ parameters for the random generation approach and record the best, median, and worst performance for each $n\in\{1,\dots,N\}$. Comparison results of three test problems (given in Section 4) are showed in Figure 1 with FDM points per dimension $K$ set to be $400$ (results with different $K$ are similar). Not surprisingly, POD is the most accurate. We note that this version of POD only serves as reference and is in general not feasible as the full solution ensemble must be generated. Our L1-ROC is one order of magnitude worse than POD, but in fact slightly better or comparable to the “Min” curve, the best possible random generation. It is roughly one order of magnitude better than the median performance of random generations.

4 Numerical results

In this section, we present the numerical results of the L1-ROC and R2-ROC methods applied to three problems, namely the fully nonlinear Poisson-Boltzmann equation (PBE), a cubic nonlinear reaction diffusion equation, and a nonlinear convection diffusion equation emulating the fluid nonlinearity.

4.1 Poisson-Boltzmann equation

We first test the following dimensionless nonlinear nonaffine Poisson-Boltzmann equation. The authors have previously designed a RBM for this equation [27]. However, due to the desire to avoid applying EIM directly, we observed limited speedup (less than one order of magnitude). Here, we show later a speedup factor of up to four orders of magnitude. Therefore, this constitutes a significant progress and underscores the power of the L1-ROC and R2-ROC approaches. Indeed, the PBE is

[TABLE]

Therefore, we are dealing with a two-dimensional parameter $\boldsymbol{\mu}:=[\sqrt{D},V]=[0.08,0.4]\times[0,5]$. The discretized training set is taken to be

[TABLE]

and a testing set

[TABLE]

which in particular does not intersect with the training set. Here, the notation $a:h:b$ denotes an equidistant discretization of the interval $[a,b]$ by elements of size $h$.

We compute the relative errors $E(n)$ over all $\boldsymbol{\mu}$ in $\Xi_{\rm test}$ of the reduced basis solution using $n$ basis functions, $\widehat{u}_{n}(\boldsymbol{\mu})$, in comparison to the high fidelity truth approximation. That is,

[TABLE]

where

[TABLE]

The error functions $E(n)$ and indicators/estimator are displayed in Figure 2 top. Clearly, both L1-ROC and R2-ROC perform similarly to the classical residual-based ROC all having stable exponential convergence.

The set of selected parameters are shown in Figure 2 middle, while the set of collocation points $X^{2N-1}$ is displayed the bottom row. We note that the distributions of chosen parameters between the traditional residual-based scheme and the nascent L1-based and R2-based schemes are very much similar which underscores the reliability of the new ROC approaches. In addition, the fact that more parameters are chosen for smaller $D$ and larger $V$ is a manifestation of the boundary layer property of the nonlinear PB equation. It is also interesting to note that more collocation points are located close or on the right boundary $x_{1}=1$. The underlying physics is that the positive voltage $V$ is applied at $x_{1}=+1$ while $u\equiv 0$ at $x_{1}=-1$.

Lastly, we showcase the vast saving of the offline time for the ROC approaches. The comparison in cumulative computation time for the residual-based, L1-ROC, R2-ROC, and the high fidelity truth approximations is shown in Figure 3. The initial nonzero start of the ROC methods is the amount of their offline time. Indeed, we record the offline time during which $N=30$ RB basis are chosen and then introduce a new test set

[TABLE]

to evaluate the cumulative time. We observe that, when $n_{\rm run}>41$, L1-ROC and R2-ROC start to save time. In comparison, the residual-based ROC is effective when $n_{\rm run}>146$ with $K=200$. The reason is that the overhead cost, devoted to calculating $\Delta_{n}^{L}$ (for L1-ROC) or $\Delta_{n}^{RR}$ (for R2-ROC), is significantly less that for $\Delta_{n}^{R}$. The latter involves (an offline-online decomposition of) the calculation of the full residual norm while the former only requires, in the L1-ROC case, obtaining an $N\times 1$ vector and evaluating its L1-norm or, in the R2-ROC case, evaluating the residual at the judiciously selected set $X^{2N-1}$. It is worth noting that the “break-even” number of runs is insensitive to $K$. Though L1-ROC and R2-ROC have a much more efficient offline procedure than the residual-based ROC, their online time for any new parameter is comparable, see Table 2. It is seen that all ROC methods accelerate the iterative truth solver by $2000\sim 50000$ times. The results also confirm that time consumption of these online ROC methods is independent of the partition numbers $K$.

4.2 Cubic reaction diffusion equation

Here, we test

[TABLE]

equipped with homogenious Dirichlet boundary. Here $f(\boldsymbol{x})=100\sin(2\pi x_{1})\cos(2\pi x_{2}),\boldsymbol{x}\in\Omega=[-,1]\times[-1,1]$,${\cal D}:=[\mu_{1},\mu_{2}]\in[0.2,5]\times[0.2,2]$ . The parameter space is discretized by a $128\times 64$ uniform tensorial grid. Denoting the step size along the $\mu_{1}$ direction by $h_{1}$, and the other by $h_{2}$, we specify the training set and test set as follows,

[TABLE]

[TABLE]

which in particular does not intersect with each other. Parameters used to compare computational time are selected from

[TABLE]

Newton’s method is used to obtain the high fidelity truth approximation. We show the $\ell+1^{\rm th}$ iteration here

[TABLE]

where $g(u;\mu_{1})=u{(u-\mu_{1})}^{2}$.

A sample high fidelity approximation $u^{\cal N}$ for $K=400$ at the center of the parameter domain is shown in Figure. 4(a). A sanity check of the iterative scheme is shown in Figure 4(b). Indeed, we set $K=200,400,800$ and $1,600$ with $\mu_{1}=2.6,\mu_{2}=1.1$ and take the solutions with $K=1,600$ as the reference. Defining

[TABLE]

with the infinity norm taken along the $x_{2}$-direction. The distribution of $E_{x}$ is shown in Fig. 4(b). This accuracy test indicates that the underlying numerical algorithm is convergent with a second-order rate. It also shows that we can terminate the offline process of the ROC methods when the relative error $E(n)$ is at the level of $\sim 10^{-4}$.

Relative errors of the RB solution at different partition numbers are displayed in Figure 5 top. Steady exponential convergence is again observed both for the L1-based and R2-based ROC matching the classical approach. The set of selected parameters are shown in Figure 5 middle, while the collocation points are shown on the bottom row. We note again that the distributions of chosen parameters between the traditional residual-based scheme and the more nascent L1/R2-based schemes are almost identical for this example underscoring the reliability of the L1-ROC and R2-ROC approaches.

Cumulative time consumption is also tested for $K=200,400,800$. With $K=200$, L1-based and R2-based ROC schemes break even when $n_{\rm run}=139$ while residual-based ROC is effective only when $n_{\rm run}>389$. It is interesting to note that the break-even point for this problem is much higher than the last one in comparison with the number of RB basis. It turns out that the reason is that the computation time for different parameter values varies dramatically for this example. Indeed, Figure 5 middle shows that many parameters with large $\mu_{1}$ and small $\mu_{2}$ are chosen. Unfortunately, the corresponding equation for these parameters need more computational time due to its nonlinear solver taking more iterations.

In order to demonstrate the time saving more intuitively, we present the calculation time at different types of given parameters by the three methods. The first kind is when $\mu_{1}$ is big and $\mu_{2}$ small, e.g. $\mu_{1}=4.55,\mu_{2}=0.42$. The second kind has the relative sizes reversed. The iteration takes about $27$ times for the first kind, while another example only takes $8$ iterations. Therefore, time consumption seems very different. However, Table 3 does indicate a speedup range of $3000\sim 17000$ when $K=400,800$.

4.3 Nonlinear convection diffusion equation

Here, we test the following nonlinear convection diffusion equation mimicking the advection terms from fluid problems

[TABLE]

equipped with zero Dirichlet boundary conditions. Here $f(\boldsymbol{x})=100\sin(2\pi x_{1})\cos(2\pi x_{2}),\boldsymbol{x}\in\Omega=[-1,1]\times[-1,1]$,${\cal D}:=[\mu_{1},\mu_{2}]\in[1,33]\times[1,5]$. We will see that the L1-ROC and R2-ROC methods handle the highly nonlinear convection term with a norm of the gradient equally well. In particular, the online cost being independent of the degrees of freedom of the truth approximation is still maintained without any direct EIM procedure.

The parameter domain is discretized by a $256\times 32$ grid. We denote the step size of $\mu_{1}$ as $h_{1}$, and the other direction by $h_{2}$ and specify the training set and test set as the following,

[TABLE]

[TABLE]

Parameters used to compare computational time are

[TABLE]

The iterative solver proceeds at the $\ell+1^{\rm th}$ iteration as

[TABLE]

with $g(u^{\ell},\mu_{1})={(||\nabla u^{\ell}||+\mu_{1})}^{1.5}$ and central finite difference is used to deal with $\nabla u^{\ell}$. A sanity check (not reported here) is performed again showing that this numerical scheme is convergent with a second-order accuracy.

Relative errors at different partition numbers, the selected parameters, and the collocation points $X^{2N-1}$ are displayed in Figure 6. Time consumption comparison shows that when $n_{\rm run}>32$, L1-ROC and R2-ROC start to save time and residual-based ROC is effective only when $n_{\rm run}>153$ with $K=200$. And the intersection points are nearly the same at different $K$. Online time of some specific parameters are displayed in Table 4. The speedup range is $10^{3}\sim 3\times 10^{4}$.

5 Conclusion

This paper proposes two novel reduced over-collocation method, dubbed L1-ROC and R2-ROC, for efficiently solving parametrized nonlinear and nonaffine PDEs. Their online computational complexity is independent of the degrees of freedom of the underlying FDM, and furthermore immune from the number of EIM expansion terms otherwise necessary to deal with the nonaffine and nonliner terms in the equation. The lack of such precomputations of nonlinear and nonaffine terms makes the method dramatically faster offline and online, and significantly simpler to implement and present than any existing RBM. An astonishing feature of the method is that the resulting break-even number of solves is comparable to the number of dimensions of the RB space. For future directions, we plan to apply these new L1-ROC and R2-ROC methods to transport equations, such as, the time dependent PNP equations [31] which plays important roles in electro chemistry and biological arenas. CFD problems involving more complicated nonlinear and nonaffine equations are also interesting and challenging directions.