# Accurate error estimation for model reduction of nonlinear dynamical systems
via data-enhanced error closure

Sridhar Chellappa<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Lihong Feng  Peter Benner

###### Abstract

Accurate error estimation is crucial in model order reduction, both to obtain
small reduced-order models and to certify their accuracy when deployed in
downstream applications such as digital twins. In existing a posteriori error
estimation approaches, knowledge about the time integration scheme is
mandatory, e.g., the residual-based error estimators proposed for the reduced
basis method. This poses a challenge when automatic ordinary differential
equation solver libraries are used to perform the time integration. To address
this, we present a data-enhanced approach for _a posteriori_ error estimation.
Our new formulation enables residual-based error estimators to be independent
of any time integration method. To achieve this, we introduce a corrected
reduced-order model which takes into account a data-driven closure term for
improved accuracy. The closure term, subject to mild assumptions, is related
to the local truncation error of the corresponding time integration scheme. We
propose efficient computational schemes for approximating the closure term, at
the cost of a modest amount of training data. Furthermore, the new error
estimator is incorporated within a greedy process to obtain parametric
reduced-order models. Numerical results on three different systems show the
accuracy of the proposed error estimation approach and its ability to produce
ROMs that generalize well.

###### keywords:

A posteriori error estimation, Model order reduction, Reduced basis method,
ODE solvers

  * •

We present a new _a posteriori_ error estimation framework for nonlinear
dynamical systems, which is independent of the ordinary differential equation
solver used.

  * •

Our approach imposes a simpler, user-defined time integration scheme on the
high-fidelity snapshots data and learns a parametrized closure term to recover
the true residual needed for error estimation.

  * •

We numerically demonstrate the approximability of the data-driven closure term
and propose efficient algorithms for its computation.

  * •

The proposed methodology eliminates one of the main intrusive elements present
in the reduced basis method, viz., the knowledge of the time integration
scheme; it further paves way for the reduced basis method to become better
integrated within existing computational software packages such as matlab® or
Python.

##  1 Introduction

Model order reduction (MOR) has become an important enabling technology that
facilitates the rapid and reliable simulations of large-scale systems in a
number of scientific disciplines. The central goal of MOR is to replace an
_expensive-to-compute_ full-order model (FOM) with a surrogate, called the
reduced-order model (ROM). The ROM has fewer degrees of freedom which enables
real-time simulation – an important requirement in recent developments such as
_digital twins_ [34, 40]. We refer to the recent books [5, 6, 7] for a
detailed background and the state-of-the-art in MOR literature.

In this work, we are interested in obtaining ROMs for parametric, nonlinear
dynamical systems that arise from the numerical discretization of partial
differential equations (PDEs). For such systems, the Reduced Basis Method
(RBM) [49, 46, 35] is a commonly used MOR approach. The RBM relies on an _a
posteriori_ error estimator to perform an iterative, greedy sampling of the
parameter space to compute solution snapshots, in order to build a linear
subspace that serves as an approximation for the solution space. The FOM
equations are then projected onto this subspace to obtain the final ROM. _A
posteriori_ , residual-based error bounds/estimators are an indispensable part
of the RBM; such estimators have been proposed for a variety of system classes
such as coercive and non-coercive elliptic systems [43, 48], linear parabolic
systems [29], non-linear or non-affine systems [52, 12, 30]. While initial
development in the RBM community focused on deriving _a posteriori_ error
estimators in the variational setting, more recent work has also focused on
deriving such estimators for PDEs that are discretized using finite volume
method or the finite difference method; see [33, 25, 55, 57, 18].

###  1.1 Motivation

For dynamical systems, all of the above works on _a posteriori_ error
estimation assume prior knowledge of the time integration scheme. The
expression of the time integration method is needed to compute the residual
vector, which is obtained by plugging an approximate solution in the residual
expression. For instance, if a Runge-Kutta method were used, the Butcher
tableau of the corresponding scheme needs to be known. However, there are
scenarios where no (or only incomplete) knowledge of the time integration
scheme is available. One such example is when automatic ordinary differential
equation (ODE) solver libraries, e.g., SUNDIALS [26, 37], ODEPACK [36], ARKODE
[47] are used to perform the time integration. Such ODE solvers are often
available readily in computational software such as matlab® (ode45, ode15s,
etc.) or Python (scipy.integrate.odeint). Another popular ODE solver library
is the TS library [1] available in PETSc [3]. Typically, automatic ODE solver
packages are computationally efficient as they employ adaptive order selection
(e.g., ode15s can switch between order 1 to order 5 methods) and/or adaptive
choice of the time step, which is beneficial in case the system being
considered exhibits _stiffness_. Moreover, much time is saved by re-using
efficient and robust software which is already available and well-maintained.
Whenever automatic ODE solvers are used, the standard RBM can no longer be
used in a straightforward manner to obtain a ROM. This is owing to the fact
that the exact expression of the time integration method used within such ODE
solvers is unknown and deriving the corresponding error estimator is an open
question. This serves as the fundamental motivation for our work.

###  1.2 Main contributions

Our work introduces a new data-enhanced paradigm for _a posteriori_ ,
residual-based error estimation for the RBM, in the absence of any knowledge
about the time integration scheme used. In addition to the discretized system
matrices, our approach only requires access to solution snapshots at different
time instances obtained from any black-box ODE solver in a library, at a small
set of system parameter samples. In our approach, we fit a (simpler) time
integration scheme of our choice to the available snapshots data. Since the
time integration method we impose on the data does not coincide with the
actual scheme used to generate the data in first place, there is a mismatch or
_defect_ at each time step. Subject to a mild assumption, this defect
corresponds to the local truncation error (LTE) of the imposed time
integration scheme. We then learn this defect term as a function of time and
the system parameters. Learning the defect using data allows us to formulate a
new _corrected_ ROM that takes into account the LTE at each time step. Taking
into account the LTE as a closure term means that the corrected ROM can
_recover_ the solution obtained from a ROM solved by a solver from an ODE
library. We use the _corrected_ ROM to compute a good approximation to the
true residual vector at each time step, leading to accurate estimation of the
true error. The computation of the defect term is a key element of our
proposed methodology. To do this efficiently, we rely on two observations.
First, we demonstrate with numerical evidence that the defect term possesses a
certain low-rank structure in space that allows us to efficiently project it
on a low-dimensional subspace. Second, if the solution to the FOM is fairly
smooth over the parameter space, then this smoothness carries over to the
defect term and it can be efficiently approximated with respect to parameter
variations using a suitable surrogate model.

###  1.3 Prior work

A number of recent works have sought to extend the RBM to a non-intrusive
framework, mainly for steady systems; see [16, 14, 15]. In [15] and [16] the
non-intrusive reduced basis (NIRB) method makes use of FEM solutions computed
on two different meshes – one fine and one coarse – to estimate the state
error due to the ROM. A different approach is proposed in [14] that involves
the empirical interpolation method (EIM) to precompute an affine decomposition
consisting of parameter-dependent and parameter-independent quantities. More
recently, the method in [31] has sought to extend the NIRB approach to time-
dependent systems. All of these methods assume no knowledge of the FOM system
matrices and rely purely on snapshots of the state variable. But, they need
knowledge of the time integration scheme used to numerically integrate the
spatially discretized system. In our new approach, we assume knowledge of the
system matrices. Therefore, our approach is an intrusive method in terms of
access to the model or the system matrices. However, our approach is non-
intrusive in terms of the time integration scheme (ODE solver) used.
Nevertheless, we do envisage an extension of our method to the case when there
is no access to the system matrices; this would be a subject for future
investigation.

To eliminate the dependence on the time integration scheme, our approach aims
to increase the accuracy of a user-imposed time integration method by addition
of a data-driven closure term. The problem of improving time integration
accuracy is an emerging topic and has received attention in the machine
learning community [51, 45, 39]. The deep Euler method (DEM) is introduced in
[51] where the authors’ motivation is to improve the first-order accuracy of
the explicit Euler time integration scheme by a factor of $\eta$, i.e., to
$\mathcal{O}(\eta\delta t)$ ($\eta\ll 1$). They achieve this by approximating
the LTE using a feed-forward neural network (FNN). Extension of the approach
to other time integration schemes beyond the explicit Euler method are also
illustrated. In the same spirit of the DEM approach, the work [45] introduces
hypersolvers, which are targeted at Neural ODEs [21]. While both DEM and
hypersolvers seek to improve the accuracy of the time integration scheme by
approximating the LTE, they differ in the assumptions made about the model.
DEM assumes the model of the PDE is known and is exact, while for
hypersolvers, the model is unknown and is approximated by a separate neural
network. Our work here differs from DEM and hypersolvers in several ways.
Firstly, our aim is to obtain a corrected ROM for better error estimation
whereas, in the aforementioned two methods the aim is simply to improve the
FOM accuracy during simulation. Secondly, our targets are parametric nonlinear
systems. Therefore, we want to learn the LTE for different time instances for
a range of system parameters. Both DEM and hypersolvers are limited to non-
parametric systems. While either method can potentially be extended to the
parametric case, it is computationally very expensive. Furthermore, both these
methods learn a closure term which is a function of the past states. Our
approach treats the closure term as a function of time and any additional
system parameters. This allows us to take advantage of the special structure
and smoothness properties of the closure term.

The rest of this paper is organized as follows. In Section 2, we present the
mathematical preliminaries of MOR and provide a short recap of the RBM along
with _a posteriori_ error estimation. We also illustrate the pitfalls of the
current error estimation approach used in RBM, through the example of the heat
equation. Section 3 contains the main contributions of this manuscript. We
introduce a user-imposed time integration scheme which incorporate the data-
driven closure term and derive the _a posteriori_ output error estimator
suited for situations where ODE solver libraries are used within the RBM. The
algorithmic aspects of the proposed approach are discussed in Section 4 while
numerical results are presented in Section 5 to support our new method. We
conclude with a summary and topics for future research in Section 6.

##  2 Mathematical background

Table 1: List of FOM variables Quantity | Variable | Dimension
---|---|---
State vector | $\mathbf{x}$ | $\mathbb{R}^{N}$
Input vector | $\mathbf{u}$ | $\mathbb{R}^{N_{I}}$
Output vector | $\mathbf{y}$ | $\mathbb{R}^{N_{O}}$
Initial condition | $\mathbf{x}_{0}$ | $\mathbb{R}^{N}$
Nonlinear vector | $\mathbf{f}$ | $\mathbb{R}^{N}$
Approximate state vector | $\widetilde{\mathbf{x}}:=\mathbf{V}\widehat{\mathbf{x}}$ | $\mathbb{R}^{N}$
System matrix | $\mathbf{A}$ | $\mathbb{R}^{N\times N}$
Input matrix | $\mathbf{B}$ | $\mathbb{R}^{N\times N_{I}}$
Output matrix | $\mathbf{C}$ | $\mathbb{R}^{N_{O}\times N}$
Table 2: List of ROM variables Quantity | Variable | Dimension
---|---|---
Left projection matrix | $\mathbf{V}$ | $\mathbb{R}^{N\times n}$
Right projection matrix | $\mathbf{W}$ | $\mathbb{R}^{N\times n}$
State vector | $\widehat{\mathbf{x}}$ | $\mathbb{R}^{n}$
Output vector | $\widehat{\mathbf{y}}$ | $\mathbb{R}^{N_{O}}$
Initial condition | $\widehat{\mathbf{x}}_{0}:=\mathbf{W}^{\textsf{T}}\mathbf{x}$ | $\mathbb{R}^{n}$
Nonlinear vector | $\widehat{\mathbf{f}}:=\mathbf{W}^{\textsf{T}}\mathbf{f}(\widetilde{\mathbf{x}})$ | $\mathbb{R}^{n}$
System matrix | $\widehat{\mathbf{A}}:=\mathbf{W}^{\textsf{T}}\mathbf{A}\mathbf{V}$ | $\mathbb{R}^{n\times n}$
Input matrix | $\widehat{\mathbf{B}}:=\mathbf{W}^{\textsf{T}}\mathbf{B}$ | $\mathbb{R}^{n\times N_{I}}$
Output matrix | $\widehat{\mathbf{C}}:=\mathbf{C}\mathbf{V}$ | $\mathbb{R}^{N_{O}\times n}$

Consider the following parametric system of ODEs:

$\displaystyle\frac{d}{dt}\mathbf{x}(\bm{\mu})$
$\displaystyle=\mathbf{A}(\bm{\mu})\mathbf{x}(\bm{\mu})+\mathbf{f}(\mathbf{x},\bm{\mu})+\mathbf{B}\mathbf{u}(t),\qquad\mathbf{x}(0)=\mathbf{x}_{0},$
(1a) $\displaystyle\mathbf{y}(\bm{\mu})$
$\displaystyle=\mathbf{C}(\bm{\mu})\mathbf{x}(\bm{\mu}).$ (1b)

Such systems often arise upon discretizing a PDE using numerical
discretization schemes. Table 1 lists each variable and its corresponding
dimension. In subsequent discussions, we refer to eq. 1 as the FOM. The state
vector $\mathbf{x}(\bm{\mu})$ can be obtained for any time instance for a
given parameter $\bm{\mu}\in\mathcal{P}\subset\mathbb{R}^{p}$ by time
integration of the FOM using any desired method like Runge-Kutta methods or
linear multi-step methods. The set $\mathcal{P}$ denotes the parameter space
and $p$ is its dimension. The number of equations $N$ in eq. 1 is often very
large to ensure a high-fidelity solution of the underlying physical process.
This poses a major challenge for its numerical solution, especially for many
instances of the parameter $\bm{\mu}$ in applications such as control and
uncertainty quantification.

###  2.1 Model order reduction

Projection-based MOR techniques offer a systematic approach to obtain a ROM
for eq. 1 in the following form:

$\displaystyle\frac{d}{dt}\widehat{\mathbf{x}}(\bm{\mu})$
$\displaystyle=\widehat{\mathbf{A}}(\bm{\mu})\widehat{\mathbf{x}}(\bm{\mu})+\widehat{\mathbf{f}}(\widetilde{\mathbf{x}},\bm{\mu})+\widehat{\mathbf{B}}\mathbf{u}(t),\qquad\widehat{\mathbf{x}}(0)=\mathbf{W}^{\textsf{T}}\mathbf{x}_{0},$
(2a) $\displaystyle\widehat{\mathbf{y}}(\bm{\mu})$
$\displaystyle=\widehat{\mathbf{C}}(\bm{\mu})\widehat{\mathbf{x}}(\bm{\mu}).$
(2b)

The above ROM is derived based on the ansatz
$\mathbf{x}\approx\widetilde{\mathbf{x}}:=\mathbf{V}\widehat{\mathbf{x}}$
applied to eq. 1. The resulting over-determined system of equations is reduced
by a Petrov-Galerkin projection, leading to eq. 2. Table 2 lists the reduced
variables and their corresponding dimensions. The number of equations $n$ in
eq. 2a is often much smaller than $N$ in eq. 1a, i.e., $n\ll N$. Therefore,
the ROM can be readily used for repeated simulations given any new values of
the parameter $\bm{\mu}$. Different MOR techniques differ in how they compute
the projection matrices $\mathbf{V},\mathbf{W}$. When $\mathbf{V}=\mathbf{W}$,
it is referred to as Galerkin projection. In the sequel, we will limit
ourselves to ROMs obtained through a Galerkin projection.

######  Remark 1.

Although eq. 2a is a ROM of dimension $n\ll N$, evaluating
$\widehat{\mathbf{f}}(\widetilde{\mathbf{x}},\bm{\mu})$ involves operations
that scale with the FOM dimension $N$ (since
$\widetilde{\mathbf{x}}=\mathbf{V}\widehat{\mathbf{x}}$ needs to be evaluated
for each time step). Thus, a direct evaluation of eq. 2a may not offer any
computational speedup over evaluating the FOM. Hyperreduction techniques [4,
17, 13] can be used to address this issue. We employ the discrete empirical
interpolation method (DEIM) for the hyperreduction in our numerical results in
Section 5.

###  2.2 Reduced basis method

The RBM is a greedy approach that builds a global projection matrix
$\mathbf{V}$ to obtain a parametric ROM for eq. 1. Since its introduction, the
RBM has been a highly successful approach to obtain ROMs for parametric
systems in a variety of applications such as process engineering [56, 18],
geosciences [23], data assimilation [24, 41], uncertainty quantification [20]
to name just a few.

When RBM is applied to dynamical systems, the POD-Greedy method [32] is
adopted. It consists of a greedy sampling in the parameter space and a
compression of the time trajectory through singular value decomposition (SVD).
Algorithm 1 sketches the pseudo-code for the standard RBM using the POD-Greedy
algorithm. In Step 1, solver denotes the time integration scheme used to solve
eqs. 1 and 2. This is determined _a priori_ by the user based on the ODE
library used, e.g. scipy.ode.integrate.

Consider a fine discretization of the parameter space $\mathcal{P}$ in the
form of a training set
$\Xi:=\\{\bm{\mu}_{1},\bm{\mu}_{2},\ldots,\bm{\mu}_{N_{p}}\\}$ consisting of
samples of the parameter $\bm{\mu}$. The POD-Greedy method starts by solving
the FOM at a randomly chosen parameter $\bm{\mu}^{*}\in\Xi$. Performing the
SVD of the solution snapshot matrix

$\displaystyle\mathbf{X}(\bm{\mu}^{*}):=\big{[}\mathbf{x}(t_{0}),\mathbf{x}(t_{1}),\ldots,\mathbf{x}(t_{K})\big{]}\in\mathbb{R}^{N\times
N_{t}}$ (3)

with $N_{t}=K+1$ yields the projection basis $\mathbf{V}$. More precisely, we
first obtain

$\mathbf{X}(\bm{\mu}^{*})=\mathbf{V}_{\bm{\mu}^{*}}\bm{\Sigma}_{\bm{\mu}^{*}}(\mathbf{W}_{\bm{\mu}^{*}})^{\textsf{T}}$

with $\mathbf{V}_{\bm{\mu}^{*}}\in\mathbb{R}^{N\times
N},\bm{\Sigma}_{\bm{\mu}^{*}}\in\mathbb{R}^{N\times N_{t}}$ and
$\mathbf{W}_{\bm{\mu}^{*}}\in\mathbb{R}^{N_{t}\times N_{t}}$. The matrix
$\bm{\Sigma}_{\bm{\mu}^{*}}$ is a rectangular diagonal matrix and contains the
singular values $\sigma_{i}$ in the locations
$\bm{\Sigma}_{\bm{\mu}^{*}}(i,\,i),i\in\\{1,2,\ldots,\min(N,N_{t})\\}$ .

We then update $\mathbf{V}$ by enriching it with the first $r_{c}$ left
singular vectors of $\mathbf{X}(\bm{\mu}^{*})$, i.e.,
$\mathbf{V}_{\bm{\mu}^{*}}(:\,,\,1:r_{c})$. In subsequent iterations,
snapshots of the FOM are similarly collected at different values of
$\bm{\mu}^{*}$ and the matrix $\mathbf{V}$ is updated with new information.
Note, however, that for all iterations after the first, we perform a SVD of
the matrix $\overline{\mathbf{X}}$ obtained after removing from
$\mathbf{X}(\bm{\mu}^{*})$ the information already represented in
$\mathbf{V}$, i.e., we set

$\overline{\mathbf{X}}:=\mathbf{X}(\bm{\mu}^{*})-\mathbf{V}(\mathbf{V}^{\textsf{T}}\mathbf{X}(\bm{\mu}^{*})).$

To ensure good conditioning, it is recommended to perform a Gram-Schmidt
orthonormalization after each update of $\mathbf{V}$. The choice of
$\bm{\mu}^{*}$ at each iteration is determined through an error estimator
$\Delta(\bm{\mu})$ as follows:

$\bm{\mu}^{*}=\operatorname*{argmin}_{\bm{\mu}\in\Xi}\Delta(\bm{\mu}).$

The error estimator serves as an upper bound for the true state error
$\|\mathbf{x}(t,\bm{\mu})-\widetilde{\mathbf{x}}(t,\bm{\mu})\|$ (or true
output error $\|\mathbf{y}(t,\bm{\mu})-\widehat{\mathbf{y}}(t,\bm{\mu})\|$).
That is,

$\displaystyle\|\mathbf{x}(t,\bm{\mu})-\widetilde{\mathbf{x}}(t,\bm{\mu})\|$
$\displaystyle\leq\Delta(\bm{\mu}),\qquad\text{or}$
$\displaystyle\|\mathbf{y}(t,\bm{\mu})-\widehat{\mathbf{y}}(t,\bm{\mu})\|$
$\displaystyle\leq\Delta(\bm{\mu}).$

All that is needed to evaluate $\Delta(\bm{\mu})$ is to solve the ROM eq. 2a
which can be formulated at each greedy iteration using $\mathbf{V}$.

Algorithm 1 POD-Greedy algorithm

0: Training set $\Xi$, tolerance (tol), Discretized model
($\mathbf{E},\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{f},\mathbf{x}_{0}$)

0: $\mathbf{V}$

1: Initialize $\mathbf{V}=[\,\,]$, $\epsilon=1+\texttt{tol}$, choose solver,
greedy parameter $\bm{\mu}^{*}$

2: while $\epsilon>\texttt{tol}$ do

3: Obtain FOM snapshots $\mathbf{X}(\bm{\mu}^{*})$ at $\bm{\mu}^{*}$ with
solver

4: Determine $\mathbf{V}_{\bm{\mu}^{*}}$ through an SVD of
$\overline{\mathbf{X}}:=\mathbf{X}(\bm{\mu}^{*})-\mathbf{V}(\mathbf{V}^{\textsf{T}}\mathbf{X}(\bm{\mu}^{*}))$

5: Update $\mathbf{V}$ as
$\mathbf{V}:=\texttt{orth}\big{(}\mathbf{V},\mathbf{V}_{\bm{\mu}^{*}}(:\,,\,1:r_{c})\big{)}$with
$\texttt{orth}\big{(}\cdot\big{)}$ denoting an orthogonalization process which
can be implemented using the modified Gram-Schmidt process, or QR algorithm

6: Get
($\widehat{\mathbf{E}},\widehat{\mathbf{A}},\widehat{\mathbf{B}},\widehat{\mathbf{C}},\widehat{\mathbf{f}},\widehat{\mathbf{x}}_{0}$)
in eq. 2 by Galerkin proj. (+ hyperreduction)

7: Solve ROM with solver

8: $\bm{\mu}^{*}:=\arg\max\limits_{\bm{\mu}\in\Xi}\Delta(\bm{\mu})$

9: Set $\epsilon=\Delta(\bm{\mu}^{*})$

10: end while

###  2.3 A posteriori error estimation for the RBM

The standard error estimation approach in the RBM literature is residual-based
[46, 35]. In order to derive the error estimator $\Delta(\bm{\mu})$, knowledge
of the time discretization scheme used to integrate the FOM and the ROM is
assumed, e.g., using implicit Euler, Crank-Nicolson method, or an implicit-
explicit (IMEX) method. Computing $\Delta(\bm{\mu})$ for a given parameter
involves determining the residual vector $\textbf{r}\in\mathbb{R}^{N}$ (or its
norm) at each time instance; some approaches to error estimation also require
the residual of a dual or adjoint system. Let us illustrate this by means of
an example.

Suppose eq. 1 is discretized in time using a first-order IMEX scheme [2]. The
linear part (involving $\mathbf{A}(\bm{\mu})$) is discretized implicitly,
while the nonlinear vector $\mathbf{f}(\mathbf{x})$ is evaluated explicitly.
The resulting discretized system reads

$\displaystyle\mathbf{E}_{\text{im}}\mathbf{x}_{\text{im}}^{k}$
$\displaystyle=\mathbf{A}_{\text{im}}\mathbf{x}_{\text{im}}^{k-1}+\delta
t\big{(}\mathbf{f}(\mathbf{x}_{\text{im}}^{k-1})+\mathbf{B}\mathbf{u}^{k}\big{)},$
(4a) $\displaystyle\mathbf{y}_{\text{im}}^{k+1}$
$\displaystyle=\mathbf{C}\mathbf{x}_{\text{im}}^{k}$ (4b)

with $\mathbf{E}_{\text{im}}:=(\mathbf{I}_{N}-\delta t\mathbf{A})$,
$\mathbf{A}_{\text{im}}:=\mathbf{I}_{N}$ where
$\mathbf{I}_{N}\in\mathbb{R}^{N\times N}$ is the identity matrix. For better
clarity, we have not shown the parameter dependence of the system matrices and
vectors. The ROM eq. 2 can be discretized in the same way and reads

$\displaystyle\widehat{\mathbf{E}}_{\text{im}}\widehat{\mathbf{x}}_{\text{im}}^{k}$
$\displaystyle=\widehat{\mathbf{A}}_{\text{im}}\widehat{\mathbf{x}}_{\text{im}}^{k-1}+\delta
t\big{(}\widehat{\mathbf{f}}(\widetilde{\mathbf{x}}_{\text{im}}^{k-1})+\widehat{\mathbf{B}}\mathbf{u}^{k}\big{)},$
(5a) $\displaystyle\widehat{\mathbf{y}}_{\text{im}}^{k+1}$
$\displaystyle=\widehat{\mathbf{C}}\widehat{\mathbf{x}}_{\text{im}}^{k}$ (5b)

with $\widehat{\mathbf{E}}_{\text{im}}:=(\mathbf{I}_{n}-\delta
t\widehat{\mathbf{A}})$, $\widehat{\mathbf{A}}_{\text{im}}:=\mathbf{I}_{n}$
with $\mathbf{I}_{n}\in\mathbb{R}^{n\times n}$ being the identity matrix. The
residual arising due to the ROM approximation can be computed by substituting
the approximate state vector $\widetilde{\mathbf{x}}_{\text{im}}^{k}$ into eq.
4a. The resulting residual at the $k$-th time step, $\textbf{r}^{k}$, reads

$\displaystyle\textbf{r}^{k}:=\mathbf{A}_{\text{im}}\widetilde{\mathbf{x}}_{\text{im}}^{k-1}+\delta
t\big{(}\mathbf{f}(\widetilde{\mathbf{x}}_{\text{im}}^{k-1})+\mathbf{B}\mathbf{u}^{k}\big{)}-\mathbf{E}_{\text{im}}\widetilde{\mathbf{x}}_{\text{im}}^{k}.$
(6)

It is clear that in order to obtain the residual $\textbf{r}^{k}$, the time
discretization scheme for the ROM should be the same as that for the FOM so
that $\widetilde{\mathbf{x}}_{\text{im}}^{k}$ in eq. 6 and
$\mathbf{x}_{\text{im}}^{k}$ in eq. 4 correspond to the same time instance
$t^{k}$. An _a posteriori_ error bound $\Delta(\bm{\mu})$ for the
approximation error
$\mathbf{e}^{k}(\bm{\mu}):=\mathbf{x}_{\text{im}}^{k}(\bm{\mu})-\widetilde{\mathbf{x}}_{\text{im}}^{k}(\bm{\mu})$
for a given parameter $\bm{\mu}$ can be computed based on the residual as
below.

######  Theorem 2.1 (Residual-based error bound).

Suppose that the nonlinear quantity $\mathbf{f}(\mathbf{x},\bm{\mu})$ is
Lipschitz continuous in the first argument for all $\bm{\mu}$ such that there
exists a constant $L_{\mathbf{f}}$ for which

$\|\mathbf{f}(\mathbf{x},\bm{\mu})-\mathbf{f}(\widetilde{\mathbf{x}}_{\textnormal{im}},\bm{\mu})\|\leq
L_{\mathbf{f}}\|\mathbf{x}-\widetilde{\mathbf{x}}_{\textnormal{im}}\|.$

Further assume that for any parameter $\bm{\mu}$ the projection error at the
first time step is
$\|\mathbf{e}^{0}(\bm{\mu})\|=\|\mathbf{x}_{\textnormal{im}}^{0}(\bm{\mu})-\widetilde{\mathbf{x}}_{\textnormal{im}}^{0}(\bm{\mu})\|=\|\mathbf{x}_{\textnormal{im}}^{0}(\bm{\mu})-\mathbf{V}\mathbf{V}^{\textsf{T}}\mathbf{x}^{0}(\bm{\mu})\|$.

The error in the state variable at the $k$-th time step,
$\|\mathbf{e}^{k}(\bm{\mu})\|=\|\mathbf{x}^{k}_{\textnormal{im}}(\bm{\mu})-\widetilde{\mathbf{x}}_{\textnormal{im}}^{k}(\bm{\mu})\|$
is given by the following expression:

$\displaystyle\|\mathbf{e}^{k}(\bm{\mu})\|\leq\Delta^{k}(\bm{\mu}):=\xi(\bm{\mu})^{k}\left\lVert\mathbf{e}^{0}(\bm{\mu})\right\rVert_{2}+\sum\limits_{i=1}^{k}\zeta(\bm{\mu})\cdot\big{(}\xi(\bm{\mu})\big{)}^{k-i}\cdot\|\mathbf{r}^{i}(\bm{\mu})\|_{2}$
(7)

where $\zeta(\bm{\mu}):=\|\big{(}\mathbf{E}(\bm{\mu})\big{)}^{-1}\|_{2}$ and
$\xi(\bm{\mu}):=\Big{(}\|\big{(}\mathbf{E}(\bm{\mu})\big{)}^{-1}\mathbf{A}(\bm{\mu})\|_{2}+\delta
t\,L_{\mathbf{f}}\|\big{(}\mathbf{E}(\bm{\mu})\big{)}^{-1}\|_{2}\Big{)}$.

###### Proof.

See Appendix A. ∎

Residual-based error bounds such as the one above are already available in the
RBM literature [32, 28] for both linear and nonlinear systems.

In many applications, only a small set of variables which are obtained as a
linear combination of the state variables are of interest. These are typically
called the output quantities of interest. Goal-oriented _a posteriori_ error
estimators for the output of interest have been discussed in several works
[57, 18, 28, 33]. They typically involve the residual of a dual or an adjoint
system. The general form of the goal-oriented error estimators in [28, 33] is

$\displaystyle\|\mathbf{y}^{k}(\bm{\mu})-\widehat{\mathbf{y}}^{k}(\bm{\mu})\|\leq\Gamma(\bm{\mu})\cdot\sum\limits_{i=1}^{k}\|\textbf{r}^{i}\|\cdot\sum\limits_{j=k}^{K}\|\textbf{r}_{\text{du}}^{j}\|$
(8)

where $\|\textbf{r}_{\text{du}}^{j}\|$ denotes the residual of an
appropriately defined dual system at the $j$-th time instance and
$\Gamma(\bm{\mu})$ is some parameter-dependent constant. Improving upon this
expression, a different goal-oriented error estimator was proposed in [57, 18]
having the following form:

$\displaystyle\|\mathbf{y}^{k}(\bm{\mu})-\widehat{\mathbf{y}}^{k}(\bm{\mu})\|\leq\breve{\Gamma}(\bm{\mu})\cdot\|\textbf{r}^{k}\|\cdot\|\textbf{r}_{\text{du}}^{k}\|$
(9)

with $\breve{\Gamma}(\bm{\mu})$ being a parameter-dependent constant. This
error estimator avoids the accumulation of the residuals over time, which is a
major drawback for the goal-oriented error estimator in eq. 8 and also the
state error estimator in eq. 7. It is also to be noted that the ROM resulting
from the use of a goal-oriented error estimator often turns out to be of a
smaller dimension.

Having reviewed the main ideas of RBM and the standard _a posteriori_ error
estimator, we next illustrate the disadvantage of this standard approach when
ODE solver libraries are used to solve the FOM/ROM within Algorithm 1.

###  2.4 RBM with ODE solver libraries

Automatic ODE solver packages are implemented and readily available in a
number of open-source and proprietary computational software. The matlab® ODE
Suite [50], for instance, implements both linear multi-step (ode15s, ode23s,
etc.) and Runge-Kutta-type solvers (ode45, ode23, etc.). In Python, using the
scipy module once can access the odeint and solve_ivp submodules, both
providing access to a variety of standard time integration schemes. In
addition to these, there are several stand-alone libraries to solve ODEs such
as SUNDIALS [37], ODEPACK [36], ARKODE [47]. All of the mentioned libraries
implement adaptive order and adaptive time-stepping which makes them highly
efficient for a variety of problems, e.g., problems exhibiting stiffness.

When adaptive ODE solver packages are used within the RBM (i.e., in Steps 3
and 7 of Algorithm 1, solver is replaced by the chosen method from the
package), the standard error estimation approach (see Section 2.3) becomes
less straightforward. One can no longer write the corresponding expression of
the residual resulting from the ROM (e.g., eq. 6). The reason for this is that
the exact expression of the time integration method used is unknown, as the
solver adaptively varies the time step and/or the order of the scheme. If a
residual expression obtained from a user-imposed time discretization scheme,
e.g., eq. 6 is used, the resulting error estimator eq. 7 no longer gives an
efficient bound for the true error, i.e., the error between
$\mathbf{V}\widehat{\mathbf{x}}^{k}(\bm{\mu})$ and $\mathbf{x}^{k}(\bm{\mu})$
obtained from simulating the ROM eq. 2 and the FOM eq. 1, respectively, using
an ODE solver from a library. We illustrate this next by obtaining a ROM for
the simple non-parametrized linear heat diffusion model and estimating its
error.

###  2.5 Example: ROM for the linear heat equation

We consider the linear heat equation in 1-D over the domain $\Omega=[0,1]$ and
time $t\in[0,1]$

$\displaystyle\frac{\partial}{\partial
t}v(z,t;\mu)-\mu\frac{\partial^{2}}{\partial z^{2}}v(z,t;\mu)=0$ (10)

where $v(z,t;\mu)$ is the state, $z\in\Omega$ is the spatial variable and the
viscosity $\mu=0.06$. We further impose Dirichlet boundary conditions
$v(0,t;\mu)=v(1,t;\mu)=0,\,\forall t\in[0,1]$. We also fix the output variable
of interest as the value of the state at the node next to the right boundary.
Employing the finite difference method, we discretize the domain $\Omega$ in
an equidistant fashion, with a grid size of $2^{-8}$. With this, the resulting
discretized ODE can be written in a form similar to eq. 1 as

$\displaystyle\frac{d}{dt}\mathbf{x}=\mathbf{A}\mathbf{x},\qquad\mathbf{x}(0)=\mathbf{x}_{0}$
(11)

where $\mathbf{x}\in\mathbb{R}^{255}$ is the discretized state vector and
$\mathbf{A}\in\mathbb{R}^{255\times 255}$. We let the initial condition
$\mathbf{x}_{0}$ to be a normal distribution, i.e.,

$\mathbf{x}_{0}:=\dfrac{1}{\sigma\sqrt{2\pi}}e^{-\dfrac{1}{2}\cdot\bigg{(}\dfrac{z-\mathrm{m}}{\sigma}\bigg{)}^{2}}$

with a mean $\mathrm{m}=0.5$ and a standard deviation of $\sigma=0.15$. Our
aim is to obtain a ROM for eq. 11 and use the state and output error
estimators in eqs. 7 and 8, respectively, to quantify the error of the ROM.
For the output error estimator, we use the one proposed in [18]. To this end,
we consider the standard POD-based ROM. This involves obtaining snapshots of
the state vector in eq. 11 for different time instances. The SVD of the
resulting snapshot matrix is used to obtain the projection matrix
$\mathbf{V}$. Figure 1 illustrates the solution to eq. 11 at $\mu=0.06$. We
also see the exponential decay of the singular values of the snapshot matrix.

Figure 1: Heat equation. Left: solution to the parametrized heat equation eq.
11 at $\mu=0.06$; Right: normalized singular values
$\sigma_{i}/\sigma_{1},i\in\\{1,2,\ldots,100\\}$.

For this example, we compute $\mathbf{V}$ consisting of the first $12$ columns
of the left singular vector matrix. For integrating eq. 11, we use the odeint
function available in the scipy package for Python. We note that odeint is a
wrapper around the LSODA solver available in the Fortran library odepack.
LSODA implements adaptive time-stepping. Moreover, it switches between methods
for non-stiff and stiff problems automatically. We use the same solver to
integrate the ROM. Once the reduced matrices are obtained through a Galerkin
projection using $\mathbf{V}$, the first step involved in estimating the error
is to integrate the ROM using odeint to evaluate and plug-in the approximate
solution $\widetilde{\mathbf{x}}$ at each time instance, into the odeint time
discretization scheme to get the residual. The residual operator at the $k$-th
time instance has the general expression:

$\displaystyle\mathcal{R}^{k}\big{[}\widetilde{\mathbf{x}}^{k},\widetilde{\mathbf{x}}^{k-1},\ldots,\widetilde{\mathbf{x}}^{k-s}\big{]}$
(12)

where the arguments for $\mathcal{R}^{k}\big{[}\cdot\big{]}$ could be the
solutions at the current and past $s$ time steps (in case the scheme used is a
linear multi-step method) or at $s$ different stage solutions (in case an
$s$-stage Runge-Kutta scheme is used). The exact form of this expression,
naturally, is dependent on the time integration method that was used to obtain
the snapshots. As we are using a solver from the odepack library, knowing the
residual operator expression is complicated and often impossible. To
circumvent this, one may choose to use a _different, but known_ time
integration method in order to compute the residual, e.g., via eq. 6. However,
this will lead to erroneous results, as we will demonstrate next. We denote by
$\textbf{r}^{k}$ the output of the residual operator
$\mathcal{R}^{k}\big{[}\cdot\big{]}$ at a given set of arguments.

Figure 2: Estimated and true errors for the heat equation eq. 11 for
$\mu=0.06$ obtained by imposing a backward Euler method on the ROM snapshots
obtained from odeint. Top: the estimated state error and the corresponding
true error (see eq. 7); Bottom: the estimated output error and the
corresponding true error (see eq. 8).

Suppose, we use a different time integration scheme, say a first-order
backward Euler method as an “approximation” to the true method used within
odeint. In this case, the approximation to the true residual reads

$\displaystyle\widetilde{\textbf{r}}^{k}=\widetilde{\mathbf{x}}^{k-1}-\bigg{(}\mathbf{I}_{N}-\delta
t\mathbf{A}\bigg{)}\widetilde{\mathbf{x}}^{k}$ (13)

with the approximate residual operator
$\widetilde{\mathcal{R}}^{k}\big{[}\ast,\star]=\mathbf{I}_{N}\ast-\bigg{(}\mathbf{I}_{N}-\delta
t\mathbf{A}\bigg{)}\star$ and $\ast,\star$ being placeholders for the
arguments of $\widetilde{\mathcal{R}}$. In general,
$\mathcal{R}^{k}\big{[}\cdot\big{]}$ and
$\widetilde{\mathcal{R}}^{k}\big{[}\cdot\big{]}$ are quite different, leading
to rather inaccurate estimation of the error when the latter is used in eq. 7.
This is illustrated in Figure 2. The top figure shows the estimated state
error obtained using eq. 7 for every time step and the corresponding true
error, both measured in the $2$-norm. The bottom figure illustrates the
estimated error for the output variable and its true error.

The estimated error measured using the wrong expression of the residual
overestimates by $2$ orders of magnitude in the best case. In this work, we
propose a scheme to suitably modify
$\widetilde{\mathcal{R}}^{k}\big{[}\cdot\big{]}$ using a closure term, such
that the resulting expression for the residual is close to the one evaluated
by $\mathcal{R}^{k}\big{[}\cdot\big{]}$. We limit our focus to the case of
output error estimation though the proposed closure technique for correcting
the residual can be straight-forwardly applied to any residual-based error
estimator.

##  3 Improving output error estimation via a data-enhanced closure approach

In this section, we introduce a data-enhanced closure strategy to ensure that
the residual resulting from the user-imposed time integration scheme, viz.
$\widetilde{\mathcal{R}}^{k}\big{[}\widetilde{\mathbf{x}}^{k},\widetilde{\mathbf{x}}^{k-1},\ldots,\widetilde{\mathbf{x}}^{k-s}\big{]}$
is close to the true residual
$\mathcal{R}^{k}\big{[}\widetilde{\mathbf{x}}^{k},\widetilde{\mathbf{x}}^{k-1},\ldots,\widetilde{\mathbf{x}}^{k-s}\big{]}$
such that the estimated output error is accurate.

###  3.1 Defect-corrected FOM and ROM

In our proposed approach, we first add a closure term to the FOM resulting
from the user-imposed time integration scheme. This closure term is derived
based on the snapshots of the true solution obtained using the ODE solver
library. More precisely, suppose we have the snapshots of the solution eq. 3
at any given parameter $\bm{\mu}$ obtained from an ODE solver or some legacy
codes:

$\displaystyle\mathbf{X}=\bigg{[}\mathbf{x}^{0}\,\,\mathbf{x}^{1}\,\,\cdots\,\,\mathbf{x}^{K}\bigg{]}\in\mathbb{R}^{N\times
N_{t}}.$ (14)

For purpose of illustration, we explain the details of the new method by
considering a first-order IMEX scheme as the user-imposed time integration
scheme. The FOM resulting from this is exactly eq. 4.

Since the user-imposed time integration scheme differs from the one used to
generate the snapshots in $\mathbf{X}$, we have a defect or a mismatch when we
insert $\mathbf{x}^{k}$ in eq. 14 into the first-order IMEX scheme. This reads

$\displaystyle\mathbf{d}^{k}:=\mathbf{E}_{\text{im}}\mathbf{x}^{k}-\Big{(}\mathbf{A}_{\text{im}}\mathbf{x}^{k-1}+\delta
t\big{(}\mathbf{f}(\mathbf{x}^{k-1})+\mathbf{B}\mathbf{u}^{k}\big{)}\Big{)}.$
(15)

######  Remark 2.

We note that this is precisely the local truncation error (LTE) [11] of the
first-order IMEX method if $\mathbf{x}^{k}$ is assumed to be the true solution
to eq. 1. This assumption is not strong as the tolerances (absolute, relative)
of the ODE solver can be set judiciously to achieve this.

We seek to modify the time-discrete FOM eq. 4 such that its solution recovers
the solution of eq. 1 computed by an ODE solver from a library. To this end,
consider the following corrected FOM (C-FOM for short) obtained by adding the
defect vector $\mathbf{d}^{k}$ as a closure term:

$\displaystyle\mathbf{E}_{\text{im}}\mathbf{x}_{\text{im,c}}^{k}$
$\displaystyle=\mathbf{A}_{\text{im}}\mathbf{x}_{\text{im,c}}^{k-1}+\delta
t\big{(}\mathbf{f}(\mathbf{x}_{\text{im,c}}^{k-1})+\mathbf{B}\mathbf{u}^{k}\big{)}+\mathbf{d}^{k},$
(16a) $\displaystyle\mathbf{y}_{\text{im,c}}^{k}$
$\displaystyle=\mathbf{C}\mathbf{x}_{\text{im,c}}^{k}.$ (16b)

In eq. 16, $\mathbf{x}_{\text{im,c}}^{k}$ is the solution obtained after
introducing the closure term and, as such, it differs from the solution
$\mathbf{x}_{\text{im}}^{k}$ to eq. 4. In fact, without adding
$\mathbf{d}^{k}$ as a closure term the local truncation error of the first-
order IMEX method is $\mathcal{O}((\delta t)^{2})$. However, the local
truncation error of the solution from the corrected IMEX method eq. 16 is of
the same order as that resulting from the ODE solver. It can be seen from
Figure 3 that the solution to the heat equation obtained using eq. 16 is
identical to the one obtained using an ODE solver.

We emphasize that the defect vector $\mathbf{d}^{k}$ represents one choice for
the closure term. In general, any other form of the closure term may be used.
Since we use the defect vector as our choice for the closure term to recover
the true residual, we use the two terms interchangeably.

Figure 3: Solution to the heat equation eq. 11; Left: solution obtained using
the ODE solver; Middle: solution using C-FOM eq. 16; Right: pointwise errors
between the solutions from the solver and C-FOM.

######  Remark 3.

For any general time integration scheme the defect vector $\mathbf{d}^{k}$ can
be shown to have the following equivalence:

$\displaystyle\mathbf{d}^{k}=-\mathcal{R}^{k}\big{[}\mathbf{x}^{k},\mathbf{x}^{k-1},\ldots,\mathbf{x}^{k-s}\big{]}.$
(17)

A corrected ROM (C-ROM) corresponding to the C-FOM can be defined by
projecting the defect vector $\mathbf{d}^{k}$ to the reduced space. The
reduced defect vector is defined as
$\widehat{\mathbf{d}}^{k}:=\mathbf{V}^{\textsf{T}}\mathbf{d}^{k}$ and the
C-ROM is:

$\displaystyle\widehat{\mathbf{E}}_{\text{im}}\widehat{\mathbf{x}}_{\text{im,c}}^{k}$
$\displaystyle=\widehat{\mathbf{A}}_{\text{im}}\widehat{\mathbf{x}}_{\text{im,c}}^{k-1}+\delta
t\big{(}\widehat{\mathbf{f}}(\widetilde{\mathbf{x}}_{\text{im,c}}^{k-1})+\widehat{\mathbf{B}}\mathbf{u}^{k}\big{)}+\widehat{\mathbf{d}}^{k},$
(18a) $\displaystyle\widehat{\mathbf{y}}_{\text{im,c}}^{k}$
$\displaystyle=\widehat{\mathbf{C}}\widehat{\mathbf{x}}_{\text{im,c}}^{k}.$
(18b)

###  3.2 An error estimator using the C-ROM

Next, we make use of the C-ROM to derive a new error estimator
$\overline{\Delta}^{k}(\bm{\mu})$ that accurately estimates the true error
$\|\mathbf{y}^{k}(\bm{\mu})-\widehat{\mathbf{y}}^{k}(\bm{\mu})\|$, where
$\mathbf{y}^{k}(\bm{\mu})$ and $\widehat{\mathbf{y}}^{k}(\bm{\mu})$ are the
output of the FOM eq. 1 and the ROM eq. 2 at time $t^{k}$, respectively. The
FOM and and the ROM can be solved using any ODE solver. As mentioned, to
derive the residual correctly, the same ODE solver must be applied to both the
FOM and ROM simulations. Our proposed error estimator makes use of a dual
system. We begin by deriving the dual system for eq. 16.

####  3.2.1 Dual system

We derive the dual system corresponding to the C-FOM eq. 16 using the method
of Lagrange multipliers. The Lagrangian can be formulated as

$\displaystyle\mathcal{L}:=\mathbf{C}\mathbf{x}_{\text{im,c}}^{k}+\big{(}\Lambda^{k}\big{)}^{\textsf{T}}\bigg{(}\mathbf{E}_{\text{im}}\mathbf{x}_{\text{im,c}}^{k}-\mathbf{A}_{\text{im}}\mathbf{x}_{\text{im,c}}^{k-1}-\delta
t\mathbf{f}(\mathbf{x}_{\text{im,c}}^{k-1})-\delta
t\mathbf{B}\mathbf{u}^{k}-\mathbf{d}^{k}\bigg{)},$ (19)

with $\Lambda\in\mathbb{R}^{N}$ being the vector of Lagrange coefficients. The
dual system can be obtained by setting
$\frac{\partial\mathcal{L}}{\partial\mathbf{x}_{\text{im,c}}^{k}}\equiv\mathbf{0}$.
This yields the system

$\displaystyle\mathbf{E}_{\text{du}}\mathbf{x}_{\text{du}}=\mathbf{C}_{\text{du}}$
(20)

with $\mathbf{E}_{\text{du}}:=\mathbf{E}_{\text{im}}^{\textsf{T}}$ and
$\mathbf{C}_{\text{du}}:=-\mathbf{C}^{\textsf{T}}$. Note that the Lagrange
multipliers $\Lambda$ are the dual state variables; for better clarity we
denote them by $\mathbf{x}_{\text{du}}$. Note that the defect vector
$\mathbf{d}^{k}$ is treated as a function of time and parameter $\bm{\mu}$.
Therefore, it does not depend on the solution $\mathbf{x}_{\text{im,c}}^{k}$
of the C-ROM.

We further define the dual ROM as

$\displaystyle\widehat{\mathbf{E}}_{\text{du}}\widehat{\mathbf{x}}_{\text{du}}=\widehat{\mathbf{C}}_{\text{du}}.$
(21)

obtained by making the ansatz
$\mathbf{x}_{\text{du}}\approx\widetilde{\mathbf{x}}_{\text{du}}=\mathbf{V}_{\text{du}}\widehat{\mathbf{x}}_{\text{du}}$.
Here, $\mathbf{V}_{\text{du}}$ is the projection matrix corresponding to the
dual system and
$\widehat{\mathbf{E}}_{\text{du}}:=\mathbf{V}_{\text{du}}^{\textsf{T}}\mathbf{E}_{\text{du}}\mathbf{V}_{\text{du}}$,
$\widehat{\mathbf{C}}_{\text{du}}:=\mathbf{V}_{\text{du}}^{\textsf{T}}\mathbf{C}_{\text{du}}$.

The residuals corresponding to the ROM eq. 18 and the dual ROM eq. 21 are,
respectively,

$\displaystyle\textbf{r}_{\text{im,c}}^{k}:=\mathbf{A}_{\text{im}}\widetilde{\mathbf{x}}_{\text{im,c}}^{k-1}+\delta
t\big{(}\mathbf{f}(\widetilde{\mathbf{x}}_{\text{im,c}}^{k-1})+\mathbf{B}\mathbf{u}^{k}\big{)}+\mathbf{d}^{k}-\mathbf{E}_{\text{im}}\widetilde{\mathbf{x}}_{\text{im,c}}^{k}$
(22)

and

$\displaystyle\textbf{r}_{\text{du}}:=\mathbf{C}_{\text{du}}-\mathbf{E}_{\text{du}}\widetilde{\mathbf{x}}_{\text{du}}.$
(23)

Following the approach in [57, 18] we define an auxiliary residual
$\breve{\textbf{r}}_{\text{im,c}}$ as

$\displaystyle\breve{\textbf{r}}_{\text{im,c}}^{k}:$
$\displaystyle=\mathbf{A}_{\text{im}}\mathbf{x}_{\text{im,c}}^{k-1}+\delta
t\big{(}\mathbf{f}(\mathbf{x}_{\text{im,c}}^{k-1})+\mathbf{B}\mathbf{u}^{k}\big{)}+\mathbf{d}^{k}-\mathbf{E}_{\text{im}}\widetilde{\mathbf{x}}_{\text{im,c}}^{k}.$
(24)

The auxiliary residual defined above will be required in the derivation of the
data-enhanced output error estimator.

####  3.2.2 Modified output term

In anticipation of the new error estimator we propose next, we introduce a
modified output variable $\overline{\mathbf{y}}_{\text{im,c}}^{k}$ defined as

$\displaystyle\overline{\mathbf{y}}_{\text{im,c}}^{k}:=\widehat{\mathbf{y}}_{\text{im,c}}^{k}-\widetilde{\mathbf{x}}_{\text{du}}^{\textsf{T}}\textbf{r}_{\text{im,c}}^{k}.$
(25)

Adding a correction term to the output in the form of a dual-weighted residual
is an established practice [48, 28, 18] and serves to derive a tighter
estimate of the true error.

####  3.2.3 Data-enhanced error estimation

Denoting the ODE solver applied to solve the FOM eq. 1 and the ROM eq. 2 as
solver, the norm of the true error we desire viz.,
$\|\mathbf{y}^{k}-\widehat{\mathbf{y}}^{k}\|$ can be written as:

$\displaystyle\|\mathbf{y}^{k}-\widehat{\mathbf{y}}^{k}\|$
$\displaystyle=\|\mathbf{y}^{k}-\overline{\mathbf{y}}_{\text{im,c}}^{k}+\overline{\mathbf{y}}_{\text{im,c}}^{k}-\widehat{\mathbf{y}}^{k}\|$
(26)
$\displaystyle\leq\|\mathbf{y}^{k}-\overline{\mathbf{y}}_{\text{im,c}}^{k}\|+\|\overline{\mathbf{y}}_{\text{im,c}}^{k}-\widehat{\mathbf{y}}^{k}\|.$
(27)

The following theorem bounds the first summand in eq. 27:

######  Theorem 3.1 (A posteriori error bound for the corrected ROM).

Given the FOM in eq. 1, the C-FOM in eq. 16 and the C-ROM eq. 18, assuming
that _$\mathbf{E}_{\text{im}}$_ is non-singular for all
$\bm{\mu}\in\mathcal{P}$, we have the following error bound for the modified
output vector in eq. 25:

$\displaystyle\|\mathbf{y}^{k}-\overline{\mathbf{y}}_{\textnormal{im,c}}^{k}\|\leq\|\mathbf{E}_{\textnormal{im}}^{-1}\|\,\|\mathbf{r}_{\textnormal{du}}\|\,\|\breve{\mathbf{r}}_{\textnormal{im,c}}^{k}\|+\|\widetilde{\mathbf{x}}_{\textnormal{du}}\|\,\|\mathbf{r}_{\textnormal{im,c}}^{k}-\breve{\mathbf{r}}_{\textnormal{im,c}}^{k}\|.$
(28)

###### Proof.

See Appendix B. ∎

Although the bound above is rigorous, it is not computable owing to the
quantity $\breve{\textbf{r}}^{k}$. Recall from eq. 24 that its computation
requires that the FOM solution $\mathbf{x}_{\text{im,c}}^{k}$ is available for
any parameter $\bm{\mu}$, which is not the case. To derive a computable error
estimator, we make use of the arguments used in [18] to get the following
error indicator:

$\displaystyle\|\mathbf{y}^{k}-\overline{\mathbf{y}}_{\textnormal{im,c}}\|\lessapprox\bigg{(}\overline{\rho}\,\|\mathbf{E}_{\textnormal{im}}^{-1}\|\,\|\mathbf{r}_{\textnormal{du}}\|+|1-\overline{\rho}|\,\|\widetilde{\mathbf{x}}_{\text{du}}\|\bigg{)}\,\|\textbf{r}_{\text{im,c}}^{k}\|.$
(29)

The quantity $\overline{\rho}$ is a measure for how close the residual
$\textbf{r}_{\text{im,c}}^{k}$ is to the auxiliary residual
$\breve{\textbf{r}}_{\text{im,c}}^{k}$ and is defined as

$\displaystyle\overline{\rho}=\frac{1}{K}\sum\limits_{k=1}^{K}\rho^{k},\qquad\rho^{k}=\frac{\|\breve{\textbf{r}}_{\text{im,c}}^{k}(\bm{\mu}^{*})\|}{\|\textbf{r}_{\text{im,c}}^{k}(\bm{\mu}^{*})\|}.$
(30)

It is evaluated only at the greedy parameter $\bm{\mu}^{*}$ for which
snapshots of the true solution $\mathbf{x}^{k}$ are available. For additional
details, we refer to [18]. Substituting eq. 29 into eq. 27 results in

$\displaystyle\|\mathbf{y}^{k}-\widehat{\mathbf{y}}^{k}\|$
$\displaystyle\lessapprox\bigg{(}\overline{\rho}\,\|\mathbf{E}_{\textnormal{im}}^{-1}\|\,\|\mathbf{r}_{\textnormal{du}}\|+|1-\overline{\rho}|\,\|\widetilde{\mathbf{x}}_{\text{du}}\|\bigg{)}\,\|\textbf{r}_{\text{im,c}}^{k}\|+\|\overline{\mathbf{y}}_{\text{im,c}}^{k}-\widehat{\mathbf{y}}^{k}\|=:\overline{\Delta}^{k}_{a}(\bm{\mu}).$
(31)

The second quantity in the above inequality is the norm of the error between
the output resulting from the user-imposed C-ROM and that obtained by solving
the ROM eq. 2 with solver. It can be obtained cheaply as only two ROMs with
small sizes need to be solved. Typically, our numerical experiments show that
this quantity is very small and lesser than the magnitude of the first
quantity. Therefore, we can safely neglect it so that only the first quantity
can act as an alternative form of the proposed data-enhanced error estimator,
i.e.,

$\displaystyle\|\mathbf{y}^{k}(\bm{\mu})-\widehat{\mathbf{y}}^{k}(\bm{\mu})\|$
$\displaystyle\lessapprox\|\mathbf{y}^{k}(\bm{\mu})-\overline{\mathbf{y}}_{\text{im,c}}^{k}(\bm{\mu})\|$
(32)
$\displaystyle\lessapprox\bigg{(}\overline{\rho}\,\|\mathbf{E}_{\textnormal{im}}^{-1}\|\,\|\mathbf{r}_{\textnormal{du}}(\bm{\mu})\|+|1-\overline{\rho}|\,\|\widetilde{\mathbf{x}}_{\text{du}}(\bm{\mu})\|\bigg{)}\,\|\textbf{r}_{\text{im,c}}^{k}(\bm{\mu})\|=:\overline{\Delta}^{k}_{b}(\bm{\mu}).$
(33)

####  3.2.4 Error estimation in presence of hyperreduction

For nonlinear systems, the efficient computation of the ROM is impeded by the
presence of the nonlinear function
$\mathbf{f}(\widetilde{\mathbf{x}}_{\text{im,c}}^{k})$ whose evaluation scales
as the dimension $N$ of the FOM. To tackle this, the DEIM approach [17] is
used in this work. Using DEIM, the residual expression in eq. 22 gets modified
as below:

$\displaystyle\textbf{r}_{\text{im,c}}^{k}$
$\displaystyle=\mathbf{A}_{\text{im}}\widetilde{\mathbf{x}}_{\text{im,c}}^{k-1}+\delta
t\big{(}\mathbf{f}(\widetilde{\mathbf{x}}_{\text{im,c}}^{k-1})+\mathscr{I}[\mathbf{f}(\widetilde{\mathbf{x}}_{\text{im,c}}^{k-1})]-\mathscr{I}[\mathbf{f}(\widetilde{\mathbf{x}}_{\text{im,c}}^{k-1})]+\mathbf{B}\mathbf{u}^{k}\big{)}+\mathbf{d}^{k}-\mathbf{E}_{\text{im}}\widetilde{\mathbf{x}}_{\text{im,c}}^{k},$
$\displaystyle=\underbrace{\bigg{(}\mathbf{A}_{\text{im}}\widetilde{\mathbf{x}}_{\text{im,c}}^{k-1}+\delta
t\big{(}\mathscr{I}[\mathbf{f}(\widetilde{\mathbf{x}}_{\text{im,c}}^{k-1})]+\mathbf{B}\mathbf{u}^{k}\big{)}+\mathbf{d}^{k}-\mathbf{E}_{\text{im}}\widetilde{\mathbf{x}}_{\text{im,c}}^{k}\bigg{)}}_{\textbf{r}_{\text{im,c},\mathscr{I}}^{k}}+\,\underbrace{\delta
t\bigg{(}\mathbf{f}(\widetilde{\mathbf{x}}_{\text{im,c}}^{k-1})-\mathscr{I}[\mathbf{f}(\widetilde{\mathbf{x}}_{\text{im,c}}^{k-1})]\bigg{)}}_{\mathbf{e}_{H}^{k}}.$
(34)

In the equation above, $\mathscr{I}[\cdot]$ denotes the hyperreduction
operator. Further, $\mathbf{e}_{H}^{k}$ refers to the error introduced by
hyperreduction at the $k$-th time step. For a detailed discussion on the
computational aspects, we refer to the work [18], where the simultaneous
adaptive construction of the RBM and the DEIM bases vectors is also discussed.
This adaptive bases construction approach is also implemented in our numerical
experiments.

We make use of the new data-enhanced error estimators
($\overline{\Delta}^{k}_{a}(\bm{\mu})$ or
$\overline{\Delta}^{k}_{b}(\bm{\mu})$) to choose the next parameter within the
greedy algorithm. For each parameter $\bm{\mu}$, we determine the average
estimated error over time as:

$\displaystyle\overline{\Delta}_{z}(\bm{\mu})=\dfrac{1}{N_{t}}\sum\limits_{k=0}^{N_{t}}\overline{\Delta}^{k}_{z}(\bm{\mu}),$
(35)

where $z$ stands for $a$ or $b$. In the numerical results, we use
$\overline{\Delta}^{k}_{b}(\bm{\mu})$ in eq. 35 in Line 11 of Algorithm 2 so
that during the greedy iterations, the ROM eq. 2 does not have to be
repeatedly solved with solver in order to evaluate the second quantity in eq.
31.

######  Remark 4.

While we have illustrated the data-enhanced error estimator using a first-
order IMEX time integration scheme, any consistent time integration scheme can
be used, including higher order ones. We will also demonstrate the use of a
second-order IMEX scheme in the numerical examples. Of course, one should be
cautious of the fact that a higher order time integration scheme comes with a
larger computational cost.

######  Remark 5.

When the defect vector $\mathbf{d}^{k}$ present in the residual (see eq. 22)
is known exactly, i.e., the C-ROM uses the same time integration method used
for the FOM, we recover the _a posteriori_ output error estimator proposed in
[18]. This shows that our newly proposed error estimator is consistent with
the case where the time integration method used is known.

######  Remark 6.

Note that the C-ROM is not the finally derived ROM, but is only used within
the greedy algorithm to derive the error estimator
$\overline{\Delta}^{k}_{z}(\bm{\mu})$ (see Steps 9-10 in Algorithm 2) which
estimates the error between the solution of the FOM eq. 1 and that of the ROM
eq. 2, both being computed using any solver. The error estimator can also be
used in the online stage where the user may wish to use any preferred solver
to compute the ROM solutions. Moreover, given a good approximation of the
defect vector, the solution of the C-ROM and the solution
$\widehat{\mathbf{x}}^{k}$ from the ODE solver are nearly the same. Therefore,
$\overline{\Delta}^{k}_{b}(\bm{\mu})$ is almost as accurate as
$\overline{\Delta}^{k}_{a}(\bm{\mu})$.

##  4 Computational aspects

While we have derived an accurate data-enhanced _a posteriori_ output error
estimator, there remain a few computational challenges. In this section we
highlight these challenges and propose efficient solutions to address them.

In eq. 33 the residual term $\textbf{r}_{\text{im,c}}^{k}$ involves the defect
vector $\mathbf{d}^{k}$ (see eq. 22). However, determining this term involves
knowing the true solution to eq. 1 obtained from the ODE solver at any
parameter $\bm{\mu}$ (see eq. 15). Without a cheap method to approximate
$\mathbf{d}^{k}$ for every parameter $\bm{\mu}$, the error estimator is not
efficient. To alleviate this, we make several observations about the function
$\mathbf{d}(t,\bm{\mu})$, which will lead to its efficient approximation.
These observations relate to

  1. 1.

a certain low-rank structure (over space) that $\mathbf{d}(t,\bm{\mu})$
possesses and

  2. 2.

a smoothness over parameter variations that it inherits from the underlying
parametric PDE.

###  4.1 Low-rank structure of the defect

For a large class of problems, there exist a low-dimensional subspace onto
which the solution snapshots can be projected, without incurring a large
error. Indeed, this forms the underlying motivation for performing model order
reduction using POD and other methods. Leveraging this fact, we assume that
the defect trajectory at a given parameter $\bm{\mu}$,
$\mathbf{d}^{k}(\bm{\mu}),\forall k\in\\{0,1,\ldots,K\\}$ can also be
efficiently approximated by a suitable low-dimensional subspace.

Figure 4: Defect/local truncation error of the heat equation. Left: space time
variation of the defect at $\mu=0.06$; Right: singular values of the matrix
$\mathbf{D}:=\\{\mathbf{d}^{k}\\}_{k=0}^{K}$. The singular values exhibit an
exponential decay illustrating the existence of a low-dimensional subspace.

Figure 4 plots the values of the defect vector $\mathbf{d}^{k}$ in eq. 15 for
the heat equation. It is evident that the defect trajectory has a certain
regularity over the spatial domain (see left figure); evidence for the low-
rank structure is also seen through the SVD performed (see right figure) on
the snapshot matrix $\mathbf{D}\in\mathbb{R}^{N\times N_{t}}$ whose columns
consist of snapshots of $\mathbf{d}^{k}$ at different values of the time
$t^{k},\forall k\in\\{0,\ldots,K\\}$. We see an exponential decay of the
relative quantity $\sigma_{i}/\sigma_{1}$, where $\sigma_{i}$ is the $i$-th
singular value. A similar exercise will be repeated in Section 5 for the three
numerical examples we consider. In each of the cases, it will be seen that the
defect trajectory is reducible spatially. We will use this fact to compute
$\mathbf{d}^{k}(\bm{\mu})$ efficiently.

The concept of Kolmogorov $n$-width is used in the RBM literature [46, 9, 44]
to quantify the approximability of the solution manifold $\mathcal{M}$
corresponding to a parametrized system of equations, with a linear subspace of
dimension $n$ denoted as $\mathcal{V}^{n}$. Consider the solution manifold
$\mathcal{M}$ for the FOM in eq. 1 defined as

$\displaystyle\mathcal{M}=\\{\mathbf{x}(t,\bm{\mu})\,:\,(t,\bm{\mu})\in\mathcal{T}\times\mathcal{P}\\}\subset\mathbb{R}^{N}.$
(36)

The Kolmogorov $n$-width of $\mathcal{M}$ using $\mathcal{V}^{n}$ can be
defined as

$\displaystyle
d_{n}(\mathcal{M}):=\inf\limits_{\text{dim}(\mathcal{V}^{n})\,=\,n}\,\,\sup\limits_{\mathbf{x}\in\mathcal{M}}\,\inf\limits_{\widetilde{\mathbf{x}}\in\mathcal{V}^{n}}\|\mathbf{x}-\widetilde{\mathbf{x}}\|.$
(37)

For the parametrized defect function $\mathbf{d}(t,\bm{\mu})$, we define the
following manifold:

$\displaystyle\mathcal{M}_{d}=\\{\mathbf{d}(t,\bm{\mu})\,:\,(t,\bm{\mu})\in\mathcal{T}\times\mathcal{P}\\}\subset\mathbb{R}^{N}.$
(38)

From our numerical examples we have observed that $\mathcal{M}_{d}$ can be
well-approximated by low-dimensional subspaces when the original parametric
problem has fast Kolmogorov $n$-width decay. This implicates that the
Kolmogorov $n$-width of $\mathcal{M}_{d}$ may inherit the behaviour of
$d_{n}(\mathcal{M})$. At present we can not strictly prove this, but from the
definition of the defect vector in eq. 15, the defect vector
$\mathbf{d}(t,\bm{\mu})$ can be seen as being in the image of the operator
$\mathcal{D}:\mathbb{R}^{n}\times\mathcal{P}\mapsto\mathbb{R}^{n}$ and
$\mathcal{D}[\mathbf{x}(t,\bm{\mu})]=\mathbf{E}_{\text{im}}\mathbf{x}(t,\bm{\mu})-\big{(}\mathbf{A}_{\text{im}}\mathbf{x}(t-\delta
t,\bm{\mu})+\delta t(\mathbf{f}(\mathbf{x}(t-\delta
t,\bm{\mu}))+\mathbf{B}\mathbf{u}(t))\big{)}$. With this observation, the
inherit property of the Kolmogorov $n$-width of $\mathcal{M}_{d}$ might be
proved based on Theorem 4.1 in [22]. We leave this as our future work.

###  4.2 Strategies to approximate the defect

As discussed in the previous section, the defect vector is assumed to possess
a certain low-rank structure in space and it inherits the smoothness of the
solution $\mathbf{x}(t,\bm{\mu})$ over parameter variations, from the
underlying parametric PDE. We use these two observations to make the
approximation of $\mathbf{d}^{k}(\bm{\mu})$ computationally efficient, such
that the output error estimator defined in eq. 33 can be used in the POD-
Greedy algorithm. To this end, we adopt a two-stage approximation strategy.

Starting from the observation about the low-rank structure of the defect, we
can approximate the defect vector at a given time and parameter using the
basis expansion

$\displaystyle\mathbf{d}(t,\bm{\mu})\approx\widetilde{\mathbf{d}}(t,\bm{\mu})=\sum_{i=1}^{n_{d}}\mathbf{v}_{d,i}\,\widehat{d}_{i}(t,\bm{\mu})$
(39)

where $\mathbf{v}_{d,i}\in\mathbb{R}^{N}$ are the expansion bases and
$\widehat{\mathbf{d}}(t,\bm{\mu}):=[\widehat{d}_{1}(t,\bm{\mu}),\ldots,\widehat{d}_{n_{d}}(t,\bm{\mu})]\in\mathbb{R}^{n_{d}}$
is the vector of expansion coefficients. We denote with
$\mathbf{V}_{d}:=\big{[}\mathbf{v}_{d,1},\mathbf{v}_{d,2},\ldots,\mathbf{v}_{d,n_{d}}\big{]}\in\mathbb{R}^{N\times
n_{d}}$ the basis matrix. If the observation regarding the rapid decay of the
singular values holds, then $n_{d}$ will be small. Given such a basis
expansion for the defect vector, we can approximate it for any given parameter
$\bm{\mu}$ and a time instance $t$ if $\widehat{\mathbf{d}}(t,\bm{\mu})$ can
be evaluated cheaply. Our two-stage approach involves:

  * •

identifying a suitable basis matrix $\mathbf{V}_{d}$ using a POD/SVD-based
approach and

  * •

learning the map $(t,\bm{\mu})\mapsto\widehat{\mathbf{d}}(t,\bm{\mu})$ for
which we propose two different approaches: one based on radial basis function
interpolation and the other using a feed-forward neural network.

####  4.2.1 SVD-based spatial reduction

In the first stage of approximation, we collect snapshots of the defect vector
$\mathbf{d}(t,\bm{\mu}),\forall
t\in\mathcal{T}_{d}:=\\{t^{1},t^{2},\ldots,t^{N_{t}}\\}$ and
$\bm{\mu}\in\Xi_{\text{defect}}$ where $\Xi_{\text{defect}}$ is a set
containing $d_{s}$ parameter samples, with $d_{s}$ typically small. Doing so
involves solving eq. 1 with solver to obtain the FOM solution snapshots.
Following this, the solution snapshots can be used to obtain the defect vector
from e.g., eq. 15111Note that the defect is obtained here for a particular
user-imposed time integration scheme.. We denote by
$\mathfrak{D}\in\mathbb{R}^{N\times N_{t}\times d_{s}}$ the third-order tensor
arranged such that each of its frontal slices corresponds to the matrix
$\mathbf{D}(\bm{\mu}):=\big{[}\mathbf{d}(t^{1},\bm{\mu}),\mathbf{d}(t^{2},\bm{\mu}),\ldots,\mathbf{d}(t^{N_{t}},\bm{\mu})\big{]}\in\mathbb{R}^{N\times
N_{t}}$ and $\bm{\mu}\in\Xi_{\text{defect}}$. We refer to the $i$-th frontal
slice as $\mathfrak{D}^{(i)}$. Next, we apply a two-step SVD reduction [53] to
$\mathfrak{D}$ which will result in
$\widehat{\mathfrak{D}}\in\mathbb{R}^{n_{d}\times N_{t}\times d_{s}}$ as
follows:

  * •

In step 1, we perform SVDs for each frontal slice separately, i.e., we perform
a SVD for each $\mathfrak{D}^{(i)}$ with $i\in\\{1,2,\ldots,d_{s}\\}$. For
every such SVD, we make use of a fixed tolerance $\texttt{tol}_{\text{SVD},t}$
to truncate the singular values and collect the first $\ell_{i}$ left singular
vectors in the matrix $\mathbf{U}_{d,i}\in\mathbb{R}^{N\times\ell_{i}}$.

  * •

In step 2, we form a matrix $\mathbf{R}$ defined as

$\mathbf{R}:=\bigg{[}\mathbf{U}_{d,1}\,\big{\lvert}\,\mathbf{U}_{d,2}\,\big{\lvert}\,\cdots\,\big{\lvert}\,\mathbf{U}_{d,d_{s}}\bigg{]}\in\mathbb{R}^{N\times(\ell_{1}+\cdots+\ell_{d_{s}})}$

whose columns consist of the truncated left singular vectors for each
parameter obtained in step 1. We then perform the SVD of $\mathbf{R}$, using a
tolerance $\texttt{tol}_{\text{SVD},\bm{\mu}}$ to obtain the projection matrix
$\mathbf{V}_{d}\in\mathbb{R}^{N\times n_{d}}$ as the first $n_{d}$ columns of
the left singular vectors of $\mathbf{R}$.

Finally, the reduced tensor $\widehat{\mathfrak{D}}$ can be obtained via a
mode-1 tensor-matrix product as

$\widehat{\mathfrak{D}}=\mathfrak{D}\times_{1}\mathbf{V}_{d}^{\textsf{T}}.$

We note that each mode-1 fiber in $\widehat{\mathfrak{D}}$ corresponds to the
reduced defect vector
$\widehat{\mathbf{d}}(t^{k},\bm{\mu})\in\mathbb{R}^{n_{d}}$.

Thus far, we have reduced the dimension of the first mode of the tensor
$\mathfrak{D}$ from $N$ to $n_{d}$, with $n_{d}\ll N$. Next, we detail the two
approaches used to approximate $\widehat{\mathbf{d}}(t,\bm{\mu})$.

######  Remark 7.

The motivation for using the two-step SVD approach is to reduce the overall
computational costs. While we have not pursued it in this work, another valid
approach to reduce the computational cost would be a randomized-SVD.

####  4.2.2 Interpolation using radial basis functions

Radial basis functions (RBFs) are a popular class of kernel methods which are
used in scattered-data approximation [10, 54]. We use RBFs to learn an
approximation $\widehat{d}_{\text{RBF}}(\bm{\mu})$ to each entry of the
reduced defect vector at each time instance, such that for any given $(j,k)$
pair, $\widehat{d}_{\text{RBF}}(\bm{\mu})$ interpolates
$\widehat{d}_{j}(t^{k},\bm{\mu})$, the $j$-th entry of the reduced defect
$\widehat{\mathbf{d}}(t^{k},\bm{\mu})$ at parameters
$\bm{\mu}\in\Xi_{\text{defect}}$.

The RBF approximation reads

$\displaystyle\widehat{d}_{\text{RBF},j}^{k}(\bm{\mu})=\sum\limits_{i=1}^{d_{s}}w_{i}\Phi(\|\bm{\mu}-\bm{\mu}_{i}\|)$
(40)

with $\\{w_{i}\\}_{i=1}^{d_{s}}$ denoting the weights and $\Phi(\cdot)$ being
the radial basis functions. The weights are obtained by imposing the
interpolation condition
$\widehat{d}_{\text{RBF},j}^{k}(\bm{\mu}_{i})=\widehat{d}_{j}(t^{k},\bm{\mu}_{i}),\,i=1,2,\dots,d_{s}$,
$\bm{\mu}_{i}\in\Xi_{\text{defect}}$. This leads to the following system of
linear equations

$\displaystyle\begin{bmatrix}\Phi(\left\lVert\bm{\mu}_{1}-\bm{\mu}_{1}\right\rVert)&\Phi(\left\lVert\bm{\mu}_{1}-\bm{\mu}_{2}\right\rVert)&\cdots&\Phi(\left\lVert\bm{\mu}_{1}-\bm{\mu}_{d_{s}}\right\rVert)\\\
\Phi(\left\lVert\bm{\mu}_{2}-\bm{\mu}_{1}\right\rVert)&\Phi(\left\lVert\bm{\mu}_{2}-\bm{\mu}_{2}\right\rVert)&\cdots&\Phi(\left\lVert\bm{\mu}_{2}-\bm{\mu}_{d_{s}}\right\rVert)\\\
\vdots&\vdots&\ddots&\vdots\\\
\Phi(\left\lVert\bm{\mu}_{d_{s}}-\bm{\mu}_{1}\right\rVert)&\Phi(\left\lVert\bm{\mu}_{d_{s}}-\bm{\mu}_{2}\right\rVert)&\cdots&\Phi(\left\lVert\bm{\mu}_{d_{s}}-\bm{\mu}_{d_{s}}\right\rVert)\end{bmatrix}\begin{bmatrix}w_{1}\\\
w_{2}\\\ \vdots\\\
w_{d_{s}}\end{bmatrix}=\begin{bmatrix}\widehat{d}_{j}(t^{k},\bm{\mu}_{1})\\\
\widehat{d}_{j}(t^{k},\bm{\mu}_{2})\\\ \vdots\\\
\widehat{d}_{j}(t^{k},\bm{\mu}_{d_{s}})\end{bmatrix}.$ (41)

Based on the RBF interpolation, the defect vector $\mathbf{d}(t^{k},\bm{\mu})$
in eq. 39 can be approximated as

$\displaystyle\mathbf{d}(t^{k},\bm{\mu})\approx\widetilde{\mathbf{d}}_{\text{RBF}}(t^{k},\bm{\mu}):=\sum_{j=1}^{n_{d}}\mathbf{v}_{d,j}\,\widehat{d}^{k}_{\text{RBF},j}(\bm{\mu})$
(42)

We obtain an RBF interpolant for each time instance $k\in\\{0,1,2,\ldots,K\\}$
and each coordinate $j\in\\{1,2,\ldots,n_{d}\\}$, resulting in a total of
$(n_{d}\cdot N_{t})$ RBF interpolants. Figure 5 graphically illustrates the
approach. Theoretically, $\widetilde{\mathbf{d}}_{\text{RBF}}(t^{k},\bm{\mu})$
in eq. 42 is valid for any $\bm{\mu}\in\mathcal{P}$. Therefore, $N_{p}$ in
Figure 5 can be arbitrarily large. We denote it by $N_{p}$ as we are only
interested in obtaining the reduced defect coefficients corresponding to the
$N_{p}$ parameter samples present in the training set $\Xi$.

In our numerical results in Section 5, we denote this method of approximating
the defect vector using a SVD spatial reduction followed by an RBF
interpolation as SVD+RBF.

######  Remark 8.

In this work, we have considered separate RBF interpolants for each time step
and each generalized spatial coordinate, leading to potentially many
interpolants. While the number of interpolants scales as
$\mathcal{O}(n_{d}N_{t})$, this can be efficiently implemented in one step by
solving a linear system with the small $d_{s}\times d_{s}$ coefficient matrix
in eq. 41, and with multiple right-hand sides (totalling $n_{d}N_{t}$). The
runtime of solving such a linear system is usually much faster than separately
solving $n_{d}N_{t}$ linear systems with the same coefficient matrix.

Figure 5: Approximation of the defect as a function of the parameter
$\bm{\mu}$. The RBF interpolant learns an approximation of the defect vector
over $N_{p}$ samples, with interpolation occurring at $d_{s}$ samples. We
construct an individual RBF interpolant for each time and generalized spatial
coordinate.

####  4.2.3 Approximation using artificial neural networks

Figure 6: Approximation of the defect coefficients as a function of the inputs
$(t,\bm{\mu})$. The neural network is trained based on data available at
$d_{s}$ parameter samples. In the inference stage, the neural network learns
the approximation of the defect for all $N_{p}$ parameter samples.

The second approach we consider to approximate the expansion coefficients
$\widehat{\mathbf{d}}(t,\bm{\mu})$ for different time and parameter values is
based on artificial neural networks (ANNs). A widely used architecture to
implement ANNs are the feed-forward neural networks (FNNs). FNNs have been
shown to be efficient for both regression and classification tasks in a
variety of applications. The basic architecture of a FNN consists of three
components: an input layer, hidden layer(s) and an output layer. The core
element of any NN in general and FNNs in particular are artificial neurons.
The hidden layer(s) in a FNN consists of neurons stacked together. Each neuron
can receive inputs from a previous layer. The output corresponds to a
nonlinear function of the weighted sum of its input signals. For a detailed
overview of ANNs and FNNs, we refer to [8].

In this work, we consider the FNN to learn an approximation to the map between
the inputs $(t,\bm{\mu})$ and the output
$\widehat{\mathbf{d}}(t,\bm{\mu})\in\mathbb{R}^{n_{d}}$. That is,

$\widehat{\mathbf{d}}_{\text{NN}}:(t,\bm{\mu})\mapsto\widehat{\mathbf{d}}_{\text{NN}}(t,\bm{\mu})\approx\widehat{\mathbf{d}}(t,\bm{\mu}).$

To train the FNN, our training data consists of the dataset
$\tau_{\text{train}}:=(\bm{\Lambda}_{\text{NN}},D_{\text{NN}})$ where
$\bm{\Lambda}_{\text{NN}}:=\\{(t_{i},\bm{\mu}_{j})\\}_{i=j=1}^{N_{t},d_{s}}$
is the input set, containing time and parameter samples and
$D_{\text{NN}}:=\\{\widehat{\mathbf{d}}(t_{i},\bm{\mu}_{j})\\}_{i=j=1}^{N_{t},d_{s}}$
is the output data set which consists of the reduced defect
$\widehat{\mathbf{d}}(t_{i},\bm{\mu}_{j})\in\mathbb{R}^{n_{d}}$ at each time
and parameter sample of the input. The neural network is implemented in
PyTorch; more details regarding the number of layers used and other
hyperparameters will be provided in the numerical section. The loss function
is the mean square loss, viz.,

$L_{\text{mse}}:=\frac{1}{2}\sum_{i=1}^{d_{s}}\sum_{k=1}^{N_{t}}\bigg{\lvert}\bigg{\lvert}\widehat{\mathbf{d}}(t^{k},\bm{\mu}_{i})-\widehat{\mathbf{d}}_{\text{NN}}(t^{k},\bm{\mu}_{i})\bigg{\rvert}\bigg{\rvert}_{2}^{2}.$

Once the neural network is trained, it can infer the values of the defect
vector $\widehat{\mathbf{d}}_{\text{NN}}(t,\bm{\mu})$ at any chosen
$(t,\bm{\mu})$. The original defect vector at a given parameter $\bm{\mu}$ and
at any time instance $\mathbf{d}(t,\bm{\mu})$ can be approximated using the
FNN-based approach as

$\displaystyle\mathbf{d}(t,\bm{\mu})\approx\widetilde{\mathbf{d}}_{\text{NN}}(t,\bm{\mu}):=\mathbf{V}_{d}\,\widehat{\mathbf{d}}_{\text{NN}}(t,\bm{\mu}).$
(43)

In our numerical results in Section 5, we denote this method of approximating
the defect vector first with a SVD spatial reduction followed by an
approximation of the coefficients $\widehat{\mathbf{d}}$ with a feedforward
neural network as SVD+FNN.

###  4.3 POD-Greedy with black-box ODE solvers

We present the POD-Greedy algorithm that supports black-box ODE solvers and
incorporates the new data-enhanced _a posteriori_ error estimator as Algorithm
2. It requires some additional inputs compared to Algorithm 1. These include a
separate training set $\Xi_{\text{defect}}\subset\Xi$ for the defect
approximation and two separate tolerances for the SVD,
$\texttt{tol}_{\text{SVD},t},\texttt{tol}_{\text{SVD},\bm{\mu}}$ which are
required for the two-step SVD method to compute $\mathbf{V}_{d}$. Before
starting of the greedy algorithm in Step 4, Steps 2-3 in Algorithm 2 are
targeted towards learning the defect vector, which is added as a closure term
to get the C-ROM. In Step 2, FOM solutions of eq. 1 are obtained for the
$d_{s}$ parameter samples in the training set $\Xi_{\text{defect}}$. Using
this data, the defect vectors (see eq. 15) induced by a user-imposed time
integration scheme (denoted $\texttt{solver}_{\text{imp}}$) are collected in
the tensor $\mathfrak{D}$. Then, in Step 3, the data in $\mathfrak{D}$ are
first compressed into a low-dimensional space and the reduced defect vectors
$\widehat{\mathbf{d}}(t,\bm{\mu})$ for all $\bm{\mu}\in\Xi$ are learned via
RBF or FNN. Afterwards, the defect vector $\mathbf{d}(t,\bm{\mu})$ is
approximated via the decoded vectors
$\widetilde{\mathbf{d}}_{\text{RBF}}(t,\bm{\mu})$ or
$\widetilde{\mathbf{d}}_{\text{NN}}(t,\bm{\mu})$. The approximation
$\widetilde{\mathbf{d}}_{\text{RBF}}(t,\bm{\mu})$ or
$\widetilde{\mathbf{d}}_{\text{NN}}(t,\bm{\mu})$ can be updated (replaced) by
the true defect vector once new FOM data is available at $\bm{\mu}^{*}$ (Step
6, Step 11). Updating the approximate defect vector with the available true
defect vector at $\bm{\mu}^{*}$ at each iteration leads to considerable
improvements in the performance of Algorithm 2, as we shall demonstrate in
Section 5. Furthermore, it reduces the amount of initial training samples
needed in $\Xi_{\text{defect}}$. Typically, the user may not know, a priori,
the number of FOM samples needed to get a good approximation of the closure
term. Therefore, $\Xi_{\text{defect}}$ can be coarsely sampled to keep the
computational cost low. Once the greedy algorithm is begun, the FOM solution
snapshots at $\bm{\mu}^{*}$ chosen at each greedy iteration are readily
available. Those snapshots can be further used to compute the true defect
vector $\mathbf{d}(t^{k},\bm{\mu}^{*})$ via eq. 15. Since the snapshots are
anyway available, the only computational costs incurred are those
corresponding to the evaluation of $\mathbf{d}(t^{k},\bm{\mu}),\forall
k\in\\{0,1,\ldots,K\\}$.

Algorithm 2 POD-Greedy algorithm for ODE solver libraries

0: Training set $\Xi$, tolerance (tol), Discretized model
($\mathbf{E},\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{f},\mathbf{x}_{0}$)

0: Training set for defect approximation $\Xi_{\text{defect}}\subset\Xi$, SVD
tolerances ($\texttt{tol}_{\text{SVD},t},\texttt{tol}_{\text{SVD},\bm{\mu}}$),
user-imposed time integration method $\texttt{solver}_{\textnormal{imp}}$

0: $\mathbf{V}$.

1: Initialize $\mathbf{V}=[\,\,]$, $\epsilon=1+\texttt{tol}$, choose solver
from a library, greedy parameter $\bm{\mu}^{*}$

2: Compute FOM snapshots for all $\bm{\mu}\in\Xi_{\text{defect}}$; determine
the defect data tensor $\mathfrak{D}$ (see Section 4.2.1)

3: Perform SVD on $\mathfrak{D}$ to get $\mathbf{V}_{d}$; obtain an
approximation $\widetilde{\mathbf{d}}_{\text{RBF}}(t,\bm{\mu})$ eq. 42 or
$\widetilde{\mathbf{d}}_{\text{NN}}(t,\bm{\mu})$ eq. 43 for
$\mathbf{d}(t,\bm{\mu})$

4: while $\epsilon>\texttt{tol}$ do

5: Obtain FOM snapshots $\mathbf{X}_{\bm{\mu}^{*}}$ at $\bm{\mu}^{*}$ with
solver

6: Compute true defect $\\{\mathbf{d}(t^{k},\bm{\mu})\\}_{k=1}^{N_{t}}$ at
$\bm{\mu}^{*}$ based on $\mathbf{X}^{*}$

7: Determine $\mathbf{V}^{*}$ through an SVD of
$\overline{\mathbf{X}}:=\mathbf{X}^{*}-\mathbf{V}(\mathbf{V}^{\textsf{T}}\mathbf{X}^{*})$

8: Update $\mathbf{V}$ as
$\mathbf{V}:=\texttt{orth}\big{(}\mathbf{V},\mathbf{V}^{*}(:\,,\,1:r_{c})\big{)}$with
$\texttt{orth}\big{(}\cdot\big{)}$ denoting an orthogonalization process which
can be implemented using the modified Gram-Schmidt process, or QR algorithm

9: Obtain the C-ROM (e.g., eq. 18) corresponding to
$\texttt{solver}_{\textnormal{imp}}$ by Galerkin proj. (+ hyperreduction)

10: Solve the C-ROM to obtain the residual and compute the error estimator
$\overline{\Delta}_{a}(\bm{\mu})$ eq. 31or $\overline{\Delta}_{b}(\bm{\mu})$
eq. 33

11:
$\bm{\mu}^{*}:=\arg\max\limits_{\bm{\mu}\in\Xi}\overline{\Delta}_{z}(\bm{\mu})$,
$z=a$ or $z=b$

12: Set $\epsilon=\overline{\Delta}_{z}(\bm{\mu}^{*})$, $z=a$ or $z=b$

13: end while

####  4.3.1 Computational cost

We now analyse the additional computational cost incurred by the proposed
algorithm Algorithm 2 when compared to the standard POD-Greedy method in
Algorithm 1.

To simplify things, we define the cost of solving the nonlinear FOM eq. 1 to
be $\mathcal{C}_{\text{FOM}}:=N_{\text{newton}}N_{L}N_{t}$; the cost of the
linear solve at each Newton iteration, viz. $N_{L}$, could range between
$N^{2}-N^{3}$ depending on the particular method being implemented in the
solver. The cost of obtaining the defect vector at each time step (see eq. 15)
is denoted by $\mathcal{C}_{\mathbf{d}}$. It depends on the user-defined time-
stepping scheme. The main contributions within $\mathcal{C}_{\mathbf{d}}$ are
matrix-vector products and evaluations of the nonlinearity. Thus, it evaluates
to $\mathcal{C}_{\mathbf{d}}:=(\mathcal{O}(N^{2})+\mathcal{O}(N))$ in the
worst-case . However, the matrix vector multiplications typically involve
sparse matrices and can be done cheaply.

We denote by
$\mathcal{C}_{\text{SVD}}:=\mathcal{O}(\text{min}(p^{2}q,pq^{2}))$ the cost of
a SVD for a matrix of dimension $\mathbb{R}^{p\times q}$. We let
$\mathcal{C}_{\text{RBF}}:=\mathcal{O}(d_{s}^{3})$ be the cost of obtaining
one RBF interpolant. The cost of training a neural network is difficult to
estimate owing to its architecture and the use of specialized hardware.
Therefore, for simplicity, we define it as
$\mathcal{C}_{\text{NN}}:=\sigma\cdot(\mathcal{C}_{\text{FP}}+\mathcal{C}_{\text{BP}})$
with $\mathcal{C}_{\text{FP}}$ denoting the cost of one forward pass,
$\mathcal{C}_{\text{BP}}$ denoting the cost of one backward pass and $\sigma$
being a constant that depends on the number of layers, epochs, the batch size
and other hyperparameters.

The overall factors contributing to the proposed algorithm are listed below:

  * •

Step 2 in Algorithm 2 requires the FOM solution at $d_{s}$ parameters. This
has cost that scales as $d_{s}\cdot\mathcal{C}_{\text{FOM}}$.

  * •

In Step 2, evaluating the defect vector at each time step and for all $d_{s}$
parameter samples incurs cost that scales as
$(d_{s}N_{t})\cdot\mathcal{C}_{\mathbf{d}}$.

  * •

The cost of the two-stage SVD in Step 3 to obtain $\mathbf{V}_{d}$ is
$(d_{s}+1)\cdot\mathcal{C}_{\text{SVD}}$.

  * •

Based on the method used for approximating the map
$(t,\bm{\mu})\mapsto\widehat{\mathbf{d}}(t,\bm{\mu})$, the costs differ:

    * -

For the RBF-based approach, the RBF coefficient matrix can be factorized once
and reused for subsequent solves. The cost is thus
$\mathcal{C}_{\text{RBF}}+(n_{d}N_{t}-1)\cdot\mathcal{O}(d_{s}^{2})$

    * -

For the NN-based approach, the cost is $\mathcal{C}_{\text{NN}}$

  * •

In the inference stage where the defect vector is approximated for all
$\bm{\mu}\in\Xi$,

    * -

the RBF-based approach has a cost $(N_{p}N_{t})\cdot d_{s}$

    * -

the NN-based approach incurs a cost scaling as
$(N_{p}N_{t})\cdot\mathcal{C}_{\text{FP}}$

The cost of the matrix tensor product to obtain $\widetilde{\mathfrak{D}}$ is
$d_{s}N_{t}\cdot n_{d}N$ for both RBF and NN-based approaches

  * •

To update the approximate defect vector in Step 6, the cost of evaluating the
true defect vector at the current greedy parameter scales as
$N_{t}\cdot\mathcal{C}_{\mathbf{d}}$

The major cost in approximating the defect will be the cost of solving the FOM
at $d_{s}$ parameter samples. From our experience, the RBF-based approach
performs better than the NN-based approach. Comparisons of corresponding run
times are detailed in Section 5.

##  5 Numerical results

To demonstrate the validity of the proposed data-enhanced output error
estimation approach, we test it on three numerical examples. These are:

  1. 1.

the viscous Burgers’ equation with one parameter

  2. 2.

the FitzHugh-Nagumo equations with two parameters and

  3. 3.

the batch chromatography equations with one parameter.

All reported numerical results were performed on a desktop computer running
Ubuntu 20.04, installed with a $12$-th generation intel®core™i5 processor,
$32$GB of RAM and a NVIDIA RTX A4000 GPU with $16$GB of memory. The
simulations were carried out in the Spyder IDE with Python 3.9.12 (miniconda).
Where required, matlab® v2019b was used to run some simulations. The radial
basis interpolation was performed using the RBF package [38].

For the two greedy algorithms (Algorithms 1 and 2), we plot the maximum
estimated errors computed using $\overline{\Delta}_{b}(\bm{\mu})$ over the
training set at every iteration. We define this as:

$\varepsilon_{\text{max}}:=\max\limits_{\bm{\mu}\in\Xi}\overline{\Delta}_{b}(\bm{\mu})$

where $\overline{\Delta}_{b}(\bm{\mu})$ is defined in eq. 33. Additionally, to
illustrate the performances of the ROM over the test set $\Xi_{\text{test}}$,
we plot the mean estimated error $\overline{\Delta}_{b}(\bm{\mu})$ for every
parameter $\bm{\mu}\in\Xi_{\text{test}}$.

###  5.1 Code availability

The companion Python code to reproduce the numerical results is available at
https://doi.org/10.5281/zenodo.8169490.

###  5.2 Burgers’ equation

##### Model description

The viscous Burgers’ equation defined in the 1-D domain $\Omega:=[0,1]$ is
given by

$\displaystyle\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial z}$
$\displaystyle=\mu\frac{\partial^{2}v}{\partial^{2}z},$ (44) $\displaystyle
v(z,0)$ $\displaystyle=\sin(2\pi z),$ $\displaystyle v(0,t)$
$\displaystyle=v(1,t)=0$

with $v:=v(z,t)\in\mathbb{R}$ denoting the state variable and $z\in\Omega$ is
the spatial variable and the time variable $t\in[0,2]$. We spatially
discretize eq. 44 with the finite difference method. The mesh size is $\Delta
z=0.001$, which results in a discretized FOM of dimension $N=1000$. As the
variable parameter, we consider the viscosity $\mu\in\mathcal{P}:=[0.005,1]$.
We sample $100$ logarithmically-spaced samples from $\mathcal{P}$ and divide
the samples randomly into a training set $\Xi$ and a testing set
$\Xi_{\text{test}}$ in the ratio $80:20$. To solve the ODE, the solver in
Algorithm 2 is the scipy library odeint. To have a uniform comparison, we
compute solutions of the FOM on a uniformly spaced time step of $\delta
t=0.01$. The output variable of interest is the value of the state at the node
just before the right boundary.

##### POD-Greedy with Algorithm 1

We first apply Algorithm 1 to the Burgers’ equation with a tolerance
$\texttt{tol}=10^{-4}$. The solver used is the odeint library in scipy. Since
the exact expression of the residual is unknown, we use a first-order IMEX
method (IMEX1) to approximate the residual and estimate the output error. As
the true residual is incorrectly approximated by the user-imposed time
integration scheme, the resulting estimated error is severely overestimated.
This results in the stagnation of the greedy algorithm as seen in Figure 7a.
For comparison, we also plot in Figure 7b the convergence of the greedy
algorithm when the exact residual is known. To remedy this, we next apply the
proposed Algorithm 2.

(a)  (b)

Figure 7: Burgers’ equation, Algorithm 1: (a) error (estimator) decay when
using a first-order IMEX method without any closure term; (b) error
(estimator) decay when the true residual is known.

##### POD-Greedy with Algorithm 2 and SVD+RBF closure approximation

First, we collect $d_{s}=16$ uniformly-spaced parameter samples from the
training set to construct $\Xi_{\text{defect}}$ and obtain the corresponding
defect vector in Step 2. The SVD tolerances
$\texttt{tol}_{\text{SVD},t},\texttt{tol}_{\text{SVD},\bm{\mu}}$ are both set
to $10^{-4}$ so that $n_{d}=47$. The user-defined solver
$\texttt{solver}_{\textnormal{imp}}$ is IMEX1. The resulting convergence of
the greedy algorithm using the RBF-based approximation of the defect vector is
shown in Figure 8a. As shown, the maximum estimated error converges
exponentially to the desired tolerance. The dimension of the ROM obtained is
$n=7$. To demonstrate the performance of the ROM, we show in Figure 8b the
mean estimated errors for the parameters in the test set $\Xi_{\text{test}}$.
The obtained errors are smaller than the desired tolerance $10^{-4}$, showing
the reliability of our error estimation approach. Furthermore, Figure 9a shows
the singular value decays of the defect vector trajectories at the $16$
parameter samples in $\Xi_{\text{defect}}$. Figure 9b shows the singular value
decay obtained from the SVD of $\mathbf{R}$ (Section 4.2.1). The singular
value decay in Figure 9b indicates the fast Kolmogorov $n$-width decay of the
defect manifold and hence good performance of Algorithm 2. In terms of
runtime, this approach requires $17$ seconds for generating the training data
for learning the defect vector and obtaining the RBF interpolants.

(a)  (b)

Figure 8: Burgers’ equation, Algorithm 2: (a) error (estimator) decay when
using the SVD+RBF method to approximate the closure term and Step 6 is not
included; (b) performance of the ROM over a test set.

(a)  (b)

Figure 9: Burgers’ equation: (a) normalized singular values of the defect
matrix $\mathbf{D}(\bm{\mu})$ for $\bm{\mu}\in\Xi_{\text{defect}}$; (b)
normalized singular values of $\mathbf{R}$.

##### POD-Greedy with Algorithm 2 and SVD+FNN closure approximation

Next, we repeat Algorithm 2 but now with the neural network-based
approximation of the defect. The feed-forward neural network has 3 hidden
layers with $16,64,64$ neurons, respectively. The activation function are the
SiLU function for the first three layers and Tanh for the last layer. We have
normalized the input data to be between $[0,1]$ and the output data is between
$[-1,1]$. The neural network is implemented in PyTorch. It is trained using
the Adam optimizer for $2000$ epochs, with the learning rate being set as
$0.005$. Initially, we do not implement Step 6 and do not update
$\widetilde{\mathbf{d}}_{\text{NN}}(t,\bm{\mu}^{*})$ with
$\mathbf{d}(t,\bm{\mu}^{*})$ at $\bm{\mu}^{*}$ (selected from the previous
iteration) when computing the error estimator
$\overline{\Delta}_{b}(\bm{\mu})$ in Step 11 at the current iteration. Using
the same $\Xi_{\text{defect}}$ with $|\Xi_{\text{defect}}|=16$ (as done for
the SVD+RBF approach), we did not obtain convergence of the greedy algorithm.
Therefore, we use a $\Xi_{\text{defect}}$ with $24$ samples. The SVD
tolerances for this case are both $0.1$ resulting in $n_{d}=6$. Using enriched
training data results in the successful convergence of Algorithm 2 as seen in
Figure 10a. However, it requires up to $14$ iterations for this convergence.
Evidently, the NN-based approach does not yield a satisfactory performance
even with more training data. Then, we implement Step 6 and update the defect
vector approximation at Step 11, which results in a significant improvement in
performance. This is shown in Figure 10b.

(a)  (b)

Figure 10: Burgers’ equation, Algorithm 2: (a) error (estimator) decay when
using the SVD+FNN method without closure updates; (b) error (estimator) decay
when using the SVD+FNN method with closure updates.

To explain the poorer performance of the SVD+FNN approach, we plot in Figure
11 the approximation of the defect vector at $\mu=0.123$ and $4$ different
time instances, viz., $t\in\\{0.01,0.2,0.5,1.5\\}$ s. We notice that the
approximation from the FNN, while qualitatively capturing the true defect,
fails to produce a very close match to the true value. This is especially the
case for latter time instances, as the magnitude of the defect vector gets
smaller. However, the RBF-based approach results in a significantly better
approximation. This might be because of the fact that each entry of the
reduced defect vector is separately learned by an individually-trained RBF
interpolation, while the whole reduced defect vector is learned by a single
and uniformly-trained FNN. However, if we use $n_{d}$ FNNs to learn the
$n_{d}$ entries of the defect vector separately, the training will become much
more expensive when $n_{d}$ is not very small.

Figure 11: Burgers’ equation: approximation of the true defect by the RBF-
based and NN-based approaches at $t\in\\{0.01,0.2.0.5,1.5\\}$ s for
$\mu=0.123$.

###  5.3 FitzHugh-Nagumo equations

##### Model description

The FitzHugh-Nagumo system models the response of an excitable neuron or cell
under an external stimulus. It finds applications in a variety of fields such
as cardiac electrophysiology and brain modeling. The nonlinear coupled system
of two partial differential equations defined in the domain $\Omega:=[0,L]$ is
given below:

$\displaystyle\epsilon\frac{\partial v_{1}(z,t)}{\partial t}$
$\displaystyle=\epsilon^{2}\frac{\partial^{2}v_{1}(z,t)}{\partial
z^{2}}+f(v_{1}(z,t))-v_{2}(z,t)+c,$ (45a) $\displaystyle\frac{\partial
v_{2}(z,t)}{\partial t}$ $\displaystyle=b\,v_{1}(z,t)-\gamma v_{2}(z,t)+c,$
(45b)

with boundary conditions

$\displaystyle\frac{\partial}{\partial z}v_{1}(0,t)$
$\displaystyle=-I_{\text{ext}}(t),\hskip 10.00002pt\frac{\partial}{\partial
z}v_{1}(L,t)=0,$ (46)

and initial conditions

$\displaystyle v_{1}(z,0)=0.001,\hskip 10.00002ptv_{2}(z,0)=0.001.$ (47)

In the above equations, $v_{1}(z,t)$ and $v_{2}(z,t)$ represent the electric
potential and the recovery rate of the potential, respectively. The spatial
variable is denoted by $z\in\Omega$ and the time $t\in[0,5]$. The nonlinear
term is represented by $f(v_{1}(z,t)):=v_{1}(v_{1}-0.1)(1-v_{1})$. The
external stimulus is $I_{\text{ext}}(t)=50000t^{3}e^{-15t}$. The system has
four free parameters $\epsilon,c,b,\gamma$. We fix $b=0.5$ and $\gamma=2$
while the two free parameters are
$\bm{\mu}=[\epsilon,c]\in\mathcal{P}:=[0.01,0.04]\times[0.025,0.075]$. A
finite difference scheme is employed to spatially discretize eqs. 45a and 45b
with $512$ nodes used for each variable leading to a FOM of dimension
$N=1024$. We sample $100$ parameters uniformly from the domain $\mathcal{P}$
and randomly divide them into the training set $\Xi$ and the test set
$\Xi_{\text{test}}$ in the ratio $70:30$. To solve the ODE, we use ode15s from
matlab®. The time discretization is done on a uniform grid with $\delta
t=0.01$. The output variables of interest are the values of the two state
variables at the node next to the leftmost boundary.

We apply Algorithms 1 and 2 to the FitzHugh-Nagumo system. This is a
particularly challenging example for both algorithms as the system exhibits a
slow singular value decay. Of particular interest is the approximation of the
limit cycle behaviour of the system for certain combinations of the two free
parameters, $(\epsilon,c)$. The RBM tolerance is set as
$\texttt{tol}=10^{-3}$.

##### POD-Greedy with Algorithm 1

Applying Algorithm 1 to this example and using a second-order IMEX scheme
(IMEX2) to compute the residual does not result in the convergence of the
greedy algorithm, as see in Figure 13. This stems from the fact that the
actual numerical scheme solver used in Algorithm 1 is the ode15s solver from
matlab®, therefore, the residual we compute using IMEX2 scheme is incorrect.

(a)  (b)

Figure 12: FitzHugh-Nagumo equations: (a) Stage 1 - normalized singular values
of the defect matrix $\mathbf{D}(\bm{\mu})$ for
$\bm{\mu}\in\Xi_{\text{defect}}$; (b) Stage 2 - normalized singular values of
$\mathbf{R}$. Figure 13: FitzHugh-Nagumo equations, Algorithm 1: error
(estimator) decay when using a second-order IMEX method without any closure
term.

(a)  (b)

Figure 14: FitzHugh-Nagumo equations, Algorithm 2: (a) error (estimator) decay
when using the SVD+RBF method with Step 6 included; (b) error (estimator)
decay when using the SVD+RBF method without Step 6.

##### POD-Greedy with Algorithm 2 and SVD+RBF closure approximation

We apply Algorithm 2 to this example using
$\widetilde{\mathbf{d}}_{\text{RBF}}(t,\bm{\mu})$ in Step 3 and with an RBM
tolerance $\texttt{tol}=10^{-3}$. The user-imposed time integration scheme
$\texttt{solver}_{\text{imp}}$ is IMEX2. Due to the challenging nature of the
problem, $d_{s}=21$ uniformly-spaced samples are chosen from the training set
to obtain $\Xi_{\text{defect}}$. The tolerance
$\texttt{tol}_{\text{SVD},t}=\texttt{tol}_{\text{SVD},\bm{\mu}}=10^{-6}$,
resulting in $n_{d}=311$. Figure 12a plots the singular value decays of
$\mathbf{D}(\bm{\mu})$ at all $\bm{\mu}\in\Xi_{\text{defect}}$ while Figure
12b shows the decay of the singular values of $\mathbf{R}$ (see Section
4.2.1). Similar to the case of the Burgers’ equation, an exponential decay of
the singular values is observed in Figure 12a. However, the second SVD shown
in Figure 12b has a relatively slower decay of the singular values. This shows
that the solution manifold of the FitzHugh-Nagumo system with respect to the
parameter variations is more difficult to be approximated by a low-dimensional
linear subspace. We incur $25$ seconds to compute the training data for the
defect trajectories and to obtain the SVD+RBF approximation of the defect
vectors. The greedy algorithm takes $11$ iterations to reach the desired
tolerance; Figure 14a shows the error convergence of Algorithm 2. The
dimension of the ROM obtained is $n=33$. Note that we have implemented Step 6
in Algorithm 2 to update the RBF approximation
$\widetilde{\mathbf{d}}_{\text{RBF}}(t,\bm{\mu}^{*})$ with
$\mathbf{d}(t,\bm{\mu}^{*})$ when we compute the error estimator
$\overline{\Delta}_{b}(\bm{\mu})$ for all $\bm{\mu}\in\Xi$ in Step 11. Without
doing so, the greedy algorithm converges nevertheless, but takes $16$
iterations (Figure 14b) and the ROM has a larger size $n=48$.

Figure 15: FitzHugh-Nagumo equations, Algorithm 2: error (estimator) decay
when using the SVD+FNN method with Step 6 included.

##### POD-Greedy with Algorithm 2 and SVD+FNN closure approximation

Next, we use $\widetilde{\mathbf{d}}_{\text{NN}}(t,\bm{\mu})$ in Step 3 of
Algorithm 2. The FNN is a $3$-layer network having, respectively, $64,64,32$
neurons in its hidden layers. The SiLU function is used for activation in all
but the last layer. In the last layer, Tanh is the activation function. The
training is carried out for $2000$ epochs using the Adam optimizer. The
learning rate is $0.002$. No special tuning was done to calibrate the
hyperparameters of the FNN. A detailed investigation on this is left for
future work. We set
$\texttt{tol}_{\text{SVD},t}=\texttt{tol}_{\text{SVD},\bm{\mu}}=10^{-3}$, such
that $n_{d}=57$. The total time for computing the training data at all the
samples in $\Xi_{\text{defect}}$ and for training the FNN is $106$ seconds,
where training the FNN dominates the total runtime. Figure 15 plots the
convergence of the greedy algorithm. It takes $10$ iterations to converge. The
resulting ROM has dimension $n=30$.

The finally derived ROM is then simulated at two parameter samples
$\bm{\mu}=(0.0267,0.0367)$ and $\bm{\mu}=(0.04,0.0472)$ taken from the test
set. Figure 16 shows the results of the ROM for the parameter
$\bm{\mu}=(0.0267,0.0367)$ obtained from Algorithm 2 using the SVD+RBF
approximation of the closure term. We see that the ROM is able to successfully
capture both the state and the output dynamics of the FOM at the test
parameter. At this parameter, the limit cycle behaviour is not very strong. At
a different test parameter ($\bm{\mu}=(0.04,0.0472)$) shown in Figure 17, the
ROM is able to successfully recover the stronger limit cycle behaviour as
well. For this case, we show results using the SVD+FNN approach. But, we note
that similar accuracy is also achieved with the SVD+RBF approach. The
corresponding results are not shown to avoid repetition.

(a)  (b)

Figure 16: FitzHugh-Nagumo equations, Algorithm 2 with SVD+RBF: performance at
the test parameter $\bm{\mu}=(0.0267,0.0367)$ (a) Limit cycle behaviour; (b)
output quantities.

(a)  (b)

Figure 17: FitzHugh-Nagumo equations, Algorithm 2 with SVD+FNN: performance at
the test parameter $\bm{\mu}=(0.04,0.0472)$ (a) Limit cycle behaviour; (b)
output quantities.

###  5.4 Batch chromatography

The last example we consider is the model of the batch chromatography
purification process. Being a coupled, nonlinear system of four PDEs, this
example poses a considerable challenge for our proposed approach.

##### Model description

Batch chromatography is an important purification process for separation of
chemicals in food and pharmaceutical industries. We consider the governing
equations for the batch chromatographic process for binary separation, i.e.,
the separation of two chemical components from a mixture. A schematic of the
complete process is shown in Figure 18.

The governing PDEs for the batch chromatography system are:

$\displaystyle\frac{\partial v_{1,z}}{\partial
t}+\frac{1-\epsilon}{\epsilon}\frac{\partial v_{2,z}}{\partial t}$
$\displaystyle=-\frac{\partial v_{1,z}}{\partial
z}+\frac{1}{\text{Pe}}\frac{\partial^{2}v_{1,z}}{\partial z^{2}},$ (48)
$\displaystyle\frac{\partial v_{2,z}}{\partial t}$
$\displaystyle=\dfrac{L}{Q/\epsilon
A_{c}}\kappa_{z}\big{(}v_{2,z}^{\text{\scriptsize{Eq}}}-v_{2,z}\big{)}$ (49)

where the state variables $v_{1,z},v_{2,z}$ refer to the concentrations of the
chemical component $z$ in the liquid and solid phase, respectively. Since we
are interested in binary separation, $z\in\\{a,b\\}$. The boundary conditions
are:

$\displaystyle\frac{\partial v_{1,z}(0,t)}{\partial
z}=\text{Pe}\big{(}v_{1,z}(0,t)-u(t)\big{)},\hskip 10.00002pt\frac{\partial
v_{1,z}(1,t)}{\partial z}=0$ (50)

and the initial conditions are:

$\displaystyle v_{1,z}(z,0)=0,\hskip 10.00002ptv_{2,z}(z,0)=0.$ (51)

The quantity $v_{2,z}^{\text{\scriptsize{Eq}}}$ in eq. 48 is the source of
nonlinearity and it denotes the adsorption equilibrium:

$\displaystyle
v_{2,z}^{\text{\scriptsize{Eq}}}=f_{z}(v_{1,a},v_{2,b}):=\dfrac{H_{z1}v_{1,z}}{1+K_{a1}v_{1,a}^{\text{f}}v_{1,a}+K_{b1}v_{2,b}^{\text{f}}v_{2,b}}+\dfrac{H_{z2}v_{1,z}}{1+K_{a2}v_{1,a}^{\text{f}}v_{1,a}+K_{b2}v_{2,b}^{\text{f}}v_{2,b}}.$

The discretization of the above PDE is performed using second-order finite
volume method. Each of the four PDEs is spatially discretized into $800$
volume elements, resulting in a FOM of dimension $N=3200$. For full details
regarding the batch chromatography PDE discretization and the terms involved,
we refer to [19]. The output variables are the concentrations of the liquid
phases ($v_{1,a},v_{1,b}$) at the rightmost node. The batch chromatography
model has two free parameters, $Q$ and $t_{\text{in}}$, which denote,
respectively, the volumetric feed flow of the solvent injected into the column
and the injection frequency of the solvent. We fix $t_{\text{in}}=0.5$ while
$Q\in[0.0667,0.1667]$. We collect $60$ uniformly-spaced samples of $Q$ and
divide them into the training set $\Xi$ and the test set $\Xi_{\text{test}}$
in the ratio $80:20$. The solver used for both algorithms is ode15s from
matlab®while the user-imposed time integration method
$\texttt{solver}_{\text{imp}}$ in Algorithm 2 is IMEX2. The time
discretization divides the time range $t\in[0,10]$ into a uniform grid with
$\delta t=0.005$.

Figure 18: Schematic of the batch chromatography process.

##### POD-Greedy with Algorithm 1

First, we show the results of applying Algorithm 1 to the batch chromatography
equations. The RBM tolerance is set to be $\texttt{tol}=10^{-3}$. Since the
exact form of the residual expression of batch chromatography equations is
unknown, (i.e., $\mathcal{R}[\cdot]$ corresponding to ode15s is not available)
we impose IMEX2 to get an approximate residual operator
$\widetilde{\mathcal{R}}[\cdot]$. However, $\widetilde{\mathcal{R}}[\cdot]$
corresponding to IMEX2 is different from $\mathcal{R}[\cdot]$ so that the
estimated error is inaccurate, resulting in stagnation of the greedy
algorithm, see Figure 19. To address this situation, we next apply the
proposed approach which involves adding a closure term.

##### POD-Greedy with Algorithm 2 and SVD+RBF closure approximation

As done for the last two examples, we first apply Algorithm 2 to the current
example and employ the RBF-based approach to approximate the closure term. The
greedy algorithm successfully converges in $10$ iterations to the desired
tolerance of $\texttt{tol}=10^{-3}$. The dimension of the ROM is $n=134$. The
convergence is shown in Figure 20. To learn the closure term, $d_{s}=10$
uniformly-spaced samples are chosen from the training set $\Xi$. The SVD
tolerances are
$\texttt{tol}_{\text{SVD},t}=\texttt{tol}_{\text{SVD},\bm{\mu}}=10^{-5}$,
leading to $n_{d}=291$. The singular value decay from the SVD of each defect
snapshot matrix $\mathbf{D}(\bm{\mu})$ for all
$\bm{\mu}\in\Xi_{\text{defect}}$ is presented in Figure 21a. It can be noticed
that the singular values decay even slower than those of the FitzHugh-Nagumo
model (see Figure 12a). In Figure 21b, the singular values of $\mathbf{R}$
also exhibit a much slower decay. This indicates that for the batch
chromatography example, both the dynamics at a given parameter and the
solution manifold with respect to the parameter variations are much more
difficult to be captured by a low-dimensional linear space. As a result, this
problem is likely to have a slow Kolmogorov $n$-width decay. The batch
chromatography equations are in fact a system of first-order hyperbolic PDEs
[42]. It is known that for problems exhibiting hyperbolic characteristics or
convection-dominance, the Kolmogorov $n$-width decay is slow [27]. The runtime
for obtaining the training data for the SVD+RBF approach and for obtaining the
RBF interpolant is $186$ seconds. To demonstrate the quality of the ROM
resulting from Algorithm 2, we compare the output obtained from the FOM and
ROM in Figure 22. The results are shown for the parameter sample $Q=0.0803$
taken from $\Xi_{\text{test}}$. We can see that the ROM is able to
successfully recover the dynamics of both output quantities. For this sample,
we obtained the mean error $\overline{\Delta}(Q)=2.082\cdot 10^{-5}$, which is
below the desired tolerance. Additionally, to show the approximation of the
state vector, we plot the space-time values and corresponding approximation
errors of the liquid phase concentration $v_{1,a},v_{1,b}$. The results are
shown in Figures 23 and 24. It can be inferred that the ROM delivers a sharp
approximation of both these state quantities for the entire duration of the
simulation at an unseen parameter during training. Note that the state vector
has error larger than the tolerance. The reason for this is that our error
estimator aims to estimate the output error rather than the whole state error.

Figure 19: Batch chromatography, Algorithm 1: error (estimator) decay when
using a second-order IMEX method without any closure term. Figure 20: Batch
chromatography, Algorithm 2: convergence when using the SVD+RBF method with
Step 6 included to approximate the closure term.

(a)  (b)

Figure 21: Batch chromatography: (a) Stage 1 - normalized singular values of
the defect matrix $\mathbf{D}(\bm{\mu})$ for $\bm{\mu}\in\Xi_{\text{defect}}$;
(b) Stage 2 - normalized singular values of $\mathbf{R}$.

The quality of the defect vector approximation using the SVD+FNN approach was
not satisfactory for this example. This owes to the particularly non-smooth
nature of the defect snapshots for different time instances and parameters.
Figure 25 shows the true reduced defect vector
$\widehat{\mathbf{d}}(t,\bm{\mu})$ and its approximation using RBF
interpolation $\widehat{\mathbf{d}}_{\text{RBF}}(t,\bm{\mu})$ corresponding to
$\bm{\mu}=0.0769$ at the time instances $t\in\\{0.05,0.25,2.5,5.0\\}$ s. While
the RBF interpolants recover a good approximation, the FNN was not successful
in capturing the entire complexity of the defect snapshots. The FNN-based
approximation was very inaccurate and we do not show those results. The neural
network struggles to capture the fast changing nature of the reduced defect
vector in the reduced coordinate space and also its wide range of magnitudes
($\pm(10^{-3}-10^{-8})$). The fact that we used separate RBF interpolants for
each coordinate and time instance led to a much better approximation than the
FNN.

Figure 22: Batch chromatography, Algorithm 2: performance of the ROM at test
parameter sample $Q=0.0803$ for the two output quantities $y_{1,a},y_{1,b}$.
Figure 23: Batch chromatography, Algorithm 2: performance of the ROM at test
parameter sample $Q=0.0803$ for the state $v_{1,a}$. Figure 24: Batch
chromatography, Algorithm 2: performance of the ROM at test parameter sample
$Q=0.0803$ for the state $v_{1,b}$.

##  6 Conclusion

In this work, we introduced a data-enhanced _a posteriori_ output error
estimator for model reduction of general parametric nonlinear dynamical
systems. The proposed error estimator does not require any knowledge of the
underlying time integration scheme used to integrate the given ODE. Applied to
the reduced basis method, the new approach enables the direct use of ODE
solver libraries, a feature that was not considered so far, to the best of our
knowledge. While it demands a modest amount of extra training data, the
proposed error estimator is efficient and can also be used in the online stage
to certify the accuracy of the ROMs. Numerical experiments performed on three
challenging examples demonstrate the benefits offered by the new approach. We
observed that the RBF-based approach performed better, compared to a NN-based
approach. An immediate extension of the proposed approach is to consider an
adaptive sampling of the training set $\Xi$. New strategies to accurately
approximate the closure term need to be considered for this. Another fruitful
line of future work could involve integrating the proposed methodology in the
Neural ODE framework [21] to make it fully non-intrusive.

Figure 25: Batch chromatography: approximation of the reduced defect vector by
the RBF-based approach at $t\in\\{0.05,0.25,2.5,5.0\\}$ s at $Q=0.0769$.

## Acknowledgments

Part of this work was performed while the first author was pursuing his
doctoral study and was supported by the International Max Planck Research
School in Process Systems Engineering (IMPRS-ProEng).

##  Appendix A Appendix A

### Proof of Theorem 2.1

###### Proof.

Consider the FOM in eq. 4 and its residual eq. 6 based on the ROM eq. 5.
Subtracting eq. 6 from eq. 4 and rearranging yields

$\displaystyle\mathbf{E}_{\text{im}}\mathbf{e}^{k}$
$\displaystyle=\mathbf{A}_{\text{im}}\mathbf{e}^{k-1}+\delta
t\bigg{(}\mathbf{f}(\mathbf{x}_{\text{im}}^{k-1})-\mathbf{f}(\widetilde{\mathbf{x}}_{\text{im}}^{k-1})\bigg{)}+\textbf{r}^{k},$
$\displaystyle\mathbf{e}^{k}$
$\displaystyle=\mathbf{E}_{\text{im}}^{-1}\mathbf{A}_{\text{im}}\mathbf{e}^{k-1}+\delta
t\mathbf{E}_{\text{im}}^{-1}\bigg{(}\mathbf{f}(\mathbf{x}_{\text{im}}^{k-1})-\mathbf{f}(\widetilde{\mathbf{x}}_{\text{im}}^{k-1})\bigg{)}+\mathbf{E}_{\text{im}}^{-1}\textbf{r}^{k}$
(52)

where
$\mathbf{e}^{k}=\|\mathbf{x}_{\text{im}}^{k}-\widetilde{\mathbf{x}}_{\text{im}}^{k}\|$
is the error in the state vector at the $k$-th time step.

Taking the norm on both sides and enforcing the Lipschitz condition results in

$\displaystyle\left\lVert\mathbf{e}^{k}\right\rVert_{2}$
$\displaystyle\leq\left\lVert\mathbf{E}_{\text{im}}^{-1}\mathbf{A}_{\text{im}}\right\rVert_{2}\left\lVert\mathbf{e}^{k-1}\right\rVert_{2}+\delta
tL_{\mathbf{f}}\left\lVert\mathbf{E}_{\text{im}}^{-1}\right\rVert_{2}\left\lVert\mathbf{e}^{k-1}\right\rVert_{2}+\left\lVert\mathbf{E}_{\text{im}}^{-1}\right\rVert_{2}\left\lVert\textbf{r}^{k}\right\rVert_{2}$
$\displaystyle\leq\bigg{(}\left\lVert\mathbf{E}_{\text{im}}^{-1}\mathbf{A}_{\text{im}}\right\rVert_{2}+\delta
tL_{\mathbf{f}}\left\lVert\mathbf{E}_{\text{im}}^{-1}\right\rVert_{2}\bigg{)}\left\lVert\mathbf{e}^{k-1}\right\rVert_{2}+\left\lVert\mathbf{E}_{\text{im}}^{-1}\right\rVert_{2}\left\lVert\textbf{r}^{k}\right\rVert_{2}.$
(53)

For notational convenience, we define
$\zeta:=\left\lVert\mathbf{E}_{\text{im}}^{-1}\right\rVert_{2}$ and
$\xi:=\big{(}\left\lVert\mathbf{E}_{\text{im}}^{-1}\mathbf{A}_{\text{im}}\right\rVert_{2}+\delta
tL_{\mathbf{f}}\left\lVert\mathbf{E}_{\text{im}}^{-1}\right\rVert_{2}\big{)}$
and rewrite appendix A as

$\displaystyle\left\lVert\mathbf{e}^{k}\right\rVert_{2}$
$\displaystyle\leq\xi\left\lVert\mathbf{e}^{k-1}\right\rVert_{2}+\zeta\left\lVert\textbf{r}^{k}\right\rVert_{2}.$
(54)

At $k=0$, the initial error is

$\displaystyle\left\lVert\mathbf{e}^{0}\right\rVert_{2}=\left\lVert\mathbf{x}_{\text{im}}^{0}-\mathbf{V}\mathbf{V}^{\textsf{T}}\mathbf{x}_{\text{im}}^{0}\right\rVert_{2}.$
(55)

Using eq. 55, the recursion in eq. 54 can be resolved for every $k$ to obtain
the expression for the error bound in eq. 7. ∎

##  Appendix B Appendix B

### Proof of Theorem 3.1

###### Proof.

The output error for the modified output in eq. 25 is

$\displaystyle\mathbf{y}^{k}-\overline{\mathbf{y}}_{\text{im,c}}^{k}=\mathbf{C}\big{(}\mathbf{x}^{k}-\widetilde{\mathbf{x}}_{\text{im,c}}^{k}\big{)}+\widetilde{\mathbf{x}}_{\text{du}}^{\textsf{T}}\textbf{r}_{\text{im,c}}^{k}.$
(56)

Multiplying by
$\big{(}\mathbf{x}^{k}-\widetilde{\mathbf{x}}_{\text{im,c}}^{k}\big{)}^{\textsf{T}}$
on both sides of eq. 20 yields

$\displaystyle\big{(}\mathbf{x}^{k}-\widetilde{\mathbf{x}}_{\text{im,c}}^{k}\big{)}^{\textsf{T}}\mathbf{E}_{\text{im}}^{\textsf{T}}\mathbf{x}_{\text{du}}$
$\displaystyle=-\big{(}\mathbf{x}^{k}-\widetilde{\mathbf{x}}_{\text{im,c}}^{k}\big{)}^{\textsf{T}}\mathbf{C}^{\textsf{T}},$

where we have made use of the fact that
$\mathbf{E}_{\text{du}}=\mathbf{E}_{\text{im}}^{\textsf{T}}$. Taking the
transpose on both sides of the above equality leads to

$\displaystyle\mathbf{x}_{\text{du}}^{\textsf{T}}\mathbf{E}_{\text{im}}\big{(}\mathbf{x}^{k}-\widetilde{\mathbf{x}}_{\text{im,c}}^{k}\big{)}$
$\displaystyle=-\mathbf{C}\big{(}\mathbf{x}^{k}-\widetilde{\mathbf{x}}_{\text{im,c}}^{k}\big{)}.$
(57)

We recall the auxiliary residual introduced in eq. 24 which can be written as

$\displaystyle\breve{\textbf{r}}_{\text{im,c}}^{k}$
$\displaystyle=\mathbf{A}_{\text{im}}\mathbf{x}_{\text{im,c}}^{k-1}+\delta
t\big{(}\mathbf{f}(\mathbf{x}_{\text{im,c}}^{k-1})+\mathbf{B}\mathbf{u}^{k}\big{)}+\mathbf{d}^{k}-\mathbf{E}_{\text{im}}\widetilde{\mathbf{x}}_{\text{im,c}}^{k},$
$\displaystyle=\mathbf{E}_{\text{im}}\big{(}\mathbf{x}_{\text{im,c}}^{k}-\widetilde{\mathbf{x}}_{\text{im,c}}^{k}\big{)},$
$\displaystyle=\mathbf{E}_{\text{im}}\big{(}\mathbf{x}^{k}-\widetilde{\mathbf{x}}_{\text{im,c}}^{k}\big{)}.$
(58)

The third equality above follows from eqs. 15 and 16 that the corrected
solution $\mathbf{x}_{\text{im,c}}^{k}$ actually recovers the FOM solution
$\mathbf{x}^{k}$ given an accurate closure term $\mathbf{d}^{k}$.

We further use the expression in appendix B to write eq. 57 as

$\displaystyle\mathbf{x}_{\text{du}}^{\textsf{T}}\breve{\textbf{r}}_{\text{im,c}}^{k}=-\mathbf{C}\big{(}\mathbf{x}^{k}-\widetilde{\mathbf{x}}_{\text{im,c}}^{k}\big{)}.$
(59)

Now, we substitute eq. 59 into eq. 56 followed by addition and subtraction of
the term
$\widetilde{\mathbf{x}}_{\text{du}}^{\textsf{T}}\breve{\textbf{r}}_{\text{im,c}}^{k}$
to get

$\displaystyle\mathbf{y}^{k}-\overline{\mathbf{y}}_{\text{im,c}}^{k}$
$\displaystyle=-\mathbf{x}_{\text{du}}^{\textsf{T}}\breve{\textbf{r}}_{\text{im,c}}^{k}+\widetilde{\mathbf{x}}_{\text{du}}^{\textsf{T}}\textbf{r}_{\text{im,c}}^{k},$
(60)
$\displaystyle=-\mathbf{x}_{\text{du}}^{\textsf{T}}\breve{\textbf{r}}_{\text{im,c}}^{k}+\widetilde{\mathbf{x}}_{\text{du}}^{\textsf{T}}\textbf{r}_{\text{im,c}}^{k}+\widetilde{\mathbf{x}}_{\text{du}}^{\textsf{T}}\breve{\textbf{r}}_{\text{im,c}}^{k}-\widetilde{\mathbf{x}}_{\text{du}}^{\textsf{T}}\breve{\textbf{r}}_{\text{im,c}}^{k},$
$\displaystyle=-\big{(}\mathbf{x}_{\text{du}}-\widetilde{\mathbf{x}}_{\text{du}}\big{)}^{\textsf{T}}\breve{\textbf{r}}_{\text{im,c}}^{k}+\widetilde{\mathbf{x}}_{\text{du}}^{\textsf{T}}\big{(}\textbf{r}_{\text{im,c}}^{k}-\breve{\textbf{r}}_{\text{im,c}}^{k}\big{)}.$

Subsequent to this, we use the expression for the dual system eq. 20 and its
residual eq. 23 we obtain

$\displaystyle\textbf{r}_{\text{du}}$
$\displaystyle=\mathbf{E}_{\text{du}}\mathbf{x}_{\text{du}}-\mathbf{E}_{\text{du}}\widetilde{\mathbf{x}}_{\text{du}}=\mathbf{E}_{\text{du}}\big{(}\mathbf{x}_{\text{du}}-\widetilde{\mathbf{x}}_{\text{du}}\big{)},$
(61) $\displaystyle\implies$
$\displaystyle\big{(}\mathbf{x}_{\text{du}}-\widetilde{\mathbf{x}}_{\text{du}}\big{)}=\mathbf{E}_{\text{du}}^{-1}\textbf{r}_{\text{du}}=\mathbf{E}_{\text{im}}^{-\textsf{T}}\textbf{r}_{\text{du}}.$

Substituting eq. 61 into eq. 60 yields

$\displaystyle\mathbf{y}^{k}-\overline{\mathbf{y}}_{\text{im,c}}^{k}=-\textbf{r}_{\text{du}}^{\textsf{T}}\mathbf{E}_{\text{im}}^{-1}\breve{\textbf{r}}_{\text{im,c}}^{k}+\widetilde{\mathbf{x}}_{\text{du}}^{\textsf{T}}\big{(}\textbf{r}_{\text{im,c}}^{k}-\breve{\textbf{r}}_{\text{im,c}}^{k}\big{)}.$
(62)

Taking the norm on either sides and using the triangle and Cauchy-Schwartz
inequalities we obtain the error bound as

$\displaystyle\left\lVert\mathbf{y}^{k}-\overline{\mathbf{y}}_{\text{im,c}}^{k}\right\rVert$
$\displaystyle=\|-\textbf{r}_{\text{du}}^{\textsf{T}}\mathbf{E}_{\text{im}}^{-1}\breve{\textbf{r}}_{\text{im,c}}^{k}\|+\|\widetilde{\mathbf{x}}_{\text{du}}^{\textsf{T}}\big{(}\textbf{r}_{\text{im,c}}^{k}-\breve{\textbf{r}}_{\text{im,c}}^{k}\big{)}\|,$
(63)
$\displaystyle\leq\|\mathbf{E}_{\text{im}}^{-1}\|\|\textbf{r}_{\text{du}}\|\|\breve{\textbf{r}}_{\text{im,c}}^{k}\|+\|\widetilde{\mathbf{x}}_{\text{du}}\|\|\textbf{r}_{\text{im,c}}^{k}-\breve{\textbf{r}}_{\text{im,c}}^{k}\big{)}\|.$
(64)

∎

## References

  * [1] Shrirang Abhyankar, Jed Brown, Emil M. Constantinescu, Debojyoti Ghosh, Barry F. Smith, and Hong Zhang.  PETSc/TS: A modern scalable ODE/DAE solver library.  Technical report, 2018.  arXiv:1806.01437.
  * [2] Uri M. Ascher, Steven J. Ruuth, and Brian T. R. Wetton.  Implicit-explicit methods for time-dependent partial differential equations.  SIAM J. Numer. Anal., 32(3):797–823, 1995.  doi:10.1137/0732037.
  * [3] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil M. Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, Jacob Faibussowitsch, William D. Gropp, Václav Hapla, Tobin Isaac, Pierre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose E. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, and Junchao Zhang.  PETSc Web page.  https://petsc.org/, 2023.  URL: https://petsc.org/.
  * [4] M. Barrault, Y. Maday, N. C. Nguyen, and A. T. Patera.  An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations.  C.R. Acad. Sci. Paris, 339(9):667–672, 2004.  doi:10.1016/j.crma.2004.08.006.
  * [5] P. Benner, S. Grivet-Talocia, A. Quarteroni, G. Rozza, and L. M. Schilder, W. Silveira, editors.  Model Order Reduction. Volume 1: System- and Data-Driven Methods and Algorithms.  De Gruyter, 2021.  doi:10.1515/9783110499001.
  * [6] P. Benner, S. Grivet-Talocia, A. Quarteroni, G. Rozza, and L. M. Schilder, W. Silveira, editors.  Model Order Reduction. Volume 2: Snapshot-Based Methods and Algorithms.  De Gruyter, 2021.  doi:10.1515/9783110671490.
  * [7] P. Benner, S. Grivet-Talocia, A. Quarteroni, G. Rozza, and L. M. Schilder, W. Silveira, editors.  Model Order Reduction. Volume 3: Applications.  De Gruyter, 2021.  doi:10.1515/9783110499001.
  * [8] Christopher M. Bishop.  Pattern recognition and machine learning.  Information Science and Statistics. Springer, New York, 2006.  doi:10.1007/978-0-387-45528-0.
  * [9] Annalisa Buffa, Yvon Maday, Anthony T. Patera, Christophe Prud’homme, and Gabriel Turinici.  A priori convergence of the greedy algorithm for the parametrized reduced basis method.  ESAIM Math. Model. Numer. Anal., 46(3):595–603, 2012.  doi:10.1051/m2an/2011056.
  * [10] M. D. Buhmann.  Radial Basis Functions: Theory and Implementations, volume 12 of Cambridge Monographs on Applied and Computational Mathematics.  Cambridge University Press, Cambridge, 2003.
  * [11] J. C. Butcher.  Numerical methods for ordinary differential equations.  John Wiley & Sons, Ltd., Chichester, Second edition, 2008.  doi:10.1002/9780470753767.
  * [12] C. Canuto, T. Tonn, and K. Urban.  A posteriori error analysis of the reduced basis method for nonaffine parametrized nonlinear PDEs.  SIAM J. Numer. Anal., 47(3):2001–2022, 2009.  doi:10.1137/080724812.
  * [13] K. Carlberg, C. Bou-Mosleh, and C. Farhat.  Efficient non-linear model reduction via a least-squares Petrov-Galerkin projection and compressive tensor approximations.  Internat. J. Numer. Methods Engrg., 86(2):155–181, 2011.  doi:10.1002/nme.3050.
  * [14] Fabien Casenave, Alexandre Ern, and Tony Lelièvre.  A nonintrusive reduced basis method applied to aeroacoustic simulations.  Adv. Comput. Math., 41(5):961–986, 2015.  doi:10.1007/s10444-014-9365-0.
  * [15] R. Chakir and J.K. Hammond.  A non-intrusive reduced basis method for elastoplasticity problems in geotechnics.  J. Comput. Appl. Math., 337:1–17, 2018.  doi:10.1016/j.cam.2017.12.044.
  * [16] Rachida Chakir and Yvon Maday.  A two-grid finite-element/reduced basis scheme for the approximation of the solution of parameter dependent PDE.  In 9e Colloque National en Calcul des Structures, Giens, France, 2009. CSMA.  URL: https://hal.archives-ouvertes.fr/hal-01420726.
  * [17] S. Chaturantabut and D. C. Sorensen.  Nonlinear model reduction via discrete empirical interpolation.  SIAM J. Sci. Comput., 32(5):2737–2764, 2010.  doi:10.1137/090766498.
  * [18] S. Chellappa, L. Feng, and P. Benner.  Adaptive basis construction and improved error estimation for parametric nonlinear dynamical systems.  Internat. J. Numer. Methods Engrg., 121(23):5320–5349, 2020.  doi:10.1002/nme.6462.
  * [19] Sridhar Chellappa.  A Posteriori Error Estimation and Adaptivity for Model Order Reduction of Large-Scale Systems.  Dissertation, Otto-von-Guericke-Universität, Magdeburg, Germany, 2023.  doi:http://dx.doi.org/10.25673/101396.
  * [20] Peng Chen, Alfio Quarteroni, and Gianluigi Rozza.  Reduced basis methods for uncertainty quantification.  SIAM/ASA J. Uncertain. Quantif., 5(1):813–869, 2017.  doi:10.1137/151004550.
  * [21] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud.  Neural ordinary differential equations.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.  URL: https://proceedings.neurips.cc/paper/2018/file/69386f6bb1dfed68692a24c8686939b9-Paper.pdf.
  * [22] Albert Cohen and Ronald DeVore.  Kolmogorov widths under holomorphic mappings.  IMA Journal of Numerical Analysis, 36(1):1–12, 03 2015.  doi:10.1093/imanum/dru066.
  * [23] Denise Degen, Karen Veroy, and Florian Wellmann.  Certified reduced basis method in geosciences: addressing the challenge of high-dimensional problems.  Comput. Geosci., 24(1):241–259, 2020.  doi:10.1007/s10596-019-09916-6.
  * [24] Markus Dihlmann and Bernard Haasdonk.  A reduced basis Kalman filter for parametrized partial differential equations.  ESAIM Control Optim. Calc. Var., 22(3):625–669, 2016.  doi:10.1051/cocv/2015019.
  * [25] M. Drohmann, B. Haasdonk, and M. Ohlberger.  Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation.  SIAM J. Sci. Comput., 34(2):A937–A969, 2012.  doi:10.1137/10081157X.
  * [26] David J Gardner, Daniel R Reynolds, Carol S Woodward, and Cody J Balos.  Enabling new flexibility in the SUNDIALS suite of nonlinear and differential/algebraic equation solvers.  ACM Transactions on Mathematical Software (TOMS), 2022.  doi:10.1145/3539801.
  * [27] Constantin Greif and Karsten Urban.  Decay of the Kolmogorov $N$-width for wave problems.  Appl. Math. Lett., 96:216–222, 2019.  doi:10.1016/j.aml.2019.05.013.
  * [28] M. Grepl.  Reduced-basis approximation a posteriori error estimation for parabolic partial differential equations.  PhD thesis, Massachussetts Institute of Technology (MIT), Cambridge, USA, 2005.  URL: http://dspace.mit.edu/handle/1721.1/7582.
  * [29] M. A. Grepl and A. T. Patera.  A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations.  ESAIM: Math. Model. Numer. Anal., 39(1):157–181, 2005.  doi:10.1051/m2an:2005006.
  * [30] Martin A. Grepl.  Certified reduced basis methods for nonaffine linear time-varying and nonlinear parabolic partial differential equations.  Math. Models Methods Appl. Sci., 22(3):1150015, 40, 2012.  doi:10.1142/S0218202511500151.
  * [31] Elise Grosjean and Yvon Maday.  Error estimate of the non-intrusive reduced basis (NIRB) two-grid method with parabolic equations.  e-prints 2211.08897, arXiv, 2022.  URL: https://arxiv.org/abs/2211.08897, doi:10.48550/arXiv.2211.08897.
  * [32] B. Haasdonk and M. Ohlberger.  Reduced basis method for finite volume approximations of parametrized linear evolution equations.  ESAIM: Math. Model. Numer. Anal., 42(2):277 – 302, 2008.  doi:10.1051/m2an:2008001.
  * [33] B. Haasdonk and M. Ohlberger.  Efficient reduced models and a posteriori error estimation for parametrized dynamical systems by offline/online decomposition.  Math. Comput. Model. Dyn. Syst., 17(2):145–161, 2011.  doi:10.1080/13873954.2010.514703.
  * [34] Dirk Hartmann, Matthias Herz, and Utz Wever.  Model Order Reduction a Key Technology for Digital Twins, pages 167–179.  Springer-Verlag, Cham, 2018.  doi:10.1007/978-3-319-75319-5_8.
  * [35] J. S. Hesthaven, G. Rozza, and B. Stamm.  Certified Reduced Basis Methods for Parametrized Partial Differential Equations.  SpringerBriefs in Mathematics. Springer International Publishing, 2016\.  doi:10.1007/978-3-319-22470-1.
  * [36] Alan C Hindmarsh.  ODEPACK, a systematized collection of ode solvers.  In R. S. Stepleman, editor, Scientific computing : applications of mathematics and computing to the physical sciences. Elsevier, 1983.
  * [37] Alan C Hindmarsh, Peter N Brown, Keith E Grant, Steven L Lee, Radu Serban, Dan E Shumaker, and Carol S Woodward.  SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers.  ACM Transactions on Mathematical Software (TOMS), 31(3):363–396, 2005.  doi:10.1145/1089014.1089020.
  * [38] Trever Hines.  Python package containing tools for radial basis function (RBF) applications.  https://github.com/treverhines/RBF, 2023.
  * [39] Zhongzhan Huang, Senwei Liang, Hong Zhang, Haizhao Yang, and Liang Lin.  Accelerating numerical solvers for large-scale simulation of dynamical system via neurvec.  e-prints 2208.03680, arXiv, 2022.  cs.CE.  URL: https://arxiv.org/abs/2208.03680.
  * [40] M. G. Kapteyn, D. J. Knezevic, D. B. P. Huynh, M. Tran, and K. E. Willcox.  Data-driven physics-based digital twins via a library of component-based reduced-order models.  Internat. J. Numer. Methods Engrg., 123(13):2986–3003, 2022.  doi:10.1002/nme.6423.
  * [41] Mark Kärcher, Sébastien Boyaval, Martin A. Grepl, and Karen Veroy.  Reduced basis approximation and a posteriori error bounds for 4D-Var data assimilation.  Optim. Eng., 19(3):663–695, 2018.  doi:10.1007/s11081-018-9389-2.
  * [42] Ju Weon Lee and Andreas Seidel-Morgenstern.  Solving hyperbolic conservation laws with active counteraction against numerical errors: Isothermal fixed-bed adsorption.  Chemical Engineering Science, 207:1309–1330, 2019.  doi:10.1016/j.ces.2019.07.053.
  * [43] L. Machiels, Y. Maday, and A. T. Patera.  Output bounds for reduced-order approximations of elliptic partial differential equations.  Comp. Meth. Appl. Mech. Eng., 190(26-27):3413–3426, 2001.  doi:10.1016/S0045-7825(00)00275-9.
  * [44] Allan Pinkus.  $n$-widths in approximation theory, volume 7 of Ergebnisse der Mathematik und ihrer Grenzgebiete (3) [Results in Mathematics and Related Areas (3)].  Springer-Verlag, Berlin, 1985.  doi:10.1007/978-3-642-69894-1.
  * [45] Michael Poli, Stefano Massaroli, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park.  Hypersolvers: Toward fast continuous-depth models.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 21105–21117. Curran Associates, Inc., 2020.  URL: https://proceedings.neurips.cc/paper/2020/file/f1686b4badcf28d33ed632036c7ab0b8-Paper.pdf.
  * [46] A. Quarteroni, A. Manzoni, and F. Negri.  Reduced Basis Methods for Partial Differential Equations, volume 92 of La Matematica per il 3+2.  Springer International Publishing, 2016.  doi:10.1007/978-3-319-15431-2.
  * [47] Daniel R Reynolds, David J Gardner, Carol S Woodward, and Rujeko Chinomona.  ARKODE: A flexible IVP solver infrastructure for one-step methods.  arXiv preprint arXiv:2205.14077, 2022.
  * [48] Dimitrios V. Rovas.  Reduced-Basis Output Bound Methods for Parametrized Partial Differential Equations.  PhD thesis, Massachussetts Institute of Technology (MIT), Cambridge, USA, 2003.  URL: https://dspace.mit.edu/handle/1721.1/16956.
  * [49] G. Rozza, D. B. P. Huynh, and A. T. Patera.  Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: application to transport and continuum mechanics.  Arch. Comput. Methods Eng., 15(3):229–275, 2008.  doi:10.1007/s11831-008-9019-9.
  * [50] Lawrence F. Shampine and Mark W. Reichelt.  The MATLAB ODE suite.  volume 18, pages 1–22. 1997.  Dedicated to C. William Gear on the occasion of his 60th birthday.  doi:10.1137/S1064827594276424.
  * [51] Xing Shen, Xiaoliang Cheng, and Kewei Liang.  Deep euler method: solving odes by approximating the local truncation error of the euler method.  e-prints 2003.09573, arXiv, 2020.  URL: https://arxiv.org/abs/2003.09573, doi:10.48550/arXiv.2003.09573.
  * [52] K. Veroy, C. Prud’Homme, D. V. Rovas, and A. T. Patera.  A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations.  In 16th AIAA Computational Fluid Dynamics Conference, Orlando, United States, 2003.  URL: https://hal.archives-ouvertes.fr/hal-01219051.
  * [53] Qian Wang, Jan S. Hesthaven, and Deep Ray.  Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem.  J. Comput. Phys., 384:289–307, 2019.  doi:10.1016/j.jcp.2019.01.031.
  * [54] Holger Wendland.  Scattered Data Approximation, volume 17 of Cambridge Monographs on Applied and Computational Mathematics.  Cambridge University Press, Cambridge, 2005.
  * [55] D. Wirtz, D. C. Sorensen, and B. Haasdonk.  A posteriori error estimation for DEIM reduced nonlinear dynamical systems.  SIAM J. Sci. Comput., 36(2):A311–A338, 2014.  doi:10.1137/120899042.
  * [56] Y. Zhang, L. Feng, S. Li, and P. Benner.  Accelerating PDE constrained optimization by the reduced basis method: application to batch chromatography.  Internat. J. Numer. Methods Engrg., 104(11):983–1007, 2015.  doi:10.1002/nme.4950.
  * [57] Y. Zhang, L. Feng, S. Li, and P. Benner.  An efficient output error estimation for model order reduction of parametrized evolution equations.  SIAM J. Sci. Comput., 37(6):B910–B936, 2015.  doi:10.1137/140998603.