# Bayesian operator inference for data-driven reduced-order modeling

Mengwu Guo Corresponding author. E-mail<EMAIL_ADDRESS>Department of
Applied Mathematics, University of Twente  Shane A. McQuarrie  Oden Institute
for Computational Engineering and Sciences, The University of Texas at Austin
Karen E. Willcox  Oden Institute for Computational Engineering and Sciences,
The University of Texas at Austin

###### Abstract

This work proposes a Bayesian inference method for the reduced-order modeling
of time-dependent systems. Informed by the structure of governing equations,
the task of learning a reduced-order model from data is posed as a Bayesian
inversion problem with Gaussian prior and likelihood. The operators defining
the reduced-order model, rather than being chosen deterministically, are
characterized probabilistically as posterior Gaussian distributions. This
embeds uncertainty into the reduced-order model, and hence the predictions
subsequently issued by the reduced-order model are endowed with uncertainty.
The learned reduced-order models are computationally efficient, which enables
Monte Carlo sampling over the posterior distributions of reduced-order
operators. Furthermore, the proposed Bayesian framework provides a statistical
interpretation of the Tikhonov regularization incorporated in the operator
inference, and the empirical Bayes approach of maximum marginal likelihood
suggests a selection algorithm for the regularization hyperparameters. The
proposed method is demonstrated by two examples: the compressible Euler
equations with noise-corrupted observations, and a single-injector combustion
process.

Keywords: data-driven reduced-order modeling; uncertainty quantification;
operator inference; Bayesian inversion; Tikhonov regularization; single-
injector combustion

##  1 Introduction

The objective of this work is to use the tools and perspectives of Bayesian
inference to enable uncertainty quantification for data-driven reduced-order
modeling. A reduced-order model (ROM) is a low-dimensional surrogate for a
high-dimensional physical system, such as the many natural and engineering
systems that are governed by partial differential equations. Ideally, the
computational cost of a ROM should be substantially reduced compared to that
of a high-fidelity model of the original system, yet without significantly
compromising the solution accuracy. Model reduction strategies typically
construct a ROM by identifying a predominant, low-dimensional latent space
that well approximates the solution manifold, then representing the reduced-
state dynamics in the latent space [4, 43]. It is crucial for the validation
of a ROM in terms of accuracy and robustness to quantify the errors and
uncertainties introduced by the modeling decisions of such a low-dimensional
approximation.

Several data-driven methods for identifying reduced-order systems have been
proposed recently, including dynamic mode decomposition [45], operator
inference that preserves the formulation structure of governing equations [15,
40], sparse regression to identify reduced-state latent dynamics [44, 8],
manifold learning using deep auto-encoders [30], data-driven approximation of
time-integration schemes [56], and Gaussian process modeling for parametric
reduced basis methods [17]. These approaches learn the low-dimensional
representations of reduced-order systems primarily from solution data, and do
not intrusively involve explicit projections of the governing physics or a
high-fidelity simulator. Different from the certification of intrusive ROMs
constructed by direct projections, which is usually regarded as an _a
posteriori_ error estimation problem, the validation of data-driven ROMs
should mainly address the various uncertainties introduced by the use of data.
Therefore, it is a natural choice to employ a probabilistic framework to
enable the uncertainty quantification of data-driven ROMs.

This work focuses on the uncertainty quantification of ROMs learned through
_operator inference_ (OpInf), a non-intrusive, projection-based model
reduction technique introduced in [40]. Operator inference learns a ROM with
polynomial structure inspired by known governing equations, inferring reduced-
order operators by a linear regression of state observations after projection
onto a low-dimensional subspace. Through data-driven discovery informed by the
basic structure of governing physics, OpInf forms a ‘glass-box’ approach to
the representation of reduced-state dynamics [15]. In the OpInf framework,
training data may be pre-processed with scaling or variable transformations to
induce an exactly or approximately polynomial structure [29, 42]. The learning
problem requires time derivatives of the projected state observations, which
are often estimated from the states with finite differences. Though OpInf has
usually been considered in a deterministic setting, the recent work [50]
considers an active sampling strategy for overcoming noise pollution in the
training data. Even without noise in the data, establishing error bounds or
uncertainty guarantees in OpInf is challenging because the learning is driven
by data, although [49] establishes probabilistic error estimates in the case
where the reduced-order operators can be ensured to be the same as those
determined by explicit projection [39].

In this work, the task of reduced-order operator inference is posed as a
Bayesian inversion problem. Bayesian inference [6] incorporates prior
knowledge and observational information to update the probabilistic
distribution of unknown system parameters, i.e., the entries of the ROM
operators in our case. The posterior distribution thus provides probabilistic
descriptions of model uncertainties. Bayesian methods are especially useful in
data-rich scenarios [9] and, in particular, have recently shown promise in
estimating the operators in differential equations, see, e.g., [21, 38, 55].
Our main contribution centers on establishing a probabilistic interpretation
of the regression problem at the heart of OpInf, which endows the resulting
ROM with intrinsic uncertainty bounds. As OpInf is formulated as a linear
regression problem, multivariate normal distributions are adopted to define
the prior and likelihood when the Bayes’ rule is applied, yielding a Bayesian
linear regression to estimate the reduced-order operators. Different from the
deterministic OpInf workflow that gives point estimates only, the proposed
Bayesian version derives posterior distributions of the reduced-order
operators, whose uncertainties are subsequently propagated into solutions
through time integration. In other words, endowing the reduced-order operators
with probabilistic descriptions results in ROM predictions with quantifiable
uncertainties. Additionally, the Bayesian interpretation of OpInf provides
insight into the numerical regularization that is often required in the
learning problem to address ill-conditioning, over-fitting, or model
inadequacy [31, 48], which also leads to a selection algorithm for the
regularization hyperparameters.

We present this work in honor of the lifetime academic achievements of
Professor J. Tinsley Oden, especially his seminal work in highlighting the
importance of _verification & validation_ [1, 2, 37] and _uncertainty
quantification_ [34, 35, 36] to achieve predictive models in computational
engineering and sciences [34, 35]. In particular, he laid the mathematical
foundation of predictive computational science [35] using Bayesian methods and
has utilized the framework in a substantial range of scientific, medical and
engineering applications [13, 20, 26, 41]. Professor Oden’s pioneering,
groundbreaking research on predictive computational science sets the stage for
this work and inspires us to enable the uncertainty quantification and model
validation of data-driven ROMs through a Bayesian approach.

The paper is structured as follows. Section 2 introduces OpInf in its
established deterministic framework. Section 3 endows OpInf with a Bayesian
perspective by establishing a posterior distribution for the operators of a
reduced-order system, and Section 4 discusses the implications for the
numerical regularization that is needed in the ROM learning formulation.
Numerical results for two examples are presented in Section 5, and Section 6
makes concluding remarks.

##  2 Deterministic operator inference

We consider problems governed by partial differential equations written, after
spatial discretization, in the form

$\displaystyle\frac{\textup{d}}{\textup{d}t}{\mathbf{q}}(t)={\mathbf{F}}({\mathbf{q}}(t),{\mathbf{u}}(t))\,,\quad\mathbf{q}(t_{0})=\mathbf{q}_{0}\,,\quad
t\in[t_{0},t_{f}]\,,$ (2.1)

in which $\mathbf{q}(t)\in\mathbb{R}^{n}$ is the state vector collecting the
values of $n$ degrees of freedom at time $t$, $\mathbf{u}(t)\in\mathbb{R}^{m}$
collects the $m$ input values at time $t$ (e.g., encoding boundary conditions
or time-dependent source terms),
$\mathbf{F}:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n}$ maps
state/input pairs to the state evolution, $t_{0}$ and $t_{f}$ are the initial
and final time, respectively, and $\mathbf{q}_{0}$ is the prescribed initial
condition. A projection-based ROM approximates the full-order system (2.1) in
the subspace spanned by basis vectors
$\mathbf{V}=[\mathbf{v}_{1}~{}\mathbf{v}_{2}~{}\cdots~{}\mathbf{v}_{r}]\in\mathbb{R}^{n\times
r}$ with reduced dimensionality $r\ll n$. One widely used approach to
construct the basis $\mathbf{V}$ is the proper orthogonal decomposition (POD)
[5, 16, 46]. For a collection of known state snapshots
$\mathbf{q}_{j}:=\mathbf{q}(t_{j})$ and inputs
$\mathbf{u}_{j}:=\mathbf{u}(t_{j})$, $j=0,1,\cdots,k-1$, at $k$ time instances
$t_{0}<t_{1}<\cdots<t_{k-1}$, organized as
$\mathbf{Q}=[~{}\mathbf{q}_{0}~{}\mathbf{q}_{1}~{}\cdots~{}\mathbf{q}_{k-1}~{}]\in\mathbb{R}^{n\times
k}$ and
$\mathbf{U}=[~{}\mathbf{u}_{0}~{}\mathbf{u}_{1}~{}\cdots~{}\mathbf{u}_{k-1}~{}]\in\mathbb{R}^{m\times
k}$, POD takes advantage of the singular value decomposition
$\mathbf{Q}={\bf\it\Phi}{\bf\it\Xi}{\bf\it\Psi}^{\mathsf{T}}$ by using the
first $r$ columns of ${\bf\it\Phi}$ as the $r$-dimensional reduced basis,
i.e., $\mathbf{V}={\bf\it\Phi}_{:,1:r}$. Projection-based methods such as POD
have the key advantage of producing ROMs that inherit structure from the full-
order system (2.1) [16]. This observation motivates us to consider, for a wide
class of problems, ROMs with the following polynomial structure:

$\frac{\textup{d}}{\textup{d}t}\hat{\mathbf{q}}(t)=\hat{\mathbf{A}}\hat{\mathbf{q}}(t)+\hat{\mathbf{H}}[\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)]+\hat{\mathbf{B}}{\mathbf{u}}(t)+\hat{\mathbf{c}}\,,\quad\hat{\mathbf{q}}(t_{0})=\mathbf{V}^{\mathsf{T}}\mathbf{q}_{0}\,,\quad
t\in[t_{0},t_{f}]\,,$ (2.2)

where $\hat{\mathbf{q}}(t)\in\mathbb{R}^{r}$ is the reduced state, $\otimes$
denotes the Kronecker product [51], and
$\hat{\mathbf{A}}\in\mathbb{R}^{r\times r}$,
$\hat{\mathbf{H}}\in\mathbb{R}^{r\times r^{2}}$,
$\hat{\mathbf{B}}\in\mathbb{R}^{r\times m}$, and
$\hat{\mathbf{c}}\in\mathbb{R}^{r}$ are reduced-order operators that together
comprise a low-dimensional representation of ${\mathbf{F}}$. The full- and
reduced-order states are related by a linear approximation
$\mathbf{q}(t)\approx\mathbf{V}\hat{\mathbf{q}}(t)$.

The data-driven operator inference (OpInf) method [40] employs the projection-
based setting of (2.2), but utilizes a regression problem to find the reduced-
order operators that yield the ROM that best matches the snapshot data in a
minimum-residual sense. Such a non-intrusive approach avoids entailing
intrusive queries and access to source code. The regression problem to infer
the reduced-order operators is written as

$\min_{\hat{\mathbf{A}},\hat{\mathbf{H}},\hat{\mathbf{B}},\hat{\mathbf{c}}}\left\\{\sum_{j=0}^{k-1}\left\|\hat{\mathbf{A}}\hat{\mathbf{q}}_{j}+\hat{\mathbf{H}}[\hat{\mathbf{q}}_{j}\otimes\hat{\mathbf{q}}_{j}]+\hat{\mathbf{B}}{\mathbf{u}}_{j}+\hat{\mathbf{c}}-\dot{\hat{\mathbf{q}}}_{j}\right\|_{2}^{2}+\mathcal{P}\left([~{}\hat{\mathbf{A}}~{}~{}\hat{\mathbf{H}}~{}~{}\hat{\mathbf{B}}~{}~{}\hat{\mathbf{c}}~{}]\right)\right\\}\,,$
(2.3)

in which the first term is a loss function defined by the residual of (2.2) at
the $k$ time instances $t_{0}<t_{1}<\cdots<t_{k-1}$, and the second term
$\mathcal{P}$ represents a regularization. Here the reduced-state dataset
$\\{\hat{\mathbf{q}}_{j}\\}_{j=0}^{k-1}$ is constructed by evaluating the
projection of known snapshots onto the reduced basis $\mathbf{V}$, i.e.,
$\hat{\mathbf{q}}_{j}=\mathbf{V}^{\mathsf{T}}\mathbf{q}_{j}$, and the time-
derivative data $\\{\dot{\hat{\mathbf{q}}}_{j}\\}_{j=0}^{k-1}$ are then
typically approximated using finite differences. The term $\mathcal{P}$ is a
regularization term that penalizes the entries of the learned quantities,
which is often necessary to induce stability in the resulting ROM [31]. The
training time instances are contained in $[t_{0},t_{k-1}]$, a subset of the
full time interval of interest $[t_{0},t_{f}]$, and the reduced-order solution
over $[t_{k},t_{f}]$ will be predictive through the inferred ROM.

Remark 1: The quadratic interaction
$\hat{\mathbf{q}}\otimes\hat{\mathbf{q}}=\text{vec}(\hat{\mathbf{q}}\hat{\mathbf{q}}^{\mathsf{T}})\in\mathbb{R}^{r^{2}}$
has only $r(r+1)/2<r^{2}$ degrees of freedom due to symmetry, i.e., the upper-
triangular elements of $\hat{\mathbf{q}}\hat{\mathbf{q}}^{\mathsf{T}}$
completely prescribe $\hat{\mathbf{q}}\otimes\hat{\mathbf{q}}$. Hence, in
practice we may learn an $r\times r(r+1)/2$ quadratic operator
$\hat{\mathbf{H}}$ as opposed to an $r\times r^{2}$ quadratic operator.

The least-squares problem (2.3) can be rewritten in a matrix form and
decoupled into $r$ least-squares problems, one for each ordinary differential
equation in the reduced system (2.2), i.e.,

$\min_{\hat{\mathbf{O}}}\left\\{\left\|\mathbf{D}\hat{\mathbf{O}}^{\mathsf{T}}-\mathbf{R}^{\mathsf{T}}\right\|_{F}^{2}+\mathcal{P}(\hat{\mathbf{O}})\right\\}=\sum_{i=1}^{r}\min_{\hat{\mathbf{o}}_{i}}\left\\{\left\|\mathbf{D}\hat{\mathbf{o}}_{i}-\mathbf{r}_{i}\right\|_{2}^{2}+\mathcal{P}_{i}(\hat{\mathbf{o}}_{i})\right\\}\,,$
(2.4)

where

$\displaystyle\begin{aligned}
\hat{\mathbf{O}}&:=[~{}\hat{\mathbf{A}}~{}~{}\hat{\mathbf{H}}~{}~{}\hat{\mathbf{B}}~{}~{}\hat{\mathbf{c}}~{}]:=[~{}\hat{\mathbf{o}}_{1}~{}~{}\cdots~{}~{}\hat{\mathbf{o}}_{r}~{}]^{\mathsf{T}}\in\mathbb{R}^{r\times
d(r,m)}\,,\\\
\mathbf{D}&:=[~{}\hat{\mathbf{Q}}^{\mathsf{T}}~{}~{}(\hat{\mathbf{Q}}\odot\hat{\mathbf{Q}})^{\mathsf{T}}~{}~{}\mathbf{U}^{\mathsf{T}}~{}~{}\mathbf{1}_{k}~{}]\in\mathbb{R}^{k\times
d(r,m)}\,,\\\
\hat{\mathbf{Q}}&:=[~{}\hat{\mathbf{q}}_{0}~{}~{}\hat{\mathbf{q}}_{1}~{}~{}\cdots~{}~{}\hat{\mathbf{q}}_{k-1}~{}]\in\mathbb{R}^{r\times
k}\,,\\\
\mathbf{R}&:=[~{}\dot{\hat{\mathbf{q}}}_{0}~{}~{}\dot{\hat{\mathbf{q}}}_{1}~{}~{}\cdots~{}~{}\dot{\hat{\mathbf{q}}}_{k-1}~{}]:=[~{}\mathbf{r}_{1}~{}~{}\cdots~{}~{}\mathbf{r}_{r}~{}]^{\mathsf{T}}\in\mathbb{R}^{r\times
k}\,,\quad\text{and}\\\
\mathbf{U}&:=[~{}\mathbf{u}_{0}~{}~{}\mathbf{u}_{1}~{}~{}\cdots~{}~{}\mathbf{u}_{k-1}~{}]\in\mathbb{R}^{m\times
k}\,,\end{aligned}$ (2.5)

in which $d(r,m)=r+r(r+1)/2+m+1$, $\mathbf{1}_{k}\in\mathbb{R}^{k}$ is a
length-$k$ vector with all entries set to 1, $\odot$ is the Khatri-Rao product
(the column-wise Kronecker product [27, 28]), and
$\hat{\mathbf{o}}_{i}\in\mathbb{R}^{d(r,m)}$ and
$\mathbf{r}_{i}\in\mathbb{R}^{k}$ are the $i$-th rows of $\hat{\mathbf{O}}$
and $\mathbf{R}$, respectively. An individual regularization $\mathcal{P}_{i}$
is defined for each $\hat{\mathbf{o}}_{i}$, i.e., it is assumed that
$\mathcal{P}(\hat{\mathbf{O}})=\sum_{i=1}^{r}\mathcal{P}_{i}(\hat{\mathbf{o}}_{i})$.
Hence OpInf comprises $r$ independent linear regressions with respect to the
unknown reduced-order operators $\hat{\mathbf{o}}_{i}$. While (2.4) is a
linear data-driven regression, the corresponding ROM (2.2) defined the
$\hat{\mathbf{o}}_{i}$ is a nonlinear dynamical system.

We usually have $k>d(r,m)$; in that case, and with
$\mathcal{P}_{i}(\hat{\mathbf{o}}_{i})=0$, the point estimate
$\tilde{\mathbf{o}}_{i}$ of $\hat{\mathbf{o}}_{i}$ through the least squares
(2.4) can be yielded in an explicit form as

$\displaystyle\tilde{\mathbf{o}}_{i}:=\mathbf{D}^{+}\mathbf{r}_{i}=(\mathbf{D}^{\mathsf{T}}\mathbf{D})^{-1}\mathbf{D}^{\mathsf{T}}\mathbf{r}_{i}\,,$
(2.6)

$\mathbf{D}^{+}$ denoting the Moore–Penrose inverse of the matrix
$\mathbf{D}$. Here we assume that $\mathbf{D}$ has full column rank; the cases
with ill-conditioned $\mathbf{D}^{\mathsf{T}}\mathbf{D}$ will be discussed in
Section 4.

##  3 Bayesian operator inference

In this section, the data-driven regression for learning the reduced-order
operators
$\hat{\mathbf{O}}=[~{}\hat{\mathbf{A}}~{}~{}\hat{\mathbf{H}}~{}~{}\hat{\mathbf{B}}~{}~{}\hat{\mathbf{c}}~{}]$
is posed as a Bayesian inference problem [6, 9]. We define a probabilistic ROM
through the resulting posterior distribution of $\hat{\mathbf{O}}$ and the
form (2.2), which enables uncertainty estimates through Monte-Carlo sampling.
We abbreviate the proposed _Bayesian operator inference_ method as BayesOpInf.

###  3.1 Bayesian inference of reduced-order operators

The quadratic structure of the ROM (2.2) introduces a model misspecification
error whenever the full-order system (2.1) includes non-quadratic dynamics.
Additionally, the quality of the least-squares estimation (2.4) may suffer
from approximation error in the estimated time derivatives of the training
snapshots and/or data noise in the snapshots themselves. To account for these
uncertainties, we use Bayesian inference to formulate the estimators of the
reduced-order operators $\hat{\mathbf{o}}_{i}$ as probabilistic distributions,
$i=1,\ldots,r$.

For the definition of a likelihood distribution, we model the residual error
in the data-driven scheme (2.4) as independent Gaussian noise, i.e., we
consider

$\displaystyle\mathbf{r}_{i}=\mathbf{D}\hat{\mathbf{o}}_{i}+\boldsymbol{\epsilon}_{i},\qquad
p(\boldsymbol{\epsilon}_{i})=\mathcal{N}(\boldsymbol{\epsilon}_{i}|\mathbf{0}_{k},\sigma_{i}^{2}\mathbf{I}_{k})\,,$
(3.1)

in which $\mathbf{0}_{k}$ is the $k$-dimensional zero vector and
$\mathbf{I}_{k}$ denotes the $k\times k$ identity matrix. Here the noise
${\bf\it\epsilon}_{i}\in\mathbb{R}^{k}$ is a random vector consisting of the
error terms in the observations at $k$ time instances. All these noise
measurements are assumed to follow independent normal distributions with zero
mean and variance $\sigma_{i}^{2}$. Hence the likelihood distribution for the
training data pair $(\mathbf{D},\mathbf{r}_{i})$ is written as

$\displaystyle
p(\mathbf{r}_{i}|~{}\mathbf{D},\hat{\mathbf{o}}_{i},\sigma_{i}^{2})=\mathcal{N}(\mathbf{r}_{i}|~{}\mathbf{D}\hat{\mathbf{o}}_{i},\sigma_{i}^{2}\mathbf{I}_{k})\,.$
(3.2)

The choice to model the entries of $\boldsymbol{\epsilon}_{i}$ as independent
(i.e., assuming the noise is not correlated across time) is convenient but not
a requirement. An alternative, kernel-based model of the residual error is
provided in Appendix 1, which accounts for correlations over time and
constructs a surrogate model of the ROM closure.

We also define the prior distribution of $\hat{\mathbf{o}}_{i}$ as a
multivariate normal distribution,

$\displaystyle
p(\hat{\mathbf{o}}_{i}|~{}\sigma_{i}^{2},{\bf\it\lambda}_{i})=\mathcal{N}(\hat{\mathbf{o}}_{i}|~{}{\bf\it\beta}_{i},\sigma_{i}^{2}\mathrm{diag}({\bf\it\lambda}_{i})^{-1})\,,$
(3.3)

where ${\bf\it\beta}_{i}\in\mathbb{R}^{d(r,m)}$ is a pre-defined mean vector
and ${\bf\it\lambda}_{i}\in\mathbb{R}^{d(r,m)}$ is a vector of positive
numbers hyperparametrizing the prior variances of $\hat{\mathbf{o}}_{i}$. Note
that the prior covariance is a diagonal matrix, implying that the prior does
not impose correlations between the components of $\hat{\mathbf{o}}_{i}$.

Conditioning on the training data, Bayes’ rule gives the posterior
distribution
$p(\hat{\mathbf{o}}_{i}|~{}\mathbf{D},\mathbf{r}_{i},\sigma_{i}^{2},{\bf\it\lambda}_{i})$
as follows:

$\begin{split}p(\hat{\mathbf{o}}_{i}|~{}\mathbf{D},\mathbf{r}_{i},\sigma_{i}^{2},{\bf\it\lambda}_{i})&\propto
p(\mathbf{r}_{i}|~{}\mathbf{D},\hat{\mathbf{o}}_{i},\sigma_{i}^{2})~{}p(\hat{\mathbf{o}}_{i}|~{}\sigma_{i}^{2},{\bf\it\lambda}_{i})\\\
&=\mathcal{N}(\mathbf{r}_{i}|~{}\mathbf{D}\hat{\mathbf{o}}_{i},\sigma_{i}^{2}\mathbf{I}_{k})~{}\mathcal{N}(\hat{\mathbf{o}}_{i}|~{}{\bf\it\beta}_{i},\sigma_{i}^{2}\mathrm{diag}({\bf\it\lambda}_{i})^{-1})=\mathcal{N}(\mathbf{D}\hat{\mathbf{o}}_{i}|~{}\mathbf{r}_{i},\sigma_{i}^{2}\mathbf{I}_{k})~{}\mathcal{N}(\hat{\mathbf{o}}_{i}|~{}{\bf\it\beta}_{i},\sigma_{i}^{2}\mathrm{diag}({\bf\it\lambda}_{i})^{-1})\\\
&=\mathcal{N}\left(\mathbf{r}_{i}|~{}\mathbf{D}{\bf\it\beta}_{i},\sigma_{i}^{2}\mathbf{D}\mathrm{diag}({\bf\it\lambda}_{i})^{-1}\mathbf{D}^{\mathsf{T}}+\sigma_{i}^{2}\mathbf{I}_{k}\right)~{}\mathcal{N}(\hat{\mathbf{o}}_{i}|~{}{\bf\it\mu}_{i},{\bf\it\Sigma}_{i})\,,\quad\text{i.e.,}\\\
p(\hat{\mathbf{o}}_{i}|~{}\mathbf{D},\mathbf{r}_{i},\sigma_{i}^{2},{\bf\it\lambda}_{i})&=\mathcal{N}(\hat{\mathbf{o}}_{i}|~{}{\bf\it\mu}_{i},{\bf\it\Sigma}_{i})\,,\end{split}$
(3.4)

in which

$\begin{split}{\bf\it\Sigma}_{i}&={\bf\it\Sigma}_{i}(\mathbf{D},\sigma_{i}^{2},{\bf\it\lambda}_{i})=\sigma_{i}^{2}\left[\mathrm{diag}({\bf\it\lambda}_{i})+\mathbf{D}^{\mathsf{T}}\mathbf{D}\right]^{-1}\\\
{\bf\it\mu}_{i}&={\bf\it\mu}_{i}(\mathbf{r}_{i},\mathbf{D},{\bf\it\lambda}_{i})=\left[\mathrm{diag}({\bf\it\lambda}_{i})+\mathbf{D}^{\mathsf{T}}\mathbf{D}\right]^{-1}\left(\mathrm{diag}({\bf\it\lambda}_{i}){\bf\it\beta}_{i}+\mathbf{D}^{\mathsf{T}}\mathbf{r}_{i}\right)\\\
&={\bf\it\beta}_{i}+\underbrace{\left[\mathrm{diag}({\bf\it\lambda}_{i})+\mathbf{D}^{\mathsf{T}}\mathbf{D}\right]^{-1}\mathbf{D}^{\mathsf{T}}(\mathbf{r}_{i}-\mathbf{D}{\bf\it\beta}_{i})}_{\delta{\bf\it\mu}_{i}(\mathbf{r}_{i},\mathbf{D},{\bf\it\lambda}_{i})}\,,\end{split}$
(3.5)

all consistent with the results of Bayesian linear regression [6].

To marginalize $(\sigma_{i}^{2},{\bf\it\lambda}_{i})$ in the distribution
$p(\hat{\mathbf{o}}_{i}|~{}\mathbf{D},\mathbf{r}_{i},\sigma_{i}^{2},{\bf\it\lambda}_{i})$
and obtain $p(\hat{\mathbf{o}}_{i}|~{}\mathbf{D},\mathbf{r}_{i})$ as the
inference of $\hat{\mathbf{o}}_{i}$, we adopt an empirical Bayes method by
maximizing the log marginal likelihood (evidence) [54] to estimate a set of
optimal values of $(\sigma_{i}^{2},{\bf\it\lambda}_{i})$. Mathematically, such
a maximization problem is written as

$\begin{split}(\sigma^{*2}_{i},{\bf\it\lambda}_{i}^{*})=\arg\max_{\sigma_{i}^{2},{\bf\it\lambda}_{i}}&~{}\log~{}p(\mathbf{r}_{i}|~{}\mathbf{D},\sigma_{i}^{2},{\bf\it\lambda}_{i})=\arg\max_{\sigma_{i}^{2},{\bf\it\lambda}_{i}}~{}\log\int
p(\mathbf{r}_{i}|~{}\mathbf{D},\hat{\mathbf{o}}_{i},\sigma_{i}^{2})~{}p(\hat{\mathbf{o}}_{i}|~{}\sigma_{i}^{2},{\bf\it\lambda}_{i})~{}\differential\hat{\mathbf{o}}_{i}\\\
\overset{\text{from~{}\eqref{eq:post}}}{=}\arg\max_{\sigma_{i}^{2},{\bf\it\lambda}_{i}}&~{}\log~{}\mathcal{N}\left(\mathbf{r}_{i}|~{}\mathbf{D}{\bf\it\beta}_{i},\sigma_{i}^{2}\left(\mathbf{D}\mathrm{diag}({\bf\it\lambda}_{i})^{-1}\mathbf{D}^{\mathsf{T}}+\mathbf{I}_{k}\right)\right)\\\
=\arg\max_{\sigma_{i}^{2},{\bf\it\lambda}_{i}}&~{}-\frac{1}{2\sigma_{i}^{2}}\|\mathbf{r}_{i}-\mathbf{D}{\bf\it\mu}_{i}(\mathbf{r}_{i},\mathbf{D},{\bf\it\lambda}_{i})\|_{2}^{2}-\frac{1}{2\sigma_{i}^{2}}\left\|\mathrm{diag}\left({\bf\it\lambda}_{i}\right)^{1/2}\delta{\bf\it\mu}_{i}(\mathbf{r}_{i},\mathbf{D},{\bf\it\lambda}_{i})\right\|_{2}^{2}\\\
&~{}-\frac{1}{2}\log~{}\left|\mathrm{diag}({\bf\it\lambda}_{i})+\mathbf{D}^{\mathsf{T}}\mathbf{D}\right|+\frac{1}{2}\log({\bf\it\lambda}_{i})^{\mathsf{T}}\mathbf{1}_{d(r,m)}-\frac{k}{2}\log(\sigma_{i}^{2})-\frac{k}{2}\log(2\pi)\,.\end{split}$
(3.6)

Thus the final posterior distribution
$p(\hat{\mathbf{o}}_{i}|~{}\mathbf{D},\mathbf{r}_{i})$ is approximated as

$p(\hat{\mathbf{o}}_{i}|~{}\mathbf{D},\mathbf{r}_{i})=p(\hat{\mathbf{o}}_{i}|~{}\mathbf{D},\mathbf{r}_{i},\sigma_{i}^{*2},{\bf\it\lambda}^{*}_{i})=\mathcal{N}\left(\hat{\mathbf{o}}_{i}|~{}{\bf\it\mu}_{i}(\mathbf{r}_{i},\mathbf{D},{\bf\it\lambda}^{*}_{i}),{\bf\it\Sigma}_{i}(\mathbf{D},\sigma_{i}^{*2},{\bf\it\lambda}^{*}_{i})\right)\,.$
(3.7)

Furthermore, taking partial derivative of the objective function in (3.6) with
respect to $\sigma_{i}^{2}$ yields that the maximum is achieved at

$\sigma_{i}^{*2}=\frac{\|\mathbf{r}_{i}-\mathbf{D}{\bf\it\mu}_{i}\|_{2}^{2}+\left\|\mathrm{diag}\left({\bf\it\lambda}_{i}^{*}\right)^{1/2}\delta{\bf\it\mu}_{i}\right\|_{2}^{2}}{k}\,.$
(3.8)

Remark 2: A pragmatic choice is to set ${\bf\it\beta}_{i}=\mathbf{0}$ and thus
$\delta{\bf\it\mu}_{i}={\bf\it\mu}_{i}$. When we assume that the prior
variance values are very large (${\bf\it\lambda}_{i}\to\mathbf{0}$), i.e., the
prior of $\hat{\mathbf{o}}_{i}$ is uninformative, it can be obtained that
${\bf\it\mu}_{i}=\tilde{\mathbf{o}}_{i}=\mathbf{D}^{+}\mathbf{r}_{i}$,
${\bf\it\Sigma}_{i}=\sigma_{i}^{*2}(\mathbf{D}^{\mathsf{T}}\mathbf{D})^{-1}$,
and
$\sigma_{i}^{*2}=\|\mathbf{r}_{i}-\mathbf{D}\tilde{\mathbf{o}}_{i}\|_{2}^{2}/k$.
In this case, the posterior mean vector of $\hat{\mathbf{o}}_{i}$ coincides
with the original least squares estimate $\tilde{\mathbf{o}}_{i}$ as in
deterministic OpInf, and $\sigma_{i}^{*2}$ is the mean squared error of the
linear regression.

Remark 3: If the pre-defined prior mean vector of $\hat{\mathbf{o}}_{i}$ is
taken as its point estimate $\tilde{\mathbf{o}}_{i}$, i.e.,
${\bf\it\beta}_{i}=\tilde{\mathbf{o}}_{i}=\mathbf{D}^{+}\mathbf{r}_{i}$, then
we have $\delta{\bf\it\mu}_{i}=\mathbf{0}$ and
${\bf\it\mu}_{i}=\tilde{\mathbf{o}}_{i}$, as well as
$\sigma_{i}^{*2}=\|\mathbf{r}_{i}-\mathbf{D}\tilde{\mathbf{o}}_{i}\|_{2}^{2}/k$.
Moreover, by simplifying (3.6) the estimate of ${\bf\it\lambda}_{i}$ is
formulated as

$\begin{split}{\bf\it\lambda}_{i}^{*}&=\arg\max_{{\bf\it\lambda}_{i}}~{}-\frac{1}{2}\log~{}|\mathrm{diag}({\bf\it\lambda}_{i})+\mathbf{D}^{\mathsf{T}}\mathbf{D}|+\frac{1}{2}\log~{}|\mathrm{diag}({\bf\it\lambda}_{i})|\\\
&=\arg\max_{{\bf\it\lambda}_{i}}~{}-\frac{1}{2}\log~{}\left|\mathbf{I}+\mathrm{diag}({\bf\it\lambda}_{i})^{-\frac{1}{2}}\mathbf{D}^{\mathsf{T}}\mathbf{D}\mathrm{diag}({\bf\it\lambda}_{i})^{-\frac{1}{2}}\right|\,.\end{split}$
(3.9)

Since
$\mathrm{diag}({\bf\it\lambda}_{i})^{-\frac{1}{2}}\mathbf{D}^{\mathsf{T}}\mathbf{D}\mathrm{diag}({\bf\it\lambda}_{i})^{-\frac{1}{2}}$
is semi-positive definite, one has
$\min\left|\mathbf{I}+\mathrm{diag}({\bf\it\lambda}_{i})^{-\frac{1}{2}}\mathbf{D}^{\mathsf{T}}\mathbf{D}\mathrm{diag}({\bf\it\lambda}_{i})^{-\frac{1}{2}}\right|=1$
when $\mathrm{diag}({\bf\it\lambda}_{i})^{-1}=\mathbf{0}$
(${\bf\it\lambda}_{i}\to\infty^{d(r,m)}$), i.e., when the prior of
$\hat{\mathbf{o}}_{i}$ is deterministic. In this case, the Bayesian
formulation reduces to deterministic OpInf, indicating that the deterministic
version essentially forms an optimal solution, provided that the Gram matrix
$\mathbf{D}^{\mathsf{T}}\mathbf{D}$ is well-conditioned.

###  3.2 A probabilistic reduced-order model

Because the OpInf regression naturally decouples along the rows of the reduced
operator matrix
$\hat{\mathbf{O}}=[~{}\hat{\mathbf{A}}~{}~{}\hat{\mathbf{H}}~{}~{}\hat{\mathbf{B}}~{}~{}\hat{\mathbf{c}}~{}]$,
we likewise treat the Bayesian inference of each row $\hat{\mathbf{o}}_{i}$
independently, $1\leq i\leq r$. Hence the inference of $\hat{\mathbf{O}}$ is
given by the posterior distribution

$p(\hat{\mathbf{O}}|~{}\mathbf{D},\mathbf{R})=\prod_{i=1}^{r}p(\hat{\mathbf{o}}_{i}|~{}\mathbf{D},\mathbf{r}_{i})\,.$
(3.10)

Using such an inferred operator $\hat{\mathbf{O}}$, we solve for the reduced
state $\hat{\mathbf{q}}(t)$ through the reduced-order system (2.2) over
$[t_{0},t_{f}]$. We may therefore view $\hat{\mathbf{q}}(t)$ as a stochastic
process because it depends on $\hat{\mathbf{O}}$, with distribution

$p(\hat{\mathbf{q}}(t)|~{}\mathbf{D},\mathbf{R})=\int
p(\hat{\mathbf{q}}(t)|~{}\hat{\mathbf{O}})p(\hat{\mathbf{O}}|~{}\mathbf{D},\mathbf{R})~{}\differential\hat{\mathbf{O}}=\int
p(\hat{\mathbf{q}}(t)|~{}\hat{\mathbf{o}}_{1},\ldots,\hat{\mathbf{o}}_{r})\prod_{i=1}^{r}p(\hat{\mathbf{o}}_{i}|~{}\mathbf{D},\mathbf{r}_{i})~{}\differential\hat{\mathbf{o}}_{i}\,,$
(3.11)

in which all the $\hat{\mathbf{o}}_{i}$’s are marginalized. Because (2.2) is
computationally inexpensive to integrate in time, we may use a Monte Carlo
sampling over the operator posterior
$p(\hat{\mathbf{O}}|~{}\mathbf{D},\mathbf{R})$ to estimate the mean function
and second-order moments of $\hat{\mathbf{q}}(t)$ (see Figure 1). In fact,
$p(\hat{\mathbf{q}}(t)|\,\hat{\mathbf{O}})$ represents a deterministic model
given by the reduced operators $\hat{\mathbf{O}}$, hence, one can write the
reduced state solution $\hat{\mathbf{q}}(t)$ determined by the reduced
operators $\hat{\mathbf{O}}$ as a function of $\hat{\mathbf{O}}$, i.e.,
$\hat{\mathbf{q}}(t;\hat{\mathbf{O}})$. Note that the solution mean and the
solution determined by the posterior operator mean are not necessarily equal
to each other, i.e.,
$\mathbb{E}[\hat{\mathbf{q}}(t;\hat{\mathbf{O}}|~{}\mathbf{D},\mathbf{R})]\neq\hat{\mathbf{q}}(t;\mathbb{E}[\hat{\mathbf{O}}|~{}\mathbf{D},\mathbf{R}])$.
We will show, however, that in practice the discrepancy between these two
quantities is often small, even in complex scenarios.

Algorithm 1 details the BayesOpInf procedure with
$\boldsymbol{\beta}_{i}=\mathbf{0}$,
$\delta\boldsymbol{\mu}_{i}=\boldsymbol{\mu}_{i}$ (see Remark 2), and a fixed
choice of the prior variance vectors $\boldsymbol{\lambda}_{i}$,
$i=1,\ldots,r$. The mean vectors $\boldsymbol{\mu}_{i}$ and covariance
matrices $\boldsymbol{\Sigma}_{i}$ are determined via (3.5) and (3.8). We
discuss the selection of the $\boldsymbol{\lambda}_{i}$ in Section 4.

Figure 1: The BayesOpInf workflow. To derive the probabilistic ROM, infer the
multivariate normal posterior distribution for the operator matrix
$\hat{\mathbf{O}}$ from snapshot and time derivative data. This is a one-time,
offline cost. To evaluate the probabilistic ROM in the online stage, draw
samples from the posterior distribution of $\hat{\mathbf{O}}$ and solve the
corresponding ROM realizations, then calculate the statistics of these ROM
solutions.

1:procedure BayesOpInf( Data matrix $\mathbf{D}\in\mathbb{R}^{k\times
d(r,m)}$, projected time derivatives
$\mathbf{R}=[~{}\mathbf{r}_{1}~{}\cdots~{}\mathbf{r}_{r}~{}]^{\mathsf{T}}\in\mathbb{R}^{r\times
k}$, regularization parameters
$\boldsymbol{\lambda}_{1},\ldots,\boldsymbol{\lambda}_{r}\in\mathbb{R}^{d(r,m)}$,
projected initial condition $\hat{\mathbf{q}}_{0}\in\mathbb{R}^{r}$, final
time $t_{f}>t_{0}$, number of Monte Carlo samples $N\in\mathbb{N}$ )

2: # Offline phase: construct operator posterior
$p(\hat{\mathbf{O}}|\>\mathbf{D},\mathbf{R})$.

3: for $i=1,\ldots,r$ do

4:
$\boldsymbol{\mu}_{i}\leftarrow\arg\min_{\boldsymbol{\mu}_{i}}\left\|\mathbf{D}\boldsymbol{\mu}_{i}-\mathbf{r}_{i}\right\|_{2}^{2}+\|\text{diag}(\boldsymbol{\lambda}_{i})^{1/2}\boldsymbol{\mu}_{i}\|_{2}^{2}$
# Mean of operator row $\hat{\mathbf{o}}_{i}$.

5:
$\sigma_{i}^{*2}\leftarrow\frac{1}{k}\left(\left\|\mathbf{D}\boldsymbol{\mu}_{i}-\mathbf{r}_{i}\right\|_{2}^{2}+\|\text{diag}(\boldsymbol{\lambda}_{i})^{1/2}\boldsymbol{\mu}_{i}\|_{2}^{2}\right)$

6:
$\boldsymbol{\Sigma}_{i}\leftarrow\sigma_{i}^{*2}\left(\mathbf{D}^{\mathsf{T}}\mathbf{D}+\textup{diag}(\boldsymbol{\lambda}_{i})\right)^{-1}$
# Covariance of operator row $\hat{\mathbf{o}}_{i}$.

7: # Online phase: sample from operator posterior
$p(\hat{\mathbf{O}}|\>\mathbf{D},\mathbf{R})$ and obtain ROM solutions.

8: for $\ell=1,\ldots,N$ do

9: for $i=1,\ldots,r$ do

10: $\hat{\mathbf{o}}_{i}^{(\ell)}\leftarrow~{}\textrm{sample
from}~{}\mathcal{N}(\boldsymbol{\mu}_{i},\boldsymbol{\Sigma}_{i})$

11: $\tilde{\mathbf{Q}}^{(\ell)}\leftarrow~{}\textrm{integrate
\eqref{eq:reducedorder}
with}~{}\hat{\mathbf{O}}=[~{}\hat{\mathbf{A}}~{}\hat{\mathbf{H}}~{}\hat{\mathbf{B}}~{}\hat{\mathbf{c}}~{}]=[~{}\hat{\mathbf{o}}_{1}^{(\ell)}~{}\cdots~{}\hat{\mathbf{o}}_{r}^{(\ell)}~{}]^{\mathsf{T}}$
from $\hat{\mathbf{q}}_{0}=\hat{\mathbf{Q}}_{:,0}$ over $[t_{0},t_{f}]$

12: return $\tilde{\mathbf{Q}}^{(1)},\ldots,\tilde{\mathbf{Q}}^{(N)}$ # ROM
solutions corresponding to the operator samples.

Algorithm 1 Bayesian Operator Inference (BayesOpInf).

Remark 4: As the posterior
$p(\hat{\mathbf{o}}_{i}|~{}\mathbf{D},\mathbf{r}_{i})$ follows a multivariate
normal distribution, the noise-free posterior prediction of
$\frac{\differential}{\differential t}\hat{q}_{i}(t)$ is a finite-rank
Gaussian process, where $\hat{q}_{i}(t)$ is the $i$th entry of
$\hat{\mathbf{q}}(t)$. Defined by the ROM form (2.2), the posterior mean and
covariance functions of $\frac{\differential}{\differential t}\hat{q}_{i}(t)$
take the form of

$\begin{split}\mathbb{E}\left[\frac{\differential}{\differential
t}\hat{q}_{i}(t)\middle|~{}\mathbf{D},\mathbf{r}_{i}\right]&=\mathbf{d}(\hat{\mathbf{q}}(t),\mathbf{u}(t))^{\mathsf{T}}{\bf\it\mu}_{i}\,,\quad\text{and}\\\
\mathbb{C}ov\left[\frac{\differential}{\differential
t}\hat{q}_{i}(t),\frac{\differential}{\differential
t}\hat{q}_{i}(t^{\prime})\middle|~{}\mathbf{D},\mathbf{r}_{i}\right]&=\mathbf{d}(\hat{\mathbf{q}}(t),\mathbf{u}(t))^{\mathsf{T}}{\bf\it\Sigma}_{i}\mathbf{d}(\hat{\mathbf{q}}(t^{\prime}),\mathbf{u}(t^{\prime}))\,,\end{split}$
(3.12)

respectively, where
$\mathbf{d}(\hat{\mathbf{q}},\mathbf{u}):=[~{}\hat{\mathbf{q}}^{\mathsf{T}}~{}~{}(\hat{\mathbf{q}}\>\otimes\>\hat{\mathbf{q}})^{\mathsf{T}}~{}~{}\mathbf{u}^{\mathsf{T}}~{}~{}1~{}]^{\mathsf{T}}\in\mathbb{R}^{d(r,m)}$.
The same result can be obtained using Gaussian process regression [54] with
the following noise-corrupted prior:

$\frac{\differential}{\differential
t}\hat{q}_{i}(t)\sim\mathcal{GP}\bigg{(}\mathbf{d}(\hat{\mathbf{q}}(t),\mathbf{u}(t))^{\mathsf{T}}{\bf\it\beta}_{i},\sigma_{i}^{2}\mathbf{d}(\hat{\mathbf{q}}(t),\mathbf{u}(t))^{\mathsf{T}}\mathrm{diag}({\bf\it\lambda}_{i})^{-1}\mathbf{d}(\hat{\mathbf{q}}(t^{\prime}),\mathbf{u}(t^{\prime}))+\sigma_{i}^{2}\delta_{\text{Dirac}}(|t-t^{\prime}|)\bigg{)}\,,$
(3.13)

in which $\delta_{\text{Dirac}}(\cdot)$ denotes the Dirac-delta function. This
prior Gaussian process is consistent with the likelihood and prior definitions
adopted in (3.2) and (3.3), respectively, for the Bayesian inference.

###  3.3 Equivalence to ridge regression

Recall the posterior mean vector of $\hat{\mathbf{o}}_{i}$ given in (3.5). The
correction term $\delta{\bf\it\mu}_{i}={\bf\it\mu}_{i}-{\bf\it\beta}_{i}$ can
be expressed alternatively as

$\delta{\bf\it\mu}_{i}=\arg\min_{{\bf\it\eta}\in\mathbb{R}^{d(r,m)}}\left\\{\|\mathbf{r}_{i}-\mathbf{D}({\bf\it\beta}_{i}+{\bf\it\eta})\|_{2}^{2}+\|\mathrm{diag}({\bf\it\lambda}_{i})^{\frac{1}{2}}{\bf\it\eta}\|_{2}^{2}\right\\}\,,$
(3.14)

implying that the posterior mean is equivalent to a point estimate through
such a regularized least squares problem, often referred to as a ridge
regression, which uses a Tikhonov regularization
$\|\mathrm{diag}({\bf\it\lambda}_{i})^{\frac{1}{2}}{\bf\it\eta}\|_{2}^{2}$.
The use of Tikhonov regularization in OpInf has been reported in [31, 48].
Thus another source of modeling uncertainty is introduced to the ROM by the
Tikhonov regularization defined with OpInf. Based on the equivalence between
the posterior mean of inferred operator and ridge regression, it can be
claimed that the uncertainty introduced by Tikhonov regularization is
quantified through the proposed Bayesian framework as well.

On the other hand, the selection of the regularization parameters
$\boldsymbol{\lambda}_{1},\ldots,\boldsymbol{\lambda}_{r}$ has been shown to
be important in addressing ill-conditioning, over-fitting, or model inadequacy
[31, 48]. The proposed Bayesian method in this work not only formulates a
framework of uncertainty quantification, but also provides a new perspective
for understanding the regularization selection process for deterministic
OpInf. We address this further in Section 4.

##  4 Regularization selection

Without regularization, the OpInf point estimate (2.6) is represented with the
Moore–Penrose inverse of matrix $\mathbf{D}$ and requires that $\mathbf{D}$
has full column rank, i.e., the Gram matrix
$\mathbf{D}^{\mathsf{T}}\mathbf{D}$ is invertible. In practice, however, the
condition number of $\mathbf{D}^{\mathsf{T}}\mathbf{D}$ may be large, which
can lead to a poor estimate of the reduced operators and compromise the
stability of the corresponding ROM. In this case, making use of a proper
regularization can be critical for stabilizing the OpInf scheme and improving
the robustness of the ROM [31, 48].

As specified in Section 3.3, the posterior mean of inferred operators is
equivalent to an estimator through ridge regression. The penalty coefficients
${\bf\it\lambda}_{i}$ in the Tikhonov regularization correspond to the prior
variances of the operators $\hat{\mathbf{o}}_{i}$. As given in (3.6), the
values of the hyperparameters ${\bf\it\lambda}_{i}$ can be determined by an
empirical Bayes method of maximizing marginal likelihood,

$\begin{split}\mathcal{L}_{i}(\sigma_{i}^{2},{\bf\it\lambda}_{i})=&~{}\log~{}\mathcal{N}\left(\mathbf{r}_{i}|~{}\mathbf{D}{\bf\it\beta}_{i},\sigma_{i}^{2}\left(\mathbf{D}\mathrm{diag}({\bf\it\lambda}_{i})^{-1}\mathbf{D}^{\mathsf{T}}+\mathbf{I}_{k}\right)\right)\\\
=&-\frac{1}{2\sigma_{i}^{2}}\|\mathbf{r}_{i}-\mathbf{D}{\bf\it\mu}_{i}\|_{2}^{2}-\frac{1}{2\sigma_{i}^{2}}\left\|\mathrm{diag}\left({\bf\it\lambda}_{i}\right)^{1/2}\delta{\bf\it\mu}_{i}\right\|_{2}^{2}\\\
&-\frac{1}{2}\log~{}\left|\mathrm{diag}({\bf\it\lambda}_{i})+\mathbf{D}^{\mathsf{T}}\mathbf{D}\right|+\frac{1}{2}\log({\bf\it\lambda}_{i})^{\mathsf{T}}\mathbf{1}_{d(r,m)}-\frac{k}{2}\log(\sigma_{i}^{2})-\frac{k}{2}\log(2\pi)\,.\end{split}$
(4.1)

In this sense, the empirical method provides a guidance for the determination
of penalty parameters ${\bf\it\lambda}_{i}$ in regularized OpInf.

Taking partial derivatives of $\mathcal{L}_{i}$ with respect to
$\sigma_{i}^{2}$ and all the entries of
${\bf\it\lambda}_{i}=[\lambda_{i,1},~{}\ldots,~{}\lambda_{i,d(r,m)}]^{\mathsf{T}}$,
the extremum conditions are

$\frac{\partial\mathcal{L}_{i}}{\partial\sigma_{i}^{2}}=0\quad\Longrightarrow\quad\sigma_{i}^{*2}=\frac{\|\mathbf{r}_{i}-\mathbf{D}{\bf\it\mu}_{i}\|_{2}^{2}+\left\|\mathrm{diag}\left({\bf\it\lambda}_{i}^{*}\right)^{1/2}\delta{\bf\it\mu}_{i}\right\|_{2}^{2}}{k}\,,$

which has already been given in (3.8), and

$\frac{\partial\mathcal{L}_{i}}{\partial\lambda_{i,j}}=0\quad\Longrightarrow\quad\lambda_{i,j}^{*}\left(\delta\mu_{i,j}^{2}+\Sigma_{i,jj}\right)=\sigma_{i}^{*2}\,,\quad
1\leq j\leq d(r,m)\,,$ (4.2)

in which $\delta\mu_{i,j}$ and $\Sigma_{i,jj}$ are the $j$-th entry of
$\delta{\bf\it\mu}_{i}$ and the $j$-th diagonal entry of ${\bf\it\Sigma}_{i}$,
respectively.111The evaluation of
$\frac{\partial\mathcal{L}_{i}}{\partial\lambda_{i,j}}$ uses the identity
$\partial_{\theta}(\log\det\mathbf{A}(\theta))=\tr(\mathbf{A}(\theta)^{-1}\partial_{\theta}\mathbf{A}(\theta))$,
in which $\mathbf{A}$ is an invertible matrix parametrized by $\theta$. Note
that (4.2) is an implicit representation of ${\bf\it\lambda}_{i}^{*}$, the
penalty coefficients of Tikhonov regularization.

For each reduced system mode $i=1,\ldots,r$, the problem of maximum marginal
likelihood in (4.1) is non-convex with $1+d(r,m)$ unknowns. Solving such a
problem using gradient-based methods is possible, but inconvenient for
moderately sized $r$ and not guaranteed to find a global optimum. To alleviate
complexity, we impose the constraint $\lambda_{i,j}=\lambda_{i}$ for all
$1\leq j\leq d(r,m)$, hence
$\boldsymbol{\lambda}_{i}=[~{}\lambda_{i}~{}\lambda_{i}~{}\cdots~{}\lambda_{i}~{}]^{\mathsf{T}}\in\mathbb{R}^{d(r,m)}$.
In this case, the regularization term in (2.4), as well as in steps 4 and 5 of
Algorithm 1, is an $L_{2}$-regularization with penalty coefficient
$\lambda_{i}>0$, i.e.,

$\displaystyle\mathcal{P}_{i}(\boldsymbol{\mu}_{i})=\mathcal{P}_{i}({\bf\it\beta}_{i}+\delta\boldsymbol{\mu}_{i})=\left\|\textrm{diag}(\boldsymbol{\lambda}_{i})^{1/2}\delta\boldsymbol{\mu}_{i}\right\|_{2}^{2}=\lambda_{i}\left\|\delta\boldsymbol{\mu}_{i}\right\|_{2}^{2}.$
(4.3)

This means that the operator entries corresponding to the reduced system mode
$\hat{q}_{i}$, i.e., the $i$th rows of $\hat{\mathbf{c}}$, $\hat{\mathbf{A}}$,
$\hat{\mathbf{H}}$, and $\hat{\mathbf{B}}$, are regularized equally, while
operator entries corresponding to different system modes are regularized
separately. Summing (4.2) over all $1\leq j\leq d(r,m)$ yields

$\lambda_{i}^{*}(\|\delta{\bf\it\mu}_{i}\|_{2}^{2}+\tr(\mathbf{\Sigma}_{i}))=d(r,m)\sigma_{i}^{*2}\,,$
(4.4)

in which $\tr(\mathbf{\Sigma}_{i})$ can alternatively be written as follows
[54]:

$\tr({\bf\it\Sigma}_{i})=\sigma_{i}^{*2}\sum_{l=1}^{d(r,m)}\frac{1}{\lambda_{i}^{*}+g_{l}}\,,$
(4.5)

where $g_{1},\ldots,g_{d(r,m)}\geq 0$ are the non-negative eigenvalues of the
Gram matrix $\mathbf{D}^{\mathsf{T}}\mathbf{D}$. This representation of
$\tr(\mathbf{\Sigma}_{i})$ has the computational advantage of avoiding
explicit matrix inversion in (3.5) to evaluate $\boldsymbol{\Sigma}_{i}$. We
then have

$\lambda^{*}_{i}=\frac{\sigma_{i}^{*2}}{\|\delta{\bf\it\mu}_{i}\|_{2}^{2}}\underbrace{\sum_{l=1}^{d(r,m)}\frac{g_{l}}{\lambda_{i}^{*}+g_{l}}}_{:=\gamma_{i}}=\frac{\gamma_{i}\sigma_{i}^{*2}}{\|\delta{\bf\it\mu}_{i}\|_{2}^{2}}\,,$
(4.6)

in which $0<\gamma_{i}<d(r,m)$. The following equality shows the connections
among the noise $\sigma_{i}^{*2}$ and the two terms in the loss function of
the ridge regression (3.14):

$\frac{\lambda_{i}^{*}\|\delta{\bf\it\mu}_{i}\|_{2}^{2}}{\gamma_{i}}=\frac{\|\mathbf{r}_{i}-\mathbf{D}{\bf\it\mu}_{i}\|_{2}^{2}}{k-\gamma_{i}}=\sigma_{i}^{*2}\,.$
(4.7)

Substituting (3.8) into (4.6), ${\lambda}_{i}^{*}$ is expressed by the
posterior mean $\delta{\bf\it\mu}_{i}$ and covariance matrix
${\bf\it\Sigma}_{i}$, both being dependent on ${\lambda}_{i}^{*}$ as well.
Thus an implicit representation of ${\lambda}_{i}^{*}$ is given by
${\lambda}_{i}^{*}=\mathcal{F}({\lambda}_{i}^{*})$, where $\mathcal{F}$
denotes the right-hand side of (4.6). In this work, we suggest a fixed-point
iterative strategy,
$({\lambda}_{i}^{*})_{s+1}=\mathcal{F}(({\lambda}_{i}^{*})_{s})$,
$s=0,1,\ldots$, as a practical approach to simultaneously satisfy (3.8) and
(4.6) and approximately maximize (4.1). The initial values
$\\{({\lambda}_{i}^{*})_{0}\\}_{i=1}^{r}$ can be set to a selection of
regularization coefficients that result in a stable ROM, but which are not
necessarily optimal. This iterative regularization update procedure is given
by Algorithm 2, which considers a zero prior mean of reduced operators, i.e.,
${\bf\it\beta}_{i}=\mathbf{0}$ and ${\bf\it\mu}_{i}=\delta{\bf\it\mu}_{i}$. We
note that while Algorithm 2 is computationally efficient, because it is a
completely data-driven selection it does not include stability information for
the reduced-order ODE system (2.2).

1:procedure BayesOpInfRegularization( Data matrix
$\mathbf{D}\in\mathbb{R}^{k\times d(r,m)}$, projected time derivatives
$\mathbf{R}=[~{}\mathbf{r}_{1}~{}\cdots~{}\mathbf{r}_{r}~{}]^{\mathsf{T}}\in\mathbb{R}^{r\times
k}$, initial regularization parameters $\lambda_{1},\ldots,\lambda_{r}\geq 0$,
convergence tolerance $\varepsilon>0$ )

2:
$g_{1},\ldots,g_{d(r,m)}\leftarrow\texttt{eig}(\mathbf{D}^{\mathsf{T}}\mathbf{D})$
# Get eigenvalues of the data Gram matrix.

3:
$\boldsymbol{\lambda}^{(0)}\leftarrow[~{}\lambda_{1}~{}\cdots~{}\lambda_{r}~{}]$

4: for $\ell=0,1,\ldots$ do

5: for $i=1,\ldots,r$ do

6:
$\boldsymbol{\mu}_{i}\leftarrow\arg\min_{\boldsymbol{\mu}_{i}}\left\|\mathbf{D}\boldsymbol{\mu}_{i}-\mathbf{r}_{i}\right\|_{2}^{2}+\lambda_{i}\left\|\boldsymbol{\mu}_{i}\right\|_{2}^{2}$
# Update operator mean.

7:
$\sigma_{i}^{*2}\leftarrow\frac{1}{k}\left(\left\|\mathbf{D}\boldsymbol{\mu}_{i}-\mathbf{r}_{i}\right\|_{2}^{2}+\lambda_{i}\left\|\boldsymbol{\mu}_{i}\right\|_{2}^{2}\right)$
# Update operator variance.

8:
$\lambda_{i}^{*}\leftarrow\displaystyle\frac{\sigma_{i}^{*2}}{\|\boldsymbol{\mu}_{i}\|_{2}^{2}}\displaystyle\sum_{l=1}^{d(r,m)}\frac{g_{l}}{\lambda_{i}+g_{l}}$
# Update regularization.

9:
$\boldsymbol{\lambda}^{(\ell+1)}\leftarrow[~{}\lambda_{1}^{*}~{}\cdots~{}\lambda_{r}^{*}~{}]$

10: if
$\|\boldsymbol{\lambda}^{(\ell+1)}-\boldsymbol{\lambda}^{(\ell)}\|_{2}\big{/}\|\boldsymbol{\lambda}^{(\ell)}\|_{2}<\varepsilon$
then # Check for convergence.

11: return $\lambda_{1}^{*},\ldots,\lambda_{r}^{*}$

12:
$\lambda_{1},\ldots,\lambda_{r}\leftarrow\lambda_{1}^{*},\ldots,\lambda_{r}^{*}$
# Continue iterating if not converged.

Algorithm 2 Iterative regularization update for BayesOpInf with
$\lambda_{i,j}=\lambda_{i}$.

##  5 Numerical examples

Section 5.1 demonstrates the proposed method on one-dimensional Euler
equations where the state observations are polluted by Gaussian noise. In
Section 5.2, we examine a combustion application in which the regularizing
terms play a crucial role for learning stable ROMs from data. The code for
these experiments can be found at https://github.com/Willcox-Research-
Group/ROM-OpInf-Combustion-2D/tree/bayes.

###  5.1 Noised Euler Equations

Consider the conservative one-dimensional compressible Euler equations for an
ideal gas,

$\displaystyle\frac{\partial}{\partial t}\left[\rho\right]$
$\displaystyle=-\frac{\partial}{\partial x}\left[\rho u\right],$
$\displaystyle\frac{\partial}{\partial t}\left[\rho u\right]$
$\displaystyle=-\frac{\partial}{\partial x}\left[\rho u^{2}+p\right],$
$\displaystyle\frac{\partial}{\partial t}\left[\rho e\right]$
$\displaystyle=-\frac{\partial}{\partial x}\left[(\rho e+p)u\right],$ (5.1)

where $u$ is velocity [m/s], $\rho$ is density [kg/m3], $p$ is pressure [Pa],
$\rho u$ is specific momentum [kg/m2s], and $\rho e$ is total energy [kgJ/m3].
The state variables are related via the ideal gas law $\rho
e=\frac{p}{\gamma-1}+\frac{1}{2}\rho u^{2}$, where $\gamma=1.4$ is the heat
capacity ratio. Our goal is to learn a reduced model of (5.1) from noisy
observations of the conservative variables $(\rho,\rho u,\rho e)$ and to
quantify the uncertainties induced by the noise. The ideal gas law enables a
change of variables $\tau:(\rho,\rho u,\rho e)\mapsto(u,p,1/\rho)$ that
transforms (5.1) into the purely quadratic system

$\displaystyle\frac{\partial u}{\partial t}$ $\displaystyle=-u\frac{\partial
u}{\partial x}-\zeta\frac{\partial p}{\partial x},$
$\displaystyle\frac{\partial p}{\partial t}$ $\displaystyle=-\gamma
p\frac{\partial u}{\partial x}-u\frac{\partial p}{\partial x},$
$\displaystyle\frac{\partial\zeta}{\partial t}$
$\displaystyle=-u\frac{\partial\zeta}{\partial x}+\zeta\frac{\partial
u}{\partial x},$ (5.2)

where $\zeta=1/\rho$ is the specific volume [m3/kg]. In the model (5.2), every
term is the product of exactly two state variables $(u,p,\zeta)$ or their
spatial derivatives, hence for this problem we seek a reduced model with
quadratic terms only:

$\displaystyle\frac{\textup{d}}{\textup{d}t}\hat{\mathbf{q}}(t)=\hat{\mathbf{H}}[\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)].$
(5.3)

The OpInf operator matrix of (2.5) in this case is
$\hat{\mathbf{O}}=\hat{\mathbf{H}}\in\mathbb{R}^{r\times d(r)}$, with column
dimension $d(r)=r(r+1)/2$.

Figure 2: One of the sixty-four training trajectories at various spatial
coordinates before noise corruption (left) and singular value decay of the
entire training set with various noise levels (right). Training data consists
of noised snapshots for $0\leq t<0.01$.

We consider the spatial domain $[0,2)$ and the time domain
$[t_{0},t_{f}]=[0,0.03]$, with periodic boundary conditions imposed at the
spatial boundaries. Training data is generated over the shorter time domain
$[0,0.01]$ for multiple initial conditions: the initial pressure is set to
$10^{5}$ Pa everywhere, and initial velocity and density profiles are
constructed via cubic spline interpolation by fixing values at the spatial
points $x\in\\{0,\frac{2}{3},\frac{4}{3}\\}$ to either $95$ or $105$ m/s for
velocity and either $20$ or $24$ kg/m3 for density. We therefore have
$2^{3\times 2}=64$ initial conditions (three interpolation nodes with two
choices of values for two variables). For each initial condition, we solve the
conservative system (5.1) by discretizing the spatial domain via first-order
finite differences with spatial resolution $\delta x=0.01$ and stepping
forward in time via the first-order explicit Euler iteration with time step
$\delta t=10^{-5}$. The result is $k=64\times 1000=64{,}000$ total training
snapshots, each of size $n=200\times 3=600$. See Figure 2 for an example
trajectory in the training set. Before processing the snapshot data for
learning, each snapshot entry is polluted with Gaussian noise proportional to
the range of the associated variable to represent observation error. For
example, to each snapshot entry corresponding to density we add error drawn
from $\mathcal{N}(0,\varsigma^{2})$ where
$\varsigma=\xi(\max_{i,j}\\{\rho_{i,j}\\}-\min_{i,j}\\{\rho_{i,j}\\})$,
$\xi=5\%$ is the noise level, and $\rho_{i,j}$ is the $i$th entry of the $j$th
density training snapshot. The noisy snapshots are then transformed from the
conservative states of (5.1) to the states of the quadratic system (5.2) via
the mapping $\tau$ and non-dimensionalized, yielding training data
$\mathbf{q}_{0},\ldots,\mathbf{q}_{k-1}\in\mathbb{R}^{n}$ to be used in OpInf.
Note that the noise in the lifted variables is not necessarily Gaussian since
$\tau$ is a nonlinear mapping.

Figure 3: The first several modes of the (noisy) projected training data,
normalized to the range $[-1,1]$, for the training trajectory shown in Figure
2. The POD projection successfully filters much, but not all, of the noise on
the first $r=9$ modes.

From the lifted, noisy snapshots we compute a POD basis
$\mathbf{V}\in\mathbb{R}^{n\times r}$ and project the training data to the
corresponding low-dimensional linear subspace as
$\hat{\mathbf{q}}_{j}=\mathbf{V}^{\mathsf{T}}\mathbf{q}_{j}$. The dimension
$r=9$ is selected based on where the decay of the singular values flattens out
due to the noise pollution as seen in Figure 2. It has been widely observed
that POD has a smoothing effect (see, e.g., [11, 12, 52]). The projected
training data, an example of which is shown in Figure 3, exhibits relatively
little noise when compared to the observed training data (see also Figure 4)
up to $r=9$ modes. However, the OpInf regression requires an estimate of the
time derivatives $\dot{\mathbf{q}}_{j}$, which will be highly inaccurate if
computed as finite differences of the noisy projected data. We therefore apply
the noise-resistant strategy of [10], which uses implicit local polynomial
regression to estimate derivatives with minimal variance, to estimate each
$\dot{\mathbf{q}}_{j}$. Therefore, uncertainty in the OpInf regression stems
from _i_) noise in the observed data, _ii_) the truncated POD representation
(model form error), and _iii_) error in the estimates of the time derivatives
of the observed data.

Figure 4: Traces of full-order and ROM solutions to the Euler problem at three
spatial locations, with the same initial condition as in Figure 2. The noisy
training observations (noise level $\xi=5\%$) is also shown as individual
points over the training interval. The BayesOpInf ROM uses $r=9$ modes, with
the regularization chosen via Algorithm 2, and 100 draws are made from the
posterior distribution to sample the mean and standard deviations.

The noise in the training data necessitates some regularization in order to
achieve stable reduced models, though in this example among stable models the
learned OpInf ROM is not highly sensitive to the choice of regularization. We
initially set $\lambda_{i,j}=50$ for all $i,j$ and use Algorithm 2 to iterate
until the relative change (as measured in step 10) drops below
$\varepsilon=0.1\%$, which occurs after only nine iterations. We then draw
$100$ samples of $\hat{\mathbf{O}}$ from the posterior operator distribution
defined in (3.7) and integrate the corresponding ROMs over the full time
domain $[t_{0},t_{f}]$ for each of the (noisy) initial conditions in the
training set. Figure 4 shows the sample mean of these draws and three standard
deviations from the sample mean at several spatial locations for the
trajectory shown in Figure 2; the results are representative for the set of
all initial conditions considered. Note that the width of three standard
deviations, which slightly increases with time, is relatively small compared
to the observed data. This is because the uncertainty band shows a credible
interval, which quantifies the accuracy of model fits, i.e., the learned
posterior distribution describes the uncertainty relative to the true
solution. Note that, over the training interval $[0,0.01]$, the credible
interval is by definition not taking into account the measurement noise in the
snapshot data222The band that takes both the propagated uncertainty from the
posterior model and the measurement errors (noise) into account is often
referred to as a prediction interval, see chapter 9 of [47]..

Figure 5: Average relative error over the set of 64 initial conditions for
increasing basis sizes, over the training regime (left) and the prediction
regime (right), for ROMs where the regularization is chosen via Algorithm 2
(BayesOpInf) and by Algorithm 3 (OpInf).

The regularization selection strategy of Algorithm 2 only approximates the
maximizers of the marginal likelihood (4.1). On the other hand, the cost of
each iteration in Algorithm 2 is dominated only by the OpInf regression
(determining each $\boldsymbol{\mu}_{i}$ and $\boldsymbol{\Sigma}_{i}$). To
compare this strategy to existing methods, Figure 5 repeats the experiments
described above with basis sizes $r=3,4,\ldots,20$ using Algorithm 2 to select
the regularization, and computes the average relative error of the ROM defined
by the posterior operator mean $\boldsymbol{\mu}_{i}$ over the set of all
initial conditions for both the training regime $t\in[0,0.01)$ (with the error
computed with respect to the true trajectories, not the noisy observations)
and the prediction regime $t\in[0.01,0.03]$. The same results are shown for
the more involved regularization selection strategy of [31] detailed in
Algorithm 3 (see Appendix 2). This strategy chooses a scalar regularization
value (i.e., $\lambda_{i,j}=\lambda$ for all $i,j$) by solving an optimization
problem based on two criteria: _i_) that the resulting ROM remains stable
throughout the full time horizon $t\in[0,0.03]$ with each of the considered
initial conditions, and _ii_) that the reconstruction error of the training
data is minimized. Note that Algorithm 3 chooses the scalar regularization
that minimizes the difference between the integrated ROM state and the _noisy_
training data. Figure 5 shows that the errors for each selection method are
highly similar, especially over the training regime, but that the BayesOpInf
approach is slightly more accurate in the prediction regime for $r>11$. This
is likley because Algorithm 2 treats the regularization of each dynamic mode
$\hat{q}_{i}$ separately, while Algorithm 3 uses a single global regularizer
for all terms of the learned operators. The key takeaway is that when the
stability of the learned model is not highly sensitive to the regularization
hyperparameters, Algorithm 2 provides a lightweight option for selecting the
regularization. The next section details an example where the choice of
regularization is critical for model stability, in which case Algorithm 2 is a
less appropriate option for choosing the $\boldsymbol{\lambda}_{i}$.

###  5.2 Single-injector Combustion Process

We now consider a combustion application on a two-dimensional domain in which
model form error, and hence the choice of regularization, plays a significant
role in the OpInf procedure. This problem presents a challenging application
for model reduction and has been used as a test problem for both intrusive
methods [22, 23, 24] and non-intrusive OpInf [25, 31, 48]. The governing
dynamics are the nonlinear conservation laws

$\displaystyle\frac{\partial\vec{q}_{\textrm{c}}}{\partial
t}+\nabla\cdot(\vec{K}-\vec{K}_{\textrm{v}})$ $\displaystyle=\vec{S},$ (5.4)

where $\vec{q}_{\textrm{c}}=(\rho,\rho v_{x},\rho v_{y},\rho e,\rho Y_{1},\rho
Y_{2},\rho Y_{3},\rho Y_{4})$ are the conservative variables, $\vec{K}$ is the
inviscid flux, $\vec{K}_{\textrm{v}}$ is the viscous flux, and $\vec{S}$
contains source terms. Here $\rho$ is the density
$[\frac{\textrm{kg}}{\textrm{m}^{3}}]$, $v_{x}$ and $v_{y}$ are the $x$ and
$y$ velocity $[\frac{\textrm{m}}{\textrm{s}}]$, $e$ is the total energy
$[\frac{\textrm{J}}{\textrm{m}^{3}}]$, and $Y_{1}$, $Y_{2}$, $Y_{3}$, $Y_{4}$
are the mass fractions for the chemical species CH4, O2, H2O, and CO2,
respectively, which follow the global one-step reaction
$\textrm{CH}_{4}+2\textrm{O}_{2}\to\textrm{CO}_{2}+2\textrm{H}_{2}\textrm{O}$
[53]. The top and bottom walls of the combustor are assigned no-slip boundary
conditions, the upstream boundary is constrained to have a constant mass flow
at the inlets, and at the downstream end a non-reflecting boundary condition
is imposed to maintain the chamber pressure via

$\displaystyle p_{\textrm{back}}(t)$
$\displaystyle=p_{\textrm{back,ref}}\left(1+0.1\sin(2\pi ft)\right)$ (5.5)

where $p_{\textrm{back,ref}}=10^{6}$ Pa and $f=5{,}000$ Hz. See [19] for more
details on the governing equations.

Figure 6: Two-dimensional domain of the single-injector combustor. The
solution is monitored pointwise at sensor locations $1$–$4$ and reported in
Figures 9–11.

As observed in [48], the governing conservative equations (5.4) take on a
nearly quadratic form when written in a specific volume form with variables
$\vec{q}=(p,v_{x},v_{y},T,\zeta,c_{1},c_{2},c_{3},c_{4}),$ where $p$ is the
pressure $[\textrm{Pa}]$, $T$ is the temperature $[\textrm{K}]$,
$\zeta=1/\rho$ is the specific volume $[\frac{\textrm{m}^{3}}{\textrm{kg}}]$,
and $c_{1},\ldots,c_{4}$ are the species molar concentrations
$[\frac{\textrm{kmol}}{\textrm{m}^{3}}]$ given by $c_{\ell}=\rho
Y_{\ell}/M_{\ell}$ with $M_{\ell}$ the molar mass of the $\ell$th species
$[\frac{\textrm{g}}{\textrm{mol}}]$. In [31], it was shown that including
specific volume, pressure, and temperature (although only two quantities are
needed to fully define the high-fidelity model, the equation of state then
defining the third quantity) led to improved ROM performance. Hence, we
construct a fully quadratic ROM of the form (2.2), with the understanding that
this structure only approximates the structure of the governing equations.

In this experiment, we consider a single set of initial conditions and the
corresponding trajectory. Our objective is to predict in time beyond the
available training data. We solve (5.4) via the finite-volume based General
Equation and Mesh Solver (GEMS), see [19] for details. The spatial domain is
discretized with $n_{x}=38{,}523$ cells (hence each snapshot has
$8n_{x}=308{,}184$ entries) and the solution is computed for $50{,}000$ time
steps with step size $\delta t=10^{-7}$ s, from $t_{0}=0.015$ s to
$t_{f}=0.020$ s. We retain the first $k=20{,}000$ snapshots for training and
reserve the remaining $30{,}000$ for testing. The training snapshots are
lifted to the specific volume variables $\vec{q}$, then scaled to $[-1,1]$ by
first subtracting out the mean profile in the temperature, pressure, and
specific volume, then normalizing each variable by its maximum absolute entry.
These processed training snapshots $\mathbf{q}_{0},\ldots,\mathbf{q}_{k-1}$
have dimension $n=9n_{x}=346{,}000$. We compute the POD basis $\mathbf{V}$ of
the training snapshots using a randomized singular-value decomposition [18]
and retain $r=38$ modes based on the associated singular value decay. As the
training data is smooth (i.e., not corrupted by observational noise), we
estimate time derivatives of the projected training data
$\hat{\mathbf{q}}_{j}=\mathbf{V}^{\mathsf{T}}\mathbf{q}_{j}$ with fourth-order
finite differences. This estimation introduces a further source of error,
though the predominant error is due to the model form discrepancy from
assuming a quadratic from in the ROM.

Figure 7: Evaluations of the first six modes (i.e., the first six entries of
$\hat{\mathbf{q}}(t)$) of $100$ draws from the posterior distribution of the
Bayesian OpInf ROM. The dashed lines show the trajectory of the ROM with
posterior mean operators ($\hat{\mathbf{o}}_{i}=\boldsymbol{\mu}_{i}$); the
shaded areas are the regions within three standard deviations of the sample
mean. Training data is observed for $t\in[0.015,0.017)$.

Previous applications of OpInf to this problem have shown that regularization
is essential to obtain stable ROMs [31, 48]. We therefore utilize the
regularization selection strategy of [31] (Algorithm 3) to choose the
regularization parameters $\lambda_{i,j}$. In this context, we select two
scalars: one to penalize the terms of the non-quadratic operators
$\hat{\mathbf{c}}$, $\hat{\mathbf{A}}$, and $\hat{\mathbf{B}}$, and one for
$\hat{\mathbf{H}}$ alone (see Appendix 2 for details). Thus the posterior mean
$\boldsymbol{\mu}_{i}$ of the operator $\hat{\mathbf{o}}_{i}$ given in (3.7)
is equivalent to the ROM inference studied in [31]. We additionally draw $100$
operator samples from (3.7) and integrate the resulting ROMs to quantify the
uncertainties in the learned model. The computational cost of integrating
$100$ ROMs is trivial compared to the high-dimensional code: computing the
$50{,}000$ unprocessed snapshots with GEMS costs approximately $1{,}000$ CPU
hours, while constructing and integrating a single ROM from the posterior
takes $\sim 0.3$ CPU seconds. Even when the mean operators
$\boldsymbol{\mu}_{i}$ of the posterior distribution (3.7) define a stable
ROM, there is no guarantee that operators drawn from the distribution will
result in stable ROMs. In this experiment, $96$ of the $100$ posterior draws
yield stable ROMs. This is an indicator of the complexity of the problem at
hand, but it highlights the constraint-free nature of the posterior
formulation.

Figure 7 plots the dominant modes of the reduced-order state as a function of
time together with the deviation defined by the posterior draws. Both the
trajectory corresponding to the mean operators $\boldsymbol{\mu}_{i}$ and the
sample mean of the $96$ stable draws are shown, but they differ very little.
We emphasize that, while the entries of the model operators are assumed to
follow a Gaussian distribution, there is no reason for the integrated
trajectories of the corresponding quadratic model to be Gaussian in nature.
Notably, the uncertainty bands remain tight until the end of the training
regime, and widen increasingly with time. The band is thinnest in the second
mode, which corresponds to the pressure signal imparted by (5.5).

Figure 8: Full-order GEMS data and Bayesian OpInf ROM reconstructions of the
pressure at the four sensor locations marked in Figure 6. The shaded areas are
the regions within three standard deviations of the sample mean, using the
same $100$ draws as in Figure 7. Training data is observed for
$t\in[0.015,0.017)$.

Figure 9: Full-order GEMS data and Bayesian OpInf ROM reconstructions of the
horizontal velocity at the four sensor locations marked in Figure 6. The
shaded areas are the regions within three standard deviations of the sample
mean, using the same $100$ draws as in Figure 7. Training data is observed for
$t\in[0.015,0.017)$.

Figures 9 and 9 show the pressure and $x$-velocity, respectively, of the
reconstructed ROM trajectory at the sensor locations, together with the data
from the high-dimensional solver. We examine these variables in particular
because they exhibit low projection error in the prediction regime with
respect to the underlying POD basis $\mathbf{V}$. The full-order trajectory
computed by GEMS is also shown for comparison. As in Figure 7, the uncertainty
bands remain quite tight in the training regime but start to noticeably
increase in width shortly after the transition to the prediction regime. The
bands are widest at the peaks and troughs of the pressure signal, correctly
reflecting that the model more accurately captures the frequency than the
amplitude. While the sample mean is sometimes point-wise inaccurate, the true
solution is generally well captured within three standard deviations of the
sample mean. Figures 11 and 11 show this more precisely for pressure and
$x$-velocity, respectively, by plotting the absolute error and three sample
standard deviation (half the width of the uncertainty bands in Figures 9 and
9) as a function of time at the sensor locations. Sharp drops in the sample
mean error indicate that the sample mean crosses the true values. The width of
the uncertainty band generally exceeds the absolute sample mean error in the
prediction regime, especially for $t>0.018$ s. The sample mean error and the
standard deviation are positively correlated in time, i.e., the uncertainty in
the prediction generally increases with the error in the model.

Figure 10: Absolute error of the sample mean and three sample standard
deviations for the pressure at the locations in Figure 6, computed from $100$
draws from the posterior distribution of the Bayesian OpInf ROM. Note the
logarithmic scale of the $y$-axis.

Figure 11: Absolute error of the sample mean and three sample standard
deviations for the horizontal velocity at the locations in Figure 6, computed
from $100$ draws from the posterior distribution of the Bayesian OpInf ROM.
Note the logarithmic scale of the $y$-axis.

##  6 Concluding remarks

The proposed BayesOpInf method non-intrusively propagates uncertainties from
training data noise and model form error to learned low-dimensional operators
of a ROM. Linear Bayesian inference with Gaussian priors, driven by projected
state observations, leads to Gaussian posterior distributions for the reduced
operators. Because constructing the posterior distributions for the reduced
operators and integrating particular realizations of the ROM are relatively
inexpensive operations, we can efficiently equip the resulting ROM solutions
with uncertainty bands through Monte Carlo sampling. Such a Bayesian approach
has an inherent connection to Gaussian processes (GPs). A further discussion
on this fact is provided in Appendix 1, wherein a GP surrogate model is used
for the likelihood definition in the Bayesian inference. One open question is
whether additional known information about the problem structure can be used
to more strongly inform the Bayesian inference, for example by imposing a non-
Gaussian distribution to model the reduced-order residual.

BayesOpInf is a natural generalization of the deterministic OpInf framework.
Since the proposed Bayesian inference is generic and carried out only using
the solution snapshots and their time-derivative data, BayesOpInf is a
promising tool for more specialized applications of OpInf, such as localized
methods [14], problems with detailed structure [3], or parametric settings
[32]. Furthermore, because the posterior mean of BayesOpInf is equivalent to
Tikhonov-regularized OpInf, the approach is applicable for any Tikhonov-type
regularization strategy. Hierarchical Bayes models may be a good candidate for
formulating an implicit, non-Tikhonov regularization in the future. For now,
we observe that the choice of the Tikhonov regularization remains an important
issue for the stability of ROMs, but coupling the state-of-the-art error-based
regularization selection of [31] with BayesOpInf leads to strong performance.
Even with carefully chosen regularization, reduced models drawn from the
posterior distribution may be unstable. Ensuring ROM stability — and, more
generally, embedding desired model structure such as symmetry or spectral
properties — through the posterior distribution remains an important topic for
future work.

In summary, BayesOpInf quantifies the modeling uncertainties stemming from the
data noise, the discretization errors of time derivatives, and the model
misspecification error when the full-order model generating the snapshots does
not match the assumed form of the ROM, as well as the uncertainties introduced
by the use of Tikhonov regularization. As demonstrated in the numerical
examples, BayesOpInf is effective in constructing accurate ROMs from noisy
snapshot data and is capable of capturing complex physics while simultaneously
providing meaningful credible intervals.

## Acknowledgments

The first author acknowledges the support from the ARPA-E DIFFERENTIATE
Program under award DE-AR0001208 during his stay at The University of Texas at
Austin and the current financial support from Sectorplan Bèta (the
Netherlands) under the focus area _Mathematics of Computational Science_. The
second and third authors aknowledge the U.S. Department of Energy, National
Nuclear Security Administration under Award Number DE-NA0003969; the U.S.
Department of Energy AEOLUS MMICC center under award DE-SC0019303; and the Air
Force Center of Excellence on Multi-Fidelity Modeling of Rocket Combustor
Dynamics (Air Force Office of Scientific Research), award FA9550-17-1-0195.

## Appendix 1: Alternative definition of the likelihood

In (3.1), an independent Gaussian noise was considered for the modeling of
residual errors. If one would like to take into account the correlation over
time and/or construct a surrogate model of the ROM closure, a Gaussian process
can be used to replace the noise term, i.e.,

$\frac{\differential}{\differential
t}\hat{q}_{i}(t)=\mathbf{d}(\hat{\mathbf{q}}(t),\mathbf{u}(t))^{\mathsf{T}}\hat{\mathbf{o}}_{i}+\epsilon_{i}(t),\quad\epsilon_{i}(t)\sim\mathcal{GP}(0,\sigma_{i}^{2}\kappa(t,t^{\prime}|{\bf\it\ell}_{i}))\,,$
(A.1)

in which $\kappa(t,t^{\prime}|{\bf\it\ell}_{i})$ is a given kernel function
with hyperparameters ${\bf\it\ell}_{i}$. In this case the likelihood is
alternatively defined as

$p(\mathbf{r}_{i}|~{}\mathbf{D},\hat{\mathbf{o}}_{i},\sigma_{i}^{2},{\bf\it\ell})=\mathcal{N}(\mathbf{r}_{i}|~{}\mathbf{D}\hat{\mathbf{o}}_{i},\sigma_{i}^{2}\mathbf{K}_{i}({\bf\it\ell}_{i}))\,,$
(A.2)

where the covariance matrix
$\mathbf{K}_{i}({\bf\it\ell}_{i})\in\mathbb{R}^{k\times k}$ is defined to be
$\mathbf{K}_{i}({\bf\it\ell}_{i})=\kappa(\mathcal{T},\mathcal{T}|{\bf\it\ell}_{i})$
and assumed to be positive definite. Here $\mathcal{T}$ denotes the set of $k$
time instances $\\{t_{0},t_{1},\cdots,t_{k-1}\\}$ corresponding to the
training data. Conditioning on the training data
$(\mathbf{D},\mathbf{r}_{i})$, the posterior distribution
$p(\hat{\mathbf{o}}_{i}|~{}\mathbf{D},\mathbf{r}_{i},\sigma_{i}^{2},{\bf\it\lambda}_{i},{\bf\it\ell}_{i})=\mathcal{N}(\hat{\mathbf{o}}_{i}|{\bf\it\mu}_{i},{\bf\it\Sigma}_{i})$
thus has new formulations of mean vector and covariance matrix as follows:

$\begin{split}{\bf\it\Sigma}_{i}&=\sigma_{i}^{2}\left[\mathrm{diag}({\bf\it\lambda}_{i})+\mathbf{D}^{\mathsf{T}}\mathbf{K}^{-1}_{i}\mathbf{D}\right]^{-1}\\\
{\bf\it\mu}_{i}&=\left[\mathrm{diag}({\bf\it\lambda}_{i})+\mathbf{D}^{\mathsf{T}}\mathbf{K}^{-1}_{i}\mathbf{D}\right]^{-1}\left(\mathrm{diag}({\bf\it\lambda}_{i}){\bf\it\beta}_{i}+\mathbf{D}^{\mathsf{T}}\mathbf{K}^{-1}_{i}\mathbf{r}_{i}\right)\\\
&={\bf\it\beta}_{i}+\underbrace{\left[\mathrm{diag}({\bf\it\lambda}_{i})+\mathbf{D}^{\mathsf{T}}\mathbf{K}^{-1}_{i}\mathbf{D}\right]^{-1}\mathbf{D}^{\mathsf{T}}\mathbf{K}^{-1}_{i}(\mathbf{r}_{i}-\mathbf{D}{\bf\it\beta}_{i})}_{:=\delta{\bf\it\mu}_{i}}\,.\end{split}$
(A.3)

Note that the correction term
$\delta{\bf\it\mu}_{i}={\bf\it\mu}_{i}-{\bf\it\beta}_{i}$ is then the solution
to the following generalized least squares problem with a Tikhonov
regularization

$\delta{\bf\it\mu}_{i}=\arg\min_{{\bf\it\eta}\in\mathbb{R}^{d(r,m)}}\left\\{\|\mathbf{r}_{i}-\mathbf{D}({\bf\it\beta}_{i}+{\bf\it\eta})\|_{\mathbf{K}^{-1}_{i}}^{2}+\|\mathrm{diag}({\bf\it\lambda}_{i})^{\frac{1}{2}}{\bf\it\eta}\|_{2}^{2}\right\\}\,.$
(A.4)

Here the norm $\|\cdot\|_{\mathbf{K}^{-1}_{i}}$ is defined as
$\|\mathbf{r}\|_{\mathbf{K}^{-1}_{i}}=\sqrt{\mathbf{r}^{\mathsf{T}}\mathbf{K}^{-1}_{i}\mathbf{r}}$,
$\mathbf{r}\in\mathbb{R}^{k}$. In the meantime, the estimate of
hyperparameters $(\sigma_{i}^{2},{\bf\it\lambda}_{i},{\bf\it\ell}_{i})$ by
maximizing marginal likelihood is reformulated as

$\begin{split}(\sigma^{*2}_{i},{\bf\it\lambda}_{i}^{*},{\bf\it\ell}^{*}_{i})=\arg\max_{\sigma_{i}^{2},{\bf\it\lambda}_{i},{\bf\it\ell}_{i}}&~{}\log~{}p(\mathbf{r}_{i}|~{}\mathbf{D},\sigma_{i}^{2},{\bf\it\lambda}_{i},{\bf\it\ell}_{i})\\\
=\arg\max_{\sigma_{i}^{2},{\bf\it\lambda}_{i},{\bf\it\ell}_{i}}&~{}\log\int
p(\mathbf{r}_{i}|~{}\mathbf{D},\hat{\mathbf{o}}_{i},\sigma_{i}^{2},{\bf\it\ell}_{i})~{}p(\hat{\mathbf{o}}_{i}|~{}\sigma_{i}^{2},{\bf\it\lambda}_{i})~{}\differential\hat{\mathbf{o}}_{i}\\\
=\arg\max_{\sigma_{i}^{2},{\bf\it\lambda}_{i},{\bf\it\ell}_{i}}&~{}\log~{}\mathcal{N}\left(\mathbf{r}_{i}|~{}\mathbf{D}{\bf\it\beta}_{i},\sigma_{i}^{2}\left(\mathbf{D}\mathrm{diag}({\bf\it\lambda}_{i})^{-1}\mathbf{D}^{\mathsf{T}}+\mathbf{K}_{i}({\bf\it\ell}_{i})\right)\right)\,.\end{split}$
(A.5)

The model (A.1) is equivalent to a universal kriging with a prior Gaussian
process as follows:

$\frac{\differential}{\differential
t}\hat{q}_{i}(t)\sim\mathcal{GP}\bigg{(}\mathbf{d}(\hat{\mathbf{q}}(t),\mathbf{u}(t))^{\mathsf{T}}{\bf\it\beta}_{i},\sigma_{i}^{2}\mathbf{d}(\hat{\mathbf{q}}(t),\mathbf{u}(t))^{\mathsf{T}}\mathrm{diag}({\bf\it\lambda}_{i})^{-1}\mathbf{d}(\hat{\mathbf{q}}(t^{\prime}),\mathbf{u}(t^{\prime}))+\sigma_{i}^{2}\kappa(t,t^{\prime}|{\bf\it\ell}_{i})\bigg{)}\,.$
(A.6)

Accordingly, the posterior prediction of $\frac{\differential}{\differential
t}\hat{q}_{i}(t)$ is again a Gaussian process whose mean and covariance
functions take the form of

$\begin{split}\mathbb{E}\left[\frac{\differential}{\differential
t}\hat{q}_{i}(t)\middle|~{}\mathbf{D},\mathbf{r}_{i}\right]=&~{}\mathbf{d}(\hat{\mathbf{q}}(t),\mathbf{u}(t))^{\mathsf{T}}{\bf\it\mu}_{i}+\kappa(t,\mathcal{T}|{\bf\it\ell}_{i})\mathbf{K}^{-1}_{i}(\mathbf{r}_{i}-\mathbf{D}{\bf\it\mu}_{i})\,,\quad\text{and}\\\
\mathbb{C}ov\left[\frac{\differential}{\differential
t}\hat{q}_{i}(t),\frac{\differential}{\differential
t}\hat{q}_{i}(t^{\prime})\middle|~{}\mathbf{D},\mathbf{r}_{i}\right]=&~{}\sigma_{i}^{2}\kappa(t,t^{\prime}|{\bf\it\ell}_{i})-\sigma_{i}^{2}\kappa(t,\mathcal{T}|{\bf\it\ell}_{i})\mathbf{K}^{-1}_{i}\kappa(\mathcal{T},t^{\prime}|{\bf\it\ell}_{i})\\\
+\big{[}\mathbf{d}(\hat{\mathbf{q}}(t),\mathbf{u}(t))^{\mathsf{T}}&-\kappa(t,\mathcal{T}|{\bf\it\ell}_{i})\mathbf{K}_{i}^{-1}\mathbf{D}\big{]}{\bf\it\Sigma}_{i}\big{[}\mathbf{d}(\hat{\mathbf{q}}(t^{\prime}),\mathbf{u}(t^{\prime}))-\mathbf{D}^{\mathsf{T}}\mathbf{K}_{i}^{-1}\kappa(\mathcal{T},t^{\prime}|{\bf\it\ell}_{i})\big{]}\,,\end{split}$
(A.7)

respectively. Note that the second term in the mean function, i.e.,
$\kappa(t,\mathcal{T}|{\bf\it\ell}_{i})\mathbf{K}^{-1}_{i}(\mathbf{r}_{i}-\mathbf{D}{\bf\it\mu}_{i})$,
formulates a surrogate of the closure in the linear model (2.2).

## Appendix 2: Error-based regularization selection

Algorithm 3 details the OpInf procedure where the regularization parameters
are selected with a search-based optimization method [31]. Regularization
parameters candidates are assessed by solving the corresponding deterministic
OpInf problem (2.4), then using the resultant ROM to reconstruct the training
data for $L\geq 1$ initial conditions and make predictions in time beyond the
training regime. The projected training snapshots are organized by initial
condition, i.e., $\hat{\mathbf{Q}}^{(\ell)}$ contains all projected snapshots
for the $\ell$th initial condition.

1:procedure OpInf( Projected training snapshots
$\hat{\mathbf{Q}}=[~{}\hat{\mathbf{Q}}^{(1)}~{}\cdots~{}\hat{\mathbf{Q}}^{(L)}~{}]\in\mathbb{R}^{r\times
k}$, training inputs $\mathbf{U}\in\mathbb{R}^{m\times k}$, projected snapshot
time derivatives
$\mathbf{R}=[~{}\mathbf{r}_{1}~{}\cdots~{}\mathbf{r}_{r}~{}]^{\mathsf{T}}\in\mathbb{R}^{r\times
k}$, final time $t_{f}>t_{0}$, bound margin $\tau\geq 1$ )

2: $\mathbf{D}\leftarrow\texttt{data\\_matrix}(\hat{\mathbf{Q}},\mathbf{U})$ #
Construct data matrix as in (2.5).

3: $B\leftarrow\tau\max_{i,j,\ell}|\hat{\mathbf{Q}}_{i,j}^{(\ell)}|$ # Set a
bound for the ROM stability criterion.

4: # Subprocedure for assessing regularization parameters.

5: procedure OpInfError( regularizer
$\boldsymbol{\lambda}\in\mathbb{R}^{d(r,m)}$ )

6: for $i=1,\ldots,r$ do # OpInf with current regularization.

7:
$\hat{\mathbf{o}}_{i}\leftarrow\arg\min_{\boldsymbol{\mu}_{i}}\left\|\mathbf{D}\boldsymbol{\mu}_{i}-\mathbf{r}_{i}\right\|_{2}^{2}+\left\|\textrm{diag}(\boldsymbol{\lambda})^{1/2}\boldsymbol{\mu}_{i}\right\|_{2}^{2}$

8: for $\ell=1,\ldots,L$ do

9: $\tilde{\mathbf{Q}}\leftarrow~{}\textrm{integrate \eqref{eq:reducedorder}
with}~{}\hat{\mathbf{O}}=[~{}\hat{\mathbf{o}}_{1}~{}\cdots~{}\hat{\mathbf{o}}_{r}~{}]^{\mathsf{T}}$
from $\hat{\mathbf{q}}_{0}=\hat{\mathbf{Q}}_{:,0}$ over $[t_{0},t_{f}]$

10: if $\max_{i,j,\ell}|\tilde{\mathbf{Q}}^{(\ell)}_{i,j}|>B$ then # Ensure
the ROM is stable on $[t_{0},t_{f}]$ for all ICs.

11: return $\infty$

12: return
$\displaystyle\frac{1}{L}\sum_{\ell=1}^{L}\left\|\hat{\mathbf{Q}}^{(\ell)}-\tilde{\mathbf{Q}}_{:,:k}^{(\ell)}\right\|$
# Compare ROM outputs to training data.

13: # Minimize the subprocedure to choose the regularization parameters.

14: $\textbf{return}\arg\min_{\boldsymbol{\lambda}}$
OpInfError$(\boldsymbol{\lambda})$

Algorithm 3 Operator Inference with error-based regularization selection
(Algorithm 1 of [31]).

In the results for Section 5.1 (Figure 5), Algorithm 3 is tailored to the
purely quadratic ROM (5.3). In this setting, there are no external inputs
$\mathbf{U}$, the data matrix of step 2 is given by
$\mathbf{D}=(\hat{\mathbf{Q}}\odot\hat{\mathbf{Q}})^{\mathsf{T}}$, and the
operator matrix of step 9 is $\hat{\mathbf{O}}=\hat{\mathbf{H}}$. The
regularization parameters $\boldsymbol{\lambda}$ are globally parameterized by
a single scalar hyperparameter $\lambda\geq 0$, i.e.,
$\boldsymbol{\lambda}=[~{}\lambda~{}\cdots~{}\lambda~{}]^{\mathsf{T}}$. Hence,
the minimization in step 14 can be carried out with a one-dimensional
optimization on $\lambda$ via, e.g., Brent’s method [7].

For the combustion problem of Section 5.2, we have a single initial condition
($L=1$) and the ROM is (2.2) with data matrix
$\mathbf{D}=[~{}\hat{\mathbf{Q}}^{\mathsf{T}}~{}~{}(\hat{\mathbf{Q}}\odot\hat{\mathbf{Q}})^{\mathsf{T}}~{}~{}\mathbf{U}^{\mathsf{T}}~{}~{}\mathbf{1}_{k}~{}]$
and reduced-order operators
$\hat{\mathbf{O}}=[~{}\hat{\mathbf{A}}~{}\hat{\mathbf{H}}~{}\hat{\mathbf{B}}~{}\hat{\mathbf{c}}~{}]$.
The regularization parameters are parameterized by two scalars
$\lambda_{1},\lambda_{2}\geq 0$:

$\boldsymbol{\lambda}=[~{}\underbrace{\lambda_{1},\ldots,\lambda_{1}}_{r~{}\textrm{times}},\underbrace{\lambda_{2},\ldots,\lambda_{2}}_{\frac{r(r+1)}{2}~{}\textrm{times}},\lambda_{1},\lambda_{1}~{}]^{\mathsf{T}}.$
(B.1)

This choice of $\boldsymbol{\lambda}$ penalizes the entries of the quadratic
operator $\hat{\mathbf{H}}$ by $\lambda_{2}$ and the entries of the remaining
operators $\hat{\mathbf{A}}$, $\hat{\mathbf{B}}$, and $\hat{\mathbf{c}}$ by
$\lambda_{1}$. More precisely, the regularization term as written in (2.3) is
given by

$\mathcal{P}\left([~{}\hat{\mathbf{A}}~{}\hat{\mathbf{H}}~{}\hat{\mathbf{B}}~{}\hat{\mathbf{c}}~{}]\right)=\lambda_{1}\left(\|\hat{\mathbf{A}}\|_{F}^{2}+\|\hat{\mathbf{B}}\|_{F}^{2}+\|\hat{\mathbf{c}}\|_{2}^{2}\right)+\lambda_{2}\|\hat{\mathbf{H}}\|_{F}^{2}.$
(B.2)

The optimization in step 14 of Algorithm 3 is carried out with a coarse grid
search in the principal quadrant of $\mathbb{R}^{2}$, followed by the Nelder-
Mead method [33].

## References

  * [1] M. Ainsworth and J. T. Oden.  A Posteriori Error Estimation in Finite Element Analysis, volume 37.  John Wiley & Sons, 2011.
  * [2] I. Babuska and J. T. Oden.  Verification and validation in computational engineering and science: basic concepts.  Computer Methods in Applied Mechanics and Engineering, 36(193):4057–4066, 2004.
  * [3] P. Benner, P. Goyal, J. Heiland, and I. P. Duff.  Operator inference and physics-informed learning of low-dimensional models for incompressible flows.  Electronic Transactions on Numerical Analysis, 56:28–51, 2022.
  * [4] P. Benner, S. Gugercin, and K. Willcox.  A survey of projection-based model reduction methods for parametric dynamical systems.  SIAM Review, 57(4):483–531, 2015.
  * [5] G. Berkooz, P. Holmes, and J. Lumley.  The proper orthogonal decomposition in the analysis of turbulent flows.  Annual Review of Fluid Mechanics, 25:539–575, 1993.
  * [6] G. E. P. Box and G. C. Tiao.  Bayesian Inference in Statistical Analysis.  John Wiley & Sons, 2011.
  * [7] R. P. Brent.  Algorithms for Minimization Without Derivatives.  Dover Publications, New York, NY, 2002.
  * [8] K. Champion, B. Lusch, J. N. Kutz, and S. L. Brunton.  Data-driven discovery of coordinates and governing equations.  Proceedings of the National Academy of Sciences, 116(45):22445–22451, 2019.
  * [9] O. A. Chkrebtii, D. A. Campbell, B. Calderhead, and M. A. Girolami.  Bayesian solution uncertainty quantification for differential equations.  Bayesian Analysis, 11(4):1239–1267, 2016.
  * [10] K. De Brabanter, J. De Brabanter, B. De Moor, and I. Gijbels.  Derivative estimation with local polynomial fitting.  Journal of Machine Learning Research, 14(1), 2013.
  * [11] D. del Castillo-Negrete, D. A. Spong, and S. P. Hirshman.  Proper orthogonal decomposition methods for noise reduction in particle-based transport calculations.  Physics of Plasmas, 15(9):092308, 2008.
  * [12] B. P. Epps and A. H. Techet.  An error threshold criterion for singular value decomposition modes extracted from PIV data.  Experiments in fluids, 48(2):355–367, 2010.
  * [13] K. Farrell, J. T. Oden, and D. Faghihi.  A Bayesian framework for adaptive selection, calibration, and validation of coarse-grained models of atomistic systems.  Journal of Computational Physics, 295:189–208, 2015.
  * [14] R. Geelen and K. Willcox.  Localized non-intrusive reduced-order modeling in the operator inference framework.  Philosophical Transactions of the Royal Society A, to appear, 2022\.
  * [15] O. Ghattas and K. Willcox.  Learning physics-based models from data: perspectives from inverse problems and model reduction.  Acta Numerica, 30:445–554, 2021.
  * [16] W. R. Graham, J. Peraire, and K. Y. Tang.  Optimal control of vortex shedding using low-order models. Part I—Open-loop model development.  International Journal for Numerical Methods in Engineering, 44(7):945–972, 1999.
  * [17] M. Guo and J. S. Hesthaven.  Data-driven reduced order modeling for time-dependent problems.  Computer Methods in Applied Mechanics and Engineering, 345:75–99, 2019.
  * [18] N. Halko, P.-G. Martinsson, and J. A. Tropp.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  SIAM Review, 53(2):217–288, 2011.
  * [19] M. E. Harvazinski, C. Huang, V. Sankaran, T. W. Feldman, W. E. Anderson, C. L. Merkle, and D. G. Talley.  Coupling between hydrodynamics, acoustics, and heat release in a self-excited unstable combustor.  Physics of Fluids, 27(4):045102, 2015.
  * [20] A. Hawkins-Daarud, S. Prudhomme, K. G. van der Zee, and J. T. Oden.  Bayesian calibration, validation, and uncertainty quantification of diffuse interface models of tumor growth.  Journal of Mathematical Biology, 67(6):1457–1485, 2013.
  * [21] S. M. Hirsh, D. A. Barajas-Solano, and J. N. Kutz.  Sparsifying priors for bayesian uncertainty quantification in model discovery.  arXiv preprint arXiv:2107.02107, 2021.
  * [22] C. Huang, K. Duraisamy, and C. Merkle.  Challenges in reduced order modeling of reacting flows.  In 2018 Joint Propulsion Conference, Cincinnati, OH, 2018.  Paper AIAA-2018-4675.
  * [23] C. Huang, K. Duraisamy, and C. L. Merkle.  Investigations and improvement of robustness of reduced-order models of reacting flow.  AIAA Journal, 57(12):5377–5389, 2019.
  * [24] C. Huang, J. Xu, K. Duraisamy, and C. Merkle.  Exploration of reduced-order models for rocket combustion applications.  In 2018 AIAA Aerospace Sciences Meeting, Orlando, FL, 2018.  Paper AIAA-2018-1183.
  * [25] P. Jain, S. McQuarrie, and B. Kramer.  Performance comparison of data-driven reduced models for a single-injector combustion process.  In AIAA Propulsion and Energy 2021 Forum, Virtual Event, 2021.  Paper AIAA-2021-3633.
  * [26] P. K. Jha, L. Cao, and J. T. Oden.  Bayesian-based predictions of COVID-19 evolution in 0texas using multispecies mixture-theoretic continuum models.  Computational Mechanics, 66(5):1055–1068, 2020.
  * [27] C. Khatri and C. R. Rao.  Solutions to some functional equations and their applications to characterization of probability distributions.  Sankhy = a: The Indian Journal of Statistics, Series A, pages 167–180, 1968.
  * [28] T. G. Kolda and B. W. Bader.  Tensor decompositions and applications.  SIAM review, 51(3):455–500, 2009.
  * [29] B. Kramer and K. Willcox.  Nonlinear model order reduction via lifting transformations and proper orthogonal decomposition.  AIAA Journal, 57(6):2297–2307, 2019.
  * [30] K. Lee and K. T. Carlberg.  Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.  Journal of Computational Physics, 404:108973, 2020.
  * [31] S. A. McQuarrie, C. Huang, and K. E. Willcox.  Data-driven reduced-order models via regularised operator inference for a single-injector combustion process.  Journal of the Royal Society of New Zealand, 51(2):194–211, 2021\.
  * [32] S. A. McQuarrie, P. Khodabakhshi, and K. E. Willcox.  Non-intrusive reduced-order models for parametric partial differential equations via data-driven operator inference.  arXiv preprint arXiv:2110.07653, 2021.
  * [33] J. A. Nelder and R. Mead.  A simplex method for function minimization.  The Computer Journal, 7(4):308–313, 1965.
  * [34] J. T. Oden.  Adaptive multiscale predictive modelling.  Acta Numerica, 27:353–450, 2018.
  * [35] J. T. Oden, I. Babuška, and D. Faghihi.  Predictive computational science: Computer predictions in the presence of uncertainty.  Encyclopedia of Computational Mechanics Second Edition, pages 1–26, 2017.
  * [36] J. T. Oden, I. Babuška, F. Nobile, Y. Feng, and R. Tempone.  Theory and methodology for estimation and control of errors due to modeling, approximation, and uncertainty.  Computer Methods in Applied Mechanics and Engineering, 194(2-5):195–204, 2005.
  * [37] J. T. Oden and S. Prudhomme.  Estimation of modeling error in computational mechanics.  Journal of Computational Physics, 182(2):496–515, 2002.
  * [38] W. Pan, Y. Yuan, J. Gonçalves, and G.-B. Stan.  A sparse bayesian approach to the identification of nonlinear state-space systems.  IEEE Transactions on Automatic Control, 61(1):182–187, 2015.
  * [39] B. Peherstorfer.  Sampling low-dimensional Markovian dynamics for pre-asymptotically recovering reduced models from data with operator inference.  SIAM Journal on Scientific Computing, 42(5):A3489–A3515, 2020.
  * [40] B. Peherstorfer and K. Willcox.  Data-driven operator inference for nonintrusive projection-based model reduction.  Computer Methods in Applied Mechanics and Engineering, 306:196–215, 2016.
  * [41] E. Prudencio, P. Bauman, D. Faghihi, K. Ravi-Chandar, and J. T. Oden.  A computational framework for dynamic data-driven material damage control, based on Bayesian inference and model selection.  International Journal for Numerical Methods in Engineering, 102(3-4):379–403, 2015.
  * [42] E. Qian, B. Kramer, B. Peherstorfer, and K. Willcox.  Lift & Learn: Physics-informed machine learning for large-scale nonlinear dynamical systems.  Physica D: Nonlinear Phenomena, 406:132401, 2020.
  * [43] A. Quarteroni, G. Rozza, et al.  Reduced order methods for modeling and computational reduction, volume 9.  Springer, 2014.
  * [44] H. Schaeffer.  Learning partial differential equations via data discovery and sparse optimization.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2197):20160446, 2017.
  * [45] P. J. Schmid.  Dynamic mode decomposition of numerical and experimental data.  Journal of Fluid Mechanics, 656:5–28, 2010.
  * [46] L. Sirovich.  Turbulence and the dynamics of coherent structures. I. Coherent structures.  Quarterly of Applied Mathematics, 45(3):561–571, 1987.
  * [47] R. C. Smith.  Uncertainty Quantification: Theory, Implementation, and Applications.  SIAM, 2013.
  * [48] R. Swischuk, B. Kramer, C. Huang, and K. Willcox.  Learning physics-based reduced-order models for a single-injector combustion process.  AIAA Journal, 58(6):2658–2672, 2020.
  * [49] W. I. T. Uy and B. Peherstorfer.  Probabilistic error estimation for non-intrusive reduced models learned from data of systems governed by linear parabolic partial differential equations.  ESAIM: Mathematical Modelling and Numerical Analysis, 55(3):735–761, 2021.
  * [50] W. I. T. Uy, Y. Wang, Y. Wen, and B. Peherstorfer.  Active operator inference for learning low-dimensional dynamical-system models from noisy data.  arXiv preprint arXiv:2107.09256, 2021.
  * [51] C. F. Van Loan.  The ubiquitous Kronecker product.  Journal of Computational and Applied Mathematics, 123(1-2):85–100, 2000.
  * [52] D. Venturi.  On proper orthogonal decomposition of randomly perturbed fields with applications to flow past a cylinder and natural convection over a horizontal plate.  Journal of Fluid Mechanics, 559:215–254, 2006.
  * [53] C. K. Westbrook and F. L. Dryer.  Simplified reaction mechanisms for the oxidation of hydrocarbon fuels in flames.  Combustion Science and Technology, 27(1-2):31–43, 1981.
  * [54] C. K. Williams and C. E. Rasmussen.  Gaussian Processes for Machine Learning.  MIT Press, 2006.
  * [55] S. Zhang and G. Lin.  Robust data-driven discovery of governing physical laws with error bars.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2217):20180305, 2018.
  * [56] Q. Zhuang, J. M. Lorenzi, H.-J. Bungartz, and D. Hartmann.  Model order reduction based on Runge-Kutta neural networks.  Data-Centric Engineering, 2, 2021.