A registration method for model order reduction: data compression and geometry reduction

TL;DR

This paper introduces a general registration method for parameterized model order reduction that aligns solution manifolds to improve linear compression, applicable across various equations and geometries.

Contribution

It proposes a novel, equation-independent registration technique that enhances model order reduction by optimizing data alignment for better compression.

Findings

Effective in reducing dimensionality of solution manifolds

Improves linear compression methods for parameterized problems

Demonstrated success on two-dimensional numerical examples

Abstract

We propose a general --- i.e., independent of the underlying equation --- registration method for parameterized Model Order Reduction. Given the spatial domain $\Omega \subset \mathbb{R}^d$ and a set of snapshots $\{ u^k \}_{k=1}^{n_{\rm train}}$ over $\Omega$ associated with $n_{\rm train}$ values of the model parameters $\mu^1,\ldots, \mu^{n_{\rm train}} \in \mathcal{P}$, the algorithm returns a parameter-dependent bijective mapping $\boldsymbol{\Phi}: \Omega \times \mathcal{P} \to \mathbb{R}^d$: the mapping is designed to make the mapped manifold $\{ u_{\mu} \circ \boldsymbol{\Phi}_{\mu}: \, \mu \in \mathcal{P} \}$ more suited for linear compression methods. We apply the registration procedure, in combination with a linear compression method, to devise low-dimensional representations of solution manifolds with slowly-decaying Kolmogorov $N$-widths; we also consider the application to problems in parameterized geometries. We present a theoretical result to show the mathematical rigor of the registration procedure. We further present numerical results for several two-dimensional problems, to empirically demonstrate the effectivity of our proposal.

Equations

$$\mathrm{find}{\mathbf{y}}_{\mu }:={\mathfrak{s}}_{\mu }(z_{\mu })s.t.{\mathcal{G}}_{\mu }(z_{\mu },v)=0\forall v\in \mathcal{X},$$

$$\mathrm{find}{\mathbf{y}}_{\mu }:={\mathfrak{s}}_{\mu }(z_{\mu })s.t.{\mathcal{G}}_{\mu }(z_{\mu },v)=0\forall v\in {\mathcal{X}}_{\mu },$$

$$\mathrm{find}{\mathbf{y}}_{\mu }={\mathfrak{s}}_{\mu ,\Phi }({\tilde{z}}_{\mu })s.t.{\mathcal{G}}_{\mu ,\Phi }({\tilde{z}}_{\mu },v)=0\forall v\in \mathcal{X},$$

$$d_N(\mathcal{M}):=\underset{{\mathcal{Z}}_N\subset \mathcal{X},\dim ({\mathcal{Z}}_N)=N}{\inf }\underset{w\in \mathcal{M}}{\sup }\Vert w-{\Pi }_{{\mathcal{Z}}_N}w\Vert ,$$

$$\mathbf{C}{\mathbf{\zeta }}_n={\lambda }_n{\mathbf{\zeta }}_n,{\zeta }_n:=\sum _{k=1}^{n_{\mathrm{train}}}{\left({\mathbf{\zeta }}_n\right)}_k{u}^k,n=1,\dots ,n_{\mathrm{train}},$$

$$N:=\min \left\{N^{\prime }:\sum _{n=1}^{N^{\prime }}{\lambda }_n\ge \left(1-to{l}_{\mathrm{pod}}\right)\sum _{i=1}^{n_{\mathrm{train}}}{\lambda }_i\right\},$$

$${\mathrm{dist}}_{\mathrm{H}}(\mathbf{\Phi }(U),V)\le \text{tol},\mathbf{\Phi }:U\to \mathbf{\Phi }(U)\mathrm{is}\mathrm{bijective}.$$

$${\mathbf{\Psi }}_{\mathbf{a}}^{\mathrm{hf}}(\mathbf{X}):=\mathbf{X}+\sum _{m=1}^{M_{\mathrm{hf}}}{a}_m{\mathbf{\phi }}_m^{\mathrm{hf}}(\mathbf{X}),$$

$$U=\mathbf{{\rm Y}}(\hat{\Omega }),V=\mathbf{\Lambda }(\hat{\Omega }),\hat{\Omega }=\{\mathbf{x}\in {\mathbb{R}}^d:f(\mathbf{x})<0\},\mathrm{where}f:{\mathbb{R}}^d\to \mathbb{R}\mathrm{is}\mathrm{convex},$$

$$\mathbf{\Phi }(\Omega )=\mathbf{\Phi }\left({\Omega }^{\mathrm{out}}\setminus \bigcup _{i=1}^Q{\Omega }_i^{\mathrm{in}}\right)=\mathbf{\Phi }\left({\Omega }^{\mathrm{out}}\right)\setminus \bigcup _{i=1}^Q\mathbf{\Phi }\left({\Omega }_i^{\mathrm{in}}\right)={\Omega }^{\mathrm{out}}\setminus \bigcup _{i=1}^Q{\Omega }_i^{\mathrm{in}}=\Omega ,$$

$$\boldsymbol{\varphi}_{m}^{\rm hf}(\mathbf{X})\cdot\mathbf{e}_{i}\,=\,0,\quad{\rm on}\;\{\mathbf{X}:\;X_{i}=0,{\rm or}\,X_{i}=1\},\quad\left\{\begin{array}[]{l}m=1,\ldots,M_{\rm hf};\\ i=1,\ldots,d.\\ \end{array}\right.$$

$$\underset{\mathbf{X}\in \overline{\Omega }}{\min }\det \left(\hat{\nabla }\mathbf{\Phi }(\mathbf{X})\right)>0.$$

$$\boldsymbol{\Phi}(\mathbf{X}=(t,1))=\left[\begin{array}[]{c}\varphi(t)\\ 1\\ \end{array}\right],\quad{\rm where}\;\varphi(t):=\mathbf{e}_{1}\cdot\boldsymbol{\Phi}(\mathbf{X}=(t,1)).$$

$$\mathbf{e}_{1}^{T}\,\widehat{\nabla}\boldsymbol{\Phi}\big{|}_{\mathbf{X}=(\bar{t},1)}=\left[\varphi^{\prime}(\bar{t}),\;\;0\right]=0,$$

$$\left\{\begin{array}[]{l}\boldsymbol{\varphi}_{m=i+(i^{\prime}-1)\bar{M}}^{\rm hf}(\mathbf{X})=\ell_{i}(X_{1})\ell_{i^{\prime}}(X_{2})\,X_{1}(1-X_{1})\,\mathbf{e}_{1}\\[5.69054pt] \boldsymbol{\varphi}_{m=\bar{M}^{2}+i+(i^{\prime}-1)\bar{M}}^{\rm hf}(\mathbf{X})=\ell_{i}(X_{1})\ell_{i^{\prime}}(X_{2})\,X_{2}(1-X_{2})\,\mathbf{e}_{2}\\ \end{array}\right.\quad i,i^{\prime}=1,\ldots,\bar{M},$$

$${\int }_{\Omega }\exp \left(\frac{ϵ-{\mathfrak{J}}_{\mathbf{a}}^{\mathrm{hf}}(\mathbf{X})}{C_{\exp }}\right)+\exp \left(\frac{{\mathfrak{J}}_{\mathbf{a}}^{\mathrm{hf}}(\mathbf{X})-1/ϵ}{C_{\exp }}\right)dX\le \delta ,$$

$$\begin{array}[]{l}\displaystyle{\int_{\Omega}\,{\rm exp}\left(\frac{\epsilon-\mathfrak{J}_{\mathbf{a}}^{\rm hf}(\mathbf{X})}{C_{\rm exp}}\right)\,+\,{\rm exp}\left(\frac{\mathfrak{J}_{\mathbf{a}}^{\rm hf}(\mathbf{X})-1/\epsilon}{C_{\rm exp}}\right)\,dX\,\geq\,\frac{\pi r^{2}}{4}\,{\rm exp}\left(\frac{\epsilon}{2C_{\rm exp}}\right)}\\[8.53581pt] \displaystyle{\geq\frac{\pi}{16}\,\left(\|\widehat{\nabla}\mathfrak{J}_{\mathbf{a}}^{\rm hf}\|_{L^{\infty}(\Omega)}\right)^{-2}\,\epsilon^{2}\,{\rm exp}\left(\frac{\epsilon}{2C_{\rm exp}}\right)=:\delta_{0},}\end{array}$$

$$\delta \ge |\Omega |\left(\exp \left(\frac{ϵ-1}{C_{\exp }}\right)+\exp \left(\frac{1-1/ϵ}{C_{\exp }}\right)\right),$$

$$ϵ=0.1,C_{\exp }=0.025ϵ,\delta =|\Omega |.$$

$$\left\{\begin{array}[]{l}\displaystyle{\min_{\mathbf{X}\in\overline{\Omega}_{\rm box}}\,\,{\rm det}\left(\widehat{\nabla}\boldsymbol{\Phi}(\mathbf{X})\right)\,>\,0,}\\[8.53581pt] \displaystyle{{\rm dist}\left(\boldsymbol{\Phi}(\mathbf{X}),\partial V\right)=0\quad\forall\,\mathbf{X}\in\partial U.}\\ \end{array}\right.$$

$$\mathrm{dist}\left(\mathbf{\Phi }(\mathbf{X}),\partial V\right)\le \text{tol}\forall \mathbf{X}\in \partial U,$$

$$\boldsymbol{\Phi}_{\mu}(\mathbf{X}):=\mathbf{X}+\sum_{m=1}^{M}\,\left(\widehat{\mathbf{a}}_{\mu}\right)_{m}\,\boldsymbol{\varphi}_{m}(\mathbf{X}),\quad{\rm where}\;\left\{\begin{array}[]{ll}\widehat{\mathbf{a}}:\mathcal{P}\to\mathbb{R}^{M},&\\ \boldsymbol{\varphi}_{m}\in{\rm span}\{\boldsymbol{\varphi}_{m}^{\rm hf}\}_{m=1}^{M_{\rm hf}}&M\leq M_{\rm hf},\\ \end{array}\right.$$

$$\begin{array}[]{l}\displaystyle{\min_{\mathbf{a}\in\mathbb{R}^{M^{\rm hf}}}\,\mathfrak{f}^{k}\left(\mathbf{a}\right)\,+\,\xi\big{|}\boldsymbol{\Psi}_{\mathbf{a}}^{\rm hf}\big{|}_{H^{2}(\Omega)}^{2},}\\[8.53581pt] \displaystyle{{\rm s.t.}\;\int_{\Omega}\,{\rm exp}\left(\frac{\epsilon-\mathfrak{J}_{\mathbf{a}}^{\rm hf}(\mathbf{X})}{C_{\rm exp}}\right)\,+\,{\rm exp}\left(\frac{\mathfrak{J}_{\mathbf{a}}^{\rm hf}(\mathbf{X})-1/\epsilon}{C_{\rm exp}}\right)\,dX\leq\delta,}\\ \end{array}$$

$$\mathfrak{f}\left(\mathbf{a},\mu,\bar{\mu}\right):=\;\int_{\Omega}\;\big{\|}u_{\mu}\circ\boldsymbol{\Psi}_{\mathbf{a}}^{\rm hf}-{u}_{\bar{\mu}}\big{\|}_{2}^{2}\,dX,$$

$${\mu }^{(k+1)}:=\arg \underset{\mu \in {\Xi }_{\mathrm{train}}\setminus \{{\mu }^{(i)}{\}}_{i=1}^k}{\min }\Vert {\mu }^{(k)}-\mu {\Vert }_2,{\mu }^{(1)}=\arg \underset{\mu \in {\Xi }_{\mathrm{train}}}{\min }\Vert \mu -\bar{\mu }{\Vert }_2,$$

$${\mathbf{a}}_0^{(1)}=0,{\mathbf{a}}_0^{(k)}={\mathbf{a}}_{\mathrm{hf}}^{({\mathrm{ne}}_k)},{\mathrm{ne}}_k:=\arg \underset{i\in \{1,\dots ,k-1\}}{\min }\Vert {\mu }^{(k)}-{\mu }^{(i)}{\Vert }_2,$$

$$-C_{\infty }\Vert {\mu }^{(k)}-{\mu }^{({\mathrm{ne}}_k)}{\Vert }_2\le {\left(\mathbf{a}\right)}_m-{\left({\mathbf{a}}_{\mathrm{hf}}^{({\mathrm{ne}}_k)}\right)}_m\le C_{\infty }\Vert {\mu }^{(k)}-{\mu }^{({\mathrm{ne}}_k)}{\Vert }_2,$$

$$\underset{\mathbf{\Phi }:\Omega \to \Omega }{\min }{\int }_{\Omega }\bar{u}(\mathbf{X})\Vert \mathbf{X}-\mathbf{\Phi }(\mathbf{X}){\Vert }_2^{2}dXs.t.\bar{u}(\mathbf{X})=u(\mathbf{\Phi }(\mathbf{X}))\mathfrak{J}(\mathbf{X})\forall \mathbf{X}\in \Omega ,$$

$$\underset{\mathbf{\Phi }:\Omega \to \Omega }{\min }{\int }_{\Omega }{\left(\bar{u}(\mathbf{X})-u(\mathbf{\Phi }(\mathbf{X}))\mathfrak{J}(\mathbf{X})\right)}^{2}dX+\frac{1}{\lambda }{\int }_{\Omega }\bar{u}(\mathbf{X})\Vert \mathbf{X}-\mathbf{\Phi }(\mathbf{X}){\Vert }_2^{2}dX,$$

$${\mathbf{a}}_{\mathrm{hf}}^k\approx {\mathbf{U}}_{\mathbf{\Phi }}{\mathbf{a}}^k,{\mathbf{U}}_{\mathbf{\Phi }}\in {\mathbb{R}}^{M_{\mathrm{hf}},M},{\mathbf{U}}_{\mathbf{\Phi }}^T{\mathbf{U}}_{\mathbf{\Phi }}=\text{mathbbm}{1}_M,M<M_{\mathrm{hf}}.$$

$$\boldsymbol{\Phi}_{\mu}(\mathbf{X})\,=\,\mathbf{X}\,+\,\sum_{m=1}^{M}\,\left(\widehat{\mathbf{a}}_{\mu}\right)_{m}\,\boldsymbol{\varphi}_{m}(\mathbf{X}),\;\left\{\begin{array}[]{l}\displaystyle{\widehat{\mathbf{a}}:=\left[\widehat{a}_{1},\ldots,\widehat{a}_{M}\right]^{T}}\\[8.53581pt] \displaystyle{\boldsymbol{\varphi}_{m}(\mathbf{X})\,=\sum_{m^{\prime}=1}^{M_{\rm hf}}\,\left(\mathbf{U}_{\boldsymbol{\Phi}}\right)_{m^{\prime},m}\,\boldsymbol{\varphi}_{m^{\prime}}^{\rm hf}(\mathbf{X})}\\ \end{array}\right.$$

$${\mathcal{Y}}_M^{\mathrm{bis},ϵ}:=\{\mathbf{\Phi }=\mathbf{X}+\mathbf{\phi }(\mathbf{X}):\mathbf{\Phi }(\Omega )=\Omega ,\det \left(\hat{\nabla }\mathbf{\Phi }\right)\ge ϵa.e.\}.$$

Full text

Given the spatial domain $\Omega\subset\mathbb{R}^{d}$ and the manifold $\mathcal{M}_{\rm u}=\{u_{\mu}:\mu\in\mathcal{P}\}$ associated with the parameter domain $\mathcal{P}\subset\mathbb{R}^{P}$ and the parametric field $\mu\mapsto u_{\mu}\in L^{2}(\Omega)$, the algorithm takes as input a set of snapshots $\{u^{k}\}_{k=1}^{n_{\rm train}}\subset\mathcal{M}_{\rm u}$ and returns a parameter-dependent bijective mapping $\boldsymbol{\Phi}:\Omega\times\mathcal{P}\to\mathbb{R}^{d}$: the mapping is designed to make the mapped manifold $\{u_{\mu}\circ\boldsymbol{\Phi}_{\mu}:\,\mu\in\mathcal{P}\}$ more suited for linear compression methods. We apply the registration procedure, in combination with a linear compression method, to devise low-dimensional representations of solution manifolds with slowly-decaying Kolmogorov $N$-widths; we also consider the application to problems in parameterized geometries. We present a theoretical result to show the mathematical rigor of the registration procedure. We further present numerical results for several two-dimensional problems, to empirically demonstrate the effectivity of our proposal.

Tommaso Taddei1

1* IMB, UMR 5251, Univ. Bordeaux; 33400, Talence, France. Inria Bordeaux Sud-Ouest, Team MEMPHIS; 33400, Talence, France, tommaso.taddei@inria.fr

Keywords: parameterized partial differential equations; model order reduction; data compression; geometry registration.

1 Introduction

1.1 Background

Numerical simulations based on mathematical models represent a valuable tool to study complex phenomena of scientific interest and/or industrial value. For several applications — including design and optimization, uncertainty quantification, active control — it is important to approximate the solution (and associated quantities of interest) to the model over a range of parameters: parameters might correspond to material properties, geometric features, operating conditions. To alleviate the computational burden associated with the evaluation of the model for many values of the parameters, parameterized model order reduction (pMOR) techniques aim to generate a reduced-order model (ROM) that approximates the original system over a prescribed parameter range. In this work, we shall develop a general — i.e., independent of the underlying equation — registration procedure for pMOR applications; we shall here focus on systems modelled by stationary partial differential equations (PDEs). In computer vision and pattern recognition, registration refers to the process of finding a spatial transformation that aligns two datasets; in this paper, registration refers to the process of finding a parametric transformation that improves the linear compressibility of a given parametric manifold.

We denote by $\mu$ the set of parameters associated with the model in the prescribed parameter domain $\mathcal{P}\subset\mathbb{R}^{P}$; we denote by $\Omega\subset\mathbb{R}^{d}$ ($d=2,3$) the spatial domain, which is assumed to be Lipschitz; we also introduce the Hilbert space $\mathcal{X}$ defined over $\Omega$, endowed with the inner product $(\cdot,\cdot)$ and the induced norm $\|\cdot\|:=\sqrt{(\cdot,\cdot)}$. Then, we introduce the problem:

[TABLE]

where $\mathbf{y}_{\mu}$ is a vector of $D$ quantities of interest, $\mathfrak{s}_{\mu}:\mathcal{X}\to\mathbb{R}^{D}$ is a continuous functional, and $\mathcal{G}_{\mu}(z_{\mu},\cdot)=0$ admits a unique solution for all $\mu\in\mathcal{P}$. We denote by $\mathcal{M}:=\{z_{\mu}:\,\mu\in\mathcal{P}\}\subset\mathcal{X}$ the solution manifold associated with the parametric problem, and we denote by $\widehat{z}_{\mu}$ the approximation to $z_{\mu}$ provided by a given ROM, for $\mu\in\mathcal{P}$.

A pMOR technique for (1) relies on three building blocks: (i) a data compression strategy for the construction of a low-dimensional approximation space of the solution manifold $\mathcal{M}$; (ii) a reduced-order statement for the rapid and reliable prediction of the solution and associated quantities of interest for any value of the parameter; and (iii) an a posteriori error estimation procedure for certification. pMOR strategies rely on an offline/online computational decomposition: during the offline stage, which is computationally expensive and performed once, a ROM for (1) is generated by exploiting several high-fidelity (finite element, finite volume,…) solutions to the mathematical model for (properly-chosen) parameter values; during the online stage, which is inexpensive and performed for any new parameter value $\mu$, the solution to the ROM is computed to rapidly obtain predictions of $z_{\mu}$ and associated quantities of interest. In the past few decades, many authors have developed pMOR strategies for a broad class of problems: we refer to the surveys [5, 22, 47] for thorough introductions to pMOR.

In (1), we assume that the domain $\Omega$ is fixed (parameter-independent); however, for several applications, the problem of interest is of the form:

[TABLE]

where $\mathcal{X}_{\mu}$ is a suitable Hilbert space defined over the parameter-dependent domain $\Omega_{\mu}\subset\mathbb{R}^{d}$. Given $\mu\in\mathcal{P}$, we remark that the solution field $z_{\mu}$ is defined over the parameter-dependent domain $\Omega_{\mu}$: since pMOR procedures rely on the definition of a low-dimensional approximation defined over a parameter-independent domain, they cannot be directly applied to (2). We should thus recast the problem in a fixed domain.

Parametric mappings are used in pMOR to recast problem (1) or (2) as

[TABLE]

where $\mathfrak{s}_{\mu,\Phi}(w):=\mathfrak{s}_{\mu}(w\circ\boldsymbol{\Phi}_{\mu}^{-1})$, and $\mathcal{G}_{\mu,\Phi}(w,v):=\mathcal{G}_{\mu}(w\circ\boldsymbol{\Phi}_{\mu}^{-1},v\circ\boldsymbol{\Phi}_{\mu}^{-1})$ for all $w,v\in\mathcal{X}$. We denote by $\widetilde{\mathcal{M}}:=\{\tilde{z}_{\mu}=z_{\mu}\circ\boldsymbol{\Phi}_{\mu}:\,\mu\in\mathcal{P}\}$ the mapped solution manifold. For (1), we require that $\boldsymbol{\Phi}_{\mu}$ is a bijection from $\Omega$ into itself, for all $\mu\in\mathcal{P}$; furthermore, we require that $v\circ\boldsymbol{\Phi}_{\mu}^{-1}\in\mathcal{X}$ for all $v\in\mathcal{X}$ —if $\mathcal{X}$ is either $L^{2}(\Omega)$ or $H^{1}(\Omega)$, it suffices to require that $\boldsymbol{\Phi}$ and its inverse are Lipschitz continuous for all $\mu\in\mathcal{P}$. For (2), we require that $\boldsymbol{\Phi}_{\mu}$ is a bijection from a reference domain ${\Omega}$ into $\Omega_{\mu}$ and that $\mathcal{G}_{\mu,\Phi}$ is well-posed in $\mathcal{X}:=\{v\circ\boldsymbol{\Phi}_{\mu}:\,v\in\mathcal{X}_{\mu}\}$. As discussed in section 1.2, the use of parameterized mappings for problems of the form (1) is motivated by approximation considerations and is not strictly necessary; on the other hand, the use of mappings for problems of the form (2) serves to restate the mathematical problem in a parameter-independent spatial domain, and is thus essential for the application of pMOR procedures.

Problem-dependent mappings are broadly used in applied mathematics to compress information and ultimately simplify the solution to PDEs. $r$-adaptivity ([64]) is used in combination with discontinuous Galerkin (DG) discretization methods to approximate discontinuous solutions to hyperbolic PDEs; similarly, problem-dependent scaling factors — in effect, mappings — are used to determine self-similar (approximate) solutions to PDEs, including the boundary layer equations (see, e.g., [41, Chapter 16.4]). Furthermore, problem-dependent mappings are employed to devise high-fidelity solvers for problems in deforming domains ([16, 44]). We describe below the two classes of pMOR tasks considered in this work.

1.2 Data compression of problems with slowly-decaying Kolmogorov widths

Most approaches to data compression for (1) rely on linear approximation spaces $\mathcal{Z}_{N}={\rm span}\{\zeta_{n}\}_{n=1}^{N}\subset\mathcal{X}$ to approximate the manifold $\mathcal{M}$: we refer to this class of methods as to linear approximation (or compression) methods; given $\mathcal{Z}_{N}$, we denote by $Z_{N}:\mathbb{R}^{N}\to\mathcal{Z}_{N}$ the linear parameter-independent operator such that $Z_{N}\boldsymbol{\alpha}:=\sum_{n=1}^{N}\alpha_{n}\,\zeta_{n}$, for all $\boldsymbol{\alpha}\in\mathbb{R}^{N}$. Two well-known linear approximation methods are the Proper Orthogonal Decomposition (POD, [7, 54, 60]) and the weak-Greedy algorithm ([53, section 7.2.2]). Success of linear methods relies on the availability of low-dimensional accurate approximation spaces for the solution manifold $\mathcal{M}$. The Kolmogorov $N$-width ([45]) provides a rigorous measure of the linear reducibility of $\mathcal{M}$: given $N>0$, the Kolmogorov $N$-width $d_{N}(\mathcal{M})$ is given by

[TABLE]

where the infimum is taken over all $N$-dimensional spaces, and $\Pi_{\mathcal{Z}_{N}}:\mathcal{X}\to\mathcal{Z}_{N}$ denotes the orthogonal projection operator onto $\mathcal{Z}_{N}\subset\mathcal{X}$. In certain engineering-relevant cases, it is possible to demonstrate that the decay of $d_{N}(\mathcal{M})$ with $N$ is extremely rapid: we refer to [53, section 8.1], [2, Theorem 3.3] and [13, 14] for further details. Recalling the (quasi-)optimality properties of POD and Greedy algorithms (see [60] and [13], respectively), rapid-decaying Kolmogorov widths justify the use of linear methods.

However, linear approximation methods are fundamentally ill-suited for several classes of relevant engineering problems. As observed in [21, Example 3.5] and [57, Example 2.5], solution fields with parameter-dependent jump discontinuities — or alternatively boundary/internal layers — cannot be adequately approximated through a low-dimensional linear expansion; similarly, linear methods are inappropriate to deal with problems with discontinuous parameter-dependent coefficients. Motivated by these considerations, several authors have proposed nonlinear approximation/compression methods to deal with these problems: we here distinguish between Eulerian and Lagrangian approaches.

Eulerian approaches consider approximations of the form $\widehat{z}_{\mu}:=Z_{N,\mu}(\widehat{\boldsymbol{\alpha}}_{\mu})$, where $Z_{N,\mu}:\mathbb{R}^{N}\to\mathcal{X}$ is a suitably-chosen operator which might depend on the parameter $\mu$ and might also be nonlinear in $\boldsymbol{\alpha}$. Eulerian techniques might rely on Grassmannian learning [1, 65], convolutional auto-encoders [25, 27], transported/transformed snapshot methods [10, 36, 48, 61], displacement interpolation [49], to determine the operator $Z_{N,\mu}$. Alternatively, they might exploit adaptive local-in-parameter and/or local-in-space enrichment strategies [3, 6, 11, 17, 18, 19, 40, 42]. Note that, since $Z_{N,\mu}$ is nonlinear and/or depends on parameter, it might be difficult to develop rapid and reliable procedures for the online computations of the coefficients $\widehat{\boldsymbol{\alpha}}_{\mu}$.

Lagrangian approaches propose to exploit a linear method $\widetilde{Z}_{N}:\mathbb{R}^{N}\to\widetilde{\mathcal{Z}}_{N}\subset\mathcal{X}$ to approximate the mapped solution $\widetilde{z}_{\mu}:=z_{\mu}\circ\boldsymbol{\Phi}_{\mu}$ where $\boldsymbol{\Phi}_{\mu}:\Omega\to\Omega$ is a suitably-chosen bijection from $\Omega$ into itself: the mapping $\boldsymbol{\Phi}_{\mu}$ should be chosen to make the mapped solution manifold $\widetilde{\mathcal{M}}$ more amenable for linear approximations. Examples of Lagrangian approaches have been proposed in111 The method of freezing proposed in [38] aims at decomposing the solution into a group component and a shape component. The approach reduces to a Lagrangian approach if the group action is induced by a mapping of the underlying spatial domain.

[24, 34, 38, 57]: in [38, 57] the construction of the map is performed separately from the construction of the solution coefficients, while in [34] the authors propose to build the mapping $\boldsymbol{\Phi}_{\mu}$ and the solution coefficients $\widehat{\boldsymbol{\alpha}}_{\mu}$ simultaneously. Simultaneous learning of mapping and solution coefficients leads to a nonlinear and non-convex optimization problem even for linear PDEs. Note that any Lagrangian method is equivalent to an Eulerian method with $Z_{N,\mu}(\boldsymbol{\alpha}):=\sum_{n=1}^{N}\alpha_{n}\,\tilde{\zeta}_{n}\circ\boldsymbol{\Phi}_{\mu}^{-1}$, while the converse is not true. On the other hand, the application of pMOR techniques to the mapped problem (3) is completely standard: this is in contrast with Eulerian approaches, which might require more involved strategies for the computation of the solution coefficients $\widehat{\boldsymbol{\alpha}}_{\mu}$ (see [27]).

1.3 Reduction of problems in parameterized domains

Given the family of parameterized domains $\{\Omega_{\mu}:\mu\in\mathcal{P}\}$, we shall here identify a reference domain $\Omega$ and a bijective mapping $\boldsymbol{\Phi}$ such that $\Omega_{\mu}=\boldsymbol{\Phi}_{\mu}({\Omega})$ for all $\mu\in\mathcal{P}$; the mapping should be computable for new values of $\mu\in\mathcal{P}$, with limited computational and memory resources. Several authors have developed geometry reduction techniques based on automatic piecewise-affine maps ([53]), Radial Basis Functions (RBFs, [31]), transfinite maps ([28, 23]) and solid extension ([30]). While the approach in [53] is restricted to a specific class of parametric deformations, the other approaches — which exploit ideas originally developed in the framework of mesh deformation and surface interpolation/approximation – are broadly applicable. Despite the many recent contributions to the field, development of rapid and reliable geometry reduction techniques for large deformations is still a challenging task.

For applications to biological systems, parametric descriptions of the boundary $\partial\Omega_{\mu}$ for all $\mu$ are unavailable. For this reason, in order to find the mapping $\boldsymbol{\Phi}$, it is important to implement strategies to systematically determine parametric descriptions of $\partial\Omega_{\mu}$. This problem, which is shared by many of the geometry reduction approaches mentioned above, is not addressed in the present paper.

1.4 Contributions and outline of the paper

We propose a registration procedure that takes as input a set of snapshots $\{u^{k}=u_{\mu^{k}}\}_{k=1}^{n_{\rm train}}\subset L^{2}(\Omega)$ with $\mu^{1},\ldots,\mu^{n_{\rm train}}\in\mathcal{P}$, and returns the parametric mapping $\boldsymbol{\Phi}:\Omega\times\mathcal{P}\to\mathbb{R}^{d}$. The key features of the approach are (i) a nonlinear non-convex optimization statement that aims at reducing the difference (in a suitable metric) between a properly-chosen reference field $\bar{u}=u_{\bar{\mu}}$ and the mapped field $\tilde{u}_{\mu^{k}}=u_{\mu^{k}}\circ\boldsymbol{\Phi}_{\mu^{k}}$ for $\mu^{1},\ldots,\mu^{n_{\rm train}}$, and (ii) a generalization procedure based on kernel regression that aims at determining a mapping for all $\mu\in\mathcal{P}$ based on the available data for $\{\mu^{k}\}_{k=1}^{n_{\rm train}}$.

For nonlinear data compression, we apply the registration procedure to the solution itself ($u_{\mu}=z_{\mu}$); for geometry reduction, we apply the registration procedure to determine a mapping from $\Omega$ to $\Omega_{\mu}$, for all $\mu\in\mathcal{P}$. To demonstrate the generality of the approach for data compression, we consider the application to a boundary layer problem and to an advection-reaction problem with a parameter-dependent sharp gradient region; on the other hand, we consider the application to two diffusion problems to demonstrate the effectivity of our proposal to geometry reduction.

We observe that the present paper is related to a number of prior works. First, as the Eulerian approach in [27], our registration algorithm is independent of the underlying equation: for this reason, we believe that the approach can be applied to a broad class of problems in science and engineering. Second, our optimization statement for the construction of $\boldsymbol{\Phi}$ is related to the recent proposal by Zahr and Persson for $r$-adaptivity in the DG framework ([64]). Third, the optimization statement is tightly linked to optimal transport ([59]), which has been employed in [24] to devise a nonlinear approximation method: in section 3.1.3, we discuss similarities and differences between the mappings obtained using our method and optimal transport maps.

Fourth, the optimization statement shares also similarities with the method of slices considered in [52, 35]. In particular, the reference field $\bar{u}$ plays the role of the template in the above mentioned papers. The key difference is that the method of slices — similarly to [38] — is associated with problems with symmetries (and more in general with problems that are invariant under a given group action) and exploits specific features of the physical system; on the other hand, our approach is exclusively driven by approximation considerations and does not require any particular structure to the elements of the parametric manifold.

The outline of the paper is as follows. In section 2, we present some necessary mathematical background that will guide us in the definition of the methodology. Then, in sections 3 and 4, we discuss the development of the registration procedure for data compression and geometry reduction; for each task, we investigate performance through the vehicle of two model problems. Finally, in section 5, we provide a short summary and we discuss several potential next steps.

1.5 Notation

By way of preliminaries, we introduce notation used throughout the paper. Given the parametric mapping $\boldsymbol{\Phi}$, we denote by $\mathbf{X}$ (resp. $\mathbf{x}$) a generic point in the reference (resp. physical) configuration; we use notation $\widehat{\nabla}=[\widehat{\partial}_{1},\ldots,\widehat{\partial}_{d}]^{T}$ and ${\nabla}=[{\partial}_{1},\ldots,{\partial}_{d}]^{T}$ to refer to reference and physical gradients and we recall that $\nabla=\left(\widehat{\nabla}\boldsymbol{\Phi}_{\mu}\right)^{-T}\,\widehat{\nabla}$; we denote by $\mathfrak{J}_{\mu}={\rm det}\left(\widehat{\nabla}\,\boldsymbol{\Phi}_{\mu}\right)$ the Jacobian determinant; finally, given the field $g\in L^{2}$, we define the corresponding mapped field $\tilde{g}_{\mu}:=g\circ\boldsymbol{\Phi}_{\mu}$.

We denote by $\mathbf{e}_{1},\ldots,\mathbf{e}_{N}$ the canonical basis in $\mathbb{R}^{N}$ and by $\|\cdot\|_{2}$ the Euclidean norm. Given the open set $A\subset\mathbb{R}^{d}$, we define the corresponding indicator function $\mathbbm{1}_{A}:\mathbb{R}^{d}\to\{0,1\}$, the distance function ${\rm dist}(\mathbf{x},A)=\inf_{\mathbf{y}\in A}\,\|\mathbf{x}-\mathbf{y}\|_{2}$, the closure $\overline{A}:=A\cup\partial A$, and given the function $w:\Omega\to\mathbb{R}^{q}$, $q\geq 1$ and $A\subset\Omega$, we define the image of $w$ in $A$ as $w(A)=\{w(\mathbf{x}):\,\mathbf{x}\in A\}\subset\mathbb{R}^{q}$. Given $\delta>0$, we denote by $A_{\delta}$ the $\delta$-neighboorhood of $A$, $A_{\delta}:=\{\mathbf{x}\in\mathbb{R}^{d}:{\rm dist}(\mathbf{x},A)<\delta\}$. Furthermore, we say that the differentiable function $\mathbf{H}:A\to\mathbb{R}^{d}$ is a diffeomorphism if it is a bijection and its inverse is differentiable; given the open set $B\subset\mathbb{R}^{d}$, we say that $A$ is diffeomorphic to $B$ (and we use notation $A\simeq B$) if there exists a diffeomorphism $\mathbf{H}$ such that $\mathbf{H}(A)=B$, and we say that $A$ is compactly embedded in $B$ (and we use notation $A\Subset B$) if $\overline{A}$ is contained in $B$. Finally, we introduce the Hausdorff distance ${\rm dist}_{\rm H}(A,B):=\max\left\{\sup_{\mathbf{x}\in A}{\rm dist}(\mathbf{x},B),\,\sup_{\mathbf{y}\in B}{\rm dist}(\mathbf{y},A)\right\}$.

Numerical results rely on a high-fidelity continuous or discontinuous Galerkin finite element (FE) discretization. We denote by $\{\texttt{D}^{k}\}_{k=1}^{n_{\rm el}}$ the elements of the mesh; we denote by $N_{\rm hf}$ the number of degrees of freedom; and we denote by $\{\mathbf{x}_{q}^{\rm qd}\}_{q=1}^{N_{\rm q}}$ the quadrature points.

In this work, we resort to POD to generate low-dimensional linear approximation spaces; we use the method of snapshots ([54]) to compute POD eigenvalues and eigenvectors. Given the snapshot set $\{u^{k}\}_{k=1}^{n_{\rm train}}\subset\mathcal{M}_{\rm u}$ and the inner product $(\cdot,\cdot)$, we define the Gramian matrix $\mathbf{C}\in\mathbb{R}^{n_{\rm train},n_{\rm train}}$, $\mathbf{C}_{k,k^{\prime}}=(u^{k},u^{k^{\prime}})$, and we define the POD eigenpairs $\{(\lambda_{n},\zeta_{n})\}_{n=1}^{n_{\rm train}}$ as

[TABLE]

with $\lambda_{1}\geq\lambda_{2}\geq\ldots\lambda_{n_{\rm train}}=0$. In our implementation, we orthonormalize the modes, that is $\|\zeta_{n}\|=1$ for $n=1,\ldots,n_{\rm train}$. To stress dependence of the POD space on the choice of the inner product, we use notation $L^{2}$-POD if $(\cdot,\cdot)=(\cdot,\cdot)_{L^{2}(\Omega)}$, $H^{1}$-POD if $(\cdot,\cdot)=(\cdot,\cdot)_{H^{1}(\Omega)}$, and $\|\cdot\|_{2}$-POD if $(\cdot,\cdot)$ is the Euclidean inner product. Finally, we shall choose the size $N$ of the POD space based on the criterion

[TABLE]

where $tol_{\rm pod}>0$ is a given tolerance.

2 Affine mappings: theoretical preliminaries

Given the diffeomorphic open sets $U$ and $V$ in $\mathbb{R}^{d}$, and the target tolerance $\texttt{tol}\geq 0$, we are interested in determining $\boldsymbol{\Phi}:U\to\mathbb{R}^{d}$ such that

[TABLE]

In data compression (cf. section 3), we have $U=V=\Omega$, while in geometry reduction (cf. section 4) $U=\Omega$ and $V=\Omega_{\mu}$. In this paper, we shall seek mappings of the form

[TABLE]

where $\{\boldsymbol{\varphi}_{m}^{\rm hf}\}_{m=1}^{M_{\rm hf}}\subset{\rm Lip}(U,\mathbb{R}^{d})$. In view of the discussion, we further define the Jacobian $\mathfrak{J}_{\mathbf{a}}^{\rm hf}:={\rm det}\left(\widehat{\nabla}\boldsymbol{\Psi}_{\mathbf{a}}^{\rm hf}\right)$.

In order to devise an effective computational methodology to determine problem-dependent parameterized mappings, we shall discuss two issues: (i) the choice of the search space ${\rm span}\{\boldsymbol{\varphi}_{m}^{\rm hf}\}_{m=1}^{M_{\rm hf}}$, and (ii) the derivation of an actionable condition for the coefficients $\mathbf{a}=\left[a_{1},\ldots,a_{M_{\rm hf}}\right]^{T}$ to enforce (6).

2.1 Main result

In Proposition 2.1, we state the main result of the paper. Proof of Proposition 2.1 is provided in section A.

Proposition 2.1**.**

Let $U,V\subset\mathbb{R}^{d}$ satisfy

[TABLE]

and $\boldsymbol{\Upsilon}$ (resp. $\boldsymbol{\Lambda}$) is a diffemorphism from $\widehat{\Omega}$ to $U$ (resp. $V$) that has a $C^{1}$ extension to $\mathbb{R}^{d}$. Given $\delta>0$, let $\boldsymbol{\Phi}:U_{\delta}\to\mathbb{R}^{d}$ be a vector-valued function that satisfies

1. (i)

$\boldsymbol{\Phi}\in C^{1}\left(U_{\delta};\mathbb{R}^{d}\right)$; 2. (ii)

$\mathfrak{J}(\mathbf{X})={\rm det}\left(\widehat{\nabla}\boldsymbol{\Phi}(\mathbf{X})\right)\geq\epsilon>0$* for all $\mathbf{X}\in U$ and a given $\epsilon>0$;* 3. (iii)

${\rm dist}\left(\boldsymbol{\Phi}(\mathbf{X}),\,\partial V\right)=0$* for all $\mathbf{X}\in\partial U$ (i.e., $\boldsymbol{\Phi}(\partial U)\subseteq\partial V$).*

Then, $\boldsymbol{\Phi}$ is a bijection that maps $U$ into $V$.

We observe that the (unit) ball and the (unit) hyper-cube satisfy (8): we have indeed that $\mathcal{B}_{1}(\mathbf{0})=\{\mathbf{x}\in\mathbb{R}^{d}:\,f(\mathbf{x})=\|\mathbf{x}\|_{2}^{2}-1<0\}$, and $(0,1)^{d}=\{\mathbf{x}\in\mathbb{R}^{d}:\,f(\mathbf{x})=2\|\mathbf{x}-[1/2,1/2]\|_{\infty}-1<0\}$. We further remark that $\partial\widehat{\Omega}=\{\mathbf{x}\in\mathbb{R}^{d}:f(\mathbf{x})=0\}$: this implies that any $U$ satisfying (8) is simply connected with connected boundary. Note that, recalling standard extension theorems in analysis, any simply connected two-dimensional smooth domain satisfies (8); the latter result is in general false for $d=3$.

Corollary 2.2 illustrates an extension of Proposition 2.1 to a broader class of domains. The result is stated for $U=V=\Omega$. The extension to the more general case $U\neq V$ is here omitted.

Corollary 2.2**.**

Let $\Omega=\Omega^{\rm out}\setminus\bigcup_{i=1}^{Q}\Omega_{i}^{\rm in}$ Let $\Omega^{\rm out},\Omega_{1}^{\rm in},\ldots,\Omega_{Q}^{\rm in}\subset\mathbb{R}^{d}$ satisfy the following assumptions.

1. (i)

$\Omega_{1}^{\rm in},\ldots,\Omega_{Q}^{\rm in},\Omega^{\rm out}$* satisfy (8).* 2. (ii)

$\Omega_{1}^{\rm in},\ldots,\Omega_{Q}^{\rm in}$* are pairwise disjoint, and $\Omega_{1}^{\rm in},\ldots,\Omega_{Q}^{\rm in}\Subset\Omega^{\rm out}$.*

Let $\boldsymbol{\Phi}:{\Omega}_{\delta}^{\rm out}\to\mathbb{R}^{d}$ be a function such that

1. (i)

$\boldsymbol{\Phi}\in C^{1}\left({\Omega}_{\delta}^{\rm out};\mathbb{R}^{d}\right)$; 2. (ii)

$\mathfrak{J}(\mathbf{X})\geq\epsilon>0$* for all $\mathbf{X}\in{\Omega}^{\rm out}$ and a given $\epsilon>0$;* 3. (iii)

$\boldsymbol{\Phi}(\partial{\Omega}^{\rm out})\subseteq\partial\Omega^{\rm out}$, $\boldsymbol{\Phi}(\partial{\Omega}_{i}^{\rm in})\subseteq\partial\Omega_{i}^{\rm in}$ for $i=1,\ldots,Q$.

Then, $\boldsymbol{\Phi}$ is a bijection that maps $\Omega$ into itself.

Proof.

We first observe that, by construction, we have $\partial\Omega=\partial\Omega^{\rm out}\cup\bigcup_{i=1}^{Q}\partial\Omega_{i}^{\rm in}$. Exploiting Proposition 2.1, we find that $\boldsymbol{\Phi}$ is a bijection that maps $\Omega^{\rm out},\Omega_{1}^{\rm in},\ldots,\Omega_{Q}^{\rm in}$ into themselves. Then, we observe that

[TABLE]

where the second identity follows from the fact that $\boldsymbol{\Phi}$ is a bijection from $\Omega^{\rm out}$ into itself and $\boldsymbol{\Phi}(\Omega_{i}^{\rm in})=\Omega_{i}^{\rm in}$ for $i=1,\ldots,Q$.

Thesis follows. ∎

2.2 Implications for $U=V=(0,1)^{d}$

Next result, which is a straightforward consequence of Proposition 2.1, provides indications for the choice of $\boldsymbol{\varphi}_{1}^{\rm hf},\ldots,\boldsymbol{\varphi}_{M_{\rm hf}}^{\rm hf}$ and of the coefficients $\mathbf{a}$ for $U=V=(0,1)^{d}$.

Proposition 2.3**.**

Let $d=2$ or $d=3$ and let $\boldsymbol{\varphi}_{1}^{\rm hf},\ldots,\boldsymbol{\varphi}_{M_{\rm hf}}^{\rm hf}$ in (7) be of class $C^{1}(\mathbb{R}^{d};\,\mathbb{R}^{d})$ and satisfy

[TABLE]

Then, for any $\bar{\mathbf{a}}\in\mathbb{R}^{M_{\rm hf}}$, $\boldsymbol{\Phi}=\boldsymbol{\Psi}_{\bar{\mathbf{a}}}^{\rm hf}$ is bijective from the unit square (or cube) $\Omega=(0,1)^{d}$ into itself if

[TABLE]

Furthermore, for any $\epsilon\in(0,1)$, there exists a ball $B=\mathcal{B}_{r_{\epsilon}}(\mathbf{0})$ of radius $r_{\epsilon}>0$ centered in $\mathbf{0}$ such that $\inf_{\mathbf{X}\in{\Omega},\mathbf{a}\in B}\mathfrak{J}_{\mathbf{a}}^{\rm hf}(\mathbf{X})\geq\epsilon$.

Proof.

We first consider the case $d=2$. Recalling Proposition 2.1, we only need to verify that $\boldsymbol{\Phi}(\partial\Omega)\subset\partial\Omega$. We here verify that $\boldsymbol{\Phi}(\mathcal{E}_{\rm top})\subset\mathcal{E}_{\rm top}$ where $\mathcal{E}_{\rm top}:=\{\mathbf{X}=(t,1):\,t\in(0,1)\}$ is the top edge of $(0,1)^{2}$: proofs for the other edges are analogous.

Exploiting (9a), we have that

[TABLE]

By contradiction, we assume that $\varphi([0,1])$ is not contained in $[0,1]$. Since $\varphi(0)=0$ and $\varphi(1)=1$, $\varphi$ has a local minimum or maximum in $(0,1)$; this implies that there exists $\bar{t}\in(0,1)$ such that $\varphi^{\prime}(\bar{t})=0$. As a result, we find

[TABLE]

which implies that $\mathfrak{J}(\mathbf{X}=(\bar{t},1))=0$. Contradiction.

To extend the result to $d=3$, we show that $\boldsymbol{\Phi}(\mathcal{E}_{\rm top})\subset\mathcal{E}_{\rm top}$ where $\mathcal{E}_{\rm top}:=\{\mathbf{X}=(t,s,1):\,t,s\in(0,1)\}$ is the top face of $(0,1)^{3}$. Towards this end, we define $\widetilde{\boldsymbol{\Phi}}(t,s)=[\Phi_{1}(t,s,1),\Phi_{2}(t,s,1)]^{T}$: clearly, $\widetilde{\boldsymbol{\Phi}}$ satisfies the hypotheses of Proposition 2.3 for $d=2$; as a result, $\widetilde{\boldsymbol{\Phi}}((0,1)^{2})=(0,1)^{2}$ and thus ${\boldsymbol{\Phi}}(\mathcal{E}_{\rm top})=\mathcal{E}_{\rm top}$. Thesis follows.

Proof of the latter statement is a direct consequence of the fact that the determinant of a matrix-valued function is continuous and that $\mathfrak{J}_{\mathbf{a}}^{\rm hf}\equiv 1$ for $\mathbf{a}=\mathbf{0}$. We omit the details. ∎

Condition (9a) imposes that each edge of the square should be mapped in itself and that each corner is mapped in itself: Figure 1 provides the geometric interpretation. In our implementation, for $d=2$, we define $\boldsymbol{\varphi}_{1}^{\rm hf},\ldots,\boldsymbol{\varphi}_{M_{\rm hf}}^{\rm hf}$ as

[TABLE]

where $\{\ell_{i}\}_{i=1}^{\bar{M}}$ are the first $\bar{M}$ Legendre polynomials and $M_{\rm hf}=2\bar{M}^{2}$; note, however, that other choices satisfying (9a) (e.g., Fourier expansions) might also be considered.

Condition (10) for the coefficients $\mathbf{a}\in\mathbb{R}^{M_{\rm hf}}$ is difficult to impose computationally; however, we might replace (10) with the approximation

[TABLE]

where $\epsilon\in(0,1)$. Provided that ${\rm exp}\left(\frac{\epsilon}{C_{\rm exp}}\right)\gg 1$ and that $\|\widehat{\nabla}\mathfrak{J}_{\mathbf{a}}^{\rm hf}\|_{L^{\infty}(\Omega)}$ is moderate, this constraint enforces that the mapping is invertible for all $\mathbf{X}\in\Omega$: more formally, for all $\epsilon>0$, there exist $C_{\rm exp},\delta,C>0$ such that if $\mathbf{a}$ satisfies (12) and $\|\widehat{\nabla}\mathfrak{J}_{\mathbf{a}}^{\rm hf}\|_{L^{\infty}(\Omega)}<C$, then $\boldsymbol{\Psi}_{\mathbf{a}}^{\rm hf}$ is globally invertible. To show this statement for $d=2$, suppose that $\mathfrak{J}_{\mathbf{a}}^{\rm hf}(\mathbf{X}^{\star})=0$ for some $\mathbf{X}^{\star}\in\Omega$; then, there exists a ball $B=\mathcal{B}_{r}(\mathbf{X}^{\star})$ of radius $r\geq\left(\|\widehat{\nabla}\mathfrak{J}_{\mathbf{a}}^{\rm hf}\|_{L^{\infty}(\Omega)}\right)^{-1}\epsilon/2$ such that $\mathfrak{J}_{\mathbf{a}}^{\rm hf}(\mathbf{X})\leq\epsilon/2$ for all $\mathbf{X}\in B$. As a result, we find222The factor $\frac{1}{4}$ follows from the identity $|\mathcal{B}_{r}(\mathbf{X})\cap\Omega|\geq\frac{\pi r^{2}}{4}$, for all $r\leq\frac{1}{2}$ and $\mathbf{X}\in\Omega$. that

[TABLE]

which implies that $\mathbf{a}$ does not satisfies (12) for $\delta<\delta_{0}$.

We can interpret $\epsilon\in(0,1)$ as the maximum allowed pointwise contraction induced by the mapping $\boldsymbol{\Psi}^{\rm hf}$ and by its inverse. We further remark that the constant $\delta$ should satisfy

[TABLE]

so that $\mathbf{a}=\mathbf{0}$ is admissible. In all our numerical examples, we choose

[TABLE]

2.3 Extension to a more general class of domains

We can exploit the results proved in sections 2.1 and 2.2 to obtain the following result. The proof is straightforward and is omitted.

Proposition 2.4**.**

Let $\boldsymbol{\varphi}_{1}^{\rm hf},\ldots,\boldsymbol{\varphi}_{M_{\rm hf}}^{\rm hf}$ in (7) be of class $C^{1}(\mathbb{R}^{d};\,\mathbb{R}^{d})$ and satisfy (9a). Let $U,V\subset\Omega_{\rm box}=(0,1)^{d}$ satisfy (8).

Then, for any $\bar{\mathbf{a}}\in\mathbb{R}^{M_{\rm hf}}$, $\boldsymbol{\Phi}=\boldsymbol{\Psi}_{\bar{\mathbf{a}}}^{\rm hf}$ is bijective from $U$ into $V$ if

[TABLE]

We remark that condition (15)2 is not practical and should in practice be replaced by

[TABLE]

for some tolerance $\texttt{tol}>0$. We are currently working on the extension of Proposition 2.4 to this more general case: as discussed in section 4, this extension would rigorously justify the approach for geometry reduction proposed in this paper.

3 Data compression

In section 3.1, we present the registration procedure for $U=V=\Omega$, for data compression. We shall here assume that $\boldsymbol{\Phi}_{\mu}(\Omega)=\Omega$ for all $\mu\in\mathcal{P}$ and $\Omega=(0,1)^{2}$: in section 2, we provide an actionable way to enforce this condition for this choice of $\Omega$. In section 3.2, we introduce the notion of Kolmogorov $N-M$ width, which is a generalization of (4). Then, in section 3.3, we discuss the integration of the proposed compression method within the standard pMOR offline/online paradigm. Finally, in section 3.4, we present numerical investigations for two two-dimensional model problems.

3.1 Registration procedure

Given the snapshot $\{u^{k}\}_{k=1}^{n_{\rm train}}\subset\mathcal{M}_{\rm u}$, we first introduce the parameterized function $\boldsymbol{\Psi}^{\rm hf}$ in (7)-(11) and a reference field $\bar{u}\in\mathcal{M}_{\rm u}$ associated with the parameter $\bar{\mu}\in\mathcal{P}$. Then, (i) for $k=1,\ldots,n_{\rm train}$, we choose $\mathbf{a}_{\rm hf}^{k}$ associated with $u^{k}$ and $\bar{u}$ by solving an optimization problem (cf. section 3.1.1); (ii) given the dataset $\{(\mu^{k},\mathbf{a}_{\rm hf}^{k})\}_{k=1}^{n_{\rm train}}$, we use a multi-target regression procedure to generate the mapping $\boldsymbol{\Phi}$ of the form

[TABLE]

for all $\mu\in\mathcal{P}$ (cf. section 3.1.4).

3.1.1 Optimization statement

We propose to choose $\mathbf{a}_{\rm hf}^{k}=\mathbf{a}_{\rm hf}(\mu^{k})$ as a solution to

[TABLE]

while the $H^{2}$ seminorm is given by $|\mathbf{v}|_{H^{2}(\Omega)}^{2}:=\sum_{i,j,k=1}^{d}\,\int_{\Omega}\,\left(\widehat{\partial}_{i,j}^{2}v_{k}\,\right)^{2}\,dX$ for all $\mathbf{v}\in H^{2}(\Omega;\mathbb{R}^{d})$. The constraint in (17a), which was introduced in (12), weakly enforces that $\mathfrak{J}_{\mathbf{a}}^{\rm hf}\in[\epsilon,1/\epsilon]$ and thus that $\boldsymbol{\Psi}^{\rm hf}$ is a bijection from $\Omega$ into itself for all admissible solutions to (17a).

Since $|v|_{H^{2}(\Omega)}=0$ for all linear polynomials, the penalty term measures deviations from linear maps: due to the condition $\boldsymbol{\Psi}_{\mathbf{a}}^{\rm hf}(\Omega)=\Omega$, we might further interpret the penalty as a measure of the deviations from the identity map. The penalty in (17a) can be further interpreted as a Tikhonov regularization, and has the effect to control the gradient of the Jacobian $\mathfrak{J}_{\mathbf{a}}^{\rm hf}$ — recalling the discussion in section 2, the latter is important to enforce bijectivity. In addition, we observe that $\mathfrak{f}^{k}(\mathbf{a})=0$ if and only if $u^{k}$ and $\bar{u}$ coincide in the mapped configuration; in particular, we have $\mathfrak{f}^{k}(\mathbf{0})=0$ if $\mu^{k}=\bar{\mu}$. The latter implies that $\mathbf{a}=\mathbf{0}$ is a solution to (17a) for $\mu^{k}=\bar{\mu}$, provided that $\delta$ satisfies (14).

The hyper-parameter $\xi$ balances accuracy — measured by $\mathfrak{f}$ — and smoothness of the mapping. We refer to the results of section 3.4 for a numerical investigation of the sensivity of the solution to the choice of $\xi$: in our experience, the choice of $\xi$ is problem-dependent; however, the same value of $\xi$ can be considered for all training points. Furthermore, since in our code $\boldsymbol{\Psi}^{\rm hf}$ is a polynomial, computation of the penalty function is straightforward: we precompute and store the matrix $\mathbf{A}^{\rm reg}\in\mathbb{R}^{M_{\rm hf},M_{\rm hf}}$ such that $A_{m,m^{\prime}}^{\rm reg}=((\boldsymbol{\varphi}_{m^{\prime}}^{\rm hf},\boldsymbol{\varphi}_{m}^{\rm hf}))_{H^{2}(\Omega)}$ for $m,m^{\prime}=1,\ldots,M_{\rm hf}$ — $((\cdot,\cdot))_{H^{2}(\Omega)}$ is the bilinear form associated with $|\cdot|_{H^{2}(\Omega)}$ — and then we compute $|\boldsymbol{\Psi}_{\mathbf{a}}^{\rm hf}|_{H^{2}(\Omega)}$ as $\mathbf{a}^{T}\,\mathbf{A}^{\rm reg}\,\mathbf{a}$.

3.1.2 Implementation

We resort to the Matlab routine fmincon [32] to solve (17): the routine relies on an interior point method ([9]) to find local minima of (17). In our implementation, we first reorder the snapshots such that

[TABLE]

In our numerical experiments, fmincon converges to local minima in $10^{2}-10^{3}$ iterations; the computational cost on a commodity laptop is $\mathcal{O}(1-5{[\rm s]})$ for all tests run. Given the snapshot set $\{u^{k}\}_{k}$, the cost per iteration is dominated by the computation of $u^{k}$ and $\nabla u^{k}$ in the mapped quadrature points $\{\boldsymbol{\Psi}_{\mathbf{a}}^{\rm hf}(\mathbf{x}_{q}^{\rm qd})\}_{q=1}^{N_{\rm q}}$ and by subsequent computation of the derivative of $\partial_{a_{m}}\mathfrak{f}^{k}$ for $m=1,\ldots,M_{\rm hf}$: for structured grids, this cost scales with $\mathcal{O}\left(N_{\rm q}\left(\log(N_{\rm hf})+M_{\rm hf}\right)\right)$.

3.1.3 Connection with optimal transport

Let us assume that $u,\bar{u}$ are probability densities over $\Omega$, that is $u,\bar{u}\geq 0$ and $\int_{\Omega}u\,dx=\int_{\Omega}\bar{u}\,dx=1$. Then, $\boldsymbol{\Phi}^{\rm opt}$ is an optimal transport map if it is a global minimizer of

[TABLE]

where $\mathfrak{J}={\rm det}\left(\widehat{\nabla}\boldsymbol{\Phi}\right)$. A barrier method to solve (19) reads as

[TABLE]

where $\lambda\gg 1$.

The first addend in (20) is closely linked to $\mathfrak{f}$ in (17b), while the second term can be interpreted as a measure of the deviation of $\boldsymbol{\Phi}$ from the identity map, exactly as $|\boldsymbol{\Phi}|_{H^{2}(\Omega)}$. Note, however, that there are important differences between (17) and (20). First, in (20), the functions $u,\bar{u}$ should be probability densities, while in (17) $u,\bar{u}$ are arbitrary real-valued functions in a Hilbert spaces $\mathcal{U}$ contained in $L^{2}(\Omega)$. Second, if $u,\bar{u}\geq 0$, solutions to (17) do not in general conserve mass. Third, if $u,\bar{u}$ are compactly supported in $\Omega$, solutions to (19) are not guaranteed to be locally invertible in $\Omega$.

3.1.4 Generalization

Given the dataset of pairs $\{(\mu^{k},\mathbf{a}_{\rm hf}^{k})\}_{k=1}^{n_{\rm train}}$, we resort to a multi-variate multi-target regression procedure to compute the mapping $\boldsymbol{\Phi}$ of the form (16). First, we resort to $\|\cdot\|_{2}$-POD to determine a low-dimensional approximation of $\{\mathbf{a}_{\rm hf}^{k}\}_{k}$:

[TABLE]

Some comments are in order. First, application of POD in (21a) leads to a (potentially substantial) reduction of the size of the mapping expansion and thus ultimately to a reduction of online costs; application of POD also reduces the importance of the choice of $M_{\rm hf}$ in (9) since it provides an automatic way of choosing the size of the expansion at the end of the offline stage. Note that, since POD is linear, the resulting expansion satisfies boundary conditions in (9a): applying Proposition 2.3, we conclude that there exists a ball $B=\mathcal{B}_{r_{\epsilon}}(\mathbf{0})$ in $\mathbb{R}^{M}$ for which $\boldsymbol{\Phi}_{\mu}$ is bijective from $\Omega$ into itself if $\widehat{\mathbf{a}}_{\mu}\in B$.

Second, any linear or nonlinear strategy for multivariate regression can be employed to construct $\widehat{\mathbf{a}}$ based on $\{(\mu^{k},\mathbf{a}_{\rm hf}^{k})\}_{k=1}^{n_{\rm train}}$: in this work, we resort to a kernel-based Ridge regression procedure based on inverse multiquadric RBFs ([63]). Third, we remark that the mapping $\boldsymbol{\Phi}_{\mu}$ is not guaranteed to be bijective for all $\mu\in\mathcal{P}$, particularly for small-to-moderate values of $n_{\rm train}$: this represents a major issue of the methodology that will be addressed in a subsequent work. We remark, nevertheless, that our approach is able to provide accurate and stable results for $n_{\rm train}=\mathcal{O}(10^{2})$ for all test cases considered in this paper. Finally, given the eigenvalues $\{\lambda_{i}\}_{i=1}^{n_{\rm train}}$ of the POD Kernel matrix $\mathbf{C}_{k,k^{\prime}}:=\mathbf{a}^{k}\cdot\mathbf{a}^{k^{\prime}}$, we choose $M$ in (21a) based on the criterion (5), for some problem-dependent tolerance $tol_{\rm POD}>0$ that will be specified in the numerical sections.

3.1.5 Review of the computational procedure

Algorithm 1 summarizes the pieces of the general approach proposed in this paper.

3.2 Two-level approximations and $N$-$M$ Kolmogorov widths

We can interpret Lagrangian approximations as two-level approximations of parametric fields, $u_{\mu}\approx\widehat{u}_{\mu}\circ\boldsymbol{\Phi}_{\mu}^{-1}$: the inner layer corresponds to the mapping process, while the outer layer is associated with the linear approximation. This observation highlights the connection between Lagrangian approaches and deep networks. As for deep vs shallow networks ([46]), the use of Lagrangian-based pMOR methods for a given class of problems should be supported by evidence of their superior approximation power compared to linear approaches. While performance of linear methods can be theoretically measured through the Kolmogorov $N$-width $d_{N}(\mathcal{M}_{\rm u})$ in (4), we here introduce the notion of Kolmogorov $N-M$ width (see also [50]) to assess performance of Lagrangian approximations.

In view of the definition of the $N-M$ width, given the $M$-dimensional space $\mathcal{Y}_{M}:={\rm span}\{\boldsymbol{\varphi}_{m}\}_{m=1}^{M}\subset{\rm Lip}(\Omega;\mathbb{R}^{d})$, and the tolerance $\epsilon\in(0,1)$, we define $\mathcal{Y}_{M}^{\rm bis,\epsilon}$ such that

[TABLE]

Note that, for any $N,M\geq 1$, $d_{N,M,\epsilon}(\mathcal{M}_{\rm u};\,\mathcal{U})\leq d_{N}\left({\mathcal{M}}_{\rm u};\,\mathcal{U}\right)$, $d_{N,0,\epsilon}(\mathcal{M}_{\rm u})=d_{N}\left({\mathcal{M}}_{\rm u};\mathcal{U}\right)$, and $d_{N,M,\epsilon=1}(\mathcal{M}_{\rm u};\mathcal{U})=d_{N}\left({\mathcal{M}}_{\rm u};\mathcal{U}\right)$.

We can provide estimates of the $N-M$ width for representative solution manifolds: we refer to section B for the proofs. First, we consider the solution $u_{\mu}$ to the one-dimensional problem

[TABLE]

Second, we consider $u_{\mu}={\rm sign}\left(x-\mu\right)$, $x\in\Omega_{\rm 1D}$, $\mu\in\mathcal{P}=[1/3,2/3]$: in this case, we find

[TABLE]

Note that the manifold in (23) presents a boundary layer at $x=0$, while the manifold in (24) is associated with an advection-reaction problem with time interpreted as parameter: therefore, these estimates suggest that Lagrangian methods might be effective for problems with boundary layers and/or travelling waves. As third and last example, we consider the parametric field

[TABLE]

3.3 Integration within the offline/online paradigm

Algorithm 2 summarizes the key steps of pMOR procedures for (1) based on Lagrangian nonlinear data compression. First, we generate a set of snapshots $u^{1},\ldots,u^{n_{\rm train}}\in\mathcal{M}_{\rm u}$ associated with the parameters $\mu^{1},\ldots,\mu^{n_{\rm train}}\in\mathcal{P}$: the field $u$ might be the solution to the PDE or a set of model coefficients. Then, we use the snapshot set $\{u^{k}\}_{k}$ to generate the problem-dependent mapping $\boldsymbol{\Phi}$, and we recast the problem in the form (3). Then, we use standard pMOR techniques to generate a ROM for (3). During the online stage, given a new parameter value $\mu^{\star}\in\mathcal{P}$, we query the ROM to estimate the output of interest and the associated prediction error.

3.3.1 A POD-Galerkin ROM for advection-diffusion-reaction problems

We discuss the application of our approach to the advection-diffusion-reaction problem:

[TABLE]

The FE discretization of (26) (we omit the superscript $(\cdot)^{\rm hf}$ to shorten notation) reads as: find $z_{\mu}\in\mathcal{X}$ such that

[TABLE]

where $\boldsymbol{\Upsilon}_{\mu}^{\rm el},\boldsymbol{\Upsilon}_{\mu}^{\rm ed},f_{\mu},f_{\mu}^{\rm ed}$ are parameter-dependent coefficients and $A^{\rm el},A^{\rm ed}$ are parameter - independent local bilinear operators.

Given the parametric mapping $\boldsymbol{\Phi}:\Omega\times\mathcal{P}\to\Omega$, provided that the data and the mapping are sufficiently smooth, we can prove that the mapped field $\tilde{z}_{\mu}=z_{\mu}\circ\boldsymbol{\Phi}_{\mu}$ solves a problem of the form (26) with coefficients:

[TABLE]

where the symbol $\widetilde{(\cdot)}$ indicates the composition with $\boldsymbol{\Phi}_{\mu}$, and $\widehat{\mathbf{n}}$ denotes the normal in the reference configuration — which coincides with $\mathbf{n}$ on $\partial\Omega$, provided that $\boldsymbol{\Phi}_{\mu}(\Omega)=\Omega$. Exploiting (27), we find that the mapping process will simply lead to different expressions for $\boldsymbol{\Upsilon}_{\mu}^{\rm el},\boldsymbol{\Upsilon}_{\mu}^{\rm ed},f_{\mu},f_{\mu}^{\rm ed}$, which can be efficiently computed.

Problem (27) with coefficients in (28) is not expected to be parametrically-affine (see, e.g., [47]); To devise an online-efficient ROM to approximate $\tilde{z}_{\mu}$, we should thus introduce the parametrically-affine approximations $\boldsymbol{\Upsilon}_{\mu}^{\rm el,eim}$, $f_{\mu}^{\rm eim}$, $\boldsymbol{\Upsilon}_{\mu}^{\rm ed,eim}$ , $f_{\mu}^{\rm ed,eim}$ of $\boldsymbol{\Upsilon}_{\mu}^{\rm el}$, $f_{\mu}$, $\boldsymbol{\Upsilon}_{\mu}^{\rm ed}$, $f_{\mu}^{\rm ed}$,

[TABLE]

We here resort to the empirical interpolation method (EIM, [4]) and to one of its extensions to vector-valued fields (cf. [56, Appendix B]): we refer to section D for further details concerning the implementation of EIM. Then we substitute these approximations in (27) to obtain a parametrically-affine surrogate of $\mathcal{G}_{\mu}$:

[TABLE]

Then, given the reduced space $\mathcal{Z}_{N}:={\rm span}\{\zeta_{n}\}_{n=1}^{N}\subset\mathcal{X}$, we define the online-efficient Galerkin ROM:

[TABLE]

We resort to POD to build the reduced space $\mathcal{Z}_{N}$. Note that the weak or strong Greedy algorithms could also be used to build the space $\mathcal{Z}_{N}$; similarly, we might also rely on Petrov-Galerkin (minimum residual) projection to estimate the solution, and we might also consider several other hyper-reduction techniques to achieve online efficiency. Since the focus of this paper is to assess the performance of the mapping procedure, we do not further discuss these aspects in the remainder.

We conclude this section with three remarks.

Remark 3.1**.**

Error estimation. As discussed in the introduction, in pMOR it is key to estimate the prediction error. For second-order elliptic problems, if we consider the norm333Similar estimates can be obtained for other norms. $\|\cdot\|:=\sqrt{\int_{\Omega}\,\|\nabla\cdot\|_{2}^{2}\,dx}$, it is straightforward to verify that

[TABLE]

where $\widehat{z}_{\mu}$ denotes a generic approximate solution to (3), $\mathfrak{c}_{\mu}=\min_{\mathbf{X}\in\overline{\Omega}}\,\lambda_{\rm min}(\mathbf{K}_{\mu})$, $\mathfrak{C}_{\mu}=\max_{\mathbf{X}\in\overline{\Omega}}\,\lambda_{\rm max}(\mathbf{K}_{\mu})$, and $\mathbf{K}_{\mu}=\,\mathfrak{J}_{\mu}\,\widehat{\nabla}\boldsymbol{\Phi}_{\mu}^{-1}\,\widehat{\nabla}\boldsymbol{\Phi}_{\mu}^{-T}$. Provided that $\mathfrak{c}_{\mu},\mathfrak{C}_{\mu}=\mathcal{O}(1)$, traditional residual-based error estimates (see, e.g., [47, Chapter 3]) in the mapped configuration can be used to sharply bound the prediction error $\|{z}_{\mu}-\widehat{z}_{\mu}\circ\boldsymbol{\Phi}_{\mu}^{-1}\|$.

Remark 3.2**.**

Online computation of $\tilde{z}_{\mu}$. The online solution to the mapped problem involves (i) the evaluation of the mapping $\boldsymbol{\Phi}_{\mu}$ and of its gradient in the $(Q_{\rm a,el}+Q_{\rm a,ed}+Q_{\rm f,el}+Q_{\rm f,ed})$ interpolation points, and (ii) the assembling of the reduced system and its solution. The second step scales with $\mathcal{O}((Q_{\rm a,el}+Q_{\rm a,ed})N^{2}+N^{3})$, while the first step depends on the supervised learning algorithm used to estimate the $M$ mapping coefficients $\widehat{\mathbf{a}}_{\mu}$: for RBF approximations, plain implementations require $\mathcal{O}(n_{\rm train}M)$ floating point operations but several acceleration techniques are now available to dramatically reduce the costs (see [51, 62]). It is difficult to establish a connection between the size of the EIM expansions in the physical and reference configuration: if some of the coefficients are non-affine, the mapping might have the effect of simplifying hyper-reduction (cf. section 4.2.1); if most or all coefficients are affine in parameter, the mapping process likely leads to much longer expansions (see, e.g., the example in section 3.4.2).

Remark 3.3**.**

Pointwise estimation of $z_{\mu}$. Computation of $\hat{z}_{\mu}\circ\boldsymbol{\Phi}_{\mu}^{-1}(\mathbf{x})$ for a given $\mathbf{x}\in\Omega$ involves the evaluation of $\boldsymbol{\Phi}_{\mu}^{-1}(\mathbf{x})$, which requires the solution to the nonlinear problem

[TABLE]

In this work, we do not address the issue of how to efficiently evaluate $\hat{z}_{\mu}\circ\boldsymbol{\Phi}_{\mu}^{-1}$.

3.4 Numerical results

3.4.1 Approximation of boundary layers

Problem statement

We consider the diffusion-reaction problem:

[TABLE]

where $\mu\in\mathcal{P}=[\mu_{\rm min},\mu_{\rm max}]$, $\mu_{\rm min}=20,\mu_{\rm max}=200$. For large values of $\mu$, the solution exhibits a boundary layer at $\Gamma_{\rm D}$. By exploiting a standard argument, we can derive the variational formulation for the lifted field $\mathring{z}_{\mu}:=z_{\mu}-1\in\mathcal{X}:=H_{0,\Gamma_{\rm D}}^{1}(\Omega)$:

[TABLE]

Given the parametric mapping $\boldsymbol{\Phi}_{\mu}:\Omega\to\Omega$ such that $\boldsymbol{\Phi}_{\mu}(\Gamma_{\rm D})=\Gamma_{\rm D}$, we find that $\widetilde{\mathring{z}}_{\mu}=\mathring{z}_{\mu}\circ\boldsymbol{\Phi}_{\mu}\in\mathcal{X}$ satisfies:

[TABLE]

where $\mathbf{K}_{\mu}^{\star}=\,\mathfrak{J}_{\mu}\,\widehat{\nabla}\boldsymbol{\Phi}_{\mu}^{-1}\,\widehat{\nabla}\boldsymbol{\Phi}_{\mu}^{-T}$ and $\sigma_{\mu}^{\star}=\mu^{2}\mathfrak{J}_{\mu}$. We here resort to a continuous P3 FE discretization with $N_{\rm hf}=11236$ degrees of freedom on a structured triangular mesh. The grid is refined close to the boundary $\Gamma_{\rm D}$, to accurately capture the boundary layer.

Construction of the mapping

We apply the registration procedure presented in section 3.1 to the solution itself ($u_{\mu}=z_{\mu}$). We set $\bar{\mu}=\sqrt{\mu_{\rm min}\,\mu_{\rm max}}$, $\xi=10^{-10}$, and we consider a polynomial expansion with $\overline{M}=8$ in (11) ($M_{\rm hf}=128$). We refer to section B for an heuristic motivation of the choice of $\bar{\mu}$. To build the regressor $\widehat{\mathbf{a}}:\mathcal{P}\to\mathbb{R}^{M}$, we consider $n_{\rm train}=70$ log-equispaced parameters in $\mathcal{P}$, and we set $tol_{\rm pod}=10^{-4}$ in (5): for this choice of the parameters, the procedure returns an affine expansion with $M=5$ terms.

Figure 2 shows the behavior of $\|\cdot\|_{2}$-POD eigenvalues associated with $\{\mathbf{a}_{\mu^{i}}^{\rm hf}\}_{i=1}^{n_{\rm train}}$, and the $L^{2}$-POD eigenvalues associated with $\{\mathbf{K}_{\mu^{i}}^{\star}\}_{i=1}^{n_{\rm train}}$.

We remark that the decay of the POD eigenvalues cannot be directly related to the behavior of the Kolmogorov $N$-width; nevertheless, POD eigenvalues provide an heuristic measure of the linear complexity of parametric manifolds and are thus shown in this paper to investigate performance of the registration algorithm.

Results

Figure 3(a) shows the behavior of the solution field $z_{\mu}$ for $\mu=20$; as anticipated above, the solution exhibits a boundary layer close to $\Gamma_{\rm D}$. In Figure 3(b) and (c), we show slices of the solution field $t\mapsto z_{\mu}(t,t)$, $t\mapsto z_{\mu}(t,0.5)$, for three parameter values along the straight lines depicted in Figure 3(a): as expected, the size of the boundary layer strongly depends on the value of $\mu$. In Figure 4, we reproduce the results of Figure 3 for the mapped field: we observe that the mapping procedure significantly reduces the sensitivity of the solution to the value of $\mu$.

Figure 5(a) shows the behavior of the relative $H^{1}$ projection error associated with the POD space in the physical (unregistered) and mapped (registered) configurations,

[TABLE]

where $\mu^{1},\ldots,\mu^{n_{\rm test}}\overset{\rm iid}{\sim}{\rm Uniform}(\mathcal{P})$ ($n_{\rm test}=200$) and the space $\{\mathcal{Z}_{N}\}_{N}$ (resp., $\{\widetilde{\mathcal{Z}}_{N}\}_{N}$) is built applying $H^{1}$-POD to the snapshot sets $\{z_{\mu^{i}}\}_{i=1}^{n_{\rm test}}$ (resp. $\{\tilde{z}_{\mu^{i}}\}_{i=1}^{n_{\rm test}}$). On the other hand, Figure 2(b) shows the $H^{1}$-POD eigenvalues associated with $\{z_{\mu^{i}}\}_{i=1}^{n_{\rm test}}$ (unregistered) and $\{\tilde{z}_{\mu^{i}}\}_{i=1}^{n_{\rm test}}$ (registered). Note that $\widetilde{E}_{\rm proj}\leq{E}_{\rm proj}$ for $N\leq 6$, while $\widetilde{E}_{\rm proj}>{E}_{\rm proj}$ for $N>6$; similar behavior can be observed for the POD eigenvalues: to reduce the sensitivity of the solution to the value of $\mu$, the mapping process introduces some small-amplitude smaller spatial scale distortions, which then ultimately control the convergence. Nevertheless, we do emphasize that for $N=1$ the relative error satisfies $\widetilde{E}_{\rm proj}\lesssim 10^{-4}$ and is thus comparable with the accuracy of the underlying truth discretization: as a result, the “complexity overhead” due to mapping-induced distortions has little practical importance for this model problem.

In Figure 6, we investigate the impact on performance of the choice of $\xi$ in (17). More in detail, we show the behavior of the proximity measure $\mathfrak{f}(\mathbf{a}_{\rm hf}(\mu),\mu,\bar{\mu})$ and of the mapping seminorm $|\boldsymbol{\Psi}_{\mathbf{a}_{\rm hf}(\mu)}^{\rm hf}|_{H^{2}(\Omega)}$ with respect to $\xi$ for three values of $\mu$ in $\mathcal{P}$, where $\mathbf{a}_{\rm hf}(\mu)$ denotes the solution to (17) for a given $\mu\in\mathcal{P}$. Interestingly, the performance is nearly independent of $\xi$ for $\xi\lesssim 10^{-8}$, for all three choices of $\mu$.

In Figure 7, we show the behavior of the proximity measure $\mathfrak{f}(\mathbf{a}_{\rm hf}(\mu),\mu,\bar{\mu})$ and of the mapping seminorm $|\boldsymbol{\Psi}_{\mathbf{a}_{\rm hf}(\mu)}^{\rm hf}|_{H^{2}(\Omega)}$ with respect to $\mu$, for two choices of the reference parameter, $\bar{\mu}=\sqrt{\mu_{\rm min}\mu_{\rm max}}$ and $\bar{\mu}=\frac{\mu_{\rm min}+\mu_{\rm max}}{2}$. For the second choice, we chose $\epsilon=0.05$. We observe that the first choice leads to superior performance in terms of accuracy and also reduces the maximum (over $\mathcal{P}$) magnitude of the mapping seminorm: similarly to [35], we empirically observe that the choice of the reference field (template) $\bar{u}$ has a significant impact on performance.

3.4.2 Approximation of advection-dominated problems

Problem statement

We consider the advection-reaction problem:

[TABLE]

Problem (35) is a parametric hyperbolic advection-reaction problem; since $\sigma_{\mu}+\frac{1}{2}\nabla\cdot\mathbf{c}_{\mu}\geq 1>0$ for all $\mathbf{x}\in\Omega$, $\mu\in\mathcal{P}$, there exists a unique solution $z_{\mu}$ to (35) for all $\mu\in\mathcal{P}$. We refer to [8, 15] for a mathematical analysis of the problem, and for the derivation of the infinite-dimensional variational form. We here resort to a DG P3 discretization on a structured triangular mesh with $N_{\rm hf}=23120$ degrees of freedom to approximate (35); we denote by $\mathcal{X}$ the broken DG space equipped with the inner product $(w,v):=\sum_{k=1}^{n_{\rm el}}\,\int_{\texttt{D}^{k}}\,w\,v\,dx,$ and the induced norm $\|\cdot\|=\sqrt{(\cdot,\cdot)}$. Note that changes in $\mu_{1}$ lead to changes in the large-gradient region associated with the solution $z_{\mu}$ (cf. Figures (8)(a) and (b)), and are thus critical for linear approximation methods.

Exploiting the reasoning in section 3.3, we can derive the variational formulation for the mapped field $\tilde{z}_{\mu}$. We provide the explicit expressions of the terms in (27) in section C. Here, we remark that the size of the expansions in (29) are chosen based on the criterion in (5), with tolerance ${tol}_{\rm eim}=5\cdot 10^{-7}$.

Construction of the mapping

We set $\bar{\mu}$ equal to the centroid of $\mathcal{P}$, and we consider the proximity measure $\mathfrak{f}$ in (17b) with $u_{\mu}=z_{\mu}$. We further consider $\xi=10^{-3}$, and we consider a polynomial expansion with $\overline{M}=6$ in (11) ($M_{\rm hf}=72$). To build the regressor $\widehat{\mathbf{a}}:\mathcal{P}\to\mathbb{R}^{M}$, we consider $n_{\rm train}=250$ uniformly-sampled parameters in $\mathcal{P}$, and we set $tol_{\rm pod}=10^{-4}$ in (5): for this choice of the parameters, the registration procedure returns an affine expansion with $M=12$ terms.

As opposed to the previous problem, we do not expect that the mapping procedure will lead to a nearly one-dimensional mapped manifold; nevertheless, we do expect that the mapping procedure will reduce the sensitivity of the solution to parametric changes, particularly in $\mu_{1}$ — which regulates the position of the high-gradient region in $\Omega$. For this reason, we here consider a larger value of the regularization parameter $\xi$ compared to the previous two examples: larger values of $\xi$ lead to smoother mappings and thus reduce the risk of overfitting and facilitate hyper-reduction. We further consider a lower value of $M_{\rm hf}$ ($M_{\rm hf}=72$ as opposed to $M_{\rm hf}=128$) compared to the previous example: in our numerical experience, the algorithm is insensitive to the choice of $M_{\rm hf}$, for $M_{\rm hf}\geq 72$.

Results

Figure 8 shows the solution to (35) for three values of the parameters $\mu^{1}=[-\pi/10,0.3,60]$, $\mu^{2}=[\pi/10,0.7,100]$ and the centroid $\bar{\mu}=[0,0.55,80]$; on the other hand, Figure 9 shows the solution to the mapped problem $\tilde{z}_{\mu}=z_{\mu}\circ\boldsymbol{\Phi}_{\mu}$ for the same three values of the parameter. As stated above, changes in $\mu_{1}$ lead to changes in the large-gradient region which corresponds to the propagation of the inflow boundary condition in the interior of $\Omega$: the mapping $\boldsymbol{\Phi}_{\mu}$ significantly reduces the sensitivity of the solution to changes in the angle $\mu_{1}$.

Figure 10(a) shows the behavior of the average relative $L^{2}$ error

[TABLE]

for the $L^{2}$-POD-Galerkin ROM associated with the unregistered configuration, and for the $L^{2}$-POD-Galerkin ROM associated with the registered configuration. Here, $\mu^{1},\ldots$, $\mu^{n_{\rm test}}$ $\overset{\rm iid}{\sim}{\rm Uniform}(\mathcal{P})$, $n_{\rm test}=20$. In order to use the same metric for both registered and unregistered ROMs, for the registered case we compute the error as

[TABLE]

Figure 10(b) shows the $L^{2}(\Omega)$-POD eigenvalues associated with $\{z_{\mu^{i}}\}_{i=1}^{n_{\rm train}}$ (unregistered) and $\{\tilde{z}_{\mu^{i}}\}_{i=1}^{n_{\rm train}}$ (registered).

We observe that the decay rates of $E_{\rm avg}$ and of the POD eigenvalues are nearly the same for both registered and unregistered configurations; however, the multiplicative constant is significantly different: $E_{\rm avg}^{\rm reg.}\approx\frac{1}{4}E_{\rm avg}^{\rm unreg.}$ for all values of $N$ considered, while $\left(\frac{\lambda_{N}}{\lambda_{1}}\right)^{\rm reg.}\approx\frac{1}{100}\left(\frac{\lambda_{N}}{\lambda_{1}}\right)^{\rm unreg.}$ for $N\geq 2$. As a result, for any given $N$, the nonlinear mapping procedure leads to a significant improvement. This empirical observation suggests a multiplicative effect between $N$ and $M$ approximations and is thus in good agreement with the estimate of the Kolmogorov $N-M$ width in (25). On the other hand, we remark that the mapping procedure leads to a significant increase in the number of EIM modes in (29). More in detail, for the registered case, we have

[TABLE]

for the unregistered case, we have

[TABLE]

Therefore, for any given $N$, the registered ROM is more expensive in terms of memory than the unregistered ROM. This behavior of EIM can be explained by observing that most coefficients in (35) are parametrically affine in the unregistered configuration.

In Figure 11, we investigate the effect of the choice of $\xi$ in (17). Figures 11(a) and (b) show the behavior of the proximity measure $\mathfrak{f}(\mathbf{a}_{\rm hf}(\mu),\mu,\bar{\mu})$ and of the mapping seminorm $|\boldsymbol{\Psi}_{\mathbf{a}_{\rm hf}(\mu)}^{\rm hf}|_{H^{2}(\Omega)}$ with respect to $\xi$ for three values of $\mu$ in $\mathcal{P}$, where $\mathbf{a}_{\rm hf}(\mu)$ denotes the solution to (17) for a given $\mu\in\mathcal{P}$. Figure 11(c) shows the decay of the mapping $\|\cdot\|_{2}$-POD eigenvalues associated with $\{\mathbf{a}_{\mu^{i}}^{\rm hf}\}_{i=1}^{n_{\rm train}}$ for two choices of $\xi$. Similarly to the previous test case, both addends of the objective function converge to finite values for $\xi\to 0^{+}$; however, their numerical values and the threshold below which curves exhibit a plateau differ significantly between the two test cases. Furthermore, we observe that reducing $\xi$ leads to an increase of the complexity of the mapping manifold and ultimately increases the complexity of the generalization (cf. section 3.1.4) step.

4 Geometry reduction

In this section, we discuss how to adapt the registration procedure introduced in section 3 to geometry reduction. In this paper, we shall assume that the reference domain $\Omega$ and the parameterized displacement field $\mathbf{d}_{\mu}:\partial\Omega\to\mathbb{R}^{d}$ are given for all $\mu\in\mathcal{P}$: given $\mathbf{X}\in\partial\Omega$, we can thus compute the corresponding material point in the physical configuration as $\mathbf{x}^{\mu}=\mathbf{X}+\mathbf{d}_{\mu}(\mathbf{X})$. As discussed in the introduction, we remark that the boundary displacement — or equivalently a parameterization of the boundary $\partial\Omega_{\mu}$ — might not be available for various classes of problems, including biological systems.

In view of the discussion, we further introduce the rectangle $\Omega_{\rm box}=(a,b)\times(c,d)$ such that $\Omega,\Omega_{\mu}\subset\Omega_{\rm box}$ for all $\mu\in\mathcal{P}$. Then, we introduce the parameterized function $\boldsymbol{\Psi}^{\rm hf}$ in (7), and we define the bases $\{\boldsymbol{\varphi}_{m}^{\rm hf}\}_{m=1}^{M_{\rm hf}}$ as in (11) — after having applied a suitable change of variables.

4.1 Registration procedure

We introduce (i) the reference points $\{\mathbf{X}_{i}\}_{i=1}^{N_{\rm bnd}}$ $\subset\partial\Omega$, (ii) the corresponding displaced points $\{\mathbf{x}_{i}^{k}=\mathbf{X}_{i}+\mathbf{d}_{\mu^{k}}(\mathbf{X}_{i})\}_{i,k}$ for some parameters $\mu^{1},\ldots,\mu^{n_{\rm train}}\in\mathcal{P}$, and (iii) the parameterized function $\boldsymbol{\Psi}^{\rm hf}$ (7) - (11). Then, for $k=1,\ldots,n_{\rm train}$, we choose $\mathbf{a}_{\rm hf}^{k}=\mathbf{a}_{\rm hf}(\mu^{k})$ as a solution to

[TABLE]

Given the dataset $\{(\mu^{k},\mathbf{a}_{\rm hf}^{k})\}_{k=1}^{n_{\rm train}}$, we proceed as in section 3.1.4 to generate the mapping $\boldsymbol{\Phi}:\Omega\times\mathcal{P}\to\mathbb{R}^{d}$.

We observe that (36) differs from (17) due to the choice of the proximity measure. Here, $\mathfrak{f}$ is an approximation of the $L^{2}(\partial\Omega)$ error over the boundary

[TABLE]

If $\mathfrak{f}^{\infty}(\mathbf{a}^{\rm hf},\mu)=0$ for some admissible $\mathbf{a}^{\rm hf}\in\mathbb{R}^{M}$ and $\Omega$ satisfies the hypotheses of Proposition 2.1, recalling Proposition 2.4, we find that $\boldsymbol{\Psi}_{\mathbf{a}}^{\rm hf}$ is a bijection from $\Omega$ to $\Omega_{\mu}$. We further remark that, if $\boldsymbol{\Psi}_{\mathbf{a}}^{\rm hf}$ is a bijection from $\Omega_{\rm box}$ into itself, it is injective in $\Omega$.

Remark 4.1**.**

Implicit surfaces. If $\partial\Omega_{\mu}$ is represented implicitly as $\partial\Omega_{\mu}=\{\mathbf{x}\in\mathbb{R}^{d}:\;\mathit{G}_{\mu}(\mathbf{x})=0\}$ for $\mathit{G}_{\mu}:\mathbb{R}^{d}\to\mathbb{R}$, then we might consider the proximity measure:

[TABLE]

where $\{\mathbf{X}_{i}\}_{i}$ is the set of control points on $\partial\Omega$. We do not consider this choice in the numerical experiments.

4.2 Model problems

We present below the two model problems considered in the numerical experiments.

4.2.1 Diffusion problem with discontinuous coefficients

We consider the diffusion problem

[TABLE]

Since $\kappa_{\mu}$ is piecewise constant with parameter-dependent jump discontinuities, pMOR techniques are not well-suited to directly tackle (37). It is thus necessary to introduce a mapping $\boldsymbol{\Phi}:\Omega_{\rm box}\times\mathcal{P}\to\Omega_{\rm box}$ such that $\kappa_{\mu}\circ\Phi_{\mu}:\Omega_{\rm box}\to\mathbb{R}$ is parameter-independent. From the geometry reduction viewpoint, given $\Omega:=\Omega_{\rm\mu=0}$, we seek $\boldsymbol{\Phi}$ such that (i) $\boldsymbol{\Phi}_{\mu}(\Omega_{\rm box})=\Omega_{\rm box}$ and (ii) $\boldsymbol{\Phi}_{\mu}(\Omega)=\Omega_{\mu}$ for all $\mu\in\mathcal{P}$. In view of the discussion, we introduce the mapped problem for a generic map $\boldsymbol{\Phi}_{\mu}:\Omega_{\rm box}\to\Omega_{\rm box}$: find $\tilde{z}_{\mu}=z_{\mu}\circ\boldsymbol{\Phi}_{\mu}\in\mathcal{X}:=H_{0}^{1}(\Omega_{\rm box})$ such that

[TABLE]

Note that (37) is a special case of the general advection-diffusion-reaction problem considered in section 3.3.

Automatic piecewise-affine maps

Since the deformation of $\Omega_{\mu}$ is rigid, we might explicitly build a piecewise linear map $\boldsymbol{\Phi}$. First, we identify a set of control points — the black dots in Figure 12(a) — and we use them to build a coarse partition of $\Omega_{\rm box}$ — shown in Figure 12(a) as well. Then, we define a mapping $\boldsymbol{\Phi}=\boldsymbol{\Phi}^{\rm aff}$ such that (i) $\Phi_{\mu}^{\rm aff}({\Omega})=\Omega_{\mu}$ for all $\mu\in\mathcal{P}$, and (ii) $\boldsymbol{\Phi}_{\mu}^{\rm aff}$ is piecewise linear in $\mathbf{X}$ in all elements of the partition. Then, we define the mapped problem (38) with $\boldsymbol{\Phi}_{\mu}=\boldsymbol{\Phi}_{\mu}^{\rm aff}$.

Note that, for this choice of the mapping, $\mathbf{K}_{\mu}^{\star}$ and $\mathfrak{J}_{\mu}$ in (38b) are piecewise-constant in each element of the partition: this implies that (38) is parametrically-affine — that is, $\mathbf{K}_{\mu}^{\star}$ and $\mathfrak{J}_{\mu}$ are linear combinations of parameter-dependent coefficients and parameter-independent spatial fields. Therefore, the solution $\widetilde{z}_{\mu}$ can be efficiently approximated using standard pMOR (e.g., Reduced Basis) techniques.

Given the user-defined control points, Rozza el al. have developed in [53] an automatic procedure to generate the partition of $\Omega_{\rm box}$ and to determine an economic piecewise-constant approximation of the form in (38). The latter relies on symbolic manipulation techniques to identify redundant terms in the affine expansions. Furthermore, for the approach to be effective, we should consider FE triangulations that are conforming to the coarse-grained partition in Figure 12(a). Finally, the choice of the control points, which is trivial in this case, might not be straightforward for more challenging problems. In conclusion, even if the approach works remarkably well in many situations, practical implementation of the procedure in [53] is rather involved and highly problem-dependent.

4.2.2 Laplace’s equation in parameterized domains

Given the parameter domain $\mathcal{P}=[0.1,0.4]^{2}\times[0,\pi/4]$, we first introduce the parametric closed curve $t\in[0,2\pi)\mapsto\boldsymbol{\gamma}_{\rm in,\mu}(t)\in\mathbb{R}^{2}$ such that

[TABLE]

We introduce the reference domain $\Omega=\Omega_{\rm box}\setminus\Omega_{\rm in}$, with $\Omega_{\rm in}=\mathcal{B}_{1}(\mathbf{0})$: note that $\Omega$ is diffeomorphic to $\Omega_{\mu}$ for all $\mu\in\mathcal{P}$. Then, given the bijection $\boldsymbol{\Phi}_{\mu}:\Omega\to\Omega_{\mu}$, and the lift $R_{\rm D}$ such that

[TABLE]

4.2.3 Application of the registration procedure

We consider the parametric function $\boldsymbol{\Psi}^{\rm hf}$ (7) - (11): thanks to this choice, control points on $\partial\Omega_{\rm box}$ are not necessary. Since, as in the examples of section 3, the conductivity matrix $\mathbf{K}_{\mu}^{\star}$ and the Jacobian determinant $\mathfrak{J}_{\mu}$ are not expected to be affine, we resort to EIM,

[TABLE]

to obtain a parametrically-affine surrogate model. As explained in section D, the basis functions are built using POD: the size of the expansions $Q_{\rm a},Q_{\rm f}$ is chosen based on the criterion in (5).

Some comments are in order. Our approach leads to a non-affine formulation —and thus requires hyper-reduction to achieve online efficiency — for both problems; on the other hand, we observe that the approach exclusively relies on the parameterization of the boundary $\partial\Omega_{\rm in,\mu}$ and can be applied for virtually any choice of the conductivity.

4.3 Numerical results

We present below results for the two model problems introduced in section 4.2. Since the focus of this section is geometry reduction, we do not discuss the construction of the ROM for the mapped problems (38) and (40).

4.3.1 Diffusion problem with discontinuous coefficients

We choose $\xi=4\cdot 10^{-4}$, $tol_{\rm POD}=10^{-5}$, $\overline{M}=6$ ($M_{\rm hf}=72$); furthermore, we consider the proximity measure (36b) where $\{\mathbf{X}_{i}\}_{i=1}^{N_{\rm bnd}}$ is an uniform grid of $\partial\Omega_{\rm in}$ with $N_{\rm bnd}=400$ and $\{\mathbf{x}_{i}^{\mu}=\mathbf{X}_{i}\,+\,\mu\}_{i}$; finally, we consider $n_{\rm train}=16^{2}$ equispaced parameters in $\mathcal{P}$. For this choice of the parameters, our procedure returns an affine expansion with $M=6$ terms.

In Figure 13, we investigate the (linear) complexity of the parametric manifolds associated with (37). In Figure 13(a), we compute the $L^{2}(\Omega_{\rm box})$-POD eigenvalues associated with $\{\kappa_{\mu^{k}}\}_{k=1}^{n_{\rm train}}$ (unregistered) and $\{\widetilde{\kappa}_{\mu^{k}}=\kappa_{\mu^{k}}\circ\boldsymbol{\Phi}_{\mu^{k}}\}_{k=1}^{n_{\rm train}}$ (registered). We observe that $\lambda_{N=2}\leq 10^{-15}\lambda_{N=1}$ for the registered case: the mapping $\boldsymbol{\Phi}_{\mu}$ is able to “fix” the position of the jump discontinuity in the reference configuration444 More precisely, if we denote by $\{\mathbf{x}_{q}^{\rm in}\}_{q}$ (resp. $\{\mathbf{x}_{q}^{\rm out}\}_{q}$) the FE quadrature points in $\Omega$ (resp. $\Omega_{\rm box}\setminus\Omega$), we have that $\boldsymbol{\Phi}_{\mu}(\{\mathbf{x}_{q}^{\rm in}\}_{q})\subset\Omega_{\mu}$ (resp. $\boldsymbol{\Phi}_{\mu}(\{\mathbf{x}_{q}^{\rm out}\}_{q})\subset\Omega_{\rm box}\setminus\Omega_{\mu}$) for all $\mu\in\mathcal{P}$.. In Figure 13(b), we show the $L^{2}(\Omega_{\rm box})$-POD eigenvalues of $\{\mathbf{K}_{\mu^{k}}^{\star}\}_{k=1}^{n_{\rm train}}$ and $\{\mathfrak{J}_{\mu^{k}}\}_{k=1}^{n_{\rm train}}$. Finally, in Figure 13(c), we show the behavior of the $H^{1}(\Omega_{\rm box})$-POD eigenvalues of $\{z_{\mu^{k}}\}_{k=1}^{n_{\rm train}}$ (unregistered) and $\{\widetilde{z}_{\mu^{k}}\}_{k=1}^{n_{\rm train}}$ (registered): even if the mapping is built based on the diffusivity coefficient, it is also effective in improving the linear reducibility of the solution manifold.

In Figure 14, we show the behavior of the mean relative error

[TABLE]

with respect to the tolerance $tol_{\rm eim}$ associated with the choice of $Q_{\rm a},Q_{\rm f}$ in (41), for $n_{\rm test}=20$. Here, $\|\cdot\|_{\star}$ denotes either the $L^{2}(\Omega)$ norm or the $H^{1}(\Omega)$ norm:

[TABLE]

In Figure 14(b), we show the behavior of $Q_{\rm a},Q_{\rm f}$ with respect to $tol_{\rm eim}$: we observe that the relative $H^{1}$ and $L^{2}$ errors are less than $10^{-3}$, for $Q_{\rm a}=19,Q_{\rm f}=9$. We remark that using piecewise-affine mappings (cf. section 4.2.1) we might obtain an exact affine form with $Q_{\rm a}=13,Q_{\rm f}=3$. Once again, we remark that our approach is fully automatic.

4.3.2 Laplace’s equation in parameterized domains

We choose $\xi=10^{-2}$, $tol_{\rm POD}=10^{-5}$, $M_{\rm hf}=288$, and we define the proximity measure $\mathfrak{f}$ (36b) based on $N_{\rm bnd}=10^{3}$ points $\{\mathbf{X}_{i}=\boldsymbol{\gamma}^{\rm ref}(2\pi i/N_{\rm bnd})\}_{i=1}^{N_{\rm bnd}}$ with $\boldsymbol{\gamma}^{\rm ref}(t)=[\cos(t),\sin(t)]$. To generate the mapping, we consider uniform grids of $\mathcal{P}$ $\{\mu^{k}\}_{k=1}^{n_{\rm train}}$ with $n_{\rm train}=4^{3},6^{3},8^{3},10^{3}$.

In Figure 15, we show boxplots of the “in-sample” and “out-of-sample” errors given by

[TABLE]

where $\{\mathbf{X}_{i}^{\rm t}=\boldsymbol{\gamma}^{\rm ref}(t^{i})\}_{i}$ and $\{\mathbf{x}_{i,j}^{\rm t}=\boldsymbol{\gamma}_{\tilde{\mu}^{j}}(t^{i})\}_{i}$, with $t^{i}=\frac{2\pi i}{N_{\rm bnd}^{\rm t}}$, $\tilde{\mu}^{1},\ldots,\tilde{\mu}^{n_{\rm test}}\overset{\rm iid}{\sim}{\rm Uniform}(\mathcal{P})$, $N_{\rm bnd}^{\rm t}=10^{4}$, $n_{\rm test}=50$. In-sample error depends on the choice of the hyper- parameters in (17) and on the choice of the discretization parameter $M_{\rm hf}$; on the other hand, out-of-sample error depends on the choice of $tol_{\rm POD}$ and on the number of training points. Interestingly, the value of $M$ is the same (and equal to $7$) for all choices of $n_{\rm train}$.

Figure 16 replicates the results of Figure 14 for the second model problem. Here, we consider $n_{\rm train}=10^{3}$, $\xi=10^{-2}$, $M_{\rm hf}=288$, $tol_{\rm POD}=10^{-5}$, and we compute $E_{\rm avg}$ (42) based on $n_{\rm test}=50$ parameters. We find that the relative $L^{2}$ and $H^{1}$ errors are below $10^{-3}$ for $Q_{\rm a}$ equal to $17$.

In the supplementary material (cf. section E), we further provide a numerical comparison of our approach with a representative mapping procedure for a given value of $\mu$.

5 Summary and discussion

In this work, we proposed a general (i.e., independent of the underlying PDE model) registration procedure, and we applied it to data compression and geometry reduction: in data compression, our registration procedure is used in combination with a linear compression method to devise a Lagrangian nonlinear compression method; in geometry reduction, our registration procedure is used to build a parametric mapping from a reference domain to a family of parameterized domains $\{\Omega_{\mu}\}_{\mu\in\mathcal{P}}$. Several numerical results empirically motivate our proposal. Although the examples considered are rather academic, we believe that our results demonstrate the applicability ìof our approach to a broad class of relevant problems in science and engineering.

In the future, we wish to extend the approach in several directions. First, in this work, we chose the reference field $\bar{u}$ (cf. section 3) and the reference domain $\Omega$ (cf. section 4) that enters in the optimization statement a priori; in the future, we wish to develop automatic procedures to adaptively choose $\bar{u}$ and $\Omega$. Second, in order to tackle problems with more complex parametric behaviors, we wish to develop partitioning techniques that leverage the use of multiple references, and we also wish to couple the approach with domain decomposition techniques. Third, as discussed in section 3.1, a major limitation of our approach is the need for several offline simulations to build the mapping $\boldsymbol{\Phi}$. To address this issue, we wish to exploit recent advances in multi-fidelity approaches to reduce the offline computational burden. Alternatively, we wish to assess performance of projection methods to simultaneously learn mapping and solution coefficients during the online stage. Fourth, the registration procedure for data compression has been developed for $\Omega=(0,1)^{2}$: we aim to combine the proximity measures (17b) and (36b) to deal with more complex domains. Furthermore, we aim to extend Proposition 2.4 to assess the sensitivity of the registration procedure to perturbations.

Finally, we wish to study the approximation properties of Lagrangian methods based on problem-dependent mappings for a range of parametric problems. Examples in section B illustrate the effectivity of Lagrangian mappings for elementary one- and two-dimensional problems; in the future, we aim to investigate whether it is possible to prove approximation results for a broader class of problems.

Acknowledgments

The author thanks Prof. Angelo Iollo (IMB, Inria), Prof. Pierre Mounoud (IMB), and Prof. Anthony Patera (MIT) for fruitful discussions.

Appendix A Proof of Proposition 2.1

A.1 Preliminaries

Given the set $U\subset\mathbb{R}^{d}$, we denote by $\mathcal{O}_{U}$ the induced (or subspace) topology on $U$, $\mathcal{O}_{U}:=\{A\cap U:\,A\,{\rm is\,open\,in}\;\mathbb{R}^{d}\}.$ It is possible to show that the ordered pair $\mathcal{T}_{U}=(U,\mathcal{O}_{U})$ is a topological space. We say that $B$ is open in $U$ if $B\in\mathcal{O}_{U}$; similarly, we say that $B$ is closed in $U$ if the complement of $B$ in $U$, $B^{\rm c}:=U\setminus B$, belongs to $\mathcal{O}_{U}$. We further say that $\mathcal{T}_{U}$ is connected if it cannot be represented as the union of two more disjoint non-empty open subsets; we say that $\mathcal{T}_{U}$ is path-connected if there exists a path joining any two points in $U$; finally, we say that $\mathcal{T}_{U}$ is simply-connected if $\mathcal{T}_{U}$ is path-connected and every path between two points can be continuously transformed into any other such path while preserving the endpoints. Next three results are key for our discussion. Theorem A.1 is a standard result in Topology that can be found in [33], while Theorem A.2 is known as Hadamard’s implicit function theorem and is proven in [26, Chapter 6]. On the other hand, we report the proof of Lemma A.1.

Theorem A.1**.**

A set $A$ in a topological space $\mathcal{T}_{U}$ is open and closed if and only if $\partial A=\emptyset$. Furthermore, the topological space $\mathcal{T}_{U}$ is connected if and only if the only open and closed sets are the empty set and $U$.

Theorem A.2**.**

Let $M_{1},M_{2}$ be smooth and connected $d$-dimensional manifolds. Suppose $\boldsymbol{\Phi}:M_{1}\to M_{2}$ is a $C^{1}$ function such that (i) $\boldsymbol{\Phi}$ is proper (i.e., for any compact set $K\subset M_{2}$, $\boldsymbol{\Phi}^{-1}(K)$ is compact in $M_{1}$), (ii) the Jacobian matrix of $\boldsymbol{\Phi}$ is everywhere invertible, and (iii) $M_{2}$ is simply connected. Then, $\boldsymbol{\Phi}$ is a homeomorphism (hence globally bijective).

Lemma A.1**.**

Let $A,B,U$ be connected, open and Lipschitz domains such that $\partial A=\partial B$, and $\overline{A\cup B}$ is strictly contained in $\overline{U}$. Then, $A=B$.

Proof.

Let $C=A\cap B$. Clearly, $C$ is open in $A$ (since it is the intersection of open sets); furthermore, $C$ is closed in $A$ since $\partial A=\partial B$. It thus follows from Theorem A.1 that $C$ is either the empty set or $A$.

By contradiction, suppose that $C=\emptyset$, and let us define $D=A\cup B$. Since $A,B$ are open, we have $\partial(\overline{A\cup B})\subset\partial A\cup\partial B=\partial A$. Let $\mathbf{x}\in\partial A$: since (i) $\partial A=\partial B$, (ii) $A\cap B=\emptyset$, and (iii) the boundary of $A$ is smooth, there exists $\epsilon>0$ such that $\mathcal{B}_{\epsilon}(\mathbf{x})\subset\overline{A\cup B}$; therefore, $\mathbf{x}\notin\partial(\overline{A\cup B})$, which implies that $\partial(\overline{A\cup B})=\emptyset$ and thus (exploiting Theorem A.1) $\overline{A\cup B}=U$. The latter contradicts the hypothesis that $\overline{A\cup B}$ is a proper subset of $U$: we can thus conclude that $C=A$.

Exploiting the same argument, we can prove that $C=B$. Thesis follows by applying the transitive property: $A=C,B=C\Rightarrow A=B$. ∎

A.2 Proof

We split the proof in four parts.

1. $\mathbf{\boldsymbol{\Phi}(\overline{{U}})\subseteq\overline{V}}$, $\mathbf{U}=\mathbf{V}=\mathbf{\widehat{\Omega}}$. Given $U=\widehat{\Omega}=\{\mathbf{x}\in\mathbb{R}^{d}:\,f(\mathbf{x})<0\}$ where $f$ is convex, we define $g(\mathbf{X})=f(\boldsymbol{\Phi}(\mathbf{X}))$. We denote by $\mathbf{X}^{\star}$ a global maximum of $g$ in $\overline{U}$, and we define $\mathbf{x}^{\star}=\boldsymbol{\Phi}(\mathbf{X}^{\star})$. Since $\boldsymbol{\Phi}$ is locally invertible at $\mathbf{X}^{\star}$, if $\mathbf{X}^{\star}$ belongs to the interior of $U$, we must have that $\mathbf{x}^{\star}$ is a local maximum of $f$: this is not possible due to the fact that $f$ is convex. As a result, we must have $\mathbf{X}^{\star}\in\partial U$: recalling (iii), we then find

[TABLE]

which implies $\boldsymbol{\Phi}(\mathbf{X})\in\overline{U}$ for all $\mathbf{X}\in\overline{U}$.

2. $\mathbf{\boldsymbol{\Phi}(\overline{U})=\overline{V}}$, $\mathbf{U=V=\widehat{\Omega}}$. Recalling Theorem A.1, we shall simply prove that $\boldsymbol{\Phi}(\overline{U})$ is open and closed in $\overline{U}$.

$\boldsymbol{\Phi}(\overline{U})$ is closed in $\overline{U}$: since $\boldsymbol{\Phi}$ is continuous, $\boldsymbol{\Phi}(\overline{U})$ is closed in $\mathbb{R}^{d}$; since $\boldsymbol{\Phi}(\overline{U})\subseteq\overline{U}$, we then find that $\boldsymbol{\Phi}(\overline{U})$ is closed with respect to the topology of $\mathcal{T}_{\overline{U}}$.

$\boldsymbol{\Phi}(\overline{U})$ is open in $\overline{U}$. To show this, it suffices to prove that, for any $\mathbf{x}\in\boldsymbol{\Phi}(\overline{U})$, there exists $B\in\mathcal{O}_{\overline{U}}$ such that $B\subset\boldsymbol{\Phi}(\overline{U})$ and $\mathbf{x}\in B$. Given $\mathbf{x}\in\boldsymbol{\Phi}(\overline{U})$, we denote by $\mathbf{X}\in\overline{U}$ a point such that $\mathbf{x}=\boldsymbol{\Phi}(\mathbf{X})$. If $\mathbf{X}$ belongs to the interior of $U$, the proof is trivial; for this reason, we focus below on the case in which $\mathbf{X}\in\partial U$ (and thus $\mathbf{x}\in\partial U$). Exploiting the local inverse-function theorem, there exists $\bar{\eta}>0$ such that for all $\eta\leq\bar{\eta}$ $\boldsymbol{\Phi}:\mathcal{B}_{\eta}(\mathbf{X})\to A_{\mathbf{x},\eta}$ is an homeomorphism and $A_{\mathbf{x},\eta}$ is an open set of $\mathbb{R}^{d}$ containing $\{\mathbf{x}\}$. Provided that

[TABLE]

we find that $\boldsymbol{\Phi}\left(\mathcal{B}_{\eta}(\mathbf{X})\cap\overline{U}\right)=A_{\mathbf{x},\eta}\cap\overline{U}$, $\boldsymbol{\Phi}\left(\mathcal{B}_{\eta}(\mathbf{X})\setminus\overline{U}\right)=A_{\mathbf{x},\eta}\setminus\overline{U}$. This implies that the set $B:=A_{\mathbf{x},\eta}\cap\overline{U}$ is contained in $\boldsymbol{\Phi}(\overline{U})$. Since $B$ is open in $\overline{U}$, we obtain the desired result.

It remains to prove (44). Since $U$ is Lipschitz, for sufficiently small values of $\eta$, we have that

[TABLE]

Furthermore, these sets are the collection of two distinct points for $U\subset\mathbb{R}^{2}$, and closed one-dimensional curves for $U\subset\mathbb{R}^{3}$. Since $\boldsymbol{\Phi}(\partial U)\subset\partial U$, we must have $\boldsymbol{\Phi}\left(\partial\mathcal{B}_{\eta}(\mathbf{X})\cap\partial U\right)\subset\partial A_{\mathbf{x},\eta}\cap\partial U$; then, since $\boldsymbol{\Phi}$ is a local homeomorphism, we conclude that $\boldsymbol{\Phi}\left(\partial\mathcal{B}_{\eta}(\mathbf{X})\cap\partial U\right)=\partial A_{\mathbf{x},\eta}\cap\partial U$.

We define $A_{1}=\boldsymbol{\Phi}\left(\mathcal{B}_{\eta}(\mathbf{X})\cap\partial U\right)$ and $A_{2}=A_{\mathbf{x},\eta}\cap\partial U$. Exploiting once again the fact that $\boldsymbol{\Phi}$ is a local homeomorphism and the fact that $\boldsymbol{\Phi}\left(\partial\mathcal{B}_{\eta}(\mathbf{X})\cap\partial U\right)=\partial A_{\mathbf{x},\eta}\cap\partial U$, we find

[TABLE]

Since $A_{1},A_{2}$ are connected, open in $\partial U$, Lipschitz, bounded sets with the same boundary, exploiting Lemma A.1, we find $A_{1}=A_{2}$, which is (44).

3. $\mathbf{\boldsymbol{\Phi}(\overline{U})=\overline{V}}$. Recalling the definitions of $\boldsymbol{\Upsilon}$ and $\boldsymbol{\Lambda}$, we find that $\widehat{\boldsymbol{\Phi}}=\boldsymbol{\Lambda}^{-1}\circ\left(\boldsymbol{\Phi}\circ\boldsymbol{\Upsilon}\right)$ satisfies the hypotheses (i)-(ii)-(iii) of Proposition 2.1 and thus $\widehat{\boldsymbol{\Phi}}(\overline{\widehat{\Omega}})=\overline{\widehat{\Omega}}$. Since $\boldsymbol{\Lambda}$ and $\boldsymbol{\Upsilon}^{-1}$ are diffeomorphisms, we obtain

[TABLE]

which is the thesis.

4. $\mathbf{\boldsymbol{\Phi}:\overline{U}\to\overline{V}}$ is bijective. From the third part of the proof, we have that $\boldsymbol{\Phi}(\overline{U})=\overline{V}$. Since $U$ is simply connected, proof follows from Theorem A.2. Note that the condition that $\boldsymbol{\Phi}$ is proper follows directly from the fact that $\boldsymbol{\Phi}$ is continuous and $U$ is bounded.

Appendix B Estimates of Kolmorogov $N$-$M$ widths

B.1 Boundary layers

We consider the one-dimensional problem:

[TABLE]

where $\mu\in\mathcal{P}:=[\mu_{\rm min},\mu_{\rm max}=\epsilon^{-2}\mu_{\rm min}]$. The solution to (45) is given by

[TABLE]

We introduce the parametrically-affine ($M=1$) mapping

[TABLE]

which is a bijection in $[0,1]$ for $0<\delta<\min\left\{1,\frac{\mu}{\bar{\mu}}\right\}$, and we set

[TABLE]

Note that the choice of $\bar{\mu}$ in (46) minimizes the maximum value attained by either the Jacobian of the mapping or by the Jacobian of the inverse ${\Phi}_{\mu}^{-1}$ over $\mathcal{P}$ in $X=0$, $\max_{\mu\in\mathcal{P}}\,\max\left\{\partial_{X}{\Phi}_{\mu}(0),\;\partial_{x}{\Phi}_{\mu}^{-1}(0)\right\}.$ Then, we find

[TABLE]

Note that in the second inequality we used the fact that $e^{-\bar{\mu}X}$ and $e^{-\mu\Phi_{\mu}(X)}$ are monotonic decreasing in $(\delta,1)$ and $e^{-\bar{\mu}\delta}=e^{-\mu\Phi_{\mu}(\delta)}=e^{-\bar{\mu}\delta}$. We observe that the right-hand side is monotonic decreasing in $\delta$. Given $\epsilon>0$, we thus choose $\delta={\rm max}\left\{\delta^{\prime}:\,\partial_{X}\Phi_{\mu}(1)\in[\epsilon,1/\epsilon],\;\forall\,\mu\in\mathcal{P}\right\}$: by tedious but straightforward calculations, we obtain $\delta=\frac{\epsilon}{1+\epsilon}$. In conclusion, we obtain

[TABLE]

and then, recalling the definition of $\delta$ and $\bar{\mu}$,

[TABLE]

which proves (23).

B.2 Shock waves

We consider the manifold ([39, section 5.1], [57, Example 2.5])

[TABLE]

which is associated with the transport problem

[TABLE]

with time interpreted as parameter. It is possible to show (see [39, section 5.1]) that the Kolmogorov $N$-width associated with $\mathcal{M}$ satisfies $d_{N}(\mathcal{M};L^{2}(\Omega))=\mathcal{O}\left(\frac{1}{\sqrt{N}}\right)$. On the other hand, if we consider the parametrically-affine ($M=1$) mapping

[TABLE]

we find $\tilde{u}(X,t)={\rm sign}(2X-1)$, which is parameter-independent. This implies that

[TABLE]

This example suggests that mapping procedures are well-suited to tackle problems with travelling fronts. On the other hand, we observe that ${\Phi}$ in (47) is not well-defined for $t\geq 1$: the reason is that the jump discontinuity exits the domain. This shows that effective applications of our approach to general parametric problems might require the partition of the time/parameter domain in several subdomains.

B.3 A two-dimensional problem

We consider the parametric field

[TABLE]

where $\mathbf{x}=[x_{1},x_{2}]\in\Omega=(0,1)^{2}$, $f_{\mu}\in{\rm Lip}([0,1])$ and $f_{\mu}([0,1])\subset[\delta,1-\delta]$, for some $\delta>0$ and for all $\mu\in\mathcal{P}\subset\mathbb{R}^{P}$. We define the manifold $\mathcal{M}_{\rm f}=\{f_{\mu}:\mu\in\mathcal{P}\}\subset{\rm Lip}([0,1])$, and we define the $M$-term approximation $\widehat{f}$ of $f$, $\widehat{f}_{\mu}(x)=\sum_{m=1}^{M}\,(\widehat{\mathbf{a}}_{\mu})_{m}\,\varphi_{m}(x)$; then, we define the mapping

[TABLE]

Exploiting the definition of $\widehat{f}$, we find that, for all $\mu\in\mathcal{P}$,

[TABLE]

where $\boldsymbol{\varphi}_{m}(\mathbf{X})\,:=\,(\varphi_{m}(X_{1})-1/2)\,(1-|2X_{2}-1|)\mathbf{e}_{2}$ for $m=1,\ldots,M$. Furthermore, we find

[TABLE]

Recalling Corollary 2.3, provided that $\widehat{f}_{\mu}([0,1])\subset[1/2\epsilon,1-1/2\epsilon]$ for some $0<\epsilon<2\delta$ and all $\mu\in\mathcal{P}$, we obtain that $\widehat{\boldsymbol{\Phi}}_{\mu}\in\mathcal{Y}_{M}^{\rm bis,\epsilon}$ (cf. (22a)) for all $\mu\in\mathcal{P}$.

By tedious but straightforward calculations, we find that $\widetilde{u}_{\mu}=u_{\mu}\circ\boldsymbol{\Phi}_{\mu}$ satisfies

[TABLE]

Note that, recalling that $f_{\mu}([0,1])\subset[\delta,1-\delta]$ and $\widehat{f}_{\mu}([0,1])\subset[1/2\epsilon,1-1/2\epsilon]$ uniformly in $\mu$,

[TABLE]

where $\Lambda_{\rm max,\mu}=\max_{X_{1}\in[0,1]}\,\Lambda_{\mu}(X_{1})$ and $\Lambda_{\rm min,\mu}=\min_{X_{1}\in[0,1]}\,\Lambda_{\mu}(X_{1})$.

To estimate the $N$-width of $\widetilde{\mathcal{M}}_{\rm u}=\{\widetilde{u}_{\mu}:\,\mu\in\mathcal{P}\}$, we partition $[0,1]$ and $[\min_{\mu}\Lambda_{\mu,\rm min},\,\max_{\mu}\Lambda_{\mu,\rm max}]$ in $n-1$ uniform intervals, $I_{i}^{x}=(\bar{x}_{i},\bar{x}_{i+1})$ $I_{j}^{y}=(\bar{y}_{j},\bar{y}_{j+1})$ $i,j=1,\ldots,n$. Then, we introduce the $N=n^{2}$ functions

[TABLE]

We observe that

[TABLE]

Since (51) holds for any choice of $\boldsymbol{\Phi}$ satisfying (49), provided that $M$ is large enough so that there exists an $M$-term approximation satisfying $\sup_{\mu\in\mathcal{P}}\,\|\widehat{f}_{\mu}-f_{\mu}\|_{{\rm Lip}([0,1])}$ $\leq\delta-1/2\epsilon$, we find that

[TABLE]

which is (25).

Appendix C DG discretization of (35)

We denote by $\mathbb{P}^{\kappa}(\widehat{\texttt{D}})$ the space of polynomials of degree at most $\kappa$ on the triangle $\widehat{\texttt{D}}$ with vertices $(0,0)-(1,0)-(0,1)$; then, we define the broken DG space

[TABLE]

where $\{\texttt{D}^{k}\}_{k=1}^{n_{\rm el}}$ are the elements of the mesh and $\boldsymbol{\Psi}_{k}^{\rm fe}:\texttt{D}^{k}\to\texttt{D}$ are the local FE mappings. For each edge $e$ of the mesh, we define the positive (resp. negative) normal $\mathbf{n}^{+}$ (resp. $\mathbf{n}^{-}$); given $w\in\mathcal{X}$ and the mesh edge $e$, we define the positive and negative limits $w^{+},w^{-}$ and the edge average and jump

[TABLE]

for all $\mathbf{x}\in e$ If $\mathbf{x}\in e\subset\partial\Omega$, we set $\{w\}:=w$ and $\mathbf{J}w:=\mathbf{n}w$.

Then, we can introduce the high-fidelity DG discretization of (35): find $z_{\mu}\in\mathcal{X}$ such that

[TABLE]

Then, we introduce

[TABLE]

where $\delta_{\partial\Omega}(\mathbf{x})=1$ if $\mathbf{x}\in\partial\Omega$ and $\delta_{\partial\Omega}(\mathbf{x})=0$ otherwise. Exploiting this notation, we can rewrite (52) as (27). Note that the associated mapped problem is also of the form (27), provided that we substitute $\mathbf{c}_{\mu},\sigma_{\mu},f_{\mu},z_{\rm D,\mu}$ with the corresponding definitions in (28).

Appendix D Empirical Interpolation Method

D.1 Review of the interpolation procedure for scalar fields

We review the Empirical Interpolation Method (EIM, [4]), and we discuss its extension to the approximation of vector-valued fields. Given the Hilbert space $\mathcal{W}$ defined over $\Omega$, the $Q$-dimensional linear space $\mathcal{W}_{Q}={\rm span}\{\psi_{q}\}_{q=1}^{Q}\subset\mathcal{W}$ and the points $\{\mathbf{x}_{q}^{\rm i}\}_{m=1}^{Q}\subset\overline{\Omega}$, we define the interpolation operator $\mathcal{I}_{Q}:\mathcal{W}\to\mathcal{W}_{Q}$ such that $\mathcal{I}_{Q}(v)(\mathbf{x}_{q}^{\rm i})=v(\mathbf{x}_{q}^{\rm i})$ for $q=1,\ldots,Q$ for all $v\in\mathcal{W}$. Given the manifold $\mathcal{F}\subset\mathcal{W}$ and an integer $Q>0$, the objective of EIM is to determine an approximation space $\mathcal{Z}_{Q}$ and $Q$ points $\{\mathbf{x}_{q}^{\rm i}\}_{q=1}^{Q}$ such that $\mathcal{I}_{Q}(f)$ accurately approximates $f$ for all $f\in\mathcal{F}$.

Algorithm 3 summarizes the EIM procedure as implemented in our code. The algorithm takes as input snapshots of the manifold $\{f^{k}\}_{k=1}^{n_{\rm train}}\subset\mathcal{F}$ and a tolerance $tol_{\rm eim}>0$, and returns the functions $\{\psi_{q}\}_{q=1}^{Q}$, the interpolation points $\{\mathbf{x}_{q}^{\rm i}\}_{q=1}^{Q}$ and the matrix $\mathbf{B}\in\mathbb{R}^{Q,Q}$ such that $\mathbf{B}_{q,q^{\prime}}=\psi_{q}(\mathbf{x}_{q^{\prime}}^{\rm i})$. It is possible to show that the matrix $\mathbf{B}$ is lower-triangular: for this reason, online computations can be performed in $\mathcal{O}(Q^{2})$ flops. Note that in [4] the authors resort to a strong Greedy procedure to generate $\mathcal{Z}_{Q}$, while here (as in several other works including [12]) we resort to POD. A thorough comparison between the two compression strategies is beyond the scope of the present work.

Remark D.1**.**

Oversampling. Several authors have proposed to consider non-interpolatory extensions of Algorithm 3: these algorithms generate a set of $Q_{\rm s}$ points $\{\mathbf{x}_{q}^{\rm i}\}_{q=1}^{Q_{\rm s}}$, a basis $\{\psi_{q}\}_{q=1}^{Q}$, and the associated system $\mathbf{B}$, with $Q_{\rm s}=\mathfrak{o}Q$, where $\mathfrak{o}>1$ is the oversampling ratio. We refer to [43] and the references therein for further details; see also [29, Algorithm 2] for a generalization to a broader class of ”measurement” functionals .

D.2 Extension to vector-valued fields

The EIM procedure can be extended to vector-valued fields. We present below the non-interpolatory extension of EIM employed in this paper; the same approach has also been employed in [56]. We refer to [58, 37] for two alternatives applicable to vector-valued fields. Given the space $\mathcal{W}_{Q}={\rm span}\{\boldsymbol{\omega}_{q}\}_{q=1}^{Q}\subset\mathcal{W}$ and the points $\{\mathbf{x}_{q}^{\rm i}\}_{q=1}^{Q}\subset\overline{\Omega}$, we define the least-squares approximation operator $\boldsymbol{\mathcal{I}}_{Q}:\mathcal{W}\to\mathcal{W}_{Q}$ such that for all $\boldsymbol{v}\in\mathcal{W}$

[TABLE]

It is possible to show that $\boldsymbol{\mathcal{I}}_{M}$ is well-defined if and only if the matrix $\mathbf{B}\in\mathbb{R}^{QD,Q}$,

[TABLE]

Algorithm 4 summarizes the procedure employed to compute $\mathcal{W}_{Q}$, $\{\mathbf{x}_{q}^{\rm i}\}_{q=1}^{Q}$ and the matrix $\mathbf{B}$. We observe that for scalar fields the procedure reduces to the one outlined in Algorithm 3. We further observe that online computational cost scales with $\mathcal{O}(DQ^{2})$, provided that $\mathbf{B}^{\dagger}$ is computed offline.

Appendix E Radial Basis Function maps for geometry reduction

E.1 RBF formulation

We illustrate how to apply RBF approximations to build the mapping $\boldsymbol{\Phi}_{\mu}$ for a given $\mu\in\mathcal{P}$. Given an even integer $N_{\rm bnd}>0$, the parameterization of $\partial\Omega_{\rm in}$ $\boldsymbol{\gamma}_{\rm in}:[0,1)\to\partial\Omega_{\rm in}$, $\boldsymbol{\gamma}_{\rm in}(t)=[\cos(2\pi t),\sin(2\pi t)]^{T}$, and the parameterization of $\partial\Omega_{\rm box}$ $\boldsymbol{\gamma}_{\rm box}:[0,1)\to\partial\Omega_{\rm box}$, we define the control points

[TABLE]

and the displaced control points

[TABLE]

where $0\leq t^{1},\ldots\leq t^{N_{\rm bnd}/2}<1$. Then, we introduce the kernel $\phi:\mathbb{R}^{+}\to\mathbb{R}$, where $\lambda>0$ is the Kernel width, and the space $\mathbb{P}^{1}$ of linear functions from $\mathbb{R}^{2}$ to $\mathbb{R}^{2}$. Finally, we define the mapping $\boldsymbol{\Phi}_{\mu}:=\mathbf{p}_{\mu}+\boldsymbol{\omega}_{\mu}$ where $(\mathbf{p}_{\mu},\boldsymbol{\omega}_{\mu})$ is the solution to the optimization problem

[TABLE]

Here, $\xi>0$ is a regularization parameter, while $\mathcal{N}_{\phi}$ is the native Hilbert space associated with the kernel $\phi$. We anticipate that in the numerical results we consider the Gaussian kernel with kernel width $\lambda>0$, $\phi(r)=e^{-\lambda r^{2}}$ and we resort to hold-out ($80\%-20\%$) validation to choose $\lambda$ and $\xi$. It is possible to show that the optimal $\boldsymbol{\omega}_{\mu}$ is of the form $\boldsymbol{\omega}_{\mu}(\mathbf{X})=\sum_{i=1}^{N_{\rm bnd}}\mathbf{a}_{i}\,\phi\left(\|\mathbf{X}-\mathbf{X}_{i}\|_{2}\right)$; as a result, solutions to (55) involve the solution to a linear system of size $d(N_{\rm bnd}+1)+d^{2}$. Note that in [31] the authors consider the pure interpolation problem, which corresponds to taking the limit $\xi\to 0^{+}$ in (55) (see [55, Proposition 2.10]).

We observe that $\boldsymbol{\Phi}_{\mu}$ defined in (55) is not guaranteed to be bijective for large deformations: this issue is shared by many approaches referenced in the introduction555 To address this issue, in the related framework of mesh deformation, several authors (see, e.g., [20]) have proposed to resort to nonlinear elasticity extensions: clearly, resorting to a nonlinear extension increases the overall computational cost. . Furthermore, computational cost scales with $\mathcal{O}(N_{\rm bnd}^{3})$: as $N_{\rm bnd}$ increases, the computational overhead associated with the mapping process might be the dominant online cost. On the other hand, for the approach presented in this paper (i) the size of the expansion does not depend on the number of control points, and (ii) the mapping is guaranteed to be globally invertible for all $\mu$ in the training set $\{\mu^{k}\}_{k=1}^{n_{\rm train}}\subset\mathcal{P}$.

E.2 Numerical results

In Figure 17, we compare performance of the RBF map obtained solving (55) with performance of the map obtained solving (17) for several values of $M_{\rm hf}$, and for a given $\mu\in\mathcal{P}$. Figure 17(a) shows the behavior of the out-of-sample error $E^{\rm geo,out}$ (cf. (43)) with respect to the number of control points $N_{\rm bnd}$; Figure 17(b) shows the behavior of the minimum Jacobian determinant over $\Omega$ with $N_{\rm bnd}$. Control points for RBF are chosen on both boundaries $\partial\Omega_{\rm box},\partial\Omega_{\rm in}$ as described above; on the other hand, in our procedure, in (17), we consider $N_{\rm bnd}=10^{3}$ points on $\partial\Omega_{\rm in}$. We choose $\epsilon=0.1$ for all tests, while we observe that higher values of $\xi$ should be considered for small $M_{\rm hf}$ to avoid overfitting: here, we set $\xi=10^{-1}$ for $M_{\rm hf}\geq 2\cdot 4^{2}$ and $\xi=10$ for $M_{\rm hf}<2\cdot 4^{2}$. We observe that convergence of the RBF mapping is relatively slow: this is due to the lack of regularity of $\partial\Omega_{\rm box}$. We further observe that for small values of $M_{\rm hf}=N_{\rm bnd}$ the RBF mapping might be singular. On the other hand, our approach is guaranteed to lead to bijective maps for all choices of $M_{\rm hf}$.