Interpretable and flexible non-intrusive reduced-order models using reproducing kernel Hilbert spaces

Abstract

This paper develops an interpretable, non-intrusive reduced-order modeling technique using regularized kernel interpolation. Existing non-intrusive approaches approximate the dynamics of a reduced-order model (ROM) by solving a data-driven least-squares regression problem for low-dimensional matrix operators. Our approach instead leverages regularized kernel interpolation, which yields an optimal approximation of the ROM dynamics from a user-defined reproducing kernel Hilbert space. We show that our kernel-based approach can produce interpretable ROMs whose structure mirrors full-order model structure by embedding judiciously chosen feature maps into the kernel. The approach is flexible and allows a combination of informed structure through feature maps and closure terms via more general nonlinear terms in the kernel. We also derive a computable a posteriori error bound that combines standard error estimates for intrusive projection-based ROMs and kernel interpolants. The approach is demonstrated in several numerical experiments that include comparisons to operator inference using both proper orthogonal decomposition and quadratic manifold dimension reduction.

Keywords: data-driven model reduction, kernel interpolation, feature maps, interpretable reduced-order model, error bounds, quadratic manifolds

1 Introduction

Large-scale numerical simulations are a crucial component of the engineering design process. For many applications, the complexity of the underlying physics and the required fidelity make such simulations highly computationally expensive, which renders many-query simulation tasks such as uncertainty quantification and design optimization infeasible. Model reduction techniques seek to mitigate high computation costs in numerical simulations by systematically extracting the relevant dynamics of a large-scale system, called the full-order model (FOM), and constructing a low-dimensional, computationally efficient reduced-order model (ROM), which can be used as a substitute for the FOM in many-query design tasks. Two appealing features of ROMs over other surrogate modeling techniques are that they aim to incorporate underlying physics from the FOM and often come equipped with rigorous error bounds. In this paper, we propose a novel model reduction framework that uses regularized kernel interpolation to compute data-driven ROMs that are interpretable, flexible, and have rigorous error estimates.

Classical projection-based model reduction techniques construct ROMs by identifying a low-dimensional linear subspace that best represents the FOM dynamics in some sense, then projecting the governing equations onto the subspace. Examples of projection-based approaches include balanced truncation [1, 7]; interpolatory projections [2, 16]; moment-matching [6, 22]; and proper orthogonal decomposition (POD) [23, 25], in which the optimal low-dimensional subspace is defined as the span of the leading left singular vectors of a representative set of state data. In recent years, several dimension reduction approaches have been proposed that aim to overcome approximation limitations of linear subspaces, including nonlinear manifolds (NMs) using autoencoders [14, 15, 28, 32, 44], quadratic manifolds (QMs) [4, 18, 19, 54], and the projection-based ROM + artificial neural network (PROM-ANN) approach [5]. These strategies are especially beneficial when applied to problems with slowly decaying Kolmogorov $n$-width, such as transport-dominated problems or problems with sharp gradients [37, 39]. In many cases, these nonlinear dimension reduction approaches can still be used to produce projection-based ROMs by inserting the state approximation into the governing FOM and projecting the residual by a test basis.

Projection-based methods have enjoyed success in a number of applications. However, a common disadvantage is that they require intrusive access to code of a given FOM. This is often an infeasible request when the FOM is defined through legacy or commercial code, and hence an intrusive projection-based ROM is unobtainable. Several so-called non-intrusive model reduction approaches have been developed recently to overcome this difficulty. These methods apply a dimension reduction technique, such as POD or an autoencoder, to project pre-computed snapshot data onto a low-dimensional latent space and learn a function that models the system dynamics within the latent space. For example, dynamic mode decomposition (DMD) [31, 46, 51, 52, 58] approximates a dynamical system by fitting a least-squares optimal linear operator to time series data. This approach has been extended to approximate nonlinear dynamical systems using Koopman operator theory, but selecting observables that yield approximately linear dynamics can be challenging [13, 35, 45, 61]. Operator inference (OpInf) [20, 29, 40] is a related method that constrains the learnable dynamics to have the same structure (e.g, polynomial) as a projection-based ROM, thereby producing interpretable nonlinear ROMs. In this method, reduced-order operators are computed by solving a linear least-squares regression that minimizes the residual of the desired reduced dynamics. Non-polynomial nonlinearities can often be incorporated by first applying a lifting transformation to the training data, then learning a polynomial ROM [43]. By contrast, neural network (NN)-based approaches [10, 33, 47, 48] typically use autoencoders for the dimension reduction and model the reduced dynamics using a NN. While these methods are very flexible in that they can model dynamics with arbitrary structure, the resulting ROMs are not interpretable. Another method, latent space dynamics identification (LaSDI) [11, 12, 17, 24, 38], can be viewed as a hybrid of OpInf and NN-based approaches that typically uses autoencoder-based dimension reduction and learns reduced-order dynamics by solving a least-squares regression problem for coefficient matrices corresponding to a library of nonlinear candidates functions. This approach is related to the SINDy algorithm [27] but does not enforce a sparsity requirement. The library of candidate functions for LaSDI is typically chosen to be polynomial, which results in solving a similar least-squares regression problem to OpInf when learning the latent dynamics. Unlike OpInf, while the resulting ROM structure is interpretable, a natural structure for the ROM dynamics cannot be deduced a priori since autoencoder-based dimension reduction does not preserve structure from the FOM. In each of these approaches, error estimates for the resulting ROMs are limited, with the exception of the recent thermodynamics-based LaSDI approach [38].

Our proposed kernel-based non-intrusive ROMs, which we call “Kernel ROMs”, share similarities with the aforementioned approaches while overcoming some noticeable drawbacks. Like other approaches, we begin by applying POD or QM dimension reduction to a set of training snapshots and learn a function that approximates the system dynamics within a latent space. However, instead of modeling the ROM dynamics as a polynomial and learning the polynomial coefficients through least-squares regression as in OpInf and LaSDI, we use regularized kernel interpolation [36, 49, 50] to model the reduced dynamics with a function belonging to a user-defined reproducing kernel Hilbert space (RKHS). The structure of the learned function depends on the positive-definite kernel that defines the RKHS. For example, if the governing FOM has a polynomial structure, we can use a kernel induced by a feature map to compute ROM dynamics that share the same polynomial structure. On the other hand, if the FOM dynamics have unknown or only partially known structure, a more generic nonlinear kernel can be used to model the unknown part of the ROM dynamics. In this sense, our proposed approach has a natural way of incorporating closure terms into the ROM dynamics. While kernel methods have been used to emulate reduced dynamics in previous work [41, 62], these approaches are intrusive in that they assume that the FOM dynamics can be sampled explicitly, and they do not demonstrate a way to inject explicit structure into the learned ROM. The authors in [3] use kernel methods to augment a DMD model, resulting in a fully data-driven surrogate model, and implicitly inject structure by modeling nonlinear terms using polynomial kernels. While this approach is similar to ours, we focus on constructing non-intrusive ROMs, and can model nonlinear terms explicitly using feature map kernels. In summary, the proposed approach is entirely data-driven, can produce interpretable and flexible ROMs, and yields computable a posteriori error bounds between the non-intrusive ROM and FOM solutions.

The outline of this paper is as follows. We first review essential aspects of regularized kernel interpolation in Section 2. We then review intrusive projection-based model reduction in Section 3, with a focus on quadratic dimension reduction and the resulting model structure. Section 4 details the application of regularization kernel interpolation in the non-intrusive model reduction setting, and a corresponding a posteriori error analysis is provided in Section 5. We demonstrate our proposed approach numerically on several examples in Section 6, including comparisons to OpInf and intrusive ROMs when possible. The results show that our proposed approach can accommodate either POD or QM dimension reduction and produces comparable results to OpInf while also yielding a computable error bound. Finally, Section 7 provides a few concluding remarks and identifies potential avenues for future development.

2 Regularized kernel interpolation

This section reviews the essentials of regularized kernel interpolation, the key ingredient for our non-intrusive model reduction approach. Section 2.1 reviews scalar-valued interpolation, which is extended to vector-valued interpolation in Section 2.2. Scenario-specific kernel design is then discussed in Section 2.3.

2.1 Scalar-valued kernel interpolation

We begin with a review of regularized kernel interpolation for scalar-valued functions. By the Moore–Aronszajn Theorem (see, e.g., [49, Theorem 3.10]), a positive-definite kernel function defines a unique Hilbert space with desirable properties. This result leads to the Representer Theorem — the key result used for computing an optimal interpolant in an RKHS — as well as a pointwise error bound on the interpolant.

Definition 2.1 (Positive-definite kernels).

A function $K:{}^{n_{x}}\times{}^{n_{x}}\to\real$ is a (real-valued) kernel function if it is symmetric, i.e., $K(\mathbf{x},\mathbf{x}^{\prime})=K(\mathbf{x}^{\prime},\mathbf{x})$ for all $\mathbf{x},\mathbf{x}^{\prime}\in{}^{n_{x}}$. A kernel function $K$ is said to be positive definite if for any matrix $\mathbf{X}=[~{}\mathbf{x}_{1}~{}~{}\cdots~{}~{}\mathbf{x}_{m}~{}]\in{}^{n_{x}\times m}$ with pairwise distinct columns, the kernel matrix $K(\mathbf{X},\mathbf{X})\in{}^{m\times m}$ with entries $K(\mathbf{X},\mathbf{X})_{ij}=K(\mathbf{x}_{i},\mathbf{x}_{j})$ is positive semi-definite.

Definition 2.2 (RKHS).

Let $K:{}^{n_{x}}\times{}^{n_{x}}\to\real$ be a positive-definite kernel function. Consider the pre-Hilbert space of functions

$${\cal H}_{K}^{0}({}^{n_{x}})=\left\{v:{}^{n_{x}}\to\real~{}\bigg{|}~{}\exists\,m\in\mathbb{N},\,{\boldsymbol{\omega}}\in{}^{m},\,\left\{\mathbf{x}\right\}_{j=1}^{m}\subset{}^{n_{x}}~{}\textup{such that}~{}v(\mathbf{x})=\sum_{j=1}^{m}\omega_{j}K(\mathbf{x}_{j},\mathbf{x})\right\}.$$

The reproducing kernel Hilbert space (RKHS) ${\cal H}_{K}({}^{n_{x}})$ induced by the kernel $K$ is the (unique) completion of ${\cal H}_{K}^{0}({}^{n_{x}})$ with respect to the norm $\left\|\cdot\right\|_{{\cal H}_{K}({}^{n_{x}})}\coloneqq\left\langle\cdot,\cdot\right\rangle_{{\cal H}_{K}({}^{n_{x}})}^{1/2}$ induced by the inner product

$$\left\langle v,v^{\prime}\right\rangle_{{\cal H}_{K}({}^{d})}\coloneqq\sum_{j=1}^{m}\sum_{k=1}^{m^{\prime}}\omega_{j}\omega_{k}^{\prime}K(\mathbf{x}_{j},\mathbf{x}_{k}^{\prime}),$$

in which $v(\mathbf{x})=\sum_{j=1}^{m}\omega_{j}K(\mathbf{x}_{j},\mathbf{x})$ and $v^{\prime}(\mathbf{x})=\sum_{j=1}^{m^{\prime}}\omega_{k}^{\prime}K(\mathbf{x}_{k}^{\prime},\mathbf{x})$.

For an ordered collection of pairwise-distinct vectors $\{\mathbf{x}_{j}\}_{j=1}^{m}\subset{}^{n_{x}}$, we use $K(\mathbf{X},\mathbf{X})$ to denote $m\times m$ kernel matrix of Definition 2.1 and define the vector $K(\mathbf{X},\mathbf{x})=[~{}K(\mathbf{x}_{1},\mathbf{x})~{}~{}\cdots~{}~{}K(\mathbf{x}_{m},\mathbf{x})~{}]^{\mathsf{T}}\in{}^{m}$. To simplify notation, we will write ${\cal H}_{K}$ for ${\cal H}_{K}({}^{n_{x}})$ when it is understood that $K$ is defined over ${}^{n_{x}}\times{}^{n_{x}}$. Importantly, for $v(\mathbf{x})=\sum_{j=1}^{m}\omega_{j}K(\mathbf{x}_{j},\mathbf{x})$, the induced RKHS norm $\left\|v\right\|_{{\cal H}_{K}}$ can be computed efficiently via the corresponding kernel matrix,

$$\left\|v\right\|_{{\cal H}_{K}}^{2}=\sum_{j=1}^{m}\sum_{k=1}^{m}\omega_{j}\omega_{k}K(\mathbf{x}_{j},\mathbf{x}_{k})={\boldsymbol{\omega}}^{\mathsf{T}}K(\mathbf{X},\mathbf{X}){\boldsymbol{\omega}}.$$ (2.1)

We now state a main result from RKHS theory that is fundamental for our method.

Definition 2.3 (Regularized kernel interpolant).

Let $v:{}^{n_{x}}\to\real$, $\left\{\mathbf{x}_{j}\right\}_{j=1}^{m}\subset{}^{n_{x}}$ be pairwise distinct, and denote $y_{j}=v(\mathbf{x}_{j})$. For a given RKHS ${\cal H}_{K}$ and a regularization parameter $\gamma\geq 0$, a regularized interpolant $s_{v}^{\gamma}\in{\cal H}_{K}$ of $v$ is a solution to the minimization problem

$$\min_{s\in{\cal H}_{K}}\;\sum_{j=1}^{m}(y_{j}-s(\mathbf{x}_{j}))^{2}+\gamma\left\|s\right\|_{{\cal H}_{K}}^{2}.$$ (2.2)

Theorem 2.1 (Representer Theorem).

The minimization problem eq. 2.2 has a solution of the form

$$s_{v}^{\gamma}(\mathbf{x})=\sum_{j=1}^{m}\omega_{j}K(\mathbf{x}_{j},\mathbf{x})={\boldsymbol{\omega}}^{\mathsf{T}}K(\mathbf{X},\mathbf{x}),$$ (2.3a)

where the coefficient vector $${\boldsymbol{\omega}}=[~{}\omega_{1}~{}~{}\cdots~{}~{}\omega_{m}~{}]^{\mathsf{T}}\in{}^{m} m\times m$$ solves the linear system

$$\big{(}K(\mathbf{X},\mathbf{X})+\gamma\mathbf{I}\big{)}{\boldsymbol{\omega}}=\left[\begin{array}[]{c}y_{1}\\ \vdots\\ y_{m}\end{array}\right].$$ (2.3e)

Moreover, if $K$ is strictly positive definite, then $s_{v}^{\gamma}$ is the unique minimizer of eq. 2.2.

See, e.g., [50, Theorem 9.3] for a proof of Theorem 2.1. A key observation from Theorem 2.1 is that a solution to the infinite-dimensional minimization problem eq. 2.2 can be obtained by solving the finite-dimensional linear system eq. 2.3e.

Without regularization ($\gamma=0$), the function $s_{v}^{0}\in{\cal H}_{K}$ exactly interpolates the data, i.e., $s_{v}(\mathbf{x}_{j})=y_{j}$ for each $j=1,\dots,m$. Moreover, $s_{v}^{0}$ satisfies the following error bound.

Theorem 2.2 (Power function error bound).

If $v\in{\cal H}_{K}$ and $s_{v}^{0}\in{\cal H}_{K}$ is an (unregularized) interpolant of $v$ corresponding to the pairwise distinct data $\left\{\mathbf{x}_{i}\right\}_{j=1}^{m}\subset{}^{n_{x}}$ and $y_{j}=v(\mathbf{x}_{j})\in\real$, then

$$|v(\mathbf{x})-s_{v}^{0}(\mathbf{x})|\leq P_{K,\mathbf{X}}(\mathbf{x})\left\|v\right\|_{{\cal H}_{K}}\qquad\forall\;\mathbf{x}\in{}^{d},$$ (2.4a)

where $$P_{K,\mathbf{X}}:{}^{n_{x}}\to\real$$ is the so-called power function defined by

$$P_{K,\mathbf{X}}(\mathbf{x})=\sqrt{K(\mathbf{x},\mathbf{x})-K(\mathbf{X},\mathbf{x})^{\mathsf{T}}K(\mathbf{X},\mathbf{X})^{-1}K(\mathbf{X},\mathbf{x})}.$$ (2.4b)

See, e.g., [49, Thm 4.9] for a proof of the bound eq. 2.4a and [49, Prop. 4.11, Prop. 4.12] for the characterization eq. 2.4b of the power function. While this error bound is for the unregularized, fully interpolatory case, it is still useful in practice for the regularized case when $\gamma>0$ is small.

2.2 Vector-valued kernel interpolation

Kernel interpolation can be readily extended to vector-valued functions. The simple extension presented here, which is sufficient for our use case, is a special case of a more general extension relying on matrix-valued kernels (see, e.g., [36, 50]).

Consider the vector-valued function $\mathbf{v}:{}^{n_{x}}\to{}^{n_{y}}$, $n_{y}>1$. As before, let $\left\{\mathbf{x}_{j}\right\}_{j=1}^{m}\subset{}^{n_{x}}$ be pairwise distinct and suppose $\mathbf{y}_{j}=\mathbf{v}(\mathbf{x}_{j})\in{}^{n_{y}}$ for $j=1,\ldots,m$. Also let $v_{i}:{}^{n_{x}}\to\real$ denote the $i$-th component of $\mathbf{v}$, $y_{j,i}$ be the $i$-th component of $\mathbf{y}_{j}$, and define the input and output data matrices

$$\mathbf{X} =[~{}\mathbf{x}_{1}~{}~{}\cdots~{}~{}\mathbf{x}_{m}~{}]\in{}^{n_{x}\times m}, \mathbf{Y} =[~{}\mathbf{y}_{1}~{}~{}\cdots~{}~{}\mathbf{y}_{m}~{}]\in{}^{n_{y}\times m}.$$ (2.5)

We construct a vector-valued regularized kernel interpolant $\mathbf{s}_{\mathbf{v}}^{\gamma}$ by fitting scalar-valued kernel interpolants to each component $v_{i}$ of $\mathbf{v}$. Consequently, $\mathbf{s}_{\mathbf{v}}^{\gamma}$ is an element of the $n_{y}$-fold Cartesian product ${\cal H}_{K}^{n_{y}}\coloneqq{\cal H}_{K}\times\dots\times{\cal H}_{K}$, which has the inner product $\left\langle\mathbf{u},\mathbf{w}\right\rangle_{{\cal H}_{K}^{n_{y}}}=\sum_{i=1}^{n_{y}}\left\langle u_{i},w_{i}\right\rangle_{{\cal H}_{K}}$ for all $\mathbf{u}=(u_{1},\ldots,u_{n_{y}})\in{\cal H}_{K}^{n_{y}}$ and $\mathbf{w}=(w_{1},\ldots,w_{n_{y}})\in{\cal H}_{K}^{n_{y}}$. The regularized kernel interpolant constructed in this manner solves the optimization problem

$$\min_{\mathbf{s}\in{\cal H}_{K}^{n_{y}}}\;\sum_{j=1}^{m}\left\|\mathbf{y}_{j}-\mathbf{s}(\mathbf{x}_{j})\right\|_{2}^{2}+\gamma\left\|\mathbf{s}\right\|_{{\cal H}_{K}^{n_{y}}}^{2},$$ (2.6)

where $\left\|\cdot\right\|_{2}$ denotes the Euclidean $2$-norm and $\left\|\cdot\right\|_{{\cal H}_{K}^{n_{y}}}^{2}=\left\langle\cdot,\cdot\right\rangle_{{\cal H}_{K}^{n_{y}}}$. To see this, note that the objective function in eq. 2.6 can be rewritten as

$$\sum_{j=1}^{m}\left\|\mathbf{y}_{j}-\mathbf{s}(\mathbf{x}_{j})\right\|_{2}^{2}+\gamma\left\|\mathbf{s}\right\|_{{\cal H}_{K}^{n_{y}}}^{2}=\sum_{i=1}^{n_{y}}\left(\sum_{j=1}^{m}\left(y_{j,i}-s_{i}(\mathbf{x}_{j})\right)^{2}+\gamma\left\|s_{i}\right\|_{{\cal H}_{K}}^{2}\right),\quad\mathbf{s}(\mathbf{x})=\left[\begin{array}[]{c}s_{1}(\mathbf{x})\\ \vdots\\ s_{n_{y}}(\mathbf{x})\end{array}\right],$$

and therefore eq. 2.6 decouples into ${n_{y}}$ independent scalar-valued regularized interpolation problems:

$$\min_{\mathbf{s}_{i}\in{\cal H}_{K}}\;\sum_{j=1}^{m}(y_{j,i}-s_{i}(\mathbf{x}_{j}))^{2}+\gamma\left\|s_{i}\right\|_{{\cal H}_{K}}^{2},\qquad i=1,\dots,n_{y}.$$ (2.7)

Theorem 2.1 can then be applied to each subproblem to yield scalar-valued interpolants $s_{v_{i}}^{\gamma}$ of the form

$$s_{v_{i}}^{\gamma}(\mathbf{x})=\sum_{j=1}^{m}\omega_{i,j}K(\mathbf{x}_{j},\mathbf{x})={\boldsymbol{\omega}}_{i}^{\mathsf{T}}K(\mathbf{X},\mathbf{x}),$$ (2.8a)

where each coefficient vector $${\boldsymbol{\omega}}_{1},\ldots,{\boldsymbol{\omega}}_{n_{y}}\in{}^{m} m\times m$$ solves an linear system,

$$\big{(}K(\mathbf{X},\mathbf{X})+\gamma\mathbf{I}\big{)}{\boldsymbol{\omega}}_{i}=\left[\begin{array}[]{c}y_{i,1}\\ \vdots\\ y_{i,m}\end{array}\right],\quad i=1,\ldots,p.$$ (2.8e)

As before, $\gamma\geq 0$ is a given regularization parameter. An interpolant of $\mathbf{v}$ can then be defined by

$$\mathbf{s}_{\mathbf{v}}^{\gamma}(\mathbf{x})=[~{}s_{v_{1}}^{\gamma}(\mathbf{x})~{}~{}\cdots~{}~{}s_{v_{n_{y}}}^{\gamma}(\mathbf{x})~{}]^{\mathsf{T}}.$$

We summarize with the following corollary of Theorem 2.1 and a straightforward extension of Theorem 2.2.

Corollary 2.1 (Vector Representer Theorem).

The minimization problem eq. 2.6 has a solution of the form

$$\mathbf{s}_{\mathbf{v}}^{\gamma}(\mathbf{x})={\boldsymbol{\Omega}}^{\mathsf{T}}K(\mathbf{X},\mathbf{x}),$$ (2.9a)

where the coefficient matrix $${\boldsymbol{\Omega}}\in{}^{m\times{n_{y}}}$$ solves the linear system

$$\big{(}K(\mathbf{X},\mathbf{X})+\gamma\mathbf{I}\big{)}{\boldsymbol{\Omega}}=\mathbf{Y}^{\mathsf{T}}.$$ (2.9b)

Moreover, if $K$ is strictly positive definite, $s_{\mathbf{v}}^{\gamma}$ is the unique minimizer.

Corollary 2.2.

Let $\mathbf{v}\in{\cal H}_{K}^{n_{y}}$ and $\mathbf{M}\in{}^{{n_{y}}\times{n_{y}}}$ be a symmetric positive definite weighting matrix with Cholesky factorization $\mathbf{M}=\mathbf{L}\mathbf{L}^{\mathsf{T}}$. If $\mathbf{s}_{\mathbf{v}}^{0}\in{\cal H}_{K}^{n_{y}}$ is an (unregularized) vector-valued interpolant of $\mathbf{v}$ of the form eq. 2.9 corresponding to the pairwise distinct data $\left\{\mathbf{x}_{i}\right\}_{i=1}^{m}\subset{}^{n_{x}}$ and $\mathbf{y}_{i}=\mathbf{v}(\mathbf{x}_{i})\in{}^{n_{y}}$, then

$$\left\|\mathbf{v}(\mathbf{x})-\mathbf{s}_{\mathbf{v}}^{0}(\mathbf{x})\right\|_{\mathbf{M}}\leq P_{K,\mathbf{X}}(\mathbf{x})\left\|\mathbf{L}\right\|_{2}\left\|\mathbf{v}\right\|_{{\cal H}_{K}^{n_{y}}}\qquad\forall\;\mathbf{x}\in{}^{n_{x}}.$$ (2.10)

Proof.

Since $\mathbf{s}_{\mathbf{v}}^{0}$ interpolates $\mathbf{v}$ component-wise using the same kernel $K$ and interpolation points $\left\{\mathbf{x}_{i}\right\}_{i=1}^{m}$, applying Theorem 2.2 yields

$$\left\|\mathbf{v}(\mathbf{x})-\mathbf{s}_{\mathbf{v}}^{0}(\mathbf{x})\right\|_{\mathbf{M}}^{2}=\left\|\mathbf{L}^{\mathsf{T}}(\mathbf{v}(\mathbf{x})-\mathbf{s}_{\mathbf{v}}^{0}(\mathbf{x}))\right\|_{2}^{2} \leq\left\|\mathbf{L}\right\|_{2}^{2}\left\|\mathbf{v}(\mathbf{x})-\mathbf{s}_{\mathbf{v}}^{0}(\mathbf{x})\right\|_{2}^{2}=\left\|\mathbf{L}\right\|_{2}^{2}\sum_{i=1}^{n_{y}}|v_{i}(\mathbf{x})-s_{v_{i}}^{0}(\mathbf{x})|^{2}$$

$$\leq\left\|\mathbf{L}\right\|_{2}^{2}\sum_{i=1}^{n_{y}}P_{K,\mathbf{X}}(\mathbf{x})^{2}\left\|v_{i}\right\|_{{\cal H}_{K}}^{2}=P_{K,\mathbf{X}}(\mathbf{x})^{2}\left\|\mathbf{L}\right\|_{2}^{2}\left\|\mathbf{v}\right\|_{{\cal H}_{K}^{p}}^{2}.$$

∎

In Section 4, we use Corollary 2.1 to develop a strategy for constructing reduced-order models (ROMs) from data; Corollary 2.2 is used in Section 5 to derive a posteriori error estimates for these ROMs.

2.3 Kernel selection

Since a positive-definite kernel $K$ uniquely defines the RKHS ${\cal H}_{K}$, the choice of kernel determines what form an interpolant can take as well as the approximation power of the optimal interpolant. We argue for the use of different types of kernels depending on how much information is available about the function $\mathbf{v}:{}^{n_{x}}\to{}^{n_{y}}$ being interpolated.

2.3.1 Unknown structure: radial basis function kernels

If the structure of $\mathbf{v}$ is unknown, one effective choice is to generate the kernel using a radial basis function (RBF). These general-purpose kernels have the form

$$K(\mathbf{x},\mathbf{x}^{\prime})=\psi(\epsilon\left\|\mathbf{x}-\mathbf{x}^{\prime}\right\|_{2}),$$ (2.11a)

where $$\psi:{}_{\geq 0}\to\real \epsilon>0$$ and . Hence, RBF kernel interpolants are given by

$$\mathbf{s}_{\mathbf{v}}^{\gamma}(\mathbf{x})={\boldsymbol{\Omega}}^{\mathsf{T}}{\boldsymbol{\psi}}_{\!\epsilon}(\mathbf{x}),\qquad{\boldsymbol{\psi}}_{\!\epsilon}(\mathbf{x})=\begin{bmatrix}\psi(\epsilon\left\|\mathbf{x}_{1}-\mathbf{x}\right\|_{2})\\ \vdots\\ \psi(\epsilon\left\|\mathbf{x}_{m}-\mathbf{x}\right\|_{2})\end{bmatrix}\in{}^{m}.$$ (2.11b)

The so-called shape parameter $\epsilon$ is a hyperparameter that should be tuned to achieve optimal performance. Table 1 provides examples of commonly used RBF generator functions $\psi$. Note that the cost of evaluating an RBF kernel interpolant is $\mathcal{O}(m(n_{x}+n_{y}))$. A thorough discussion of the use of RBFs in kernel interpolation can be found in, e.g., [64].

Name	$\psi(x)$
Gaussian	$\exp(-x^{2})$
Basic Matérn	$\exp(-x)$
Inverse Quadratic	$(1+x^{2})^{-1}$
Inverse Multiquadric	$(1+x^{2})^{-1/2}$
Thin Plate Spline	$x^{2}\log(x)$

Table 1: Examples of common RBF kernel-generating functions.

2.3.2 Known structure: feature map kernels

If the structure of $\mathbf{v}$ is known, kernels induced by feature maps can often be used to endow the interpolant with matching structure, which can result in more accurate and interpretable approximations than when using general-purpose kernels. A feature map kernel can be written as

$$K(\mathbf{x},\mathbf{x}^{\prime})={\boldsymbol{\phi}}(\mathbf{x})^{\mathsf{T}}\mathbf{G}{\boldsymbol{\phi}}(\mathbf{x}^{\prime}),$$ (2.12a)

where $${\boldsymbol{\phi}}:{}^{n_{x}}\to{}^{n_{\phi}} \mathbf{G}\in{}^{n_{\phi}\times n_{\phi}}$$ is called the feature map and is a symmetric positive definite weighting matrix. It can be easily verified that feature map kernels are positive definite kernels (see, e.g., [49]). A feature map kernel results in a kernel interpolant of the form

$$\mathbf{s}_{\mathbf{v}}^{\gamma}(\mathbf{x}) ={\boldsymbol{\Omega}}^{\mathsf{T}}K(\mathbf{X},\mathbf{x})={\boldsymbol{\Omega}}^{\mathsf{T}}{\boldsymbol{\phi}}(\mathbf{X})^{\mathsf{T}}\mathbf{G}{\boldsymbol{\phi}}(\mathbf{x})=\underbrace{{\boldsymbol{\Omega}}^{\mathsf{T}}{\boldsymbol{\phi}}(\mathbf{X})^{\mathsf{T}}\mathbf{G}}_{\mathbf{C}}{\boldsymbol{\phi}}(\mathbf{x})=\mathbf{C}{\boldsymbol{\phi}}(\mathbf{x}),$$ (2.12b)

where ${\boldsymbol{\phi}}(\mathbf{X})\coloneqq[~{}{\boldsymbol{\phi}}(\mathbf{x}_{1})~{}~{}\cdots~{}~{}{\boldsymbol{\phi}}(\mathbf{x}_{m})~{}]\in{}^{n_{\phi}\times m}$. Importantly, the matrix $\mathbf{C}\in{}^{n_{y}\times n_{\phi}}$ can be computed once and reused repeatedly for online kernel evaluations. After constructing $\mathbf{C}$, the cost of evaluating a feature map kernel interpolant is therefore $\mathcal{O}(n_{\phi}n_{y})$, plus the expense of evaluating ${\boldsymbol{\phi}}$ once.

The advantage of feature map kernels is that one can imbue $\mathbf{s}_{\mathbf{v}}^{\gamma}$ with specific structure by designing the feature map ${\boldsymbol{\phi}}$ accordingly. For example, if

$${\boldsymbol{\phi}}(\mathbf{x})=\begin{bmatrix}\mathbf{x}\\ \mathbf{x}\otimes\mathbf{x}\end{bmatrix}\in{}^{n_{x}+n_{x}^{2}},$$ (2.13)

where $\otimes$ denotes the Kronecker product [59], then the associated kernel interpolant can be written as

$$\mathbf{s}_{\mathbf{v}}^{\gamma}(\mathbf{x})=\mathbf{C}_{1}\mathbf{x}+\mathbf{C}_{2}[\mathbf{x}\otimes\mathbf{x}],\qquad\mathbf{C}=[~{}\mathbf{C}_{1}~{}~{}\mathbf{C}_{2}~{}]\in{}^{n_{y}\times(n_{x}+n_{x}^{2})}.$$ (2.14)

Therefore, if it is known that $\mathbf{v}$ has linear-quadratic structure, then using a kernel induced by the feature map eq. 2.13 results in a kernel interpolant that has the same linear-quadratic structure.

2.3.3 Hybrid approach

For the purposes of model reduction, it is critical to keep the cost of evaluating the kernel interpolant low. The cost of evaluating an RBF kernel interpolant eq. 2.11b scales with the number of training samples $m$; by contrast, the cost of evaluating a feature map kernel interpolant eq. 2.12b is independent of $m$, but depends on the feature dimension $n_{\phi}$. If a feature map that fully specifies the desired structure requires a large $n_{\phi}$, one alternative is to define a new kernel that sums a less aggressive feature map kernel with an RBF kernel:

$$K(\mathbf{x},\mathbf{x}^{\prime})=c_{\phi}{\boldsymbol{\phi}}(\mathbf{x})^{\mathsf{T}}\mathbf{G}{\boldsymbol{\phi}}(\mathbf{x}^{\prime})+c_{\psi}\psi(\epsilon\left\|\mathbf{x}-\mathbf{x}^{\prime}\right\|_{2}),$$ (2.15a)

where $$c_{\phi},c_{\psi}\in{}_{>0} {\boldsymbol{\phi}} n_{\phi}$$ are positive weighting coefficients and is chosen to keep from being too large. The resulting kernel interpolant then has the form

$$\mathbf{s}_{\mathbf{v}}^{\gamma}(\mathbf{x}) ={\boldsymbol{\Omega}}^{\mathsf{T}}K(\mathbf{X},\mathbf{x})={\boldsymbol{\Omega}}^{\mathsf{T}}\big{(}c_{\phi}{\boldsymbol{\phi}}(\mathbf{X})^{\mathsf{T}}\mathbf{G}{\boldsymbol{\phi}}(\mathbf{x})+c_{\psi}{\boldsymbol{\psi}}_{\!\epsilon}(\mathbf{x})\big{)}=\mathbf{C}{\boldsymbol{\phi}}(\mathbf{x})+c_{\psi}{\boldsymbol{\Omega}}^{\mathsf{T}}{\boldsymbol{\psi}}_{\!\epsilon}(\mathbf{x}),$$ (2.15b)

where $\mathbf{C}=c_{\phi}{\boldsymbol{\Omega}}^{\mathsf{T}}{\boldsymbol{\phi}}(\mathbf{X})^{\mathsf{T}}\mathbf{G}$ now incorporates the weighting coefficient $c_{\phi}$. The idea is to use the feature map to incorporate dominant structure while relying on the RBF to approximate additional, potentially expensive terms. Note that this framework also applies to scenarios where the structure of $\mathbf{v}$ is only partially known.

As an example, consider the case where $\mathbf{v}$ is a quartic polynomial, i.e.,

$$\mathbf{v}(\mathbf{x})=\mathbf{A}_{1}\mathbf{x}+\mathbf{A}_{2}[\mathbf{x}\otimes\mathbf{x}]+\mathbf{A}_{3}[\mathbf{x}\otimes\mathbf{x}\otimes\mathbf{x}]+\mathbf{A}_{4}[\mathbf{x}\otimes\mathbf{x}\otimes\mathbf{x}\otimes\mathbf{x}],$$ (2.16)

where $\mathbf{A}_{1}\in{}^{n_{x}\times n_{x}}$, $\mathbf{A}_{2}\in{}^{n_{x}\times n_{x}^{2}}$, $\mathbf{A}_{3}\in{}^{n_{x}\times n_{x}^{3}}$, and $\mathbf{A}_{4}\in{}^{n_{x}\times n_{x}^{4}}$. One option is to fully capture the structure using a quartic feature map,

$${\boldsymbol{\phi}}(\mathbf{x})=\begin{bmatrix}\mathbf{x}\\ \mathbf{x}\otimes\mathbf{x}\\ \mathbf{x}\otimes\mathbf{x}\otimes\mathbf{x}\\ \mathbf{x}\otimes\mathbf{x}\otimes\mathbf{x}\otimes\mathbf{x}\end{bmatrix}\in{}^{n_{x}+n_{x}^{2}+n_{x}^{3}+n_{x}^{4}}.$$ (2.17)

However, evaluating the associated kernel interpolant costs $\mathcal{O}(n_{x}^{4}n_{y})$ operations, which is quite large for moderate $n_{x}$. Using the linear-quadratic feature map eq. 2.13 decreases $n_{\phi}$ from $n_{x}^{4}$ to $n_{x}^{2}$, and supplementing with an RBF kernel results in a kernel interpolant of the form

$$\mathbf{s}_{\mathbf{v}}^{\gamma}(\mathbf{x})=\mathbf{C}_{1}\mathbf{x}+\mathbf{C}_{2}[\mathbf{x}\otimes\mathbf{x}]+c_{\psi}{\boldsymbol{\Omega}}^{\mathsf{T}}{\boldsymbol{\psi}}_{\!\epsilon}(\mathbf{x}).$$ (2.18)

This interpolant does not fully represent the quartic structure of eq. 2.16, but it can be evaluated with only $\mathcal{O}((n_{x}^{2}+m)n_{y})$ operations. In this case, the RBF term acts as a type of closure term for structure that is not accounted for by the feature map.

Remark 2.1 (Input normalization).

In some cases, in particular when using high-order polynomial feature maps, the kernel matrix $K(\mathbf{X},\mathbf{X})$ used for determining ${\boldsymbol{\Omega}}$ may be poorly conditioned. Increasing the regularization constant $\gamma$ can improve the conditioning of the system eq. 2.9b, but this can also degrade the accuracy of the resulting kernel interpolant. Applying a normalization to the inputs can help remedy the situation: for any injective ${\boldsymbol{\nu}}:{}^{n_{x}}\to{}^{n_{x}}$, if $K$ is positive definite, then the function $K_{\boldsymbol{\nu}}:{}^{n_{x}}\times{}^{n_{x}}\to\real$ defined by

$$K_{\boldsymbol{\nu}}(\mathbf{x},\mathbf{x}^{\prime})=K({\boldsymbol{\nu}}(\mathbf{x}),{\boldsymbol{\nu}}(\mathbf{x}^{\prime}))$$ (2.19)

is also a positive-definite kernel function [49], and choosing ${\boldsymbol{\nu}}$ judiciously can improve the conditioning of $K_{\boldsymbol{\nu}}(\mathbf{X},\mathbf{X})$ compared to $K(\mathbf{X},\mathbf{X})$. A common choice is ${\boldsymbol{\nu}}(\mathbf{x})={\boldsymbol{\Sigma}}^{-1}(\mathbf{x}-\bar{\mathbf{x}})$, where ${\boldsymbol{\Sigma}}=\operatorname{diag}({\boldsymbol{\sigma}})\in{}^{n_{x}\times n_{x}}$ and $\bar{\mathbf{x}}\in{}^{n_{x}}$ with components

$$\sigma_{i}=\max_{j}(\mathbf{X}_{ij})-\min_{j}(\mathbf{X}_{ij}),\qquad\bar{x}_{i}=\min_{j}(\mathbf{X}_{ij}),\qquad i=1,\dots,d,$$ (2.20)

which maps the entries of each row of inputs to the interval $[0,1]$. In this case, an effective choice for the weighting matrix $\mathbf{G}$ in feature map kernels is $\mathbf{G}=(1/n_{\phi})\mathbf{I}$, where $n_{\phi}$ is the feature map dimension.

3 Intrusive projection-based model reduction

We now return to the model reduction setting and give a brief overview of intrusive projection-based ROMs, which inherit certain structure from the systems they emulate. Section 4 presents a non-intrusive alternative to intrusive model reduction for which kernel interpolation is the key ingredient and which can be designed to mimic the structure inheritance enjoyed by projection-based ROMs.

3.1 Generic projection-based reduced-order models

We consider high-dimensional systems of ordinary differential equations (ODEs) of the form

$$\frac{\textrm{d}}{\textrm{d}t}\mathbf{q}(t)=\mathbf{f}(\mathbf{q}(t)),\qquad\mathbf{q}(0)=\mathbf{q}_{0}({\boldsymbol{\mu}}),$$ (3.1)

where $\mathbf{q}:[0,T]\to{}^{n_{q}}$ is the state, $\mathbf{f}:{}^{n_{q}}\to{}^{n_{q}}$ governs the state evolution, $\mathbf{q}_{0}({\boldsymbol{\mu}})\in{}^{n_{q}}$ is the initial condition parameterized by ${\boldsymbol{\mu}}\in{\cal D}\subset{}^{n_{\mu}}$, and $T>0$ is the final desired simulation time. Models of this form often arise from semi-discretizations of time-dependent partial differential equations (PDEs), in which case the large state dimension $n_{q}$ corresponds to the fidelity of the underlying mesh. We call eq. 3.1 the full-order model (FOM).

A ROM for eq. 3.1 is a low-dimensional system of ODEs whose solution can be used to approximate the FOM state $\mathbf{q}(t)$. To that end, we consider a low-dimensional state approximation,

$$\mathbf{q}(t)\approx\mathbf{g}(\tilde{\mathbf{q}}(t)),$$ (3.2)

where $\mathbf{g}:{}^{r}\to{}^{n_{q}}$ and $\tilde{\mathbf{q}}:[0,T]\to{}^{r}$ is the reduced-order state, with $r\ll n_{q}$. The function $\mathbf{g}$ represents a decompression operation, mapping from reduced coordinates to the original high-dimensional space. We assume the existence of a corresponding compression map $\mathbf{h}:{}^{n_{q}}\to{}^{r}$, mapping high-dimensional states to reduced coordinates, such that $\mathbf{h}\circ\mathbf{g}$ is the identity. Importantly, $(\mathbf{g}\circ\mathbf{h})^{2}=\mathbf{g}\circ(\mathbf{h}\circ\mathbf{g})\circ\mathbf{h}=\mathbf{g}\circ\mathbf{h}$, i.e., $\mathbf{g}\circ\mathbf{h}$ is a projection. The evolution for the reduced state $\tilde{\mathbf{q}}(t)$ is then given by

$$\frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)=\frac{\textrm{d}}{\textrm{d}t}\mathbf{h}(\mathbf{g}(\tilde{\mathbf{q}}(t)))=\mathbf{h}^{\prime}(\mathbf{g}(\tilde{\mathbf{q}}(t)))\frac{\textrm{d}}{\textrm{d}t}\mathbf{g}(\tilde{\mathbf{q}}(t))\approx\mathbf{h}^{\prime}(\mathbf{g}(\tilde{\mathbf{q}}(t)))\mathbf{f}(\mathbf{g}(\tilde{\mathbf{q}}(t))),$$ (3.3)

in which $\mathbf{h}^{\prime}:{}^{n_{q}}\to{}^{n_{q}\times n_{q}}$ is the Jacobian of $\mathbf{h}$ and where the final step comes from inserting the approximation eq. 3.2 into the FOM eq. 3.1. The resulting system

$$\frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)=\mathbf{h}^{\prime}(\mathbf{g}(\tilde{\mathbf{q}}(t)))\mathbf{f}(\mathbf{g}(\tilde{\mathbf{q}}(t))),\qquad\tilde{\mathbf{q}}(0)=\mathbf{h}(\mathbf{q}_{0}({\boldsymbol{\mu}}))$$ (3.4)

is the projection-based ROM for eq. 3.1 corresponding to $\mathbf{g}$ and $\mathbf{h}$.

As written, eq. 3.4 is not highly practical because it involves mapping up to the high-dimensional state space, performing computations in that space, then compressing the results. However, for many common choices of $\mathbf{f}$, $\mathbf{g}$, and $\mathbf{h}$, eq. 3.4 simplifies in such a way that all computations can be performed in the reduced space, as we will demonstrate shortly.

3.2 Linear and quadratic dimension reduction

Classical model reduction methods typically define $\mathbf{g}$ and $\mathbf{h}$ as affine functions. In this work, we consider a slightly generalized approximation introduced in [26] and leveraged in [4, 18, 19, 54]: let

$$\mathbf{g}(\tilde{\mathbf{q}})=\bar{\mathbf{q}}+\mathbf{V}\tilde{\mathbf{q}}+\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}],$$ (3.5)

where $\bar{\mathbf{q}}\in{}^{n_{q}}$ is a fixed reference vector, $\mathbf{V}\in{}^{n_{q}\times r}$ has orthonormal columns, and $\mathbf{W}\in{}^{n_{q}\times r^{2}}$ satisfies $\mathbf{V}^{\mathsf{T}}\mathbf{W}=\bf 0$. This approximation defines an $r$-dimensional quadratic manifold embedded in ${}^{n_{q}}$. An appropriate compression map corresponding to eq. 3.5 is given by

$$\mathbf{h}(\mathbf{q}) =\mathbf{V}^{\mathsf{T}}(\mathbf{q}-\bar{\mathbf{q}}),$$ (3.6)

which has Jacobian $\mathbf{h}^{\prime}(\mathbf{q})=\mathbf{V}^{\mathsf{T}}$ and satisfies

$$\begin{aligned} (\mathbf{h}\circ\mathbf{g})(\tilde{\mathbf{q}})&=\mathbf{V}^{\mathsf{T}}\big{(}(\bar{\mathbf{q}}+\mathbf{V}\tilde{\mathbf{q}}+\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}])-\bar{\mathbf{q}}\big{)}=\mathbf{V}^{\mathsf{T}}\bar{\mathbf{q}}+\tilde{\mathbf{q}}-\mathbf{V}^{\mathsf{T}}\bar{\mathbf{q}}=\tilde{\mathbf{q}},\end{aligned}$$ (3.7)

since $\mathbf{V}^{\mathsf{T}}\mathbf{V}$ is the identity and $\mathbf{V}^{\mathsf{T}}$ annihilates $\mathbf{W}$. With $\mathbf{g}$ and $\mathbf{h}$ thus defined, the ROM eq. 3.4 becomes

$$\frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)=\tilde{\mathbf{f}}(\tilde{\mathbf{q}}(t))\coloneqq\mathbf{V}^{\mathsf{T}}\mathbf{f}\big{(}\bar{\mathbf{q}}+\mathbf{V}\tilde{\mathbf{q}}(t)+\mathbf{W}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)]\big{)},\qquad\tilde{\mathbf{q}}(0)=\mathbf{V}^{\mathsf{T}}(\mathbf{q}_{0}({\boldsymbol{\mu}})-\bar{\mathbf{q}}),$$ (3.8)

a system of $r\ll n_{q}$ ODEs defined by the function $\tilde{\mathbf{f}}:{}^{r}\to{}^{r}$.

The choices of $\bar{\mathbf{q}}$, $\mathbf{V}$, and $\mathbf{W}$ dictate the quality of the approximation eq. 3.5 and of the resulting ROM eq. 3.8. To make an informed selection, we assume access to a limited set of training data: given a set of training parameters ${\boldsymbol{\mu}}_{1},\ldots,{\boldsymbol{\mu}}_{M}\subset{\cal D}$ and observation times $t_{0},t_{1},\ldots,t_{n_{t}}$, let

$$\mathbf{q}_{k}^{(\ell)}\coloneqq\mathbf{q}(t_{k};{\boldsymbol{\mu}}_{\ell})\in{}^{n_{q}},\qquad\begin{aligned} \ell&=1,\ldots,M,\\ k&=0,1,\ldots,n_{t},\end{aligned}$$ (3.9)

which are snapshots of the full-order state solution to the FOM eq. 3.1. The reference vector $\bar{\mathbf{q}}$ is usually set to zero, the initial condition at a fixed training parameter value, or the average snapshot, i.e.,

$$\bar{\mathbf{q}}=\frac{1}{Mn_{t}}\sum_{\ell=1}^{M}\sum_{k=0}^{n_{t}}\mathbf{q}_{k}^{(\ell)}.$$ (3.10)

The model reduction framework developed in Section 4 applies for any $\bar{\mathbf{q}}$, $\mathbf{V}$ and $\mathbf{W}$ such that $\mathbf{V}^{\mathsf{T}}\mathbf{V}=\mathbf{I}$ and $\mathbf{V}^{\mathsf{T}}\mathbf{W}=\bf 0$, but we focus on two best-practice cases.

First, if $\mathbf{W}=\bf 0$, the manifold defined by $\mathbf{g}$ has no curvature and reduces to an affine subspace (or a linear subspace if $\bar{\mathbf{q}}=\bf 0$) of ${}^{n_{q}}$. In this case, we select $\mathbf{V}$ using proper orthogonal decomposition (POD) [9, 21, 56]. Define

$$\mathbf{Q}\coloneqq\begin{bmatrix}(\mathbf{q}_{0}^{(1)}-\bar{\mathbf{q}})&\cdots&(\mathbf{q}_{n_{t}}^{(1)}-\bar{\mathbf{q}})&(\mathbf{q}_{0}^{(2)}-\bar{\mathbf{q}})&\cdots&(\mathbf{q}_{n_{t}}^{(M)}-\bar{\mathbf{q}})\end{bmatrix}\in{}^{n_{q}\times M(n_{t}+1)},$$ (3.11)

the matrix of snapshots stacked column-wise and shifted by the reference snapshot. The rank-$r$ POD basis matrix $\mathbf{V}$ is given by the first $r$ left singular vectors of $\mathbf{Q}$. With this choice, $\mathbf{g}\circ\mathbf{h}$ is the optimal $r$-dimensional approximator for the (shifted) training snapshots in an $L^{2}$ sense.

Second, to construct a nonzero $\mathbf{W}$, we use the greedy-optimal quadratic manifold (QM) approach of [54]. This method iteratively selects the columns of $\mathbf{V}$ from the left singular vectors of $\mathbf{Q}$ and solves a least-squares problem to determine $\mathbf{W}$,

$$\min_{\mathbf{v}_{i},\mathbf{W}}\;\left\|(\mathbf{I}-\mathbf{V}\mathbf{V}^{\mathsf{T}})\mathbf{Q}-\mathbf{W}[\mathbf{V}^{\mathsf{T}}\mathbf{Q}\odot\mathbf{V}^{\mathsf{T}}\mathbf{Q}]\right\|_{F}^{2}+\rho\left\|\mathbf{W}\right\|_{F}^{2},$$ (3.12)

where $\mathbf{v}_{i}$ is the final column of $\mathbf{V}$ and all other columns are fixed from previous iterations. Here, $\odot$ indicates the Khatri–Rao (column-wise Kronecker) product, and $\rho\geq 0$ is a scalar regularization parameter. Traditional POD always sets $\mathbf{v}_{i}$ to the $i$-th left singular vector, but here each $\mathbf{v}_{i}$ can be chosen from among any of the left singular vectors that have not yet been selected, which can lead to substantial accuracy gains.

Remark 3.1 (Kronecker redundancy).

The product $\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}$ contains redundant terms, i.e., $\tilde{q}_{i}\tilde{q}_{j}$ appears twice for each $i\neq j$, which means two columns of $\mathbf{W}\in{}^{r\times r^{2}}$ act on the same quadratic state interaction in the product $\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]$. As a consequence, the learning problem eq. 3.12 has infinitely many solutions. In practice, this issue is avoided by replacing $\otimes$ in eq. 3.5 with a compressed Kronecker product $\tilde{\otimes}$, defined by

$$\tilde{\mathbf{q}}\,\tilde{\otimes}\,\tilde{\mathbf{q}}\coloneqq\left[~{}\tilde{q}_{1}^{2}~{}~{}\tilde{q}_{1}\tilde{q}_{2}~{}~{}\tilde{q}_{2}^{2}~{}~{}\dots~{}~{}\tilde{q}_{r-1}\tilde{q}_{r}~{}~{}\tilde{q}_{r}^{2}~{}\right]^{\mathsf{T}}\in{}^{r(r+1)/2},$$ (3.13)

which leads to a matrix $\tilde{\mathbf{W}}\in{}^{r\times r(r+1)/2}$ such that $\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]=\tilde{\mathbf{W}}[\tilde{\mathbf{q}}\,\tilde{\otimes}\,\tilde{\mathbf{q}}]$ for all $\tilde{\mathbf{q}}\in{}^{r}$. Then, if $\odot$ applies $\tilde{\otimes}$ column-wise, the optimization eq. 3.12 has a unique solution. Similar adjustments can be made for higher-order Kronecker products.

3.3 Intrusive reduced-order models for quadratic systems

The key observation in projection-based model reduction is that projection preserves certain structure. Suppose that the function $\mathbf{f}$ defining the dynamics of the FOM eq. 3.1 has linear-quadratic structure, i.e.,

$$\frac{\textrm{d}}{\textrm{d}t}\mathbf{q}(t)=\mathbf{f}(\mathbf{q}(t))\coloneqq\mathbf{A}\mathbf{q}(t)+\mathbf{H}[\mathbf{q}(t)\otimes\mathbf{q}(t)],$$ (3.14)

where $\mathbf{A}\in{}^{n_{q}\times n_{q}}$, $\mathbf{H}\in{}^{n_{q}\times n_{q}^{2}}$. It is assumed that $\mathbf{H}$ is symmetric in the sense that $\mathbf{H}[\mathbf{q}\otimes\mathbf{p}]=\mathbf{H}[\mathbf{p}\otimes\mathbf{q}]$ for all $\mathbf{q},\mathbf{p}\in{}^{n_{q}}$. Models with quadratic structure arise from quadratic PDEs, but can also result from applying lifting transformations to models with other structure [30, 43]. With a linear state approximation ($\bar{\mathbf{q}}=\bf 0$ and $\mathbf{W}=\bf 0$), the ROM eq. 3.8 can be written as

$$\frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)=\tilde{\mathbf{A}}\tilde{\mathbf{q}}(t)+\tilde{\mathbf{H}}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)],\qquad\tilde{\mathbf{q}}(0)=\mathbf{V}^{\mathsf{T}}\mathbf{q}_{0}({\boldsymbol{\mu}}),$$ (3.15)

in which $\tilde{\mathbf{A}}=\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{V}\in{}^{r\times r}$ and $\tilde{\mathbf{H}}=\mathbf{V}^{\mathsf{T}}\mathbf{H}[\mathbf{V}\otimes\mathbf{V}]\in{}^{r\times r^{2}}$. Constructing eq. 3.15 is an intrusive process because $\tilde{\mathbf{A}}$ and $\tilde{\mathbf{H}}$ depend explicitly on $\mathbf{A}$ and $\mathbf{H}$; however, we need not have access to $\mathbf{A}$ and $\mathbf{H}$ to observe that the quadratic structure is preserved.

In the QM case ($\mathbf{W}\neq\bf 0$, but still with $\bar{\mathbf{q}}=\bf 0$ for convenience), the ROM eq. 3.8 has quartic dynamics,

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)&=\tilde{\mathbf{A}}\tilde{\mathbf{q}}(t)+\tilde{\mathbf{H}}_{2}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)]+\tilde{\mathbf{H}}_{3}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)]+\tilde{\mathbf{H}}_{4}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)],\\ \tilde{\mathbf{q}}(0)&=\mathbf{V}^{\mathsf{T}}\mathbf{q}_{0}({\boldsymbol{\mu}}),\end{aligned}$$ (3.16)

where $\tilde{\mathbf{A}}=\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{V}$, $\tilde{\mathbf{H}}_{2}=\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{W}+\mathbf{V}^{\mathsf{T}}\mathbf{H}(\mathbf{V}\otimes\mathbf{V})$, $\tilde{\mathbf{H}}_{3}=\mathbf{V}^{\mathsf{T}}\mathbf{H}(\mathbf{V}\otimes\mathbf{W}+\mathbf{W}\otimes\mathbf{V})$, and $\tilde{\mathbf{H}}_{4}=\mathbf{V}^{\mathsf{T}}\mathbf{H}(\mathbf{W}\otimes\mathbf{W})$. Again, this process is intrusive, but the key result is that if one knows the structure of the FOM dynamics, one can also deduce the structure of the projection-based ROM. See Appendix A for the case when $\bar{\mathbf{q}}\neq\bf 0$, in which a constant term appears in the reduced dynamics.

4 Non-intrusive model reduction via kernel interpolation

This section leverages regularized kernel interpolation to construct ROMs akin to eq. 3.8, denoted

$$\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t)=\hat{\mathbf{f}}(\hat{\mathbf{q}}(t)),\qquad\hat{\mathbf{q}}(0)=\mathbf{V}^{\mathsf{T}}(\mathbf{q}_{0}({\boldsymbol{\mu}})-\bar{\mathbf{q}}),$$ (4.1)

where $\hat{\mathbf{q}}:[0,T]\to{}^{r}$ and $\hat{\mathbf{f}}:{}^{r}\to{}^{r}$. The structure of $\hat{\mathbf{f}}$ can be informed by intrusive projection, but unlike projection, defining $\hat{\mathbf{f}}$ through kernel interpolation does not require access to FOM operators such as $\mathbf{A}$ or $\mathbf{H}$ in eq. 3.14. We use the notation $\hat{\cdot}$ to mark non-intrusive objects and differentiate from intrusive objects, which are marked with $\tilde{\cdot}$.

4.1 Kernel reduced-order models

We pose the problem of learning an appropriate $\hat{\mathbf{f}}$ for the ROM eq. 4.1 as a regression, which requires data for the state $\hat{\mathbf{q}}(t)$ and its time derivative. For the former, we reduce the FOM snapshots eq. 3.9 using the compression map $\mathbf{h}$, that is,

$$\hat{\mathbf{q}}_{k}^{(\ell)}\coloneqq\mathbf{h}(\mathbf{q}_{k}^{(\ell)})=\mathbf{V}^{\mathsf{T}}(\mathbf{q}_{k}^{(\ell)}-\bar{\mathbf{q}})\in{}^{r}.$$ (4.2)

If the time step between observations is sufficiently small, an accurate approximation for the time derivatives of the state can be computed from finite differences of the reduced states, for example,

$$\dot{\hat{\mathbf{q}}}_{k}^{(\ell)}\coloneqq\frac{\hat{\mathbf{q}}_{k}^{\ell}-\hat{\mathbf{q}}_{k-1}^{\ell}}{t_{k}-t_{k-1}}\approx\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t)\big{|}_{t=t_{k}}.$$ (4.3)

The ROM function $\hat{\mathbf{f}}$ can then be defined as the solution to a minimization problem,

$$\hat{\mathbf{f}}=\underset{\mathbf{s}\in S}{\arg\min}\;\sum_{\ell=1}^{M}\sum_{k=0}^{n_{t}}\left\|\dot{\hat{\mathbf{q}}}_{k}^{(\ell)}-\mathbf{s}(\hat{\mathbf{q}}_{k}^{(\ell)})\right\|_{2}^{2}+R(\mathbf{s}),$$ (4.4)

where $S$ is some set of functions and $R:S\to{}_{\geq 0}$ is a regularization function.

The generic minimization eq. 4.4 encompasses several data-driven approaches which each use different choices for the space $S$ and the regularizer $R$. By defining a kernel $K$ and an associated RKHS $S=\mathcal{H}_{K}^{r}$, and setting $R(\mathbf{s})=\gamma\|\mathbf{s}\|_{\mathcal{H}_{K}^{r}}$, we obtain a vector regularized kernel interpolation problem,

$$\hat{\mathbf{f}}=\underset{\mathbf{s}\in{\cal H}_{K}^{r}}{\arg\min}\;\sum_{\ell=1}^{M}\sum_{k=0}^{n_{t}}\left\|\dot{\hat{\mathbf{q}}}_{k}^{(\ell)}-\mathbf{s}(\hat{\mathbf{q}}_{k}^{(\ell)})\right\|_{2}^{2}+\gamma\left\|\mathbf{s}\right\|_{{\cal H}_{K}^{r}}^{2},$$ (4.5)

which is eq. 2.6 with $\mathbf{x}_{j}=\hat{\mathbf{q}}_{k}^{(\ell)}$, $\mathbf{y}_{j}=\dot{\hat{\mathbf{q}}}_{k}^{(\ell)}$ and $m=M(n_{t}+1)$ after some minor reindexing for $k$ and $\ell$. Corollary 2.1 gives an explicit representation for $\hat{\mathbf{f}}$, resulting in the ROM

$$\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t)=\hat{\mathbf{f}}(\hat{\mathbf{q}}(t))\coloneqq{\boldsymbol{\Omega}}^{\mathsf{T}}K(\hat{\mathbf{Q}},\hat{\mathbf{q}}(t)),\qquad\hat{\mathbf{q}}(0)=\mathbf{V}^{\mathsf{T}}(\mathbf{q}_{0}({\boldsymbol{\mu}})-\bar{\mathbf{q}}),$$ (4.6a)

where $${\boldsymbol{\Omega}}\in{}^{M(n_{t}+1)\times r}$$ solves the linear system

$$\big{(}K(\hat{\mathbf{Q}},\hat{\mathbf{Q}})+\gamma\mathbf{I}\big{)}{\boldsymbol{\Omega}}=\hat{\mathbf{Z}}^{\mathsf{T}},$$ (4.6b)

with interpolation input and output matrices

$$\begin{aligned} \hat{\mathbf{Q}}&=\begin{bmatrix}\hat{\mathbf{q}}_{0}^{(1)}&\cdots&\hat{\mathbf{q}}_{n_{t}}^{(1)}&\hat{\mathbf{q}}_{0}^{(2)}&\cdots&\hat{\mathbf{q}}_{n_{t}}^{(M)}\end{bmatrix}\in{}^{r\times M(n_{t}+1)},\\ \hat{\mathbf{Z}}&=\begin{bmatrix}\dot{\hat{\mathbf{q}}}_{0}^{(1)}&\cdots&\dot{\hat{\mathbf{q}}}_{n_{t}}^{(1)}&\dot{\hat{\mathbf{q}}}_{0}^{(2)}&\cdots&\dot{\hat{\mathbf{q}}}_{n_{t}}^{(M)}\end{bmatrix}\in{}^{r\times M(n_{t}+1)}.\end{aligned}$$ (4.6c)

Note that the cost of evaluating $\hat{\mathbf{f}}$ is $\mathcal{O}(rMn_{t})$, plus the cost of evaluating the kernel term $K(\hat{\mathbf{Q}},\hat{\mathbf{q}}(t))$.

Remark 4.1.

If the time derivatives of the FOM snapshots $\frac{\textrm{d}}{\textrm{d}t}\mathbf{q}_{k}^{(\ell)}$ are available, the time derivatives of the reduced state can instead be computed as

$$\dot{\hat{\mathbf{q}}}_{k}^{(\ell)}=\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}_{k}^{(\ell)} \approx\frac{\textrm{d}}{\textrm{d}t}(\mathbf{h}(\mathbf{q}_{k}^{(\ell)}))=\mathbf{h}^{\prime}(\mathbf{q}_{k}^{(\ell)})\frac{\textrm{d}}{\textrm{d}t}\mathbf{q}_{k}^{(\ell)}=\mathbf{V}^{\mathsf{T}}\frac{\textrm{d}}{\textrm{d}t}\mathbf{q}_{k}^{(\ell)}.$$ (4.7)

4.2 Specifying structure through kernel design

We now employ the observations of Section 2.3 to endow Kernel ROMs with structure. If the structure of the FOM function $\mathbf{f}$ is unknown, an RBF kernel is a reasonable general-purpose choice for $K$. However, if the structure of $\mathbf{f}$ is known, a feature map kernel can be employed so that the resulting $\hat{\mathbf{f}}$ has the same structure of the intrusive projection-based ROM function $\tilde{\mathbf{f}}$. This is best shown by example.

Consider again the quartic QM ROM eq. 3.16. Using the quartic feature map ${\boldsymbol{\phi}}$ of eq. 2.17 (with $\mathbf{x}=\hat{\mathbf{q}}$ and $n_{x}=r$) to define a feature map kernel $K(\hat{\mathbf{q}},\hat{\mathbf{q}}^{\prime})={\boldsymbol{\phi}}(\hat{\mathbf{q}})^{\mathsf{T}}\mathbf{G}{\boldsymbol{\phi}}(\hat{\mathbf{q}}^{\prime})$, the Kernel ROM eq. 4.6 takes the form

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t)&=\mathbf{C}{\boldsymbol{\phi}}(\hat{\mathbf{q}}(t))\\ &=\hat{\mathbf{A}}\hat{\mathbf{q}}(t)+\hat{\mathbf{H}}_{2}[\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)]+\hat{\mathbf{H}}_{3}[\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)]+\hat{\mathbf{H}}_{4}[\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)],\end{aligned}$$ (4.8)

in which $\mathbf{C}={\boldsymbol{\Omega}}^{\mathsf{T}}{\boldsymbol{\phi}}(\hat{\mathbf{Q}})^{\mathsf{T}}\mathbf{G}=[~{}\hat{\mathbf{A}}~{}~{}\hat{\mathbf{H}}_{2}~{}~{}\hat{\mathbf{H}}_{3}~{}~{}\hat{\mathbf{H}}_{4}~{}]$. This ROM has the same dynamical structure as eq. 3.16 but can be constructed non-intrusively. The structure can be tailored by adjusting the feature map: if the FOM eq. 3.14 is linear ($\mathbf{H}=\bf 0$), then the intrusive QM ROM eq. 3.16 simplifies to a quadratic form,

$$\frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)=\tilde{\mathbf{A}}\tilde{\mathbf{q}}(t)+\tilde{\mathbf{H}}_{1}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)],\qquad\tilde{\mathbf{H}}_{1}=\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{W},$$ (4.9)

which can be mimicked by a Kernel ROM by employing a linear-quadratic feature map as in eq. 2.13.

Remark 4.2 (Input terms).

Kernel ROMs can be designed to account for known input terms by including them in the feature map. Suppose we wish to construct a ROM with the structure

$$\frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)=\tilde{\mathbf{f}}(\tilde{\mathbf{q}}(t))+\tilde{\mathbf{b}}(\mathbf{u}(t)),$$ (4.10)

where $\tilde{\mathbf{b}}:{}^{n_{u}}\to{}^{r}$ and $\mathbf{u}:[0,T]\to{}^{n_{u}}$ model, for example, time-varying boundary conditions or forcing terms. In this case, we can construct feature maps ${\boldsymbol{\phi}}_{q}$ and ${\boldsymbol{\phi}}_{u}$ which aim to emulate the structures of $\tilde{\mathbf{f}}$ and $\tilde{\mathbf{b}}$, respectively, and define a concatenated feature map

$${\boldsymbol{\phi}}(\hat{\mathbf{q}};\mathbf{u})=\begin{bmatrix}{\boldsymbol{\phi}}_{q}(\hat{\mathbf{q}})\\ {\boldsymbol{\phi}}_{u}(\mathbf{u})\end{bmatrix}.$$ (4.11)

The resulting Kernel ROM has the form

$$\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t)=\mathbf{C}_{q}{\boldsymbol{\phi}}_{q}(\hat{\mathbf{q}}(t))+\mathbf{C}_{u}{\boldsymbol{\phi}}_{u}(\mathbf{u}(t)),$$ (4.12)

whose structure can be tailored to that of eq. 4.10 by designing ${\boldsymbol{\phi}}_{q}$ and ${\boldsymbol{\phi}}_{u}$ appropriately.

As discussed in Section 2.3, feature map kernels can lead to cost savings over generic kernels. Let $n_{\phi}$ be the dimension of the feature map, i.e., ${\boldsymbol{\phi}}(\hat{\mathbf{q}})\in{}^{n_{\phi}}$. Because the matrix $\mathbf{C}\in{}^{r\times n_{\phi}}$ can be computed once and reused, the cost of evaluating the ROM function $\hat{\mathbf{f}}$ online is $\mathcal{O}(rn_{\phi})$. Hence, if $n_{\phi}<Mn_{t}$, a feature map kernel is less expensive to evaluate than a generic kernel. If $n_{\phi}>Mn_{t}$ (e.g., due to a moderate reduced state dimension $r$), it can be beneficial to reduce $n_{\phi}$ and add a more generic element to the kernel to compensate. For instance, in place of the quartic ROM eq. 4.8, we may choose a quadratic feature map and add an RBF term to account for the cubic and quartic nonlinearities, resulting in a ROM of the form

$$\frac{\textrm{d}}{\textrm{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)]+c_{\psi}{\boldsymbol{\Omega}}^{\mathsf{T}}{\boldsymbol{\psi}}_{\!\epsilon}(\hat{\mathbf{q}}(t)),$$ (4.13)

where ${\boldsymbol{\psi}}_{\!\epsilon}:{}^{r}\to{}^{M(n_{t}+1)}$ is as in eq. 2.11 and $c_{\psi}>0$ is a weighting coefficient as in eq. 2.15b. We test ROMs with this hybrid structure in Section 6. Note that this strategy can also apply to cases where the desired ROM structure is only partially known or representable by a feature map kernel.

4.3 Comparison to operator inference

Our kernel-based method is philosophically similar to the operator inference (OpInf) framework pioneered in [40], with a few key differences. Like our method, OpInf stipulates the form of a ROM based on structure that arises from intrusive projection, and the objects defining the ROM are learned from a regression problem of reduced states and corresponding time derivatives. However, the learning problems in each approach use different candidate function spaces and regularizers, resulting in different ROMs even when the same training data and model structure are used for both procedures.

Generally speaking, OpInf constructs ROMs of the form

$$\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t)=\hat{\mathbf{O}}{\boldsymbol{\phi}}(\hat{\mathbf{q}}(t))$$ (4.14a)

for a specified feature map $${\boldsymbol{\phi}}:{}^{r}\to{}^{n_{\phi}}$$ by solving the regularized residual minimization problem

$$\min_{\hat{\mathbf{O}}\in{}^{r\times n_{\phi}}}\;\sum_{\ell=1}^{M}\sum_{k=0}^{n_{t}}\left\|\dot{\hat{\mathbf{q}}}_{k}^{(\ell)}-\hat{\mathbf{O}}{\boldsymbol{\phi}}(\hat{\mathbf{q}}_{k}^{(\ell)})\right\|_{2}^{2}+\|{\boldsymbol{\Gamma}}\hat{\mathbf{O}}^{\mathsf{T}}\|_{F}^{2},$$ (4.14b)

where ${\boldsymbol{\Gamma}}\in{}^{n_{\phi}\times n_{\phi}}$. This is the generic learning problem eq. 4.4 with the function space $S$ given by

$$S_{{\boldsymbol{\phi}}}=\left\{\mathbf{s}:{}^{r}\to{}^{r}\;:\;\mathbf{s}(\hat{\mathbf{q}})=\hat{\mathbf{O}}{\boldsymbol{\phi}}(\hat{\mathbf{q}})\quad\textup{for some}\quad\hat{\mathbf{O}}\in{}^{r\times n_{\phi}}\right\},$$ (4.15)

and where $R$ is a Tikhonov regularizer. The so-called operator matrix $\hat{\mathbf{O}}$ satisfies the linear system

$$(\mathbf{D}^{\mathsf{T}}\mathbf{D}+{\boldsymbol{\Gamma}}^{\mathsf{T}}{\boldsymbol{\Gamma}})\hat{\mathbf{O}}^{\mathsf{T}}=\mathbf{D}^{\mathsf{T}}\hat{\mathbf{Z}}^{\mathsf{T}},$$ (4.16)

where $\mathbf{D}={\boldsymbol{\phi}}(\hat{\mathbf{Q}})^{\mathsf{T}}$ and $\hat{\mathbf{Q}}$ and $\hat{\mathbf{Z}}$ are the training data matrices in eq. 4.6c. As with our kernel-based approach, the feature map is chosen to emulate the structure of a projection-based ROM. For example, the OpInf regression to learn a linear-quadratic ROM of the form eq. 3.15 is given by

$$\min_{\hat{\mathbf{A}}\in{}^{r\times r},\hat{\mathbf{H}}\in{}^{r\times r^{2}}}\;\sum_{\ell=1}^{M}\sum_{k=0}^{n_{t}}\left\|\dot{\hat{\mathbf{q}}}_{k}^{(\ell)}-\left(\hat{\mathbf{A}}\hat{\mathbf{q}}_{k}^{(\ell)}+\hat{\mathbf{H}}[\hat{\mathbf{q}}_{k}^{(\ell)}\otimes\hat{\mathbf{q}}_{k}^{(\ell)}]\right)\right\|_{2}^{2}+\left\|{\boldsymbol{\Gamma}}[~{}\hat{\mathbf{A}}~{}~{}\hat{\mathbf{H}}~{}]^{\mathsf{T}}\right\|_{F}^{2},$$ (4.17)

and the solution $\hat{\mathbf{O}}=[~{}\hat{\mathbf{A}}~{}~{}\hat{\mathbf{H}}~{}]$ satisfies eq. 4.16 with $\mathbf{D}=[~{}\hat{\mathbf{Q}}^{\mathsf{T}}~{}~{}(\hat{\mathbf{Q}}\odot\hat{\mathbf{Q}})^{\mathsf{T}}~{}]\in{}^{M(n_{t}+1)\times(r+r^{2})}$. The underlying feature map in this case is the linear-quadratic map eq. 2.13. In practice, the compressed Kronecker product of Remark 3.1 is used so that eq. 4.17 has a unique solution.

For a kernel ROM with the kernel specified entirely by a feature map, the resulting ROM can be expressed in terms of the training data and the feature map as

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t)&=\hat{\mathbf{Z}}\big{(}K(\hat{\mathbf{Q}},\hat{\mathbf{Q}})+\gamma\mathbf{I}\big{)}^{-1}K(\hat{\mathbf{Q}},\hat{\mathbf{q}}(t))\\ &=\underbrace{\hat{\mathbf{Z}}\big{(}{\boldsymbol{\phi}}(\hat{\mathbf{Q}})^{\mathsf{T}}\mathbf{G}{\boldsymbol{\phi}}(\hat{\mathbf{Q}})+\gamma\mathbf{I}\big{)}^{-1}{\boldsymbol{\phi}}(\hat{\mathbf{Q}})^{\mathsf{T}}\mathbf{G}}_{\mathbf{C}}{\boldsymbol{\phi}}(\hat{\mathbf{q}}(t)),\end{aligned}$$ (4.18)

whereas the OpInf ROM with the same training data and feature map is given by

$$\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t) =\underbrace{\hat{\mathbf{Z}}{\boldsymbol{\phi}}(\hat{\mathbf{Q}})^{\mathsf{T}}\big{(}{\boldsymbol{\phi}}(\hat{\mathbf{Q}}){\boldsymbol{\phi}}(\hat{\mathbf{Q}})^{\mathsf{T}}+{\boldsymbol{\Gamma}}^{\mathsf{T}}{\boldsymbol{\Gamma}}\big{)}^{-1}}_{\hat{\mathbf{O}}}{\boldsymbol{\phi}}(\hat{\mathbf{q}}(t)).$$ (4.19)

These models share the same nonlinear structure due to the final term ${\boldsymbol{\phi}}(\hat{\mathbf{q}}(t))$, but the coefficients on the feature map differ: the Kernel ROM coefficient matrix $\mathbf{C}$ solves the $M(n_{t}+1)\times M(n_{t}+1)$ linear system eq. 4.6b, while the solution $\hat{\mathbf{O}}$ to the OpInf regression satisfies an $n_{\phi}\times n_{\phi}$ linear system eq. 4.16. Furthermore, OpInf is in general restricted to the feature map formulation eq. 4.14, though it has in some cases been augmented with additional nonlinear terms through, e.g., the discrete empirical interpolation method [8]; by contrast, Kernel ROMs can be designed to have general nonlinear (RBF) structure or hybrid structure such as in eq. 4.13, depending on the choice of kernel. Finally, establishing error bounds is an open problem for OpInf ROMs, whereas Kernel ROMs inherit properties from the underlying RKHS which lead to error estimates.

5 Error estimates

We now derive several a posteriori error estimates for Kernel ROMs that relate the FOM solution $\mathbf{q}(t)$, the intrusive ROM solution $\tilde{\mathbf{q}}(t)$, and the Kernel ROM solution $\hat{\mathbf{q}}(t)$. These results require three main ingredients: the so-called local logarithmic Lipschitz constant, a Grönwall-type inequality, and standard error results for kernel interpolants. In this section, $\mathbf{M}\in{}^{n_{x}\times n_{x}}$ denotes a symmetric positive definite weighting matrix with Cholesky factorization $\mathbf{M}=\mathbf{L}\mathbf{L}^{\mathsf{T}}$. The $\mathbf{M}$-weighted inner product and norm are denoted with $\left\langle\mathbf{x},\mathbf{y}\right\rangle_{\mathbf{M}}\coloneqq\mathbf{x}^{\mathsf{T}}\mathbf{M}\mathbf{y}=\left\langle\mathbf{L}^{\mathsf{T}}\mathbf{x},\mathbf{L}^{\mathsf{T}}\mathbf{y}\right\rangle$ and $\left\|\mathbf{x}\right\|_{\mathbf{M}}\coloneqq\left\langle\mathbf{x},\mathbf{x}\right\rangle_{\mathbf{M}}^{1/2}=\|\mathbf{L}^{\mathsf{T}}\mathbf{x}\|_{2}$, respectively.

5.1 Preliminaries

We begin with the definition of the local logarithmic Lipschitz constant. The reader is directed to, e.g., [57, 63] for a more complete overview.

Definition 5.1.

For a function $\mathbf{b}:{}^{n_{x}}\to{}^{n_{x}}$, the local logarithmic Lipschitz constant at $\mathbf{x}\in{}^{n_{x}}$ with respect to $\mathbf{M}$ is defined as

$$\Lambda_{\mathbf{M}}[\mathbf{b}](\mathbf{x})\coloneqq\sup_{\mathbf{z}\in{}^{d}}\frac{\left\langle\mathbf{z}-\mathbf{x},\mathbf{b}(\mathbf{z})-\mathbf{b}(\mathbf{x})\right\rangle_{\mathbf{M}}}{\left\|\mathbf{z}-\mathbf{x}\right\|_{\mathbf{M}}^{2}}.$$ (5.1)

The local logarithmic Lipschitz constant can be seen as a nonlinear generalization of the logarithmic norm of a matrix.

Definition 5.2 (Logarithmic norm).

The logarithmic norm of a matrix $\mathbf{B}\in{}^{n_{x}\times n_{x}}$ with respect to $\mathbf{M}$ is defined as

$$\lambda_{\mathbf{M}}(\mathbf{B})\coloneqq\sup_{\mathbf{x}\in{}^{n_{x}}}\frac{\left\langle\mathbf{x},\mathbf{B}\mathbf{x}\right\rangle_{\mathbf{M}}}{\left\|\mathbf{x}\right\|_{\mathbf{M}}^{2}}=\max\sigma\left(\frac{1}{2}\left(\tilde{\mathbf{B}}+\tilde{\mathbf{B}}^{\mathsf{T}}\right)\right),$$ (5.2)

where $\sigma(\frac{1}{2}(\tilde{\mathbf{B}}+\tilde{\mathbf{B}}^{\mathsf{T}}))$ is the spectrum of $\frac{1}{2}(\tilde{\mathbf{B}}+\tilde{\mathbf{B}}^{\mathsf{T}})$ and $\tilde{\mathbf{B}}=\mathbf{L}^{\mathsf{T}}\mathbf{B}\mathbf{L}^{-\mathsf{T}}.$

If $\mathbf{b}$ is an affine function, i.e., $\mathbf{b}(\mathbf{x})=\mathbf{B}\mathbf{x}+\mathbf{d}$ for some $\mathbf{B}\in{}^{n_{x}\times n_{x}}$ and $\mathbf{d}\in{}^{n_{x}}$, then $\Lambda_{\mathbf{M}}[\mathbf{b}](\mathbf{x})=\lambda_{\mathbf{M}}(\mathbf{B})$. Note that the local logarithmic Lipschitz constant and the logarithmic norm can be negative, unlike a standard Lipschitz constant. We also note that if $\mathbf{b}$ is differentiable, then $\Lambda_{\mathbf{M}}[\mathbf{b}](\mathbf{x})$ can be approximated by the logarithmic norm of the Jacobian $\mathbf{b}^{\prime}(\mathbf{x})$:

$$\frac{\left\langle\mathbf{z}-\mathbf{x},\mathbf{b}(\mathbf{z})-\mathbf{b}(\mathbf{x})\right\rangle_{\mathbf{M}}}{\left\|\mathbf{z}-\mathbf{x}\right\|_{\mathbf{M}}^{2}} =\frac{\left\langle\mathbf{z}-\mathbf{x},\mathbf{b}^{\prime}(\mathbf{x})(\mathbf{z}-\mathbf{x})\right\rangle_{\mathbf{M}}}{\left\|\mathbf{z}-\mathbf{x}\right\|_{\mathbf{M}}^{2}}+\mathcal{O}(\left\|\mathbf{z}-\mathbf{x}\right\|_{\mathbf{M}})\approx\frac{\left\langle\mathbf{z}-\mathbf{x},\mathbf{b}^{\prime}(\mathbf{x})(\mathbf{z}-\mathbf{x})\right\rangle_{\mathbf{M}}}{\left\|\mathbf{z}-\mathbf{x}\right\|_{\mathbf{M}}^{2}}\leq\lambda_{\mathbf{M}}(\mathbf{b}^{\prime}(\mathbf{x})).$$

We also need the following Grönwall-type inequality.

Lemma 5.1 (Grönwall inequality).

Let $T>0$ and $\alpha,\beta:[0,T]\to\real$ be integrable functions. If $u:[0,T]\to\real$ is differentiable and satisfies $u^{\prime}(t)\leq\beta(t)u(t)+\alpha(t)$ for all $t\in(0,T)$, then

$$u(t)\leq\int_{0}^{t}\alpha(s)e^{\int_{s}^{t}\beta(\tau)d\tau}ds+e^{\int_{0}^{t}\beta(\tau)d\tau}u(0)$$

for any $0\leq t\leq T$.

See, e.g., [63, Lemma 2.6] for a proof.

5.2 Error bounds

We now present an a posteriori error analysis for Kernel ROMs, which follows the approach detailed in [63]. The strategy is to view the Kernel ROM function $\hat{\mathbf{f}}$ as a regularized kernel interpolant of the intrusive projection-based ROM function $\tilde{\mathbf{f}}$, plus a discrepancy term ${\boldsymbol{\delta}}$ that accounts for the approximation error between $\tilde{\mathbf{f}}(\hat{\mathbf{q}}_{k}^{(\ell)})$ and the time derivative estimates $\dot{\hat{\mathbf{q}}}_{k}^{(\ell)}$ used to train the interpolant.

First, define the Kernel ROM reconstruction error

$$\mathbf{e}(t)\coloneqq\mathbf{q}(t)-\mathbf{g}(\hat{\mathbf{q}}(t)),$$ (5.3)

where $\mathbf{q}(t)$ is the solution to the FOM eq. 3.1, $\hat{\mathbf{q}}(t)$ is the solution to the Kernel ROM eq. 4.1, and $\mathbf{g}$ is the decompression map eq. 3.5. The reconstruction error evolves according to the system

$$\frac{\textrm{d}}{\textrm{d}t}\mathbf{e}(t)=\mathbf{f}(\mathbf{q}(t))-\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\hat{\mathbf{f}}(\hat{\mathbf{q}}(t)),\qquad\mathbf{e}(0)=\mathbf{q}_{0}-\mathbf{g}(\mathbf{V}^{\mathsf{T}}(\mathbf{q}_{0}-\bar{\mathbf{q}})),$$ (5.4)

where $\mathbf{g}^{\prime}(\hat{\mathbf{q}})=\mathbf{V}+\mathbf{W}[\mathbf{I}\otimes\hat{\mathbf{q}}+\hat{\mathbf{q}}\otimes\mathbf{I}]\in{}^{n_{q}\times r}$ is the Jacobian of $\mathbf{g}$. Although we use a QM to define the reconstruction mapping $\mathbf{g}$, the following error analysis holds for any reconstruction mapping of the same structure, namely with $\mathbf{g}$ taken to be the sum of an affine part and a nonlinear part.

Theorem 5.1 (A posteriori error).

If $\hat{\mathbf{f}}$ is an unregularized kernel interpolant of $\tilde{\mathbf{f}}+{\boldsymbol{\delta}}\in{\cal H}_{K}^{r}$ where $\left\|{\boldsymbol{\delta}}(\hat{\mathbf{q}}(s))\right\|_{\mathbf{M}}<\delta(s)$, then

$$\left\|\mathbf{e}(t)\right\|_{\mathbf{M}} \leq\int_{0}^{t}(\alpha_{P}(s)+\alpha_{K}(s))e^{\int_{s}^{t}\beta(\tau)d\tau}ds+e^{\int_{0}^{t}\beta(\tau)d\tau}\left\|\mathbf{e}_{N}(0)\right\|_{\mathbf{M}}, \forall \;t\in[0,T],$$ (5.5)

where

$$\alpha_{P}(s) =\left\|(\mathbf{I}-\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\mathbf{V}^{\mathsf{T}})\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(s)))\right\|_{\mathbf{M}},$$ (5.6a)

$$\alpha_{K}(s) =\left\|\mathbf{g}^{\prime}(\hat{\mathbf{q}}(s))\right\|_{\mathbf{M}}\left(P_{K,\tilde{\mathbf{Q}}}(\hat{\mathbf{q}}(s))\left\|\mathbf{L}\right\|_{2}\|\tilde{\mathbf{f}}+{\boldsymbol{\delta}}\|_{{\cal H}_{K}^{r}}+\delta(s)\right),\quad\text{and}$$ (5.6b)

$$\beta(s) =\Lambda_{\mathbf{M}}[\mathbf{f}](\mathbf{g}(\hat{\mathbf{q}}(s))).$$ (5.6c)

Proof.

Notice that the evolution equations in eq. 5.4 can be rewritten as

$$\frac{\textrm{d}}{\textrm{d}t}\mathbf{e}(t) =\mathbf{f}(\mathbf{q}(t))-\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\hat{\mathbf{f}}(\hat{\mathbf{q}}(t))-\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(t)))+\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(t)))$$

$$\quad-\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\mathbf{V}^{\mathsf{T}}\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(t)))+\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\mathbf{V}^{\mathsf{T}}\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(t)))-\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t)){\boldsymbol{\delta}}(\hat{\mathbf{q}}(t))+\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t)){\boldsymbol{\delta}}(\hat{\mathbf{q}}(t))$$

$$=\mathbf{f}(\mathbf{q}(t))-\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(t)))+(\mathbf{I}-\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\mathbf{V}^{\mathsf{T}})\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(t)))$$

$$\quad+\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\left(\mathbf{V}^{\mathsf{T}}\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}))+{\boldsymbol{\delta}}(\hat{\mathbf{q}}(t))-\hat{\mathbf{f}}(\hat{\mathbf{q}}(t))\right)-\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t)){\boldsymbol{\delta}}(\hat{\mathbf{q}}(t)).$$

Taking the $\mathbf{M}$-weighted inner product with $\mathbf{e}(t)$ and using the definition of the logarithmic Lipschitz constant and Corollary 2.2 yields

$$\left\langle\mathbf{e}(t),\frac{\textrm{d}}{\textrm{d}t}\mathbf{e}(t)\right\rangle_{\mathbf{M}}$$

$$=\left\langle\mathbf{e}(t),\mathbf{f}(\mathbf{q}(t))-\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(t)))\right\rangle_{\mathbf{M}}+\left\langle\mathbf{e}(t),(\mathbf{I}-\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\mathbf{V}^{\mathsf{T}})\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(t)))\right\rangle_{\mathbf{M}}$$

$$\quad+\left\langle\mathbf{e}(t),\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\left(\mathbf{V}^{\mathsf{T}}\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}))+{\boldsymbol{\delta}}(\hat{\mathbf{q}}(t))-\hat{\mathbf{f}}(\hat{\mathbf{q}}(t))\right)\right\rangle_{\mathbf{M}}-\left\langle\mathbf{e}(t),\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t)){\boldsymbol{\delta}}(\hat{\mathbf{q}}(t))\right\rangle_{\mathbf{M}}$$

$$\leq\Lambda_{\mathbf{M}}[\mathbf{f}](\mathbf{g}(\hat{\mathbf{q}}(t)))\left\|\mathbf{e}(t)\right\|_{\mathbf{M}}^{2}+\left\|(\mathbf{I}-\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\mathbf{V}^{\mathsf{T}})\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(t)))\right\|_{\mathbf{M}}\left\|\mathbf{e}(t)\right\|_{\mathbf{M}}$$

$$\quad+\left\|\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\right\|_{\mathbf{M}}\|\underbrace{\mathbf{V}^{\mathsf{T}}\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(t)))}_{\tilde{\mathbf{f}}(\hat{\mathbf{q}}(t))}+{\boldsymbol{\delta}}(\hat{\mathbf{q}}(t))-\hat{\mathbf{f}}(\hat{\mathbf{q}}(t))\|_{\mathbf{M}}\left\|\mathbf{e}(t)\right\|_{\mathbf{M}}+\left\|\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\right\|_{\mathbf{M}}\underbrace{\left\|{\boldsymbol{\delta}}(\hat{\mathbf{q}}(t))\right\|_{\mathbf{M}}}_{\leq\delta(t)}\left\|\mathbf{e}(t)\right\|_{\mathbf{M}}$$

$$\leq\Lambda_{\mathbf{M}}[\mathbf{f}](\mathbf{g}(\hat{\mathbf{q}}(t)))\left\|\mathbf{e}(t)\right\|_{\mathbf{M}}^{2}+\left\|(\mathbf{I}-\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\mathbf{V}^{\mathsf{T}})\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(t)))\right\|_{\mathbf{M}}\left\|\mathbf{e}(t)\right\|_{\mathbf{M}}$$

$$\quad+\left\|\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\right\|_{\mathbf{M}}P_{K,\tilde{\mathbf{Q}}}(\hat{\mathbf{q}}(t))\left\|\mathbf{L}\right\|_{2}\|\tilde{\mathbf{f}}+{\boldsymbol{\delta}}\|_{{\cal H}_{K}^{r}}\left\|\mathbf{e}(t)\right\|_{\mathbf{M}}+\delta(t)\left\|\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\right\|_{\mathbf{M}}\left\|\mathbf{e}(t)\right\|_{\mathbf{M}}.$$

Therefore,

$$\frac{\textrm{d}}{\textrm{d}t}\left\|\mathbf{e}(t)\right\|_{\mathbf{M}}=\frac{\left\langle\mathbf{e}(t),\frac{\textrm{d}}{\textrm{d}t}\mathbf{e}(t)\right\rangle_{\mathbf{M}}}{\left\|\mathbf{e}(t)\right\|_{\mathbf{M}}} \leq\Lambda_{\mathbf{M}}[\mathbf{f}](\mathbf{g}(\hat{\mathbf{q}}(t)))\left\|\mathbf{e}(t)\right\|_{\mathbf{M}}+\left\|(\mathbf{I}-\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\mathbf{V}^{\mathsf{T}})\mathbf{f}(\mathbf{g}(\hat{\mathbf{q}}(t)))\right\|_{\mathbf{M}}$$

$$\quad+\left\|\mathbf{g}^{\prime}(\hat{\mathbf{q}}(t))\right\|_{\mathbf{M}}\left(P_{K,\hat{\mathbf{Q}}}(\hat{\mathbf{q}}(t))\left\|\mathbf{L}\right\|_{2}\|\tilde{\mathbf{f}}+{\boldsymbol{\delta}}\|_{{\cal H}_{K}^{r}}+\delta(t)\right).$$

Applying Lemma 5.1 yields the result. ∎

A caveat to the result in Theorem 5.1 is that it relies on Corollary 2.2, which requires zero regularization. However, as we demonstrate empirically in Section 6, the error bound eq. 5.5 still holds when the regularization hyperparameter $\gamma$ is small. Secondly, computing the local logarithmic Lipschitz constant $\lambda_{\mathbf{M}}[\mathbf{f}]$ is difficult to do in general. In practice, we instead approximate it using the logarithmic norm of $\mathbf{f}^{\prime}(\mathbf{g}(\hat{\mathbf{q}}(t)))$. Lastly, we note that the estimate eq. 5.5 requires evaluating the FOM right-hand side $\mathbf{f}$, and therefore is a code-intrusive error bound. We leave the non-intrusive estimation of the bound eq. 5.5 to future work.

We also obtain the following a posteriori error result for intrusive projection-based ROMs by examining the special case $\hat{\mathbf{f}}=\tilde{\mathbf{f}}$.

Corollary 5.1.

The following error estimate holds for all $t\in[0,T]$:

$$\left\|\mathbf{q}(t)-\mathbf{g}(\tilde{\mathbf{q}}(t))\right\|_{\mathbf{M}} \leq\int_{0}^{t}\left\|(\mathbf{I}-\mathbf{g}^{\prime}(\tilde{\mathbf{q}}(s))\mathbf{V}^{\mathsf{T}})\mathbf{f}(\mathbf{g}(\tilde{\mathbf{q}}(s)))\right\|_{\mathbf{M}}e^{\int_{s}^{t}\Lambda_{\mathbf{M}}[\mathbf{f}](\mathbf{g}(\hat{\mathbf{q}}_{I}(\tau)))d\tau}ds$$ (5.7)

$$\quad+e^{\int_{0}^{t}\Lambda_{\mathbf{M}}[\mathbf{f}](\mathbf{g}(\tilde{\mathbf{q}}(\tau)))d\tau}\left\|\mathbf{q}(0)-\mathbf{g}(\hat{\mathbf{q}}(0))\right\|_{\mathbf{M}}.$$

We conclude with an error result comparing the intrusive projection-based ROM solution $\tilde{\mathbf{q}}(t)$ and the Kernel ROM solution $\hat{\mathbf{q}}(t)$. Let $\hat{\mathbf{e}}(t)=\tilde{\mathbf{q}}(t)-\hat{\mathbf{q}}(t)$, which satisfies the ODE

$$\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{e}}(t)=\tilde{\mathbf{f}}(\tilde{\mathbf{q}}(t))-\hat{\mathbf{f}}(\hat{\mathbf{q}}(t)),\qquad\hat{\mathbf{e}}(0)=\mathbf{0}.$$ (5.8)

We then have the following.

Proposition 5.1.

Let $\hat{\mathbf{M}}\in{}^{r\times r}$ be a symmetric positive definite weighting matrix with Cholesky factorization $\hat{\mathbf{M}}=\hat{\mathbf{L}}\hat{\mathbf{L}}^{\mathsf{T}}$. If $\hat{\mathbf{f}}$ is an unregularized kernel interpolant of $\tilde{\mathbf{f}}+{\boldsymbol{\delta}}\in{\cal H}_{K}^{r}$ where $\left\|{\boldsymbol{\delta}}(\hat{\mathbf{q}}(s))\right\|_{\hat{\mathbf{M}}}<\delta(s)$, then

$$\left\|\hat{\mathbf{e}}(t)\right\|_{\hat{\mathbf{M}}} \leq\int_{0}^{t}\left(P_{K,\tilde{\mathbf{Q}}}(\hat{\mathbf{q}}(s))\|\hat{\mathbf{L}}\|_{2}\|\tilde{\mathbf{f}}+{\boldsymbol{\delta}}\|_{{\cal H}_{K}^{r}}+\delta(s)\right)e^{\int_{s}^{\mathsf{T}}\Lambda_{\hat{\mathbf{M}}}[\hat{\mathbf{f}}](\hat{\mathbf{q}}(\tau))d\tau}ds, \forall\;t \in(0,T).$$ (5.9)

Proof.

The dynamics in eq. 5.8 can be rewritten as

$$\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{e}}(t) =\tilde{\mathbf{f}}(\tilde{\mathbf{q}}(t))-\tilde{\mathbf{f}}(\hat{\mathbf{q}}(t))+\tilde{\mathbf{f}}(\hat{\mathbf{q}}(t))+{\boldsymbol{\delta}}(\hat{\mathbf{q}}(t))-\hat{\mathbf{f}}(\hat{\mathbf{q}}(t))-{\boldsymbol{\delta}}(\hat{\mathbf{q}}(t)).$$

Taking the $\hat{\mathbf{M}}$-weighted inner product with $\hat{\mathbf{e}}(t)$ yields

$$\left\langle\hat{\mathbf{e}}(t),\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{e}}(t)\right\rangle_{\hat{\mathbf{M}}}$$

$$=\left\langle\hat{\mathbf{e}}(t),\tilde{\mathbf{f}}(\tilde{\mathbf{q}}(t))-\tilde{\mathbf{f}}(\hat{\mathbf{q}}(t))\right\rangle_{\hat{\mathbf{M}}}+\left\langle\hat{\mathbf{e}}(t),\tilde{\mathbf{f}}(\hat{\mathbf{q}}(t))+{\boldsymbol{\delta}}(\hat{\mathbf{q}}(t))-\hat{\mathbf{f}}(\hat{\mathbf{q}}(t))\right\rangle_{\hat{\mathbf{M}}}-\left\langle\hat{\mathbf{e}}(t),{\boldsymbol{\delta}}(\hat{\mathbf{q}}(t))\right\rangle_{\hat{\mathbf{M}}}$$

$$\leq\Lambda_{\hat{\mathbf{M}}}[\tilde{\mathbf{f}}](\hat{\mathbf{q}}(t))\left\|\hat{\mathbf{e}}(t)\right\|_{\hat{\mathbf{M}}}^{2}+\left\|\tilde{\mathbf{f}}(\hat{\mathbf{q}}(t))+{\boldsymbol{\delta}}(\hat{\mathbf{q}}(t))-\hat{\mathbf{f}}(\hat{\mathbf{q}}(t))\right\|_{\hat{\mathbf{M}}}\left\|\hat{\mathbf{e}}(t)\right\|_{\hat{\mathbf{M}}}+\delta(t)\left\|\hat{\mathbf{e}}(t)\right\|_{\hat{\mathbf{M}}}$$

$$\leq\Lambda_{\hat{\mathbf{M}}}[\tilde{\mathbf{f}}](\hat{\mathbf{q}}(t))\left\|\hat{\mathbf{e}}(t)\right\|_{\hat{\mathbf{M}}}^{2}+P_{K,\hat{\mathbf{Q}}}(\hat{\mathbf{q}}(t))\|\hat{\mathbf{L}}\|_{2}\|\tilde{\mathbf{f}}+{\boldsymbol{\delta}}\|_{{\cal H}_{K}^{r}}\left\|\hat{\mathbf{e}}(t)\right\|_{\hat{\mathbf{M}}}+\delta(t)\left\|\hat{\mathbf{e}}(t)\right\|_{\hat{\mathbf{M}}}.$$

Therefore

$$\frac{\textrm{d}}{\textrm{d}t}\left\|\hat{\mathbf{e}}(t)\right\|_{\hat{\mathbf{M}}} =\frac{\left\langle\hat{\mathbf{e}}(t),\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{e}}(t)\right\rangle_{\hat{\mathbf{M}}}}{\left\|\hat{\mathbf{e}}(t)\right\|_{\hat{\mathbf{M}}}}\leq\Lambda_{\hat{\mathbf{M}}}[\tilde{\mathbf{f}}](\hat{\mathbf{q}}(t))\left\|\hat{\mathbf{e}}(t)\right\|_{\hat{\mathbf{M}}}+P_{K,\hat{\mathbf{Q}}}(\hat{\mathbf{q}}_{N}(t))\|\hat{\mathbf{L}}\|_{2}\|\tilde{\mathbf{f}}+{\boldsymbol{\delta}}\|_{{\cal H}_{K}^{r}}+\delta(t).$$

Applying Lemma 5.1 yields the result. ∎

6 Numerical results

In this section, we test Kernel ROMs on several numerical examples using both POD and QM for dimension reduction. In each experiment, we construct Kernel ROMs with three kernel designs: 1) a feature map kernel encoding the full structure of the projection-based ROM, abbreviated “FM”; 2) an RBF kernel, marked “RBF”; and 3) a feature map-RBF hybrid kernel, labeled “Hybrid”. We also compare to the performance of intrusive projection-based ROMs in the first two examples and to OpInf in all three examples.

6.1 1D Advection-diffusion equation

We first consider a linear PDE, the advection-diffusion equation in one spatial dimension with periodic boundary conditions:

$$\frac{\partial}{\partial t}q(x,t)-\kappa\frac{\partial^{2}}{\partial x^{2}}q(x,t)+\beta\frac{\partial}{\partial x}q(x,t)=0, \qquad x\in(0,1),\quad t\in(0,T),$$ (6.1a)

$$q(0,t)=q(1,t),\quad\frac{\partial}{\partial x}q(0,t)=\frac{\partial}{\partial x}q(1,t), \qquad t\in(0,T),$$ (6.1b)

$$q(x,0)=q_{0}(x;{\boldsymbol{\mu}})\coloneqq e^{-({x-\mu_{1}})^{2}/\mu_{2}^{2}}, \qquad x\in(0,1).$$ (6.1c)

Here, $\kappa>0$ is the diffusion parameter, $\beta\geq 0$ is the advection parameter, $T>0$ is the final time, and ${\boldsymbol{\mu}}=(\mu_{1},\mu_{2})$ parameterizes the initial condition. For this experiment, we set $\kappa=10^{-2}$, $\beta=1$, and $T=1$. The initial condition is a Gaussian pulse with center $\mu_{1}\in[0.25,0.35]$ and width $\mu_{2}\in[0.05,0.15]$. The dynamics of eq. 6.1 are linear, but advective phenomena can be difficult to capture with linear dimension reduction methods such as POD.

Figure 1: Solutions of the full-order advection diffusion model eq. 6.2 with initial conditions eq. 6.1c for various choices of ${\boldsymbol{\mu}}$.

Spatially discretizing eq. 6.1 with an upwind finite difference scheme over a grid of $n_{q}+1$ uniformly spaced points in the spatial domain $[0,1]$ results in a linear FOM of the form

$$\frac{\textrm{d}}{\textrm{d}t}\mathbf{q}(t) =\mathbf{A}\mathbf{q}(t),\qquad\mathbf{q}(0)=\mathbf{q}_{0}({\boldsymbol{\mu}}),$$ (6.2)

where $\mathbf{q}(t),\mathbf{q}_{0}({\boldsymbol{\mu}})\in{}^{n_{q}}$ and $\mathbf{A}\in{}^{n_{q}\times n_{q}}$. We use $n_{q}=256$ spatial degrees of freedom in this experiment. To collect training data, we sample $M=10$ initial conditions corresponding to $10$ Latin hypercube samples from the parameter domain ${\cal D}=[0.25,0.35]\times[0.05,0.15]$ and integrate the FOM eq. 6.2 using a fully implicit variable-order backwards difference formula (BDF) time stepper with quasi-constant step size, executed with scipy.interpolate.solve_ivp() in Python [60, 55]. The solution is recorded at $n_{t}=256$ equally spaced time instances after the initial condition, resulting in $M(n_{t}+1)=2570$ total training snapshots. We also solve the FOM at the testing parameter value $\bar{{\boldsymbol{\mu}}}=(0.3,0.1)$, which is not included in the training set. Figure 1 plots the FOM states for two training parameter values and the testing parameter value.

POD	${\boldsymbol{\phi}}(\hat{\mathbf{q}})=\begin{bmatrix}1\\ \hat{\mathbf{q}}\end{bmatrix},\qquad\mathbf{G}=\frac{1}{1+r}\mathbf{I}_{1+r}$
QM	${\boldsymbol{\phi}}(\hat{\mathbf{q}})=\begin{bmatrix}1\\ \hat{\mathbf{q}}\\ \hat{\mathbf{q}}\otimes\hat{\mathbf{q}}\end{bmatrix},\qquad\mathbf{G}=\begin{bmatrix}\mathbf{I}_{1+r}&\mathbf{0}\\ \mathbf{0}&\left\|\mathbf{W}\right\|_{F}\mathbf{I}_{r^{2}}\end{bmatrix}$

Table 2: Feature maps and weighting matrices in POD and QM Kernel ROMs for the 1D advection-diffusion example.

The training snapshots are used to compute POD and QM state approximations with the reference vector $\bar{\mathbf{q}}$ set to the average training snapshot. Since the FOM eq. 6.2 is linear and $\bar{\mathbf{q}}\neq\bf 0$, the intrusive projection-based POD ROM of dimension $r$ has affine structure,

$$\frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)=\tilde{\mathbf{c}}+\tilde{\mathbf{A}}\tilde{\mathbf{q}}(t),$$ (6.3)

where $\tilde{\mathbf{c}}\in{}^{r}$ and $\tilde{\mathbf{A}}\in{}^{r\times r}$, whereas the intrusive QM ROM has the form

$$\frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)=\tilde{\mathbf{c}}+\tilde{\mathbf{A}}\tilde{\mathbf{q}}(t)+\tilde{\mathbf{H}}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)],$$ (6.4)

with $\tilde{\mathbf{c}}\in{}^{r}$, $\tilde{\mathbf{A}}\in{}^{r\times r}$, and $\tilde{\mathbf{H}}\in{}^{r\times r^{2}}$. For both POD and QM, we construct feature map Kernel ROMs and OpInf ROMs with the corresponding intrusive ROM structure. The underlying feature maps ${\boldsymbol{\phi}}$ and weighting matrices $\mathbf{G}$ are listed in Table 2. Note that the second diagonal block in the weight $\mathbf{G}$ for QM is scaled by $\left\|\mathbf{W}\right\|_{F}$ to account for the fact that $\tilde{\mathbf{H}}$ in the intrusive QM ROM eq. 6.4 also depends on $\mathbf{W}$. We also construct an RBF Kernel ROM using a Gaussian kernel-generating RBF $\psi$ (see Table 1) with fixed shape parameter $\epsilon=10^{-1}$. This ROM has the same evolution equations in the POD and QM cases, since the compression map $\mathbf{h}$ is the same in both instances, but we report results for both POD and QM decompression maps $\mathbf{g}$. Finally, we construct hybrid Kernel ROMs using the POD feature map from Table 2 with weighting coefficient $c_{\phi}=1$ and a Gaussian RBF kernel with $\epsilon=0.1$ and weighting coefficient $c_{\psi}=10^{-3}$, yielding ROMs with the following structure:

$$\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t)=\hat{\mathbf{c}}+\hat{\mathbf{A}}\hat{\mathbf{q}}(t)+10^{-3}{\boldsymbol{\Omega}}^{\mathsf{T}}{\boldsymbol{\psi}}_{\!\epsilon}(\hat{\mathbf{q}}(t)).$$ (6.5)

For QM, the RBF term takes the place of the quadratic nonlinearity $\tilde{\mathbf{H}}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)]$, but for POD, the RBF term is purely supplementary. Kernel input normalization as in Remark 2.1 is not needed in this problem. Performance is measured with a relative $L^{\infty}$-$L^{2}$ error between the FOM and reconstructed ROM states,

$$\mathbf{e}(\mathbf{q},\hat{\mathbf{q}})=\frac{\max_{k}\;\left\|\mathbf{q}(t_{k})-\mathbf{g}(\hat{\mathbf{q}}(t_{k}))\right\|_{2}}{\max_{k}\;\left\|\mathbf{q}(t_{k})\right\|_{2}},$$ (6.6)

where the ROMs are integrated with the same BDF time stepper as the FOM and the maxima are taken over time indices $k\in\{0,1,\ldots,n_{t}\}$. The ROM error is bounded from below by the projection error $\mathbf{e}(\mathbf{q},\mathbf{h}(\mathbf{q}))$.

Figure 2: Relative ROM error at the test parameter $\bar{{\boldsymbol{\mu}}}=(0.3,0.1)$ as a function of number of basis vectors in linear POD (left) and quadratic manifold (right) reduced state approximations for the advection-diffusion problem eq. 6.1.

Results are reported in Figure 2, which compares ROM and projection errors at the testing parameter value $\bar{{\boldsymbol{\mu}}}$ for both POD and QM as a function of the reduced dimension $r$. For each Kernel ROM, the regularization hyperparameter $\gamma$ for the learning problem eq. 4.5 is selected to minimize the ROM error over the training data, i.e.,

$$\gamma=\underset{\gamma}{\arg\min}\sum_{\ell=1}^{M}\sum_{k=0}^{n_{t}}\big{\|}\hat{\mathbf{q}}_{k}^{(\ell)}-\hat{\mathbf{q}}(t_{k};{\boldsymbol{\mu}}_{\ell},\gamma)\big{\|}_{2},$$ (6.7)

where $\hat{\mathbf{q}}_{k}^{(\ell)}$ are the training snapshots eq. 4.2 and $\hat{\mathbf{q}}(t;{\boldsymbol{\mu}}_{\ell},\gamma)$ denotes the solution to the Kernel ROM with regularization $\gamma$ evaluated for training parameter ${\boldsymbol{\mu}}_{\ell}$. In this experiment, we do this via a grid search over $\gamma\in\{10^{-14},10^{-13},\ldots,10^{2}\}$ for each Kernel ROM. This procedure is adapted from best practices for OpInf [34, 42]; a similar selection is carried out for OpInf ROMs with the regularization matrix ${\boldsymbol{\Gamma}}$ parameterized so that

$$\big{\|}{\boldsymbol{\Gamma}}\hat{\mathbf{O}}^{\mathsf{T}}\big{\|}_{F}^{2}=\gamma_{1}^{2}(\|\hat{\mathbf{c}}\|_{2}^{2}+\|\hat{\mathbf{A}}\|_{F}^{2})+\gamma_{2}^{2}\|\hat{\mathbf{H}}\|_{F}^{2},$$ (6.8)

where $\gamma_{1},\gamma_{2}\geq 0$. This is the state-of-the-art procedure for OpInf and results in accurate ROMs. Indeed, Figure 2 shows that each of the POD-based ROMs yield errors that are nearly identical to the POD projection error for $r\leq 15$. The POD RBF Kernel ROM error plateaus for $r>15$, possibly due to the RBF shape parameter being fixed independent of $r$. The POD hybrid Kernel ROM error begins to plateau for $r>17$, again possibly due to the fixed RBF shape parameter and fixed weighting coefficients $c_{\phi}$ and $c_{\psi}$. The OpInf ROMs and feature map Kernel ROMs match the projection error for $r\leq 20$, but deviate slightly from the projection and intrusive ROM errors for some values of $r>20$.

Figure 3: Effect of the quadratic manifold regularization parameter $\rho$ on the projection error and intrusive ROM error for two different reduced state dimensions $r$, using data from the advection-diffusion problem eq. 6.1 at the test parameter $\bar{{\boldsymbol{\mu}}}=(0.3,0.1)$.

The QM regularization parameter $\rho\geq 0$ in eq. 3.12 plays an important role in the stability and accuracy of QM ROMs, see Appendix B for a stability analysis of the intrusive QM ROM for a linear FOM. Figure 3 plots the value of $\rho$ versus the projection error and the intrusive ROM error for two choices of the reduced dimension $r$. As is evident from eq. 3.12, $\left\|\mathbf{W}\right\|_{F}\to 0$ as $\rho$ increases, which is why the QM projection and QM ROM errors approach their POD counterparts for large enough $\rho$. Note that the optimal $\rho$ varies with the reduced state dimension $r$. Furthermore, at least for $r=12$, the best $\rho$ for the reconstruction error is not necessarily the best $\rho$ for the intrusive QM ROM error. To account for this, the QM results in Figure 2 report only the best results for each ROM after testing each of the QM regularization values $\rho\in\{10^{-3},10^{-2},\ldots,10^{8}\}$. In other words, Figure 2 shows a best-case scenario comparison. The QM OpInf ROMs and QM feature map Kernel ROMs again show highly similar performance, while the QM RBF and QM hybrid Kernel ROM errors plateau for $r>17$. Note that the POD and QM projection errors are close for $r>15$, indicating that in this particular problem QM results in diminishing returns over POD for large enough $r$.

Figure 4: Approximate error bounds for POD and QM Kernel ROMs for the advection-diffusion problem eq. 6.1.

Next, we compute the error bound from Theorem 5.1 for the feature map Kernel ROMs for $r\in\{6,12\}$. Although the computed Kernel ROMs use a nonzero regularization $\gamma\neq 0$, the computed error bounds still hold. We estimate the norm $\|\tilde{\mathbf{f}}+{\boldsymbol{\delta}}\|_{{\cal H}_{K}^{r}}$ with the norm of the interpolant $\|\hat{\mathbf{f}}\|_{{\cal H}_{K}^{r}}$, which can be computed quickly and explicitly using equation eq. 2.1. The local logarithmic Lipschitz constant $\Lambda_{\mathbf{M}}[\mathbf{f}](\mathbf{g}(\hat{\mathbf{q}}(s)))$ is estimated using the logarithmic norm $\lambda_{\mathbf{M}}(\mathbf{f}(\hat{\mathbf{q}}(s)))$, and the weighting matrix $\mathbf{M}$ is taken to be $\frac{1}{r}\mathbf{I}_{r}$. We also examine feature map Kernel ROMs where the chosen feature map does not match the true projection-based ROM form, i.e. POD with a quadratic feature map and QM with a linear feature map. The results are displayed in Figure 4, which shows that the computed error estimates indeed bound the true error without dramatically overestimating it. In the POD cases with linear ROMs, the $\alpha_{P}$ term, which is related to the POD projection errors, is what dominates the error bound computation, while the $\alpha_{K}$ term, which corresponds to the pointwise kernel error bound from Corollary 2.2, is negligible. For the QM Quadratic ROM with $r=6$, $\alpha_{P}$ again dominates the error bound and the $\alpha_{K}$ is negligible. However, for the QM Quadratic ROM with $r=12$, the $\alpha_{K}$ term is much larger. This may indicate that the chosen quadratic feature map may yield a non-optimal model form for the Kernel ROM. Indeed, since POD with $r=12$ already yields small ROM errors, one may expect that a QM is unnecessary for $r=12$, and thus the quadratic term in the Kernel ROM may be extraneous. To test this, we remove the quadratic term, which comes from the quadratic component of $\mathbf{g}$, and compute the error bound for a linear QM Kernel ROM with $r=12$. We observe that the $\alpha_{K}$ term is once again negligible in this case. On the other hand, adding a quadratic term to the POD ROM with $r=12$ also substantially increases $\alpha_{K}$. Therefore, we can infer that a larger $\alpha_{K}$ contribution may indicate that a non-optimal model form (i.e., feature map) was used for the Kernel ROM.

6.2 1D Burgers’ equation

Figure 5: Solutions of the full-order Burgers’ model eq. 6.10 with initial conditions eq. 6.9c for various choices of ${\boldsymbol{\mu}}$.

We now consider the 1D viscous Burgers’ equation with homogeneous Dirichlet boundary conditions, which is nonlinear with respect to the state:

$$\frac{\partial}{\partial t}q(x,t)-\nu\frac{\partial^{2}}{\partial x^{2}}q(x,t)+q(x,t)\frac{\partial}{\partial x}q(x,t)=0, \qquad x\in(0,1),\quad t\in(0,T),$$ (6.9a)

$$q(0,t)=0,\quad q(1,t)=0, \qquad t\in(0,T),$$ (6.9b)

$$q(x,0)=q_{0}(x;{\boldsymbol{\mu}})\coloneqq e^{-(x-\mu_{1})^{2}/\mu_{2}^{2}}, \qquad x\in(0,1).$$ (6.9c)

Here, $\nu>0$ is the viscosity, which we set to $\nu=10^{-2}$ for our experiments. Solutions to this system are characterized by sharp gradients along an advection front. Just as in the previous problem, we consider parameterized Gaussian initial conditions with center $\mu_{1}\in[0.25,0.35]$ and $\mu_{2}\in[0.05,0.15]$, set the final time to $T=1$, use $n_{q}=n_{t}=256$ spatial degrees of freedom and temporal observations, and draw $M=10$ latin hypercube samples of the parameters ${\boldsymbol{\mu}}=[~{}\mu_{1}~{}~{}\mu_{2}~{}]^{\mathsf{T}}$ to use for generating training data. The spatial discretization uses uniform centered finite differences, yielding a quadratic FOM of the form

$$\frac{\textrm{d}}{\textrm{d}t}\mathbf{q}(t)=\mathbf{A}\mathbf{q}(t)+\mathbf{H}[\mathbf{q}(t)\otimes\mathbf{q}(t)],\qquad\mathbf{q}(0)=\mathbf{q}_{0}({\boldsymbol{\mu}}),$$ (6.10)

where $\mathbf{q}(t),\mathbf{q}_{0}({\boldsymbol{\mu}})\in{}^{n_{q}}$, $\mathbf{A}\in{}^{n_{q}\times n_{q}}$, and $\mathbf{H}\in{}^{n_{q}\times n_{q}^{2}}$. We again use a BDF time integrator to solve the FOM (and constructed ROMs) at the parameter samples, resulting in $M=10$ trajectories of $n_{t}+1=257$ snapshots each. The FOM states for a few parameter values are displayed in Figure 5.

POD

${\boldsymbol{\phi}}(\hat{\mathbf{q}})=\begin{bmatrix}1\\ \hat{\mathbf{q}}\\ \hat{\mathbf{q}}\otimes\hat{\mathbf{q}}\end{bmatrix},\qquad\mathbf{G}=\frac{1}{1+r+r^{2}}\mathbf{I}_{1+r+r^{2}}$

QM

${\boldsymbol{\phi}}(\hat{\mathbf{q}})=\begin{bmatrix}1\\ \hat{\mathbf{q}}\\ \hat{\mathbf{q}}\otimes\hat{\mathbf{q}}\\ \hat{\mathbf{q}}\otimes\hat{\mathbf{q}}\otimes\hat{\mathbf{q}}\\ \hat{\mathbf{q}}\otimes\hat{\mathbf{q}}\otimes\hat{\mathbf{q}}\otimes\hat{\mathbf{q}}\end{bmatrix},\qquad\mathbf{G}=\begin{bmatrix}\mathbf{I}_{1+r+r^{2}}&\mathbf{0}&\mathbf{0}\\ \mathbf{0}&\left\|\mathbf{W}\right\|_{F}\mathbf{I}_{r^{3}}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}&\left\|\mathbf{W}\right\|_{F}^{2}\mathbf{I}_{r^{4}}\\ \end{bmatrix}$

Table 3: Feature maps and weighting matrices for POD and QM Kernel ROMs for the 1D Burgers’ example.

For both POD and QM, we use $\bar{\mathbf{q}}=\bf 0$, hence the intrusive POD ROM takes the quadratic form eq. 3.15, whereas the intrusive QM ROM has the quartic form eq. 3.16. For this problem, we apply the kernel input normalization discussed in Remark 2.1, which is helpful for balancing the contribution of higher-order terms. We therefore construct feature map Kernel ROMs to mirror the structure of the intrusive models, with the addition of a constant term that arises due to the input scaling (see Appendix A), by using the feature maps and weighting matrices listed in Table 3. Similar to before, we learn Gaussian RBF Kernel ROMs with fixed shape parameter $\epsilon=10^{-1}$ and hybrid Kernel ROMs using the POD feature map from Table 3 with weighting coefficient $c_{\phi}=1$ and a Gaussian RBF kernel with $\epsilon=0.1$ and weighting coefficient $c_{\psi}=10^{-3}$, which result in ROMs of the form

$$\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t)=\hat{\mathbf{c}}+\hat{\mathbf{A}}\hat{\mathbf{q}}(t)+\tilde{\mathbf{H}}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)]+10^{-3}{\boldsymbol{\Omega}}^{\mathsf{T}}{\boldsymbol{\psi}}_{\!\epsilon}(\tilde{\mathbf{q}}(t)).$$ (6.11)

We also learn OpInf ROMs with the intrusive ROM structure, with the regularization designed so

$$\big{\|}{\boldsymbol{\Gamma}}\hat{\mathbf{O}}^{\mathsf{T}}\big{\|}_{F}^{2}=\gamma_{1}^{2}\|\hat{\mathbf{A}}\|_{F}^{2}+\gamma_{2}^{2}\|\hat{\mathbf{H}}\|_{F}^{2}$$ (6.12)

for the POD OpInf ROM, and

$$\big{\|}{\boldsymbol{\Gamma}}\hat{\mathbf{O}}^{\mathsf{T}}\big{\|}_{F}^{2}=\gamma_{3}^{2}(\|\hat{\mathbf{H}}_{2}\|_{F}^{2}+\|\hat{\mathbf{H}}_{3}\|_{F}^{2})+\gamma_{4}^{2}\|\hat{\mathbf{H}}_{4}\|_{F}^{2}$$ (6.13)

for the QM OpInf ROM, performing a grid search for $\gamma_{1},\ldots,\gamma_{4}>0$. The relative $L^{\infty}$-$L^{2}$ error eq. 6.6 is used to evaluate ROM performance at the testing parameter value $\bar{{\boldsymbol{\mu}}}=(0.3,0.1)$.

Figure 6: Relative ROM error at the test parameter $\bar{{\boldsymbol{\mu}}}=(0.3,0.1)$ as a function of number of basis vectors in linear POD (left) and quadratic manifold (right) state approximations.

Figure 6 reports results for various reduced dimensions $r$. All POD ROM errors are nearly identical to the POD projection error. For the QM ROMs, the OpInf and Kernel ROMs have very similar performance for $r\leq 14$. The feature map Kernel ROM errors plateau for $r>14$, while the OpInf and RBF Kernel ROMs plateau for $16<r<18$ and increase slightly for $r=19,20$. Notably, the hybrid Kernel ROM continues to match the projection and intrusive ROM error as $r$ increases, indicating that the RBF term in eq. 6.11 acts as a more accurate closure term for the ROM dynamics at larger values of $r$ compared to the cubic and quartic nonlinearities of the OpInf and FM Kernel ROMs. Unlike the advection-diffusion case, the QM projection and intrusive ROM errors are notably lower than the corresponding POD errors, and thus QM dimension reduction may be beneficial for this problem.

Figure 7: Approximate error bound for POD and QM kernel ROMs for different values of $r$.

We next compute the error bound from Theorem 5.1 for the FM Kernel ROMs for $r=6,12$. The quantities $\|\tilde{\mathbf{f}}+{\boldsymbol{\delta}}\|_{{\cal H}_{K}^{r}}$, $\delta(s)$, $\Lambda_{\mathbf{M}}[\mathbf{f}](\mathbf{g}(\hat{\mathbf{q}}(s)))$ are estimated in the same way as in the advection-diffusion case. We use feature map Kernel ROMs corresponding to the quadratic and quartic feature maps in Table 3 and examine the cases when the chosen feature map does not match the true projection-based ROM form, i.e. POD with a quartic feature map and QM with a quadratic feature map. Figure 7 displays the results and shows that the computed error estimate again bounds the true error without dramatically overestimating it. In the POD cases with quadratic ROMs, the $\alpha_{P}$ term dominates the error bound contribution, while the $\alpha_{K}$ term is negligible in the $r=6$ case, but less negligible in the $r=12$ case. For the QM quartic ROMs, the $\alpha_{P}$ and $\alpha_{K}$ terms contribute similarly to the error bound evaluation. The $r=6$ case contrasts with the advection-diffusion QM ROM case in that the $\alpha_{K}$ term is non-negligible despite having a model form that should reproduce the projection-based ROM model form. We again compute the error bound for a QM ROM with the cubic and quartic terms removed, which come from the quadratic part of $\mathbf{g}$, resulting in a QM quadratic ROM. As in the advection-diffusion example, the $\alpha_{K}$ term decreases significantly, which may indicate that a quadratic model form may be the better choice for a QM Burgers ROM. To again test if an incorrect model form significantly increases $\alpha_{K}$, we compute a POD quartic ROM and observe that $\alpha_{K}$ is much larger than for the POD quadratic ROM, as expected. This further evidences that a larger $\alpha_{K}$ contribution may indicate that a non-optimal model form is being used for the Kernel ROM.

6.3 2D Euler–Riemann problem

Our last numerical example uses the 2D conservative Euler equations

$$\frac{\partial}{\partial t}\begin{bmatrix}\rho\\ \rho u\\ \rho v\\ \rho E\end{bmatrix}+\frac{\partial}{\partial x}\begin{bmatrix}\rho u\\ \rho u^{2}+p\\ \rho uv\\ (E+p)u\end{bmatrix}+\frac{\partial}{\partial x}\begin{bmatrix}\rho v\\ \rho uv\\ \rho v^{2}+p\\ (E+p)v\end{bmatrix}=0,$$ (6.14)

where $u$ is the $x$-velocity, $v$ is the $y$-velocity, $\rho$ is the fluid density, $p$ is the pressure, and $E$ is the energy. The system is closed by the state equation

$$p=(\gamma-1)\left(\rho E-\frac{1}{2}\rho(u^{2}+v^{2})\right),$$ (6.15)

where $\gamma=1.4$ is the specific heat ratio. The spatial domain is the unit square $\Omega=(0,1)\times(0,1)$ with homogeneous Neumann boundary conditions on each side, and the time domain is $(0,0.8)$.

The initial condition is given by a classical Riemann problem as follows. The spatial domain is divided into four quadrants with a vertical dividing line at $x=0.8$ and a horizontal dividing line at $y=0.8$. The initial pressure is set to $p_{BL}=0.029$ in the bottom left quadrant; in the top right quadrant, the initial velocities are fixed at $u_{TR}=v_{TR}=0$, and the initial density is $\rho_{TR}=1.5$. We parameterize the initial condition by setting the upper-right quadrant pressure to $p_{TR}\in\left\{0.5,0.75,1.0,1.25,1.5\right\}$ and compute remaining quantities following the relations in [53, Configuration 3]. For testing, we consider the initial upper-right quadrant pressure to $\bar{p}_{TR}=1.125$. In every case, the discontinuities of the initial condition propagate through the domain, a highly challenging scenario for projection-based model reduction.

(a) Training, $p_{TR}=0.5$

(b) Training, $p_{TR}=1.5$

(c) Testing, $p_{TR}=1.125$

Figure 8: Pressure snapshots for the 2D Euler equations corresponding to different initial values of $p_{TR}$.

We collect FOM snapshots using the open-source Python library pressio-demoapps to simulate eq. 6.14, which uses a cell-centered finite volume scheme. For this example, we use a $256\times 256$ uniform Cartesian mesh, resulting in a FOM with state dimension $n_{q}=256\times 256\times 4=262,144$, and a Weno5 scheme for inviscid flux reconstruction. The FOM time stepping is done using pressio-demoapps’ SSP3 scheme for times $t\in(0,0.8)$ with time step $\Delta t=0.001$, while the ROM is integrated with BDF time stepping. The first $2000$ normalized POD singular values are plotted in Figure 9; the slow decay indicates the high difficulty of the problem for POD-based methods.

Figure 9: First $2000$ normalized singular values $\sigma_{k}/\sigma_{1}$ for the 2D Euler example.

Before computing ROMs, the FOM state variables are first transformed via the map

$$\begin{bmatrix}\rho\\ \rho u\\ \rho v\\ \rho E\end{bmatrix}\mapsto\begin{bmatrix}u\\ v\\ p\\ \zeta\end{bmatrix},$$ (6.16)

where $\zeta=1/\rho$ is the specific volume. A discretized FOM using the specific volume formulation is purely quadratic,

$$\frac{\textrm{d}}{\textrm{d}t}\mathbf{q}(t)=\mathbf{H}[\mathbf{q}(t)\otimes\mathbf{q}(t)],\qquad\mathbf{q}(0)=\mathbf{q}_{0}(p_{TR}).$$ (6.17)

This FOM is not formed explicitly, but it motivates an appropriate structure for feature map Kernel ROMs using POD or QM. In both cases, we set $\bar{\mathbf{q}}$ to the average training snapshot and apply the kernel input normalization from Remark 2.1, leading to a POD ROM structure

$$\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t)=\hat{\mathbf{c}}+\hat{\mathbf{A}}\hat{\mathbf{q}}(t)+\hat{\mathbf{H}}[\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)],$$ (6.18)

whereas the QM ROMs have the quartic form

$$\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t)=\hat{\mathbf{c}}+\hat{\mathbf{A}}\hat{\mathbf{q}}(t)+\hat{\mathbf{H}}_{2}[\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)]+\hat{\mathbf{H}}_{3}[\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)]+\hat{\mathbf{H}}_{4}[\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)].$$ (6.19)

Since we use pressio-demoapps to collect FOM data, this example only considers the purely non-intrusive cases. That is, we do not compute intrusive ROMs for this problem and do not evaluate the a posteriori error bound as in the previous examples.

The POD and QM OpInf ROMs are constructed to have the same structure as eq. 6.18 and eq. 6.19, respectively. Notice that this is the same structure as for Burgers’ equation. Consequently, the feature map Kernel ROMs use the same feature maps ${\boldsymbol{\phi}}$ and weighting matrices $\mathbf{G}$ as in Table 3. As in both previous examples, the RBF Kernel ROMs use a Gaussian RBF kernel with fixed shape parameter $\epsilon=10^{-1}$. The hybrid Kernel ROMs use the sum of the kernel induced by the POD feature map from Table 3 with weighting coefficient $c_{\phi}=1$ and the same Gaussian RBF kernel with $\epsilon=0.1$ and weighting coefficient $c_{\psi}=10^{-3}$, resulting in a right-hand side of the form eq. 6.11. The error metric that we consider is the relative $L^{\infty}$-$L^{1}$ norm

$$\mathbf{e}(\mathbf{q},\hat{\mathbf{q}})=\frac{\max_{k}\;\left\|\mathbf{q}(t_{k})-\mathbf{g}(\hat{\mathbf{q}}(t_{k}))\right\|_{1}}{\max_{k}\;\left\|\mathbf{q}(t_{k})\right\|_{1}},$$ (6.20)

The $L^{1}$ norm is more appropriate than $L^{2}$ for this problem due to the discontinuities in the solution.

We plot the error eq. 6.20 versus the reduced dimension $r$ for the POD OpInf, feature map Kernel, RBF Kernel, and hybrid Kernel ROMs in Figure 10. For $r=5,10,15$, each of the ROMs obtain nearly identical performance. For $r>15$, the projection error and the Kernel ROM errors plateau, with the Kernel ROMs yielding a $<2\%$ difference in error compared to the projection error. The Hybrid and FM Kernel ROMs have nearly identical errors, while the RBF yields slightly different but very similar errors. The OpInf ROM increases slightly in error for $r>15$, yet still obtains errors within a few percent of the projection error. We note that the plateauing of the ROM and projection errors for the tested ROM sizes is expected since the singular value decay is slow, as shown in Figure 9.

Figure 10: Relative errors for several non-intrusive POD ROMs for the Euler–Riemann problem eq. 6.14.

We omit a similar comparison for the QM ROMs for this problem because the resulting ROMs are highly dependent on the QM regularization $\rho$, and require very large values of $\rho$ to obtain a stable ROM. To illustrate this, we compute QM Kernel FM ROMs for $r=10,20$ for QM regularizations $\rho\in\left\{10^{0},10^{1},\dots,10^{12}\right\}$ and plot the resulting errors, see Figure 11. For $r=10$, we observe that the QM Kernel ROM errors are very large for $\rho<10^{10}$, whereas the corresponding QM projection errors are relatively small. The QM ROM errors do not approach the QM projection errors until $\rho=10^{11}$, where a slightly better error compared to POD is achieved. For $r=20$, the QM ROMs for $\rho<10^{8}$ are unstable and do not finish the time integration, while for $\rho=10^{8},10^{9},10^{10}$, the QM ROM errors still exceed the POD errors. The QM ROMs for $\rho=10^{11},10^{12}$ obtain yield the best errors, but because the QM regularization $\rho$ is so large, the resulting ROM errors are no better than POD.

Figure 11: QM regularization $\rho$ versus relative $L^{\infty}$-$L^{1}$ error.

7 Conclusion

This paper develops a novel non-intrusive model reduction technique grounded in regularized kernel interpolation. While previous approaches approximate the ROM dynamics by solving a data-driven polynomial regression problem, our approach yields an optimal approximant to the ROM dynamics from an RKHS, which is determined by the choice of kernel. In particular, using kernels induced by feature maps allows one to imbue interpretable structure into the resulting ROM. Furthermore, using an RBF kernel or a hybrid approach using the sum of a feature map and an RBF kernel allows one to compute effective non-intrusive ROMs that incorporate no structure or partial structure. The hybrid approach also provides a natural way of incorporating closure terms into our ROM formulation, and this approach was demonstrated to be effective in each of the numerical examples. Since the approximant lives in an RKHS, we can leverage the pointwise error bound from Theorem 2.2, a standard result from RKHS theory, as well as standard intrusive ROM error estimates to derive an a posteriori error estimate for our Kernel ROMs in Theorem 5.1. This error estimate, as well as the added flexibility afforded by arbitrary choices of kernel, are key innovations of our approach.

Future work will focus on expanding the applicability and efficiency of Kernel ROMs. In particular, we will extend our approach to problems where the FOM right-hand side $\mathbf{f}$ is parametrized, which is the case in many engineering applications of interest. Second, we will implement a greedy sampling procedure to build a minimal training set for the kernel interpolants. This is particularly relevant when using an RBF interpolant, since the computation cost of evaluating the RBF interpolant is proportional to the amount of training data whenever the kernel is not entirely prescribed by feature maps. Third, we will develop a method for non-intrusively approximating the a posteriori error bound in Theorem 5.1. As mentioned in Section 6, evaluating the bound eq. 5.5 requires access to the FOM right-hand side $\mathbf{f}$, which we assume that we cannot access in the fully non-intrusive setting. Therefore, in future work, it will be necessary to develop an accurate estimator for the quantities in eq. 5.6.

Appendix A Quadratic systems with QM approximations

This appendix considers a linear-quadratic FOM,

$$\frac{\textrm{d}}{\textrm{d}t}\mathbf{q}(t)=\mathbf{f}(\mathbf{q}(t))\coloneqq\mathbf{A}\mathbf{q}(t)+\mathbf{H}[\mathbf{q}(t)\otimes\mathbf{q}(t)],$$ (3.14)

and derives the structure of the corresponding intrusive projection-based ROM with a QM approximation,

$$\mathbf{g}(\tilde{\mathbf{q}})=\bar{\mathbf{q}}+\mathbf{V}\tilde{\mathbf{q}}+\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}],$$ (3.5)

for nonzero $\bar{\mathbf{q}}\in{}^{n_{q}}$ and $\mathbf{W}\in{}^{n_{q}\times r^{2}}$. Specifically, we show that a nonzero reference vector $\bar{\mathbf{q}}$ causes a constant term to appear in the ROM dynamics.

Using the product rule $(\mathbf{X}\otimes\mathbf{Y})(\mathbf{Z}\otimes\mathbf{U})=(\mathbf{X}\mathbf{Z})\otimes(\mathbf{Y}\mathbf{U})$, we have

$$\mathbf{g}(\tilde{\mathbf{q}})\otimes\mathbf{g}(\tilde{\mathbf{q}}) =\big{(}\bar{\mathbf{q}}+\mathbf{V}\tilde{\mathbf{q}}+\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]\big{)}\otimes\big{(}\bar{\mathbf{q}}+\mathbf{V}\tilde{\mathbf{q}}+\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]\big{)}$$

$$=\bar{\mathbf{q}}\otimes\bar{\mathbf{q}}+\bar{\mathbf{q}}\otimes(\mathbf{V}\tilde{\mathbf{q}})+\bar{\mathbf{q}}\otimes(\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}])$$

$$\qquad+(\mathbf{V}\tilde{\mathbf{q}})\otimes\bar{\mathbf{q}}+(\mathbf{V}\tilde{\mathbf{q}})\otimes(\mathbf{V}\tilde{\mathbf{q}})+(\mathbf{V}\tilde{\mathbf{q}})\otimes(\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}])$$

$$\qquad+(\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}])\otimes\bar{\mathbf{q}}+(\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}])\otimes(\mathbf{V}\tilde{\mathbf{q}})+(\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}])\otimes(\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}])$$

$$=\bar{\mathbf{q}}\otimes\bar{\mathbf{q}}+(\bar{\mathbf{q}}\otimes\mathbf{V})\tilde{\mathbf{q}}+(\bar{\mathbf{q}}\otimes\mathbf{W})[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]$$

$$\qquad+(\mathbf{V}\otimes\bar{\mathbf{q}})\tilde{\mathbf{q}}+(\mathbf{V}\otimes\mathbf{V})[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]+(\mathbf{V}\otimes\mathbf{W})[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]$$

$$\qquad+(\mathbf{W}\otimes\bar{\mathbf{q}})[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]+(\mathbf{W}\otimes\mathbf{V})[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]+(\mathbf{W}\otimes\mathbf{W})[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]$$

$$=\bar{\mathbf{q}}\otimes\bar{\mathbf{q}}+(\bar{\mathbf{q}}\otimes\mathbf{V}+\mathbf{V}\otimes\bar{\mathbf{q}})\tilde{\mathbf{q}}+(\bar{\mathbf{q}}\otimes\mathbf{W}+\mathbf{V}\otimes\mathbf{V}+\mathbf{W}\otimes\bar{\mathbf{q}})[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]$$

$$\qquad+(\mathbf{V}\otimes\mathbf{W}+\mathbf{W}\otimes\mathbf{V})[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]+(\mathbf{W}\otimes\mathbf{W})[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}].$$

Therefore,

$$\mathbf{f}(\mathbf{g}(\tilde{\mathbf{q}})) =\mathbf{A}\big{(}\bar{\mathbf{q}}+\mathbf{V}\tilde{\mathbf{q}}+\mathbf{W}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]\big{)}+\mathbf{H}[\mathbf{g}(\tilde{\mathbf{q}})\otimes\mathbf{g}(\tilde{\mathbf{q}})]$$

$$=\mathbf{A}\bar{\mathbf{q}}+\mathbf{H}[\bar{\mathbf{q}}\otimes\bar{\mathbf{q}}]+\big{(}\mathbf{A}\mathbf{V}+\mathbf{H}(\bar{\mathbf{q}}\otimes\mathbf{V}+\mathbf{V}\otimes\bar{\mathbf{q}})\big{)}\tilde{\mathbf{q}}$$

$$\qquad+\big{(}\mathbf{A}\mathbf{W}+\mathbf{H}(\bar{\mathbf{q}}\otimes\mathbf{W}+\mathbf{V}\otimes\mathbf{V}+\mathbf{W}\otimes\bar{\mathbf{q}})\big{)}[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]$$

$$\qquad+\mathbf{H}(\mathbf{V}\otimes\mathbf{W}+\mathbf{W}\otimes\mathbf{V})[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}]+\mathbf{H}(\mathbf{W}\otimes\mathbf{W})[\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}\otimes\tilde{\mathbf{q}}],$$

so that the intrusive projection-based ROM eq. 3.8 can be written as

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)&=\mathbf{V}^{\mathsf{T}}\mathbf{f}(\mathbf{g}(\tilde{\mathbf{q}}(t)))\\ &=\tilde{\mathbf{c}}+\tilde{\mathbf{A}}\tilde{\mathbf{q}}(t)+\tilde{\mathbf{H}}_{2}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)]+\tilde{\mathbf{H}}_{3}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)]+\tilde{\mathbf{H}}_{4}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)],\end{aligned}$$ (A.1a)

where

$$\begin{aligned} \tilde{\mathbf{c}}&=\mathbf{V}^{\mathsf{T}}(\mathbf{A}\bar{\mathbf{q}}+\mathbf{H}[\bar{\mathbf{q}}\otimes\bar{\mathbf{q}}])\in{}^{r},\\ \tilde{\mathbf{A}}&=\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{V}+\mathbf{V}^{\mathsf{T}}\mathbf{H}(\bar{\mathbf{q}}\otimes\mathbf{V}+\mathbf{V}\otimes\bar{\mathbf{q}})\in{}^{r\times r},\\ \tilde{\mathbf{H}}_{2}&=\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{W}+\mathbf{V}^{\mathsf{T}}\mathbf{H}(\bar{\mathbf{q}}\otimes\mathbf{W}+\mathbf{V}\otimes\mathbf{V}+\mathbf{W}\otimes\bar{\mathbf{q}})\in{}^{r\times r^{2}},\\ \tilde{\mathbf{H}}_{3}&=\mathbf{V}^{\mathsf{T}}\mathbf{H}(\mathbf{V}\otimes\mathbf{W}+\mathbf{W}\otimes\mathbf{V})\in{}^{r\times r^{3}},\\ \tilde{\mathbf{H}}_{4}&=\mathbf{V}^{\mathsf{T}}\mathbf{H}(\mathbf{W}\otimes\mathbf{W})\in{}^{r\times r^{4}}.\end{aligned}$$ (A.1b)

The quartic polynomial structure of eq. A.1 also arises when $\bar{\mathbf{q}}=\bf 0$ but a Kernel ROM is constructed with the input scaling preprocessing step of Remark 2.1. In that case, the matrices in eq. A.1 reduce to

$$\tilde{\mathbf{A}} =\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{V}, \tilde{\mathbf{H}}_{2} =\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{W}+\mathbf{V}^{\mathsf{T}}\mathbf{H}(\mathbf{V}\otimes\mathbf{V}),$$

$$\tilde{\mathbf{H}}_{3} =\mathbf{V}^{\mathsf{T}}\mathbf{H}(\mathbf{V}\otimes\mathbf{W}+\mathbf{W}\otimes\mathbf{V}), \tilde{\mathbf{H}}_{4} =\mathbf{V}^{\mathsf{T}}\mathbf{H}(\mathbf{W}\otimes\mathbf{W}),$$

with $\tilde{\mathbf{c}}=\bf 0$. However, the Kernel ROM targets a shifted and scaled reduced state $\hat{\mathbf{q}}(t)={\boldsymbol{\Sigma}}^{-1}(\tilde{\mathbf{q}}(t)-\bar{\mathbf{x}})$ for some ${\boldsymbol{\Sigma}}\in{}^{r\times r}$ and $\bar{\mathbf{x}}\in{}^{r}$, which evolves according to

$$\frac{\textrm{d}}{\textrm{d}t}\hat{\mathbf{q}}(t) =\frac{\textrm{d}}{\textrm{d}t}\left[{\boldsymbol{\Sigma}}^{-1}\tilde{\mathbf{q}}(t)-{\boldsymbol{\Sigma}}^{-1}\bar{\mathbf{x}}\right]$$

$$={\boldsymbol{\Sigma}}^{-1}\big{(}\tilde{\mathbf{A}}\tilde{\mathbf{q}}(t)+\tilde{\mathbf{H}}_{2}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)]+\tilde{\mathbf{H}}_{3}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)]+\tilde{\mathbf{H}}_{4}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)]\big{)}$$

$$={\boldsymbol{\Sigma}}^{-1}\big{(}\tilde{\mathbf{A}}({\boldsymbol{\Sigma}}\hat{\mathbf{q}}(t)+\bar{\mathbf{x}})+\tilde{\mathbf{H}}_{2}[({\boldsymbol{\Sigma}}\hat{\mathbf{q}}(t)+\bar{\mathbf{x}})\otimes({\boldsymbol{\Sigma}}\hat{\mathbf{q}}(t)+\bar{\mathbf{x}})]$$

$$\qquad\qquad+\tilde{\mathbf{H}}_{3}[({\boldsymbol{\Sigma}}\hat{\mathbf{q}}(t)+\bar{\mathbf{x}})\otimes({\boldsymbol{\Sigma}}\hat{\mathbf{q}}(t)+\bar{\mathbf{x}})\otimes({\boldsymbol{\Sigma}}\hat{\mathbf{q}}(t)+\bar{\mathbf{x}})]$$

$$\qquad\qquad+\tilde{\mathbf{H}}_{4}[({\boldsymbol{\Sigma}}\hat{\mathbf{q}}(t)+\bar{\mathbf{x}})\otimes({\boldsymbol{\Sigma}}\hat{\mathbf{q}}(t)+\bar{\mathbf{x}})\otimes({\boldsymbol{\Sigma}}\hat{\mathbf{q}}(t)+\bar{\mathbf{x}})\otimes({\boldsymbol{\Sigma}}\hat{\mathbf{q}}(t)+\bar{\mathbf{x}})]\big{)}$$

$$=\hat{\mathbf{c}}+\hat{\mathbf{A}}\hat{\mathbf{q}}(t)+\hat{\mathbf{H}}_{2}[\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)]+\hat{\mathbf{H}}_{3}[\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)]+\hat{\mathbf{H}}_{4}[\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)\otimes\hat{\mathbf{q}}(t)],$$

where

$$\hat{\mathbf{c}} ={\boldsymbol{\Sigma}}^{-1}\big{(}\tilde{\mathbf{A}}\bar{\mathbf{x}}+\tilde{\mathbf{H}}_{2}[\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}]+\tilde{\mathbf{H}}_{3}[\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}]+\tilde{\mathbf{H}}_{4}[\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}]\big{)}\in{}^{r},$$

$$\hat{\mathbf{A}} ={\boldsymbol{\Sigma}}^{-1}\big{(}\tilde{\mathbf{A}}{\boldsymbol{\Sigma}}+\tilde{\mathbf{H}}_{2}({\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}+\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}})+\tilde{\mathbf{H}}_{3}({\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}+\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}+\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}})$$

$$\qquad\qquad+\tilde{\mathbf{H}}_{4}({\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}+\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}+\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}+\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}})\big{)}\in{}^{r\times r},$$

$$\hat{\mathbf{H}}_{2} ={\boldsymbol{\Sigma}}^{-1}\big{(}\tilde{\mathbf{H}}_{2}({\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}})+\tilde{\mathbf{H}}_{3}({\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}+{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}+\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}})$$

$$\qquad\qquad+\tilde{\mathbf{H}}_{4}({\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}+{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}+{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}$$

$$\qquad\qquad\qquad+\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}+\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}+\bar{\mathbf{x}}\otimes\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}})\big{)}\in{}^{r\times r^{2}},$$

$$\hat{\mathbf{H}}_{3} ={\boldsymbol{\Sigma}}^{-1}\big{(}\tilde{\mathbf{H}}_{3}({\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}})$$

$$\qquad\qquad+\tilde{\mathbf{H}}_{4}({\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}+{\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}+{\boldsymbol{\Sigma}}\otimes\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}}+\bar{\mathbf{x}}\otimes{\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}})\big{)}\in{}^{r\times r^{3}},$$

$$\hat{\mathbf{H}}_{4} ={\boldsymbol{\Sigma}}^{-1}\tilde{\mathbf{H}}_{4}({\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}}\otimes{\boldsymbol{\Sigma}})\in{}^{r\times r^{4}}.$$

The salient point is that none of these matrices need to be constructed explicitly when using a non-intrusive model reduction method: only the desired structure is needed to design the non-intrusive ROM.

Appendix B Stability for linear systems

The following stability result illustrates the importance of the regularization hyperparameter $\rho\geq 0$ when solving the minimization problem eq. 3.12 for computing $\mathbf{W}$. Applying the QM approach with reference state $\bar{\mathbf{q}}=\mathbf{0}$ to a linear FOM

$$\frac{\textrm{d}}{\textrm{d}t}\mathbf{q}(t)=\mathbf{A}\mathbf{q}(t),\qquad\mathbf{q}(0)=\mathbf{q}_{0}({\boldsymbol{\mu}}),$$ (B.1)

results in a ROM with quadratic dynamics

$$\frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)=\tilde{\mathbf{A}}\tilde{\mathbf{q}}(t)+\tilde{\mathbf{H}}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)],\qquad\tilde{\mathbf{q}}(0) =\mathbf{V}^{\mathsf{T}}\mathbf{q}_{0}({\boldsymbol{\mu}}),$$ (B.2)

where $\tilde{\mathbf{A}}=\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{V}$ and $\tilde{\mathbf{H}}=\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{W}$. We then have a stability estimate for the ROM solution

Proposition B.1.

Let $\lambda$ denote the maximum eigenvalue of $\mathbf{A}_{sym}=\frac{1}{2}(\mathbf{A}+\mathbf{A}^{\mathsf{T}})$ (the symmetric part of $\mathbf{A}$). Then the following stability estimate for the QM ROM eq. B.2 holds for all $t\in[0,T]$:

$$\left\|\tilde{\mathbf{q}}(t)\right\|_{2}\leq\left\|\mathbf{A}\right\|_{2}\left\|\mathbf{W}\right\|_{2}\int_{0}^{t}\|\tilde{\mathbf{q}}(s)\otimes\tilde{\mathbf{q}}(s)\|_{2}e^{\lambda(t-s)}ds+e^{\lambda t}\left\|\mathbf{V}^{\mathsf{T}}\mathbf{q}_{0}\right\|_{2}.$$ (B.3)

Proof.

Observe that

$$\begin{aligned} \tilde{\mathbf{q}}(t)^{\mathsf{T}}\frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)&=\tilde{\mathbf{q}}(t)^{\mathsf{T}}\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{V}\tilde{\mathbf{q}}(t)+\tilde{\mathbf{q}}(t)^{\mathsf{T}}\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{W}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)]\\ &=\tilde{\mathbf{q}}(t)^{\mathsf{T}}\mathbf{V}^{\mathsf{T}}\mathbf{A}_{sym}\mathbf{V}\tilde{\mathbf{q}}(t)+\tilde{\mathbf{q}}(t)^{\mathsf{T}}\mathbf{V}^{\mathsf{T}}\mathbf{A}\mathbf{W}[\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)]\\ &\leq\lambda\left\|\mathbf{V}\tilde{\mathbf{q}}(t)\right\|_{2}^{2}+\left\|\mathbf{A}\right\|_{2}\left\|\mathbf{W}\right\|_{2}\left\|\mathbf{V}\tilde{\mathbf{q}}(t)\right\|_{2}\|\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\|_{2}\\ &=\lambda\left\|\tilde{\mathbf{q}}(t)\right\|_{2}^{2}+\left\|\mathbf{A}\right\|_{2}\left\|\mathbf{W}\right\|_{2}\left\|\tilde{\mathbf{q}}(t)\right\|_{2}\|\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\|_{2},\end{aligned}$$ (B.4)

where the last line follows from the orthonormality of $\mathbf{V}$. The bound eq. B.4 implies that

$$\frac{\textrm{d}}{\textrm{d}t}\left\|\tilde{\mathbf{q}}(t)\right\|_{2}=\frac{\tilde{\mathbf{q}}(t)^{\mathsf{T}}\frac{\textrm{d}}{\textrm{d}t}\tilde{\mathbf{q}}(t)}{\left\|\tilde{\mathbf{q}}(t)\right\|_{2}}\leq\lambda\left\|\tilde{\mathbf{q}}(t)\right\|_{2}+\left\|\mathbf{A}\right\|_{2}\left\|\mathbf{W}\right\|_{2}\|\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\|_{2}$$

Applying Lemma 5.1 with $u(t)=\left\|\tilde{\mathbf{q}}(t)\right\|_{2}$, $\alpha(t)=\left\|\mathbf{A}\right\|_{2}\left\|\mathbf{W}\right\|_{2}\|\tilde{\mathbf{q}}(t)\otimes\tilde{\mathbf{q}}(t)\|_{2}$, and $\beta(t)=\lambda$ yields the result. ∎

Proposition B.1 indicates that the magnitude of $\left\|\mathbf{W}\right\|_{2}$ has a crucial impact on the stability of the resulting QM ROM. Consequently, it is important to apply sufficient regularization (i.e., choose $\rho$ large enough) when computing $\mathbf{W}$ to ensure that $\left\|\mathbf{W}\right\|_{2}$ remains small.

Appendix C Main nomenclature

Notation

Kernel interpolation
$K:{}^{n_{x}}\times{}^{n_{x}}\to\real$	symmetric kernel function
${\cal H}_{K}$	reproducing kernel Hilbert space
$\mathbf{x}_{j}\in{}^{n_{x}}$	inputs for kernel interpolation
$\mathbf{y}_{j}\in{}^{n_{y}}$	outputs for kernel interpolation
$\mathbf{v}:{}^{n_{x}}\to{}^{n_{y}}$	function to interpolate: $\mathbf{y}_{j}=\mathbf{v}(\mathbf{x}_{j})$
$\gamma\geq 0$	kernel regularization parameter
${\boldsymbol{\Omega}}\in{}^{m\times n_{y}}$	coefficient matrix for kernel interpolation
$\mathbf{s}_{\mathbf{v}}^{\gamma}\in{\cal H}_{K}^{n_{y}}$	kernel interpolant of $\mathbf{v}$ with regularization $\gamma$
${\boldsymbol{\psi}}_{\!\epsilon}:{}^{n_{x}}\to{}^{m}$	RBF kernel evaluation function
${\boldsymbol{\phi}}:{}^{n_{x}}\to{}^{n_{\phi}}$	feature map
$\mathbf{G}\in{}^{n_{\phi}\times n_{\phi}}$	weighting matrix for feature map kernels
$\mathbf{C}\in{}^{n_{y}\times n_{\phi}}$	post-feature map kernel coefficients
Full-order models
$\mathbf{q}(t)\in{}^{n_{q}}$	full-order model state
$\mathbf{f}:{}^{n_{q}}\to{}^{n_{q}}$	full-order model dynamics function
${\boldsymbol{\mu}}\in{}^{n_{\mu}}$	parameters for the initial condition
$\mathbf{Q}\in{}^{n_{q}\times M(n_{t}+1)}$	shifted state snapshot matrix (all trajectories)
Reduced-order models
$\tilde{\mathbf{q}}(t)\in{}^{r}$	intrusive reduced-order model state
$\hat{\mathbf{q}}(t)\in{}^{r}$	non-intrusive reduced-order model state
$\mathbf{V}\in{}^{n_{q}\times r}$	proper orthogonal decomposition (POD) basis matrix
$\mathbf{W}\in{}^{r\times r(r+1)/2}$	quadratic manifold (QM) weight matrix
$\mathbf{g}:{}^{r}\to{}^{n_{q}}$	decompression map
$\mathbf{h}:{}^{n_{q}}\to{}^{r}$	compression map
$\mathbf{e},\hat{\mathbf{e}}$	error quantities