Kernel manifolds: nonlinear-augmentation dimensionality reduction using reproducing kernel Hilbert spaces

TL;DR

This paper introduces a kernel-based nonlinear augmentation method for dimensionality reduction that leverages reproducing kernel Hilbert spaces to improve approximation accuracy with low training cost.

Contribution

It generalizes quadratic manifold methods by learning an optimal nonlinear correction in a reproducing kernel Hilbert space, allowing flexible nonlinear structures and improved accuracy.

Findings

Monotonically decreasing error with increasing latent space dimension

Lower training cost compared to existing methods

Outperforms proper orthogonal decomposition and recent quadratic manifold approaches

Abstract

This paper generalizes recent advances on quadratic manifold (QM) dimensionality reduction by developing kernel methods-based nonlinear-augmentation dimensionality reduction. QMs, and more generally feature map-based nonlinear corrections, augment linear dimensionality reduction with a nonlinear correction term in the reconstruction map to overcome approximation accuracy limitations of purely linear approaches. While feature map-based approaches typically learn a least-squares optimal polynomial correction term, we generalize this approach by learning an optimal nonlinear correction from a user-defined reproducing kernel Hilbert space. Our approach allows one to impose arbitrary nonlinear structure on the correction term, including polynomial structure, and includes feature map and radial basis function-based corrections as special cases. Furthermore, our method has relatively low training cost and has monotonically decreasing error as the latent space dimension increases. We compare our approach to proper orthogonal decomposition and several recent QM approaches on data from several example problems.

Figures

Tables

Table 1. Table 1 : Examples of RBF kernels.

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

Table 2. Table 2 : Fixed regularization values for different dimensionality reduction approaches for the 2D advection-diffusion-reaction example.

Manifold type	regularization $\lambda$
Kernel QM	${10}^{-2}$
Kernel RBF	${10}^{-11}$
Alternating QM	${10}^0$
Greedy QM	${10}^0$

Table 3. Table 3 : Fixed regularization values for different dimensionality reduction approaches for the surface heating example.

Manifold type	regularization $\lambda$
Kernel QM	${10}^{-9}$
Kernel RBF	${10}^{-9}$
Alternating QM	${10}^{-5}$
Greedy QM	${10}^{-4}$

Table 4. Table 4 : Fixed regularization λ \lambda values for different dimensionality reduction approaches in 3D solid mechanics flexible bracket example for the displacement, velocity, and acceleration fields.

Manifold type	Displacement reg.	Velocity reg.	Acceleration reg.
Kernel QM	${10}^{-7}$	${10}^{-7}$	${10}^{-7}$
Kernel RBF	${10}^{-7}$	${10}^{-7}$	${10}^{-7}$
Alternating QM	${10}^{-11}$	${10}^6$	${10}^8$
Greedy QM	${10}^{-11}$	NA	NA

Table 5. Table 5 : Fixed regularization values for different dimensionality reduction approaches for the 2D Euler double Mach reflection example.

Manifold type	regularization $\lambda$
Kernel QM	${10}^6$
Kernel RBF	${10}^{-3}$
Greedy QM	${10}^{10}$
Alternating QM	${10}^8$

Equations

$$\mathcal{O}:\mathbf{\mu }\mapsto \mathbf{q}(\mathbf{\mu }){\in }^N.$$

$$\mathcal{Q}=\mathcal{O}(\mathcal{D}){\subset }^N.$$

$$\mathbf{g}(\mathbf{h}(\mathbf{q}))\approx \mathbf{q}\forall \mathbf{q}\in \mathcal{Q}.$$

$$\mathbf{e}(\mathbf{q})=\mathbf{q}-\mathbf{g}(\mathbf{h}(\mathbf{q}))$$

$$\tilde{\mathbf{Q}}=\mathbf{Q}-\bar{\mathbf{q}}{1}^{\mathsf{T}},$$

$$\bar{\mathbf{q}}=\frac{1}{M}\sum _{j=1}^M{\mathbf{q}}_j,$$

$$\mathbf{U}\mathbf{\Sigma }{\mathbf{Y}}^{\mathsf{T}}=\tilde{\mathbf{Q}},$$

$$\mathbf{V}=\mathbf{U}(:,1:r).$$

$$\underset{\mathbf{V}{\in }^{N\times r},{\mathbf{V}}^{\mathsf{T}}\mathbf{V}={\mathbf{I}}_r}{\min }{\Vert \tilde{\mathbf{Q}}-\mathbf{V}{\mathbf{V}}^{\mathsf{T}}\tilde{\mathbf{Q}}\Vert }_F^2,$$

$$\mathbf{h}(\mathbf{q})$$

$$\mathbf{g}(\hat{\mathbf{q}})=\bar{\mathbf{q}}+\mathbf{V}\hat{\mathbf{q}}+\overline{\mathbf{V}}\mathbf{\Xi }\mathbf{\varphi }(\hat{\mathbf{q}}),$$

$$\underset{\mathbf{\Xi }{\in }^{m\times n_{\varphi }}}{\min }\sum _{j=1}^M{\Vert {\mathbf{q}}_j-\left(\bar{\mathbf{q}}+\mathbf{V}{\mathbf{V}}^{\mathsf{T}}({\mathbf{q}}_j-\bar{\mathbf{q}})+\overline{\mathbf{V}}\mathbf{\Xi }\mathbf{\varphi }({\mathbf{V}}^{\mathsf{T}}({\mathbf{q}}_j-\bar{\mathbf{q}}))\right)\Vert }_2^2+\lambda {\Vert \mathbf{\Xi }\Vert }_F^2,$$

$${\mathbf{V}}^{\mathsf{T}}\mathbf{V}$$

$${\hat{\mathbf{q}}}_j$$

$$\underset{\mathbf{\Xi }{\in }^{m\times n_{\varphi }}}{\min }\sum _{j=1}^M{\Vert {\hat{\mathbf{p}}}_j-\mathbf{\Xi }\mathbf{\varphi }({\hat{\mathbf{q}}}_j)\Vert }_2^2+\lambda {\Vert \mathbf{\Xi }\Vert }_F^2.$$

$$\hat{\mathbf{Q}}$$

$$\left(\mathbf{\varphi }(\hat{\mathbf{Q}})\mathbf{\varphi }(\hat{\mathbf{Q}}{)}^{\mathsf{T}}+\lambda {\mathbf{I}}_{n_{\varphi }}\right){\mathbf{\Xi }}^{\mathsf{T}}=\mathbf{\varphi }(\hat{\mathbf{Q}}){\hat{\mathbf{P}}}^{\mathsf{T}},$$

$$\mathbf{V}$$

$$\mathbf{n}$$

$$\mathbf{\varphi }(\hat{\mathbf{q}})={\left[\begin{array}{cccccc}{\hat{\mathbf{q}}}_1^2& {\hat{\mathbf{q}}}_2{\hat{\mathbf{q}}}_1& {\hat{\mathbf{q}}}_2^2& \dots & {\hat{\mathbf{q}}}_n{\hat{\mathbf{q}}}_{n-1}& {\hat{\mathbf{q}}}_n^2\end{array}\right]}^{\mathsf{T}}{\in }^{n(n+1)/2}.$$

$$\hat{\mathbf{q}}$$

$${\hat{\mathbf{q}}}^{\prime }$$

$$\mathbf{g}(\hat{\mathbf{q}})=\bar{\mathbf{q}}+\mathbf{V}\hat{\mathbf{q}}+\overline{\mathbf{V}}\mathbf{n}(\hat{\mathbf{q}}),$$

$$\displaystyle{\cal H}_{K}^{0}(^{d})=\left\{f:^{d}\to~\bigg{|}~\exists\,n\in\mathbb{N},\,\mathbf{w}\in^{n},\,\left\{\mathbf{x}_{j}\right\}_{j=1}^{n}\subset^{d}~\textup{such that}~f(\mathbf{x})=\sum_{j=1}^{n}w_{j}K(\mathbf{x}_{j},\mathbf{x})\right\}$$

$${⟨f,f^{\prime }⟩}_{{\mathcal{H}}_K{(}^d)}:=\sum _{j=1}^n\sum _{k=1}^{n^{\prime }}{w}_j{w}_k^{\prime }K({\mathbf{x}}_j,{\mathbf{x}}_k^{\prime }),$$

$$\underset{\mathbf{s}\in {\mathcal{H}}_K^p}{\min }\sum _{j=1}^n{\Vert {\mathbf{y}}_j-\mathbf{s}({\mathbf{x}}_j)\Vert }_2^2+\lambda {\Vert \mathbf{s}\Vert }_{{\mathcal{H}}_K^p}^2.$$

$${\mathbf{s}}_{\mathbf{f}}^{\lambda }(\mathbf{x})={\mathbf{\Omega }}^{\mathsf{T}}K(\mathbf{X},\mathbf{x}),$$

$$(K(\mathbf{X},\mathbf{X})+\lambda \mathbf{I})\mathbf{\Omega }={\mathbf{Y}}^{\mathsf{T}}.$$

$$\underset{\mathbf{n}\in {\mathcal{H}}_K^m}{\min }\sum _{j=1}^M{\Vert {\mathbf{q}}_j-\left(\bar{\mathbf{q}}+\mathbf{V}{\mathbf{V}}^{\mathsf{T}}({\mathbf{q}}_j-\bar{\mathbf{q}})+\overline{\mathbf{V}}\mathbf{n}({\mathbf{V}}^{\mathsf{T}}({\mathbf{q}}_j-\bar{\mathbf{q}}))\right)\Vert }_2^2+\lambda {\Vert \mathbf{n}\Vert }_{{\mathcal{H}}_K^m}^2,$$

$${\hat{\mathbf{q}}}_j$$

Full text

Keywords: nonlinear dimensionality reduction, quadratic manifolds, kernel methods, surrogate modeling

1 Introduction

In recent years, data-driven surrogate modeling of complex physical systems has been an active research area in scientific machine learning (SciML). A common challenge of constructing surrogate models is that data from simulations of partial differential equations (PDEs) are often high dimensional, with the dimensionality dictated by the discretization methods used and the fidelity of the corresponding mesh. Consequently, surrogate modeling techniques typically involve a dimensionality reduction component to extract low-dimensional features from high-dimensional simulation data, which can then be used to construct a low-dimensional, computationally efficient surrogate model. Crucially, the dimensionality reduction technique used fundamentally limits the achievable accuracy of the resulting surrogate model.

Classical linear dimensionality reduction techniques, such as proper orthogonal decomposition (POD) [19, 21], seek to identify a low-dimensional affine subspace that best captures the high-dimensional data. POD identifies this subspace as the span of the leading left singular vectors of a representative set of high-dimensional snapshot data. The POD basis induces an affine encoder map that projects high-dimensional snapshots onto generalized coordinates that lie in a low-dimensional latent space, which can then be mapped to an approximate reconstruction of the high-dimensional snapshot via an affine decoder map. While POD is effective in numerous applications, it is also well-known to poorly approximate data from problems with slowly decaying Kolmogorov $n$-width [32], such as transport-dominated problems or problems with sharp gradients. In these cases, a large POD basis is required to accurately approximate the given data. However, despite lying in a high-dimensional affine subspace, these data may still possess an inherent low-dimensional nonlinear structure. This assumption motivates using nonlinear dimensionality reduction for problems with slowly decaying Kolmogorov $n$-width.

A number of nonlinear dimensionality reduction approaches have recently been proposed to circumvent the Kolmogorov $n$-width barrier. These techniques are often collectively referred to as nonlinear manifolds to distinguish them from linear dimensionality reduction, and include autoencoders (AEs) [28, 24, 36, 9, 11, 17, 8, 7, 13], quadratic manifolds (QMs) [3, 6, 15, 14, 43, 42, 22, 37], the projection-based reduced-order model + artificial neural network (PROM-ANN) technique [4], and the recent Gaussian process (GP) and radial basis function (RBF) interpolation-based nonlinear closure approach [10]. AEs consist of fully nonlinear encoder and decoder maps, which can offer substantial accuracy improvements over linear dimensionality reduction. However, they can be difficult to train due to their dependence on the neural network architecture and large data requirements. Furthermore, they often lack interpretable structure and can be slow to evaluate, which is detrimental in use cases requiring fast online evaluations, such as reduced-order modeling. To address these issues, the QM, PROM-ANN, and GP/RBF-interpolation approaches instead augment linear dimensionality reduction with a nonlinear correction term in the decoder. We characterize this class of approaches as nonlinear-augmentation manifolds. QMs model the nonlinear augmentation term as a quadratic polynomial, while PROM-ANN models it as an ANN and the GP/RBF-interpolation approach models it as the mean function of a GP or as an RBF interpolant. These approaches have interpretable structure and, with the exception of PROM-ANN, comparatively low training cost relative to AEs, while still providing improved accuracy over linear dimensionality reduction approaches. However, QMs, and more generally, feature-map based nonlinear-augmentation manifolds, are typically limited to polynomial corrections, which can restrict the achievable accuracy of these methods, while becoming increasingly difficult to train as the polynomial order increases. Although using GPs and RBF interpolants as in [10] can address these difficulties, this approach can be further generalized to improve the flexibility of nonlinear-augmentation manifolds. In this paper we develop a nonlinear-augmentation dimensionality reduction technique using regularized kernel interpolation, which generalizes recent work on QMs and GP/RBF interpolation-based nonlinear closures, yielding greater expressivity in the nonlinear approximation, reduced training cost compared to QMs, and improved accuracy on out-of-sample data.

Regularized kernel interpolation is a generalization of RBF interpolation relying on mature mathematical theory from reproducing kernel Hilbert spaces (RKHS) [39, 40, 30, 46]. While kernel methods are a common tool in classical data science and machine learning, they have become increasingly popular in SciML in recent years, with notable applications in operator and equation learning [5, 29, 23], tensor decomposition [27], and reduced-order and surrogate modeling [12, 2]. In our dimensionality reduction method, which we refer to as kernel manifold dimensionality reduction, one can model the nonlinear augmentation term with arbitrary nonlinearities by choosing an appropriate kernel, and includes QMs and GP/RBF-based closures as special cases. Furthermore, unlike feature map-based nonlinear augmentations, regularized kernel interpolation allows one to model the nonlinear term using polynomials of arbitrarily high order without increasing the number of degrees of freedom in the learning problem. We also note that our approach is distinct from kernel principal component analysis (KPCA) [41] and kernel POD (KPOD) [38] in that it augments classical POD with a nonlinear correction term, whereas KPCA/KPOD returns a nonlinear kernel-based encoder without an explicit decoder, thus requiring one to learn an approximate decoder function through, e.g., kernel interpolation or neural networks. While we purely focus on dimensionality reduction in this paper, we emphasize that our method is general-purpose, and can be used for standalone dimensionality reduction, intrusive model reduction as in [10], or non-intrusive model reduction through integration with operator inference (OpInf) [16, 25, 34, 35, 14] or the recent kernel methods-based approach of [12].

The outline of this paper is as follows. We first describe the problem setup and review POD and feature map-based nonlinear-augmentation dimensionality reduction in Section 2. Next, we describe our new kernel manifold approach in Section 3, beginning with a brief review of regularized kernel interpolation in Section 3.1. We then compare our method with POD and the recently proposed Alternating QM [14] and Greedy QM [43] approaches on data from several example problems in Section 4. The results show that our approach can yield improved accuracy over POD and QM in the small latent dimension size regime in several cases while requiring reduced training cost compared to Greedy QM. Finally, Section 5 provides several concluding remarks and directions for future work.

2 Review of feature map-based nonlinear-augmentation manifolds

This section reviews a class of nonlinear manifold dimensionality reduction that augments linear dimensionality reduction with a nonlinear correction term in the decoder map. These approaches are typically applied to data from large-scale physics simulations. We first discuss a unified problem setting in Section 2.1. We then briefly review the proper orthogonal decomposition (POD) in Section 2.2, since nonlinear-augmentation manifolds approaches use POD as a starting point. Lastly, we review the feature map-based nonlinear augmentation approach in Section 2.3.

2.1 Problem setting

Let ${\cal D}\subset^{p}$ denote a set of parameters for a given physical system and let ${\cal O}:{\cal D}\to^{N}$ denote an operator that maps a parameter ${\boldsymbol{\mu}}\in{\cal D}$ to a discretized state of the physical system on a mesh of size $N$:

[TABLE]

For example, ${\cal O}$ can represent the solution operator of a discretized PDE. We seek a low-dimensional representation of the solution space of the operator ${\cal O}$,

[TABLE]

In particular, we wish to construct an encoder map $\mathbf{h}:{\cal Q}\to{\widehat{\cal Q}}\subset^{r}$ and a decoder map $\mathbf{g}:{\widehat{\cal Q}}\to{\cal Q},$ where $r\ll N$, such that

[TABLE]

In this setting, ${\cal Q}\subset^{N}$ is referred to as the full space and ${\widehat{\cal Q}}\subset^{n}$ is referred to as the latent space. We will commonly use the projection error

[TABLE]

to assess the quality of the encoder and decoder. In the following sections, we review methods for constructing $\mathbf{h}$ and $\mathbf{g}$ using data samples from the full space ${\cal Q}$.

2.2 Proper orthogonal decomposition

POD computes an affine encoder $\mathbf{h}$ and affine decoder $\mathbf{g}$ optimal affine subspace $\bar{\mathbf{q}}+{\cal V}\subset^{N}$ that best approximates ${\cal Q}$, where $\bar{\mathbf{q}}\in^{N}$ is a fixed offset vector and ${\cal V}\subset^{N}$ is a linear subspace of dimension $r\ll N$. We review the essential details here and refer the reader to, e.g.,  [21, 19] for a more complete overview.

Let $\mathbf{Q}=\begin{bmatrix}\mathbf{q}_{1},\dots,\mathbf{q}_{M}\end{bmatrix}\in^{N\times M}$ be a data matrix consisting of columns $\left\{\mathbf{q}_{j}\right\}_{j=1}^{M}\subset{\cal Q}$. We then shift the data $\mathbf{Q}$ by the offset vector $\bar{\mathbf{q}}$,

[TABLE]

where $\mathbf{1}\in^{N}$ is the vector of all ones. The offset $\bar{\mathbf{q}}$ can be chosen by the user, but is often taken to be the mean of the snapshots:

[TABLE]

or set to zero: $\bar{\mathbf{q}}=\mathbf{0}$. It can also be used to strongly enforce Dirichlet boundary conditions in a projection-based reduced-order model (ROM) [20]. We then compute the singular value decomposition (SVD) of $\widetilde{\mathbf{Q}}$,

[TABLE]

and define the POD basis matrix $\mathbf{V}\in^{N\times r}$ to be the first $r$ leading left singular vectors of $\widetilde{\mathbf{Q}}$,

[TABLE]

The subspace ${\cal V}$ is taken to be the span of the columns of $\mathbf{V}$, ${\cal V}={\rm span}(\mathbf{V})$. Furthermore, the basis matrix $\mathbf{V}$ minimizes

[TABLE]

and hence $\mathbf{V}$ optimally reconstructs the shifted data $\widetilde{\mathbf{Q}}$. Lastly, the encoder $\mathbf{h}$ and decoder $\mathbf{g}$ are defined as

[TABLE]

While POD can successfully reduce the dimensionality of a given dataset $\mathbf{Q}$ in a number of applications with small projection errors, it is well-known that it struggles to approximate data coming from transport-dominated problems or from problems with sharp gradients. In these cases, a large basis size $r$ is needed to reconstruct the data to achieve a desired accuracy. This issue is known as the Kolmogorov $n$-width barrier [32]. To overcome this barrier, a number of recent works have proposed that augment the affine decoder $\mathbf{g}$ from equation eq. 2.8 with a nonlinear correction term. In particular, these include quadratic manifolds (QMs) [3, 6, 15, 14, 43, 42, 22, 37], the projection-based ROM + artificial neural network (PROM-ANN) approach [4], and recent work using Gaussian processes and radial basis function (RBF) interpolation to fit nonlinear closure terms in the latent space [10]. We review a class of these approaches in the following section.

2.3 Feature map-based nonlinear-augmentation manifolds

The feature map-based nonlinear-augmentation manifold approach, which we will abbreviate to “feature map manifold” (FM), takes an affine encoder $\mathbf{h}$ of the form eq. 2.8 and augments the affine decoder map $\mathbf{g}$ with a nonlinear correction term, resulting in a decoder of the form

[TABLE]

where $\mathbf{V}\in^{N\times r}$ is a reduced basis matrix, $\overline{\mathbf{V}}\in^{N\times m}$ is an augmenting basis matrix, ${\boldsymbol{\phi}}:^{r}\to^{n_{\phi}}$ is a nonlinear feature map, and ${\boldsymbol{\Xi}}\in^{m\times n_{\phi}}$ is a coefficient matrix. The basis matrices $\mathbf{V}$ and $\overline{\mathbf{V}}$ are typically chosen beforehand, while the coefficient matrix ${\boldsymbol{\Xi}}$ is learned by minimizing the norms of the projection errors $\left\|\mathbf{q}_{j}-\mathbf{g}(\mathbf{h}(\mathbf{q}_{j}))\right\|_{2}$ of the given data $\left\{\mathbf{q}_{j}\right\}_{j=1}^{M}\subset{\cal Q}$:

[TABLE]

where $\lambda\geq 0$ is a regularization parameter.

If $\mathbf{V}$ and $\overline{\mathbf{V}}$ are chosen to be orthonormal matrices satisfying

[TABLE]

and the projections of $\mathbf{q}_{j}-\bar{\mathbf{q}}$ by $\mathbf{V}$ and $\overline{\mathbf{V}}$, respectively, are denoted as

[TABLE]

then one can show that eq. 2.10 is equivalent to

[TABLE]

With the data matrices $\widehat{\mathbf{Q}}$ and $\widehat{\mathbf{P}}$ defined as

[TABLE]

a minimizer ${\boldsymbol{\Xi}}\in^{m\times n_{\phi}}$ can be computed by solving the normal equations

[TABLE]

where ${\boldsymbol{\phi}}$ applied to a matrix is understood to operate column-wise.

A typical choice for $\mathbf{V}\in^{N\times r}$ is the same POD basis from Section 2.2 consisting of the first $r$ leading POD modes, with the augmenting basis $\overline{\mathbf{V}}\in^{N\times m}$ chosen to be the subsequent $m$ modes:

[TABLE]

see, e.g., [14, Algorithm 1]. With this choice of $\mathbf{V}$ and $\overline{\mathbf{V}}$, $\widehat{\mathbf{q}}_{j}$ from equation eq. 2.12 corresponds to the projection of $\mathbf{q}_{j}-\bar{\mathbf{q}}$ onto the first $r$ leading modes and $\widehat{\mathbf{p}}_{j}$ corresponds to the projection onto the subsequent $m$ modes. Moreover, solving eq. 2.13 is equivalent to constructing a function

[TABLE]

that approximates a map from the lower order modal coefficients $\widehat{\mathbf{q}}$ to higher order modal coefficients $\widehat{\mathbf{p}}$. We will leverage this interpretation in the kernel-based approach in Section 3.

Remark 2.1**.**

The feature map ${\boldsymbol{\phi}}:^{n}\to^{n_{\phi}}$ is often chosen to be a quadratic feature map that omits duplicated entries,

[TABLE]

Using a decoder of the form eq. 2.9 with the quadratic feature map eq. 2.18 is referred to as quadratic manifold (QM) dimensionality reduction.

Remark 2.2**.**

The approaches in [3, 15, 43] use $\overline{\mathbf{V}}=\mathbf{I}_{N}$ and construct $\mathbf{V}$ using POD modes. Since $\mathbf{V}^{\mathsf{T}}\overline{\mathbf{V}}=\mathbf{V}^{\mathsf{T}}\neq\mathbf{0}$ in this case, the simplified minimization problem defined by Equations eq. 2.12 and eq. 2.13 no longer holds, and thus eq. 2.10 must be solved directly. While [3, 15] use the $r$ leading POD modes to construct $\mathbf{V}$, the greedy QM approach in [43] greedily selects (potentially high-order) POD modes, which can lead to substantial accuracy gains; see [43] for further details.

Remark 2.3**.**

The nonlinear correction term $\overline{\mathbf{V}}{\boldsymbol{\Xi}}{\boldsymbol{\phi}}(\widehat{\mathbf{q}})$ in the decoder eq. 2.9 can also be interpreted as a closure modeling term, since $\overline{\mathbf{V}}{\boldsymbol{\Xi}}{\boldsymbol{\phi}}(\widehat{\mathbf{q}})$ approximates the effect of discarded POD modes without introducing additional degrees of freedom. However, it is important to note that an exact mapping $\mathbf{n}_{*}:\widehat{\mathbf{q}}\mapsto\widehat{\mathbf{p}}$ from low-order modes to high-order modes is not well-defined in general. Indeed, consider the case $r=m=1$ and $\bar{\mathbf{q}}=\mathbf{0}$. Then $\mathbf{V}=\mathbf{v}_{1}=\mathbf{U}(:,1)$ and $\overline{\mathbf{V}}=\mathbf{v}_{2}=\mathbf{U}(:,2).$ Suppose that $\mathbf{q}:=\mathbf{v}_{1}\in{\cal Q}$ and $\mathbf{q}^{\prime}:=\mathbf{v}_{1}+\mathbf{v}_{2}\in{\cal Q}$. Notice that

[TABLE]

However,

[TABLE]

Therefore $\mathbf{n}_{*}$ is not well-defined. Nevertheless, fitting an approximate map of the form eq. 2.17 is useful for dimensionality reduction.

While the FM approach outlined in this section is relatively general, typically ${\boldsymbol{\phi}}$ is chosen to be a low-order polynomial. Extending the approach to higher-order polynomials results in an increasingly large linear system eq. 2.15 that must be solved to compute the coefficients ${\boldsymbol{\Xi}}$. Moreover, using a high-order polynomial for ${\boldsymbol{\phi}}$ would increase the evaluation cost of the decoder $\mathbf{g}$, which can be detrimental in use-cases where $\mathbf{g}$ is evaluated many times in an online stage, such as reduced-order modeling. In the next section, we further generalize the FM approach to a kernel-based nonlinear manifold to address these issues.

3 Kernel-based nonlinear-augmentation manifolds

In the kernel-based nonlinear-augmentation manifold approach, which we will abbreviate to “kernel manifold” (KM), we again use an affine encoder $\mathbf{h}$ and augment the affine decoder with a nonlinear correction term. We express the decoder as

[TABLE]

where $\mathbf{V}\in^{N\times r}$ is a reduced basis matrix and $\overline{\mathbf{V}}\in^{N\times m}$ is an augmenting basis matrix as before, and $\mathbf{n}:^{r}\to^{m}$ is an optimal nonlinear function that is learned through regularized kernel interpolation. We briefly review regularized kernel interpolation in Section 3.1 and present the details of our new approach in Section 3.2.

3.1 Review of regularized kernel interpolation

Definition 3.1** (Positive-definite kernels).**

A function $K:^{d}\times^{d}\to$ 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^{d}$. A kernel function $K$ is said to be positive-definite if for any $\mathbf{X}=[~\mathbf{x}_{1}~~\cdots~~\mathbf{x}_{n}~]\in^{n}$ with pairwise distinct columns, the kernel matrix $K(\mathbf{X},\mathbf{X})\in^{n\times n}$ with entries $K(\mathbf{X},\mathbf{X})_{ij}=K(\mathbf{x}_{i},\mathbf{x}_{j})$ is positive semi-definite. Moreover, if $K(\mathbf{X},\mathbf{X})$ is strictly positive-definite, then $K$ is also said to be strictly positive-definite.

Definition 3.2** (Reproducing kernel Hilbert space (RKHS)).**

Consider the pre-Hilbert space

[TABLE]

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

[TABLE]

where $f(\mathbf{x})=\sum_{j=1}^{n}w_{j}K(\mathbf{x}_{j},\mathbf{x})$ and $f^{\prime}(\mathbf{x})=\sum_{j=1}^{n^{\prime}}w_{k}^{\prime}K(\mathbf{x}_{k}^{\prime},\mathbf{x})$.

For a set of pairwise-distinct vectors $\{\mathbf{x}_{j}\}_{j=1}^{n}\subset^{d}$, we use the notation $K(\mathbf{X},\mathbf{X})\in^{n\times n}$ to denote the kernel matrix of Definition 3.1 and define the vector $K(\mathbf{X},\mathbf{x})=[~K(\mathbf{x}_{1},\mathbf{x})~~\cdots~~K(\mathbf{x}_{n},\mathbf{x})~]^{\mathsf{T}}\in^{n}$. To simplify notation, we write ${\cal H}_{K}$ instead of ${\cal H}_{K}(^{d})$ when it is clear that $K$ is defined over ${}^{d}\times^{d}$.

While the RKHS ${\cal H}_{K}$ is a Hilbert space of scalar-valued functions, one can easily construct an RKHS of vector-valued functions by taking the $p$-fold Cartesian product of ${\cal H}_{K}$, i.e. ${\cal H}_{K}^{p}={\cal H}_{K}\times\dots\times{\cal H}_{K}$, with inner product $\left\langle\mathbf{f},\mathbf{f}^{\prime}\right\rangle_{{\cal H}_{K}^{p}}=\sum_{i=1}^{p}\left\langle f_{i},f_{i}^{\prime}\right\rangle_{{\cal H}_{K}}$ for all $\mathbf{f}=(f_{1},\ldots,f_{p}),\mathbf{f}^{\prime}=(f_{1}^{\prime},\ldots,f_{p}^{\prime})\in{\cal H}_{K}^{p}$. We now state the Representer Theorem, a key result from the RKHS literature underlying our methodology.

Definition 3.3** (Regularized kernel interpolant).**

Let $\mathbf{f}:^{d}\to^{p}$, $\left\{\mathbf{x}_{j}\right\}_{j=1}^{n}\subset^{d}$ be pairwise distinct, and denote $\mathbf{y}_{j}=\mathbf{f}(\mathbf{x}_{j})$. For an RKHS ${\cal H}_{K}^{p}$ and regularization parameter $\lambda\geq 0$, a regularized interpolant $\mathbf{s}_{\mathbf{f}}^{\lambda}\in{\cal H}_{K}^{p}$ of $\mathbf{f}$ solves the minimization problem

[TABLE]

Theorem 3.1** (Representer Theorem).**

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

[TABLE]

where the coefficient matrix ${\boldsymbol{\Omega}}\in^{n\times p}$ solves the linear system

[TABLE]

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

This result is a straightforward extension of the Representer Theorem for scalar-valued functions (see, e.g., [40, Theorem 9.3]). A key observation from Theorem 3.1 is that a solution to the infinite-dimensional minimization problem eq. 3.2 can be obtained by solving the finite-dimensional linear system eq. 3.4.

3.2 Learning kernel-based nonlinear augmentations

For our KM approach, we form the decoder eq. 3.1 where the nonlinear function $\mathbf{n}:^{r}\to^{m}$ belongs to the RKHS ${\cal H}_{K}^{m}={\cal H}_{K}^{m}(^{r})$ for a given positive-definite kernel $K:^{r}\times^{r}\to$. In particular, we use the optimal $\mathbf{n}\in{\cal H}_{K}^{m}$ that minimizes the norms of the projection errors $\left\|\mathbf{q}_{j}-\mathbf{g}(\mathbf{h}(\mathbf{q}_{j}))\right\|_{2}$ of the given data $\left\{\mathbf{q}_{j}\right\}_{j=1}^{M}\subset{\cal Q}$:

[TABLE]

where $\lambda\geq 0$ is a regularization parameter.

For our approach, we select $\mathbf{V}\in^{N\times r}$ to be the POD basis from Section 2.2 consisting of the first $r$ leading POD modes and the augmenting basis $\overline{\mathbf{V}}\in^{N\times m}$ to be the subsequent $m$ modes as in eq. 2.16. Thus $\mathbf{V}$ and $\overline{\mathbf{V}}$ satisfy eq. 2.11. As in eq. 2.12, again let

[TABLE]

i.e., $\widehat{\mathbf{q}}_{j}\in^{r}$ denotes the first $r$ leading POD modal coefficients and $\widehat{\mathbf{p}}_{j}\in^{m}$ the subsequent $m$ modal coefficients. One can again show that eq. 3.5 is equivalent to

[TABLE]

With the data matrices $\widehat{\mathbf{Q}}$ and $\widehat{\mathbf{P}}$ from eq. 2.14, we can apply Theorem 3.1 to eq. 3.7 and obtain the solution

[TABLE]

where ${\boldsymbol{\Omega}}\in^{M\times m}$ solves

[TABLE]

We then obtain the desired nonlinear decoder

[TABLE]

with $\mathbf{n}$ given by eq. 3.8.

An important distinction between our KM approach and the FM approach in Section 2.3 is that the nonlinear function $\mathbf{n}$ is optimal in the RKHS ${\cal H}_{K}^{m}$, whereas the coefficient matrix ${\boldsymbol{\Xi}}\in^{m\times n_{\phi}}$ from Section 2.3 is optimal in the corresponding vector space of matrices. Secondly, the linear system eq. 3.9 scales with the amount of data $M$, whereas the linear system eq. 2.15 scales with the dimension of the feature map ${\boldsymbol{\phi}}$. This can be problematic when using high order polynomial feature maps ${\boldsymbol{\phi}}$, where $n_{\phi}$ can be much larger than the amount of data $M$. Moreover, in the FM approach, the structure of the feature map ${\boldsymbol{\phi}}$ must be explicitly defined, whereas the nonlinear term $\mathbf{n}$ can have implicit structure enforced by different types of kernels. Lastly, our KM allows one to construct a feature map-based reconstruction of the form eq. 2.9 as a special case.

Remark 3.1**.**

The projection-based ROM + artificial neural network (PROM-ANN) approach of [4] also considers decoders of the form eq. 3.10 with $\mathbf{V}$ and $\overline{\mathbf{V}}$ consisting of POD modes, but instead models the nonlinear function $\mathbf{n}$ as an ANN. Furthermore, [4] and [10] suggest selecting $r$ and $m$ by first identifying $n_{\rm tot}:=r+m$ satisfying the energy criterion

[TABLE]

for some $\nu>0$ where $\sigma_{j}$ is the $j$th singular value of $\widetilde{\mathbf{Q}}$, then setting $r$ to some small fraction of $n_{\rm tot}$ and $m$ to $m=n_{\rm tot}-r$.

3.3 Kernel selection and equivalence to feature map approach

By the Moore–Aronszajn Theorem (see, e.g., [39, Theorem 3.10]), a positive-definite kernel $K$ uniquely defines an RKHS ${\cal H}_{K}$. Consequently, the choice of kernel determines what form $\mathbf{n}\in{\cal H}_{K}^{m}$ can take and how accurately it can be used to approximate the data $\mathbf{Q}$. We summarize several effective choices for kernels and show that specific types of kernels, i.e. feature map kernels, can result in nonlinear augmentations of the same form as eq. 2.9.

An effective general-purpose choice of kernel are radial basis function (RBF) kernels. These kernels have the form

[TABLE]

where $\psi:_{\geq 0}\to$ for all $\mathbf{x},\mathbf{x}^{\prime}\in^{n}$ and $\epsilon>0$ is a shape parameter, and hence an RBF interpolant of the data eq. 2.14 has the form

[TABLE]

Table 1 provides some examples of commonly used RBF kernels. A thorough discussion of RBFs in kernel interpolation can be found in, e.g., [48].

Another type of kernel is the polynomial kernel $K:^{d}\times^{d}\to$

[TABLE]

where $c,\rho\in$, $c\geq 0$, $\rho>0$, are hyper-parameters and $\ell\in\mathbb{N}$ is the order of the polynomial. Typical values of $c$ and $\rho$ are $c=1$ and $\rho=1/d$. Using a polynomial kernel results in a polynomial nonlinear term $\mathbf{n}$. Therefore, the resulting decoder eq. 3.10 defines a polynomial manifold. However, using a polynomial kernel only implicitly defines the polynomial structure of eq. 3.10. In some use cases, such as reduced-order modeling, it may be beneficial to know the explicit polynomial structure. In these cases, using a feature map kernel may be preferable.

A feature map ${\boldsymbol{\phi}}:^{d}\to^{n_{\phi}}$ and a strictly positive definite matrix $\mathbf{G}\in^{n_{\phi}\times n_{\phi}}$ induce a postive-definite kernel $K:^{d}\times^{d}\to$ defined as the $\mathbf{G}$-weighted inner product of two feature vectors:

[TABLE]

The resulting interpolant $\mathbf{n}$ can then be expressed as

[TABLE]

and thus the decoder $\mathbf{g}$ is of the form

[TABLE]

Therefore, feature map-based nonlinear decoders of the form eq. 2.9 can be obtained as a special case of the KM approach. In particular, taking ${\boldsymbol{\phi}}$ to be the quadratic feature map eq. 2.18 results in a QM.

Remark 3.2**.**

Although the GP-based closure of [10, Section 3.2] relies on a probabilistic framework, the authors’ use of the GP mean for $\mathbf{n}$ is computationally equivalent to our KM framework with an RBF kernel and regularization $\lambda$ equal to the GP prior covariance $\sigma^{2}$. Similarly, the RBF interpolation-based closure of [10, Section 3.3] is a special case of our KM approach using an RBF kernel.

Remark 3.3** (Input normalization).**

Depending on the input data $\mathbf{X}$ and choice of kernel $K$, the kernel matrix $K(\mathbf{X},\mathbf{X})$ may have poor conditioning, which can adversely affect the solve for the coefficient matrix ${\boldsymbol{\Omega}}.$ Although increasing the regularization parameter $\lambda$ may ameliorate the conditioning of the linear system eq. 3.4, it may also increase the error of the resulting kernel interpolant at the given training data. Alternatively, one can improve the conditioning of $K(\mathbf{X},\mathbf{X})$ by normalizing the inputs $\mathbf{X}$. This relies on the following observation: for an injective ${\boldsymbol{\nu}}:^{d}\to^{d}$, if $K$ is a positive definite kernel, then the function $K_{\boldsymbol{\nu}}:^{d}\times^{d}\to$ defined by

[TABLE]

is also a positive-definite kernel [39, Prop. 2.7]. Therefore, choosing an appropriate ${\boldsymbol{\nu}}$ can improve the conditioning of $K_{\boldsymbol{\nu}}(\mathbf{X},\mathbf{X})$ compared to $K(\mathbf{X},\mathbf{X})$. An effective choice is the function ${\boldsymbol{\nu}}(\mathbf{x})=\mathbf{M}^{-1}(\mathbf{x}-\bar{\mathbf{x}})$, where $\mathbf{M}=\operatorname{diag}(\mathbf{m})\in^{d\times d}$ and $\bar{\mathbf{x}}\in^{d}$ are given by

[TABLE]

which maps the components of each row of $\mathbf{X}$ to the interval $[0,1]$. With this choice of ${\boldsymbol{\nu}}$, we set the weighting matrix $\mathbf{G}$ in a kernel induced by a feature map ${\boldsymbol{\phi}}:^{d}\to^{n_{\phi}}$ to $\mathbf{G}=(1/n_{\phi})\mathbf{I}$.

4 Numerical results

We compare our KM approach with POD and two recent QM approaches: the Greedy QM approach from [43] and the Alternating minimization QM approach from [14]. For the Alternating QM approach, we use the same stopping tolerance of $10^{-3}$ for the alternating minimization procedure as in [14]. For the KM approach, we compare several RBF kernels with varying shape parameter $\epsilon$ and a feature map kernel induced by the quadratic feature map eq. 2.18. Where appropriate, the RBF KM will be denoted “Kernel RBF” and the KM with kernel induced by the quadratic feature map will be denoted “Kernel QM”. Recall from Section 3.3 that a feature map kernel with a quadratic feature map results in a QM of the form eq. 2.9.

4.1 2D Advection-diffusion-reaction – boundary layer

The first problem we consider is given by state snapshot data from the 2D advection-diffusion-reaction equation with homogeneous Dirichlet boundary conditions, adapted from one of the examples in [33]. The 2D advection-diffusion-reaction equation on the unit square $\Omega=(0,1)^{2}$ is given by

[TABLE]

where $\alpha>0$ is the diffusion parameter, ${\boldsymbol{\beta}}=[\beta_{1},\beta_{2}]\in^{2}$, $\beta_{1},\beta_{2}\geq 0$ is the advection parameter, $\gamma\geq 0$ is the reaction parameter, $T>0$ is the final time, and $f:(0,T)\times\Omega\to$ is a forcing function. We discretize the PDE eq. 4.1 using a centered finite difference scheme for the second derivatives and a first-order upwind backward difference scheme for the first derivatives, resulting in a linear system of ODEs of the form

[TABLE]

For this experiment, we take $258$ uniformly spaced grid points in both the $x$ and $y$ directions, resulting in a full state dimension of $N=256^{2}=65,536$. The time stepping is done using scipy.interpolate.solve_ivp() in Python with a BDF scheme [45, 44], time step $\Delta t=0.01$, and final time $T=5$. We take an advection parameter of ${\boldsymbol{\beta}}=(1/2)[\cos(\pi/3),\sin(\pi/3)]$, a reaction parameter of $\gamma=1$, forcing $f=1$, and vary the diffusion parameter $\alpha$. This combination of parameters yields a sharp boundary layer for smaller values of $\alpha$. A similar test case was used in [33, Section 7.2] using a fixed value of $\alpha=10^{-3}$ and a finite element discretization. Figure 1 plots profiles of the solution at $x=0.75$ for different values of $\alpha$ to illustrate how varying $\alpha$ affects the sharpness of the boundary layer.

To collect training data $\mathbf{Q}$, we simulate eq. 4.1 for $\ell=10$ logarithmically spaced values of $\alpha\in[10^{-6},10^{-1}]$. This results in a training dataset

[TABLE]

where $N=65,536$ and $M=\ell(N_{t}+1)=5,010$ and where $\mathbf{q}_{j}^{(k)}$ corresponds to the state of eq. 4.2 at time step $j$ and the $k$th diffusion training diffusion parameter. The offset $\bar{\mathbf{q}}$ is taken to be the mean of the training data eq. 2.6. We then test different dimensionality reduction approaches for $\alpha=10^{-3}$, which is outside of training set. The singular value decay for the first $500$ normalized singular values is plotted in Figure 2. The error metric that we use for this example is the relative $\ell^{2}$ projection error

[TABLE]

where $\mathbf{q}_{j}^{*}$ is test data for $\alpha=10^{-3}$ at time step $j$.

For this example, the Kernel RBF approach uses the quadratic Matérn kernel with shape parameter $\epsilon=10^{-2}.$ We do not use the input normalization discussed in Remark 3.3 for any of the KMs computed. The Kernel QM, Kernel RBF, Alternating QM, and Greedy QM approaches use the regularization values in Table 2. For further details on the selection of the RBF kernel, shape parameter, and regularization values, see Section A.1.

We examine the effect of the number of augmenting modes on the projection error by varying the number of augmenting modes $m\in[10,100]$ for the Kernel QM, Kernel RBF, and Alternating QM approaches for $r\in\left\{10,20,30\right\}$. The Greedy QM approach is excluded from this test because $\overline{\mathbf{V}}=\mathbf{I}$ in the greedy case, and hence does not depend on a choice of number of augmenting modes. Figure 3 indicates that for a fixed value of $r$, the accuracy improvement by increasing the number of augmenting modes $m$ plateaus as $m$ exceeds a small integer multiple of $r$. For $r=10$, the Kernel QM and Alternating QM errors plateau for $m>2r=20$, while the Kernel RBF error decreases for $m\leq 5r=50$, then increases slightly and plateaus for $m>50$. For $r=20$, each method plateaus and achieves similar error for $m>2r=40$, while for $r=30$, each method’s error plateaus for $m\geq 30$. These results further indicate that the benefit of nonlinear augmentation is more substantial for smaller values of $r$. Since Kernel RBF at $r=10$ is the slowest method to plateau for $m\geq 5r$, we fix $m=5r$ for each approach for the remaining tests on the current numerical example.

Next, we compare the different approaches by varying the latent dimension size $r$. As indicated above, we take the number of augmenting modes to be $m=5r$ for the Kernel QM, Kernel RBF, and Alternating QM approaches. Figure 4 shows that the Kernel RBF approach achieves the lowest errors for $r\leq 12$, while each method plateaus with similar error for $r\geq 14$. Notably, Figure 4 also shows that the Kernel QM, Kernel RBF, and Alternating QM approaches have significantly smaller training times than Greedy QM, especially for larger values of $r$. For this numerical example, we can conclude that the Kernel RBF approach achieves the best performance for small values of $r$, while each method performs similarly for larger values of $r$. To visualize the error decrease of the Kernel RBF approach over POD on this numerical example, we plot the full state $\mathbf{q}$ and the relative projection errors $|\mathbf{q}_{j}^{*}-\mathbf{g}(\mathbf{h}(\mathbf{q}_{j}^{*}))|/\sqrt{\sum_{j=0}^{N_{t}}\left\|\mathbf{q}_{j}^{*}\right\|^{2}}$ using POD, Alternating QM, and Kernel RBF with $r=6$ and $m=5r=30$ in Figure 5.

4.2 High-speed aerodynamics – analytic model of surface heating

The next problem considered is an analytic function developed to parametrically emulate prominent features of high-speed aerodynamic flowfields. The physical interpretation for this model problem is the variation in aerothermal quantities of interest over conical geometries at varying flight conditions.

Specifically, this model emulates the variation in heat flux over a flared-cone geometry, similar to the HIFiRE-1 test vehicle [26], though with significant simplification of the underlying physical models and geometry. The test function is defined so that a non-dimensional surface quantity, $s$, corresponding to the heat flux magnitude can be computed on a quasi-three-dimensional domain, defined as a function in cylindrical coordinates $(z,r,\theta)$, where $z\in[0,2]$ is the axial coordinate, $r\in[0,R]$, $R>0$, is a radial coordinate, and $\theta\in[0,2\pi]$ is a circumferential angle measure. We consider variation of the field in response to ${\boldsymbol{\mu}}=[\mu_{1}\;\mu_{2}]^{\mathsf{T}}\in{\cal D}=[-10^{\circ},10^{\circ}]\times[0.4,0.8]$ consisting of two parameters: the angle of attack, $\mu_{1}\in[-10^{\circ},10^{\circ}]$, dictating the pitch angle orientation relative to the freestream, and a laminar-to-turbulent transition location parameter, $\mu_{2}\in[0.4,0.8]$, controlling the location of turbulent transition (and associated increase in heat flux) on the cone surface, where the dimensionless field quantity $s:[0,2]\times[0,2\pi]\times{\cal D}\to$ can be expressed as

[TABLE]

Here, $\sigma(x)$ is a scaled sigmoid function modified to approximate a sharp discontinuity as

[TABLE]

and $\delta$ is a small, nonzero value of $10^{-6}$.

The first term in eq. 4.5 constitutes a rapid quadratic decay of $s$ over the surface modulated by the angle of attack parameter. The second and third terms in the equation represent the onset and end of a laminar-to-turbulent transition front on the conical portion of the vehicle. An assumption is made that the transition is pulled forward on the windward surface of the geometry, and in this parameter field the transition front is constrained to the conical region. The fourth term represents heating augmentation associated with shock-boundary layer at the collar compression corner, resulting in a large jump in $s$ modulated by the angle of attack with a reciprocal decay. It is important to note that this model is designed to qualitatively capture features of high-speed heating profiles, but is not based on first-principles physics. We plot the field quantity $s$ on a flared-cone geometry for several different parameters ${\boldsymbol{\mu}}$ in Figure 6 to demonstrate the effect of the parameters ${\boldsymbol{\mu}}$ on $s$.

For a given parameter ${\boldsymbol{\mu}}$, we compute $s$ on a uniform grid for 100 equispaced $z$ values in $[0,2]$ and 100 equispaced $\theta$ values in $[0,2\pi]$, resulting in snapshots $\mathbf{q}({\boldsymbol{\mu}})\in^{N}$, $N=10,000$. To collect training and testing data for this problem, we compute $\mathbf{q}({\boldsymbol{\mu}})$ over a uniform grid of 50 equispaced $\mu_{1}$ values in $[-10^{\circ},10^{\circ}]$ and 50 equispaced $\mu_{2}$ values in $[0.4,0.8]$, resulting in $M=2500$ total snapshots, ordered as

[TABLE]

We then scale $\mathbf{Q}^{0}$ as

[TABLE]

and separate every other column of $\mathbf{Q}$ into training and testing sets, $\mathbf{Q}_{\rm train}$ and $\mathbf{Q}_{\rm test}$, respectively. The offset $\bar{\mathbf{q}}$ is taken to be the mean of the training data $\mathbf{Q}_{\rm train}$, and the error metric used is the mean relative $\ell_{2}$ projection error

[TABLE]

where $M_{\rm test}=1250$. Figure 7 plots the decay of the normalized singular values of the training data. Notice that the singular values for this example decay more slowly than for the boundary layer example, indicating that the current example is a more demanding test case for our KM approach.

For this example, the Kernel RBF approach uses the Gaussian kernel with shape parameter $\epsilon=10^{-1}.$ The Kernel QM and Kernel RBF approaches both use the input normalization discussed in Remark 3.3, and the Kernel QM, Kernel RBF, Alternating QM, and Greedy QM approaches use the regularization values in Table 3. For further details on the selection of the RBF kernel, shape parameter, and regularization values for this example, see Section A.2.

First, we vary the number of augmenting modes $m\in[10,150]$ for the Kernel QM, Kernel RBF, and Alternating QM approaches for $r\in\left\{10,20,30\right\}$ to examine the effect of the number of augmenting modes on the mean relative projection error. As in the advection-diffusion-reaction case, the Greedy QM approach is excluded from this test because it does not include an augmenting basis. Like the previous example, Figure 8 shows that for a fixed $r$, the error plateaus as $m$ exceeds a small integer multiple of $r$ for each method. Interestingly, the errors begin to plateau for larger values of $m$ for this problem than for the advection-diffusion-reaction problem. For $r=10$, the Alternating QM and Kernel QM errors plateau for $m\geq 40$ augmenting modes, while for $r=20$, the errors plateau for $m\geq 100$ modes, and for $r=30$, they plateau for $m\geq 130$ modes. On the other hand, the Kernel RBF approach plateaus for $m\geq 100$ for $r=10$, but does not plateau in the tested range of $m$ for $r=20$ and $r=30$. Thus, in contrast with the advection-diffusion-reaction problem, the additional accuracy benefit of using an RBF kernel over a QM is not only restricted to smaller values of $r$. For the remaining comparisons on this example, we fix $m=10r$ for each of the Kernel RBF, Kernel QM, and Alternating QM approaches.

We now compare the different approaches by computing the error as a function of the latent space dimension $r$. From Figure 9, we see that each nonlinear augmentation approach outperforms POD. The Kernel QM and Alternating QM have nearly identical errors for each $r$. The Kernel RBF approach attains the lowest errors for $r\leq 30$, while Greedy QM has roughly the same error as Kernel RBF for $r=35$ and slightly lower error for $r=40$. Greedy QM yields similar errors to the other QM approaches for $r\leq 15$, but decreases more rapidly for $r>15$. The KM approaches have smaller training times than both the Alternating QM and Greedy QM approaches, and the training time for the Greedy QM approach increases substantially with increasing $r$. Since the Kernel RBF approach has the smallest errors for $r\leq 30$ and has nearly the same error as Greedy QM for $r=35,40$ with substantially lower training times, we can conclude that Kernel RBF has the best performance for this example. Furthermore, as noted above, unlike the advection-diffusion-reaction example, the Kernel RBF has accuracy benefits over QMs for each value of $r$ tested rather than for only small values of $r$.

We plot the full state $\mathbf{q}({\boldsymbol{\mu}})$ for ${\boldsymbol{\mu}}=[-10^{\circ},\,0.4082]^{\mathsf{T}}$ and the relative projection errors $|\mathbf{q}({\boldsymbol{\mu}})-\mathbf{g}(\mathbf{h}(\mathbf{q}({\boldsymbol{\mu}})))|/\left\|\mathbf{q}({\boldsymbol{\mu}})\right\|_{2}$ for POD, Alternating QM, and Kernel RBF in Figure 10 for $r=10$ and $m=100$. Notice that the POD projection errors are largest towards the tip of the conical geometry. These errors are decreased in the Alternating QM projection error, and decreased further in the Kernel RBF projection error.

4.3 Solid mechanics – 3D flexible bracket

The third example is a moderately sized 3D hyper-elastic solid dynamics problem, referred to herein as the flexible bracket. The governing PDEs are the following equations of motion, written in terms of displacement vector $\mathbf{d}\in\mathbb{R}^{3}$, material density $\rho>0$, and the first Piola-Kirchhoff stress $\mathbf{P}\in\mathbb{R}^{3\times 3}$:

[TABLE]

In (4.8), $\ddot{\mathbf{d}}:=\frac{\partial^{2}\mathbf{d}}{\partial t^{2}}$ represents the acceleration. Encoded into $\mathbf{P}$ is a constitutive model for the underlying material. Herein, we specify a Neohookean-type material model extended to the compressible regime, with Young’s modulus $E=200\times 10^{9}$ GPa, Poisson’s ratio $\nu=0.25$ and density $\rho=7800$ kg/m3. These properties correspond to a real material, namely steel. For more details about the Neohookean-type constitutive model used to generate the results herein, the reader is referred to [31]. We emphasize that the governing PDE (4.8) is nonlinear and contains generic (non-polynomial) nonlinearities, as the Helmholtz free-energy density corresponding to the Neohookean constitutive model is nonlinear.

The flexible bracket geometry, mesh and initial conditions for our test case are shown in Figure 11. In the $xy-$plane is a square plate that is 0.205 m wide and features a circular whole with diameter 0.075 m. This plate is connected to a vertical plate of the same dimensions in the $yz-$plane. On the left side, at $x=-0.105$ m, the bracket is clamped, meaning we set a boundary condition of $\mathbf{d}=\mathbf{0}$. The problem is initialized by setting a constant initial condition \dot{\mathbf{d}}(0;\mu):=\frac{\partial\mathbf{d}}{\partial t}(0;\mu)=\left(\begin{array}[]{ccc}0,&0,&\mu\end{array}\right)^{\intercal}\in\mathbb{R}^{3} m/s for a specified value of the $z$–velocity $\mu\in\mathbb{R}$ on the right-most side of the bracket, corresponding to $x=0.1$ m (see Figure 11(b)). As a result of the velocity initial condition, the flexible bracket begins to vibrate and flap. The presence of the circular hole weakens the bracket as well as complicates the propagation and interaction of 3D waves moving throughout the bracket, leading to truly nonlinear 3D dynamics. A linear version of this problem was considered in [18].

To generate training data for our KM methods, we first discretize and solve eq. 4.8 using a Julia finite element code known as Norma.jl [1]. The flexible bracket geometry is discretized using a unstructured four-node tetrahedral mesh consisting of 1174 elements having a total of 481 nodes, as shown in Figure 11(a). The resulting problem has 1443 degrees of freedom. The semi-discrete version of eq. 4.8 is then advanced forward in time from time 0 to $20.0\times 10^{-3}$ seconds using an implicit Newmark time-integration scheme with parameters $\beta=0.25$ and $\gamma=0.5$. For the cases with initial condition $\mu\leq 150$ m/s, a time-step of $\Delta t=1.0\times 10^{-5}$ seconds is utilized. For the harder $\mu=200$ m/s initial condition case, we employ a smaller time-step of $\Delta t=1.0\times 10^{-6}$ seconds. Each of the solutions are then uniformly downsampled to retain $1001$ snapshots per sample. To create the training data for our methods, we simulate the flexible bracket for four different initial conditions on the $z$–velocity $\mu$: 10 m/s, 50 m/s, 150 m/s and 200 m/s. We then test the different dimensionality reduction approaches discussed herein on the displacement, velocity, and acceleration fields for the $\mu=100$ m/s case, which is outside our training set. Figure 12 shows the norm of the displacement vector, $||\mathbf{d}||_{2}$, amplified by a factor of five, at the final simulation time for various initial conditions $\mu$. The reader can observe that the set of $\mu$ considered in our study gives rise to nontrivial differences in the solution, making this an excellent test case for evaluating the predictive power of the various KM approaches. Figure 13 shows the singular value decay for the first $1500$ normalized singular values for the flexible bracket test case. The fact that there is an incredibly slow singular value decay up to index 1000 makes the flexible bracket problem a good candidate for the KM approaches described and proposed herein.

Results comparing the performance of the various dimension reduction approaches are summarized in Figure 14–Figure 16. We employ input normalization as described in Remark 3.3 for each KM approach tested. The error metric that we use for this example is the relative $\ell^{2}$ projection error eq. 4.4, where $\mathbf{q}_{j}^{*}$ is test data for the displacement solution corresponding to a given value of the parameter $\mu$ (the $z$–velocity initial condition) at time step $j$.

For this example, we compute separate nonlinear-augmentation manifolds for the displacement, velocity, and acceleration fields. For each field, the Kernel RBF approach uses the linear Matérn kernel with shape parameter $\epsilon=10^{-1}.$ Each Kernel QM and Kernel RBF computed uses the input normalization as discussed in Remark 3.3. Table 3 lists the regularization values used for Kernel QM, Kernel RBF, Alternating QM, and Greedy QM for the displacement, velocity, and acceleration fields. Note that the regularization values needed for the velocity and acceleration Alternating QMs are quite large. This is likely due to the complexity of the underlying structures in the data, which in turn increases the values of $r$ and $m$ necessary to obtain low errors, thus resulting in increasingly poor conditioning of the system eq. 2.15 as $r$ and $m$ increase. This poor conditioning is mitigated in the KM approaches because the size of eq. 3.4 is independent of $r$ and because normalizing the inputs also improves the conditioning of eq. 3.4. For further details on the selection of the RBF kernel, shape parameter, and regularization values for this example, see Section A.3.

We examine the effect of the number of augmenting modes $m$ on the displacement, velocity, and acceleration projection errors for the Kernel QM, Kernel RBF, and Alternating QM approaches for $r=20,40,60$. Unlike the previous examples, Figure 14 shows that Kernel QM and Alternating QM show little improvement in error as $m$ increases for $r=20,40$ across each of the displacement, velocity, and acceleration fields. For $r=60$, the Kernel QM error improves compared to Alternating QM and decreases as $m$ increases, whereas $r=60$ Alternating QM only decreases noticeably as $m$ increases for the velocity field. In contrast, the Kernel RBF substantially improves in error as $m$ increases for each of the fields and for each value of $r$ tested. Only for the $r=10$ displacement case does Kernel RBF plateau; in the other cases, the Kernel RBF error decreases monotonically. These results indicate that the quadratic nonlinear augmentations of the QM approaches poorly model the corrections needed to improve the accuracy for this problem, while RBF-based augmentations are a better model for the nonlinear correction terms.

Since Kernel RBF has monotonically decreasing error as $m$ increases for this example, we set $m=5r$ for each of the Kernel QM, Kernel RBF, and Alternating QM methods in the following comparisons. Lastly, we note that, although the trends among the different nonlinear augmentation methods are similar for displacement, velocity, and acceleration, the errors obtained for velocity and acceleration are significantly larger than for displacement. This result is expected for solid dynamic problems: generally, the velocity error is one order of magnitude higher than the displacement error, and the acceleration error is one order of magnitude higher than the velocity error.

Next, in Figure 15, we compare the relative errors and train times of the different approaches as a function of latent dimension size for the displacement, velocity, and acceleration fields. As noted above, the velocity and acceleration fields yield much larger errors than displacement, thus requiring much larger values of $r$ to yield errors similar to displacement. Consequently, we omit Greedy QM from the velocity and acceleration comparisons due to the excessive training times required for larger $r$ values. Figure 15 shows that for the displacement field, each of the QMs yield nearly identical errors for $r\leq 30$, while Kernel QM yields slightly smaller errors than Alternating QM and Greedy QM for $r\geq 35$. Each of the QMs only yield slight error improvements over POD. Kernel RBF yields slighly larger error than POD at $r=5$, but rapidly decreases in error as $r$ increases, with nearly an order of magnitude improvement in error over the other approaches. For the velocity field, Kernel QM and Alternating QM have nearly identical errors, and yield nearly an order of magnitude improvement over POD for $r>50$. Kernel QM performs similarly for the acceleration field, while Alternating QM sharply increases in error for $r>50$. We note that since we use the fixed regularization parameters in Table 4 for each value of $r$, the Alternating QM error for $r>50$ may be improved by fine-tuning the regularization for different values of $r$. Kernel RBF still yields the smallest errors among the different methods for the velocity and acceleration fields. Lastly, we note that Greedy QM still requires the largest train times for the displacement field, while Alternating QM requires the largest train times for the velocity and acceleration fields among the compared methods.

Finally, in Figure 16, we plot the full solution states and the relative projection errors for the displacement, velocity, and acceleration fields at the final time $t=20.0\times 10^{-3}$ seconds. Projection errors are shown for the POD, Alternating QM and Kernel RBF KMs. We utilize a value of $r=25$ and $m=125$ when predicting the displacement, a value of $r=50$ and $r=250$ when predicting the velocity and a value of $r=100$ and $m=500$ when predicting the acceleration. The reader can observe that the results in Figure 16 are consistent with the previously-presented results in Figure 14 and Figure 15. The Kernel RBF KM is the most accurate by at least an order of magnitude. Generally, POD is the least accurate, with the exception of the acceleration field. While it may initially seem surprising that the Alternating QM is the least accurate when it comes to the acceleration field, this result is consistent with the spike in the acceleration projection error that occurs around $r=100$, as shown in Figure 15. The noisiness in the velocity and acceleration solutions in Figure 16 is expected, and can be attributed to the complex nonlinear interaction of propagating waves within the bracket geometry.

4.4 2D Euler – double Mach reflection

Our last numerical example is the double Mach reflection, a commonly used verification and validation test case for a compressible Euler equation solver originally proposed by Woodward and Colella [47]. The double Mach reflection case models a shocktube experiment in which a planar shock is driven down a tube containing a wedge, as shown in Figure 17.

When the shock interacts with the wedge, a complicated shock interaction occurs, as seen in Figure 18. This shock structure is a challenge to the shock capturing schemes utilized by Euler solvers. As in the original case proposed by Woodward and Colella, we consider a Mach 10 shock, a specific heat ratio $\gamma=1.4$ and an upstream non-dimensionalized pressure of $p=1.0$.

To parametrize the problem, we consider the pressure fields $p(x,y)$ computed for different wedge half-angles $\alpha$. This results in a range of different shock structures that are difficult to capture accurately with a low-dimensional linear subspace.

The double Mach reflection is modeled by the 2D conservative Euler equations

[TABLE]

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

[TABLE]

where $\gamma=1.4$ is the specific heat ratio. The spatial domain is the unit square $\Omega=(0,4)\times(0,1)$ and the time domain is $(0,0.25)$.

We collect snapshots using the open-source Python library pressio-demoapps111pressio.github.io/pressio-demoapps to simulate eq. 4.9, which uses a cell-centered finite volume scheme. For this example, we use a $600\times 150$ uniform Cartesian mesh, resulting in a discretization with state dimension $n_{q}=600\times 150\times 4=360,000$, and a WENO5 scheme for inviscid flux reconstruction. The time stepping is done using pressio-demoapps’ SSP3 scheme for times $t\in(0,0.25)$ with time step $\Delta t=2.5\times 10^{-4}$. We plot decay of the normalized singular values in Figure 20. It can be observed that the singular value decay for this example is exceptionally slow, indicating that this is a particularly difficult test case. The error metric that we use for this example is the relative $\ell^{1}$ projection error

[TABLE]

where $\mathbf{q}_{j}^{*}$ is test data at time step $j$.

Figure 19 plots the initial pressure for varying wedge angles. Notice that each initial pressure contains a sharp front. The training set is comprised of $\alpha=5^{\circ},15^{\circ},25^{\circ},35^{\circ},45^{\circ}$, and all testing is done with $\alpha=30^{\circ}$.

For this example, the Kernel RBF approach uses the inverse quadratic RBF kernel with shape parameter $\epsilon=10^{-3}.$ The Kernel QM and Kernel RBF approaches do use the input normalization discussed in Remark 3.3, and the Kernel QM, Kernel RBF, Alternating QM, and Greedy QM approaches use the regularization values in Table 5. Note that similar to the 3D bracket example, the optimal regularization values for this problem are especially large. This is likely due to the sharp shocks that appear in the training data, which require large POD basis sizes to resolve, thus worsening the conditioning of the systems eq. 2.15 and eq. 3.4 for training the Alternating QM and KMs, respectively. For further details on the selection of the RBF kernel, shape parameter, and regularization values for this example, see Section A.4.

Figure 21 plots the relative projection error as a function of the number of augmenting modes for $r=10,20,30$ for Kernel QM, Kernel RBF, and Alternating QM. For $r=10$, Alternating QM outperforms both Kernel QM and Kernel RBF, yielding smaller relative errors and plateaus in error for $r\geq 30$. Kernel RBF yields smaller errors than Kernel QM and plateaus in error for $r\geq 40$, while Kernel QM with $r=10$ plateaus in error for $r\geq 40$. For $r=20$, Alternating QM again outperforms both Kernel QM and Kernel RBF and plateaus in error for $r\geq 60$. In contrast with the $r=10$ case, Kernel QM outperforms Kernel RBF and plateaus in error for $r\geq 60$. Interestingly, $r=20$ Kernel RBF has nearly identical performance as $r=10$. For $r=30$, Alternating QM and Kernel QM have nearly identical performance as $m$ increases, and plateau in error for $r\geq 50$, while Kernel RBF has the largest errors for $r=30$. However, note that all of the relative errors computed for this comparison are in the range $[0.07,0.09]$, so each dimensionality reduction approach has overall similar performance. For the following comparison, we fix $m=4r$ for Kernel QM, Kernel RBF, and Alternating QM.

Figure 22 compares the different approaches by plotting the relative projection error as a function of the latent dimension size $r$ for $r\in[5,40].$ Notice that Kernel QM and Alternating QM have nearly identical errors for $r\leq 20$, and Kernel QM slightly increases in error for $r\geq 25$ while Alternating QM continues to decrease in error. Kernel RBF has slightly smaller errors than Kernel QM and Alternating QM for $r=5$, but has larger errors than Kernel QM and Alternating QM for $r\geq 10$. Greedy QM has errors nearly identical to POD for each value of $r$, with $r=10$ yielding slightly higher error than POD. This is because Greedy QM requires an especially large regularization value for this particular problem, thus effectively negating the benefits of the additional quadratic term; see Figure 31. Unlike the previous examples, Alternating QM overall outperforms the KM and Greedy QM approaches for this particular example. Thus we cannot conclude that the RBF KM approach will always outperform other nonlinear-augmentation manifolds. However, as noted above, the difference in performance among Kernel QM, Kernel RBF, and Alternating QM is slight; the error obtained for each method and for each value of $r$ lies in the range $[0.06,0.09]$. Furthermore, the relative error for POD with $r\geq 20$ also lies in this range, indicating that the nonlinear-augmentation approaches provide limited accuracy gains over POD for this particular example.

Lastly, Figure 23 plots the full state pressure for $\alpha=30^{\circ}$ and times $t=0,0.125,0.25$ and the corresponding projection errors for POD, Kernel RBF, and Alternating QM for $r=10$ and $m=40$. For each approach, the error is largest near the shock fronts. We also see that Kernel RBF and Alternating QM have smaller errors than POD in regions further away from the shock fronts.

5 Discussion and conclusion

This paper generalizes recent work on QMs and the GP/RBF-based nonlinear closure approach of [10] by developing nonlinear-augmentation dimensionality reduction relying on regularized kernel interpolation. Like previous approaches, our method augments POD with a nonlinear correction term in the decoder map, which we model as an optimal nonlinear function lying in a given RKHS. We leverage the Representer Theorem to fit a nonlinear mapping from low-order POD coefficients to high-order POD coefficients, which enables one to fit nonlinear augmentation terms with arbitrary nonlinearities through the choice of kernel. We show that feature-map based nonlinear augmentations, specifically QMs, as well as the GP/RBF approach, can be obtained as special cases of our method.

Our numerical examples compare the KM approach to POD and several recent state-of-the-art QM approaches on data from several example problems and show that KMs using RBF kernels can yield substantial accuracy gains in the small latent dimension size regime. These results suggest that mappings from low-order POD modal coefficients to an optimal nonlinear correction often cannot be sufficiently approximated using quadratic nonlinearities, making RBF kernels an attractive alternative. Additionally, the results on the double Mach reflection test case show that nonlinear-augmentation manifolds may only provide slight improvements over POD for especially difficult test cases. This is because nonlinear-augmentation approaches still rely on an underlying POD, and therefore will at best match the projection error of linear dimensionality reduction with a large POD basis. Thus if a large POD basis is still insufficient for approximating the given problem, a nonlinear-augmenation manifold will be as well. Therefore, for sufficiently difficult problems, an alternative nonlinear dimensionality reduction technique, such as autencoders, may be preferable.

Future work will focus on improving the attainable accuracy in the KM framework. In particular, we will explore other kernels, such as high-order polynomial kernels or sums and products of kernels, and examine their effect on the KM’s accuracy. We will also investigate if a natural nonlinear form for the map $\mathbf{n}$ mapping low-order POD coefficients to high-order POD coefficients can be deduced, thus informing the selection of an “optimal” kernel for use in the KM approach. Lastly, we will incorporate KMs as the dimensionality reduction component in a ROM, such as the kernel-based ROMs of [12].

Acknowledgements

Support for this work was received through Sandia National Laboratories’ Laboratory Directed Research and Development (LDRD) program and through the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Mathematical Multifaceted Integrated Capability Centers (MMICCs) program, under Field Work Proposal 22025291 and the Multifaceted Mathematics for Predictive Digital Twins (M2dt) project. A.N. Diaz acknowledges support from Sandia National Laboratories’ S. Scott Collis Fellowship in Data Science. Additionally, the writing of this manuscript was funded in part by I.K. Tezaur’s Presidential Early Career Award for Scientists and Engineers (PECASE). The authors wish to thank Alejandro Mota for creating the Norma.jl finite element code and for helping to set up the 3D flexible bracket problem within this code.

Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC (NTESS), a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration (DOE/NNSA) under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. SAND2025-10914O

Appendix A RBF and regularization selection

Each of the nonlinear-correction manifolds that we compare is highly dependent on the regularization parameter $\lambda$. Additionally, the KM approach using an RBF kernel also depends on the choice of RBF kernel and shape parameter $\epsilon$. While in practice, one can use a hyper-parameter optimization package to tune these choices, as is done in [10], we instead investigate the dependence of the compared approaches on the regularization, RBF kernel, and shape parameters for each of the numerical experiments in Section 4.

A.1 2D Advection-diffusion-reaction – boundary layer

We first examine the effect of using different RBF kernels with varying shape parameters $\epsilon$ and regularization $\lambda$ on the testing error achieved by the KM approach. For this test, we fix the latent size to $r=10$ and the number of augmenting modes to $m=20$. Figure 24 shows that increasing the shape parameter shifts both the optimal regularization value and the lowest relative projection error. Among the combinations tested, the quadratic Matérn kernel with shape $\epsilon=10^{-2}$ achieved the lowest error of $E=5.07\times 10^{-4}$ at regularization $\lambda=10^{-10}$. Therefore, we use the quadratic Matérn kernel with shape parameter $\epsilon=10^{-2}$ in the following comparisons with the varying QM approaches.

Next we examine the effect of the regularization parameter $\lambda$ on the different KM and QM approaches for different latent space dimensions. See Figure 25. For the Kernel RBF, Kernel QM, and Alternating QM approaches, the number of augmenting modes $m$ was set to $m=2r$. From Figure 25, we see that for the three latent sizes examined, the Kernel QM and Alternating QM approaches benefit from larger regularization values, whereas the Kernel RBF achieves smaller errors for smaller regularization values. Interestingly, the Greedy QM approach has a different optimal regularization range for each value of $r$. In the following tests where we vary the latent dimension, we fix the regularization for each dimension reduction approach to the values listed in Table 2. However, for optimal results, tuning the regularization for each latent space dimension would be necessary.

A.2 High-speed aerodynamics – analytic model of surface heating

As in the previous numerical example, we first compare KMs with different RBF kernels and varying shape parameters $\epsilon$ and varying regularization $\lambda$. For this example, each KM utilizes the input normalization discussed in Remark 3.3. We fix the latent size to $r=10$ and the number of augmenting modes to $m=100$. Figure 26 shows that over the regularization range tested each KM attains a similar minimum error for each RBF and shape parameter combination. For each shape value, as the regularization increases, the error also increases from a minimum value to the POD projection error. The regularization value for which this increase begins increases as the shape parameter increases. Among the combinations tested, the Gaussian kernel with shape $\epsilon=10^{-1}$ achieved the lowest error of $E=3.75\times 10^{-2}$ at regularization $\lambda=10^{-12}$. Therefore, we use the Gaussian kernel with shape parameter $\epsilon=10^{-1}$ in the following comparisons with the varying QM approaches.

Next we examine the effect of the regularization parameter $\lambda$ on the different KM and QM approaches for different latent space dimensions. See Figure 27. For both KM approaches and the Alternating QM approach, the number of augmenting modes $m$ was set to $m=10r$. Figure 27 shows that for $r=10$ and $r=20$, the Kernel QM, Alternating QM, and Greedy QM are relatively insensitive to changes in regularization for smaller regularization values, while for $r=30$, there is a clear optimal regularization value in the tested range for each method. The Kernel RBF approach increases in error as regularization increases for $r=10$ and $r=20$, but has an optimal regularization value at $\lambda=10^{-10}$ for $r=30$. Based on the errors shown in Figure 27, we fix the regularization for each dimension reduction approach to the values listed in Table 3 for the remaining experiments.

A.3 Solid mechanics - 3D flexible bracket

As in the previous examples, we first examine the effect of different RBFs and shape parameters on the KM projection error for varying regularization values $\lambda$. We conduct this comparison for each of the displacement, velocity, and acceleration fields. For displacement, we use $r=10$ and $m=50$, while for velocity and acceleration we use $r=30$ and $m=150$. Figure 28 shows across each field and shape parameter, the basic and linear Matérn kernels are the least sensitive to the regularization value $\lambda$, and the linear Matérn kernel yields the lowest projection error for each field. We also observe that increasing the shape parameter $\epsilon$ also increases the optimal regularization value for the Gaussian, quadratic Matérn, inverse quadratic, and inverse multiquadric kernels. Due to the insensitivity to regularization and yielding the lowest error, we use the linear Matérn kernel with shape $\epsilon=10^{-1}$ for the remaining comparisons on the flexible bracket example.

Next we examine the effect of the regularization value $\lambda$ on the Kernel QM, Kernel RBF, Alternating QM, and Greedy QM approaches for the displacement, velocity, and acceleration fields. Since the velocity and acceleration fields require a significantly larger value of $r$ to obtain similar errors to displacement for smaller values of $r$, Greedy QM is omitted from the comparison for the velocity and acceleration fields. Figure 29 shows that for displacement, each of the QM approaches are relatively insensitive to changes in the regularization and yield errors nearly identical to POD. For the velocity field, Kernel QM and Alternating QM yield lower errors than POD for $r=40$ and $r=60$. Interestingly, for velocity at $r=60$, a large regularization value of $\lambda=10^{6}$ for Alternating QM yields the optimal error. For the acceleration field, Alternating QM yields identical errors for each regularization value tested. Kernel QM yields identical acceleration errors as Alternating QM for $r=20$ for each regularization value, while for $r=40$, the two QMs have the same error for smaller regularization values and the Kernel QM error increases to match the POD error for larger regularization values. For $r=60$, Kernel QM yields optimal acceleration errors at $\lambda=10^{-7}$, slightly outperforming Alternating QM. The Kernel RBF approach yields the lowest projection errors for each field, regularization value, and latent space dimension.

A.4 2D Euler – double Mach reflection

Figure 30 plots relative projection error as a function of regularization $\lambda$ for the Kernel RBF approach using several different RBF kernels and shape parameters $\epsilon$. For this comparison, we fix $r=10$ and $m=20$. Notice that using a shape parameter $\epsilon=10^{-3}$ yields the lowest errors. Increasing the shape parameter to $\epsilon=10^{-2}$ increases the optimal regularization values for each RBF kernel except the basic Matérn kernel, which is relatively insensitive to changes in regularization. Taking the shape parameter $\epsilon=10^{-1}$ degrades the accuracy of each of the KMs, but also reduces the sensitivity of the error with respect to regularization. For the remaining comparisons on the double Mach reflection example, we use the inverse quadratic kernel with shape parameter $10^{-3}$.

Figure 31 plots relative projection error as a function of regularization $\lambda$ for the Kernel QM, Kernel RBF, Alternating QM, and Greedy QM approaches for several different values of $r$. For this comparison, we take $m=2r$ for the Kernel QM, Kernel RBF, and Alternating QM approaches. Observe that as $r$ increases, the best error achieved by each method gets closer to the POD error. Furthermore, as $r$ increases, the QM approaches require significantly larger regularization values; otherwise, the error greatly exceeds the POD error. The Kernel RBF approach is the least sensitive to regularization, but also yields larger errors compared to the Kernel QM and Alternating QM methods at their optimal regularization values.