GFN for reduced operator learning in multifidelity applications]GFN: A graph feedforward network for resolution-invariant reduced operator learning in multifidelity applications

Oisín M. Morrison^1 , Federico Pichi^1,2 , and Jan S. Hesthaven^1
$^1$ Chair of Computational Mathematics and Simulation Science, École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland
$^2$ mathLab, Mathematics Area, SISSA, via Bonomea 265, I-34136 Trieste, Italy.

This work presents a novel resolution-invariant model order reduction strategy for multifidelity applications. We base our architecture on a novel neural network layer developed in this work, the graph feedforward network, which extends the concept of feedforward networks to graph-structured data by creating a direct link between the weights of a neural network and the nodes of a mesh, enhancing the interpretability of the network. We exploit the method's capability of training and testing on different mesh sizes in an autoencoder-based reduction strategy for parametrised partial differential equations. We show that this extension comes with provable guarantees on the performance via error bounds. The capabilities of the proposed methodology are tested on three challenging benchmarks, including advection-dominated phenomena and problems with a high-dimensional parameter space. The method results in a more lightweight and highly flexible strategy when compared to state-of-the-art models, while showing excellent generalisation performance in both single fidelity and multifidelity scenarios.

Keywords: graph neural networks, model order reduction, operator learning, resolution invariance, multifidelity surrogate modelling, parametrised PDEs.

Code availability: <https://github.com/Oisin-M/GFN>

§ INTRODUCTION AND MOTIVATION

Traditional numerical solvers for the high-fidelity approximation of partial differential equations (PDEs) are computationally prohibitive in real-time and many-query contexts, underpinning the need for quicker evaluations of the numerical solution of PDEs. Reduced order models (ROMs) have arisen as a means to address this issue for parametrised PDEs, accelerating the process via the creation of efficient computational models [10], with fewer degrees of freedom.

Model order reduction (MOR) is therefore concerned with the creation of offline-online surrogate models that are cheaper whilst maintaining high levels of accuracy. One such popular approach is the reduced basis method, characterised by the creation of a reduced space, typically obtained via proper orthogonal decomposition (POD) or greedy methods [79, 46]. Such methods are generally linear, rendering them inefficient for problems exhibiting a slow decay in the Kolmogorov $n$-width [36, 73]. To combat this, a range of nonlinear ROM techniques have been developed, with a notable subclass of such methods being machine learning-based approaches. In particular, much promise has been shown in incorporating autoencoder-based architectures into the ROM context [62, 33, 77, 88] due to their feasibility as a nonlinear extension of principal component analysis [57], displaying superior compression capabilities compared to linear methods especially when dealing with problems with slow Kolmogorov $n$-width decay [36, 73].

Such autoencoder-based MOR strategies typically benefit from being non-intrusive approaches, removing the requirement for any knowledge of the generation procedure for the training data, therefore leading to possible exploitation of experimental data [31, 30, 59].
Moreover, contrarily to physics-informed neural networks (PINNs) [83] or sparse identification of nonlinear dynamics (SINDy) [14], these are usually not physics-based methods, and hence render it possible to learn the solution's behaviour even without prior knowledge.

In this direction, some approaches seek to recover the operator mappings, also known as operator inference, on some fixed discretised domain [75, 78, 56]. Instead of approximating the solution at some fixed locations, advanced neural operators (NOs) go further by directly approximating the solution operator itself [55, 13]. As a result, NOs benefit from the desirable property of being independent from a fixed discretisation, possibly achieving super- and sub-resolution [64, 86, 94, 55].

The development of so-called resolution-invariant methods, i.e. methods without a dependence on a fixed discretisation, can be of the utmost importance when dealing with experimental data. This is especially the case in climate modelling, where information depends on both weather stations and satellites, the number and locations of which can vary greatly [26, 25, 72]. These considerations have lead to similar and complementary approaches including the aforementioned NOs, graph neural networks (GNNs) and multifidelity reduction strategies.

Multifidelity ROMs are typically focused on leveraging cheaper computational data, known as low-fidelity data, in addition to more expensive high-fidelity data, to build an efficient offline MOR strategy. Low-fidelity may be interpreted as simplification of the governing equations, low-order approximations, or scarce information, but here we focus on the case of coarser discretisations, aiming to obtain a ROM capable of learning from data on multiple fidelities. State-of-the-art approaches for ROM [39, 52, 70, 19], however, are not generally capable of learning from any arbitrary resolution, but rather from a prespecified fixed set of possible discretisations.

On the other hand, GNNs are capable of leveraging data from graphs of arbitrary cardinality, while preserving the structure of the problem by enforcing geometric priors into the model. A large number of such works (reviewed in [101, 102, 66]) have been focused on the generalisation of pooling, unpooling and convolutional layers to graphs, achieving excellent results across a variety of domains including drug development [93], traffic prediction [49], finance [97], climate [60] and fluid mechanics [44, 9].
Moreover, being designed for unstructured data, GNNs have an advantage over many current state-of-the-art NOs such as Fourier neural operators [64] and convolutional neural operators [86] which require data on structured grids.

GFN-ROM is a nonlinear non-intrusive multifidelity ROM capable of dealing with unstructured data, interplaying between MOR and resolution-invariant techniques.

Contributions. In this work, we make the following contributions:

* We introduce a novel neural network layer called the graph feedforward network (GFN) which generalises the feedforward layer for data on arbitrary discretisations.

* We introduce GFN-ROM, a new autoencoder-based ROM approach for multifidelity data using GFNs.

* We provide bounds on the super- and sub-resolution errors using GFNs and GFN-ROM.

* We demonstrate the remarkable generalisation properties of GFN-ROM through a collection of challenging problems.

To the best of our knowledge, GFN-ROM is the only graph-based resolution-invariant ROM. By associating weights and biases to nodes in a mesh, GFN-ROM also allows for better interpretability compared to approaches such as DL-ROM or GCA-ROM. Furthermore, we show that GFN-ROM is closely related to the standard definitions of NOs and GNNs, by retaining important properties of both concepts (see <ref>).

The remainder of the paper is organised as follows. [sec:ROM]Section <ref>sec:ROM presents the standard intrusive and non-intrusive framework for MOR. [sec:MROM]Section <ref>sec:MROM defines the challenges in multifidelity applications, while [sec:resinv_method]Section <ref>sec:resinv_method introduces resolution invariant methods. The core of the manuscript is presented in [sec:GFN]Section <ref>sec:GFN where the GFN is introduced to leverage data on any discretisation. [sec:GFN-ROM]Section <ref>sec:GFN-ROM and [sec:adaptive_training]Section <ref>sec:adaptive_training present the exploitation of GFN in the ROM context, and how to perform adaptive multifidelity training. Finally, [sec:results]Section <ref>sec:results shows the performance of the proposed methodology on complex MOR benchmarks, and [sec:conc]Section <ref>sec:conc gives some conclusions and highlights future perspectives.

§ REDUCED ORDER MODELS

Full order models are designed to solve PDEs via high-fidelity systems of equations. Such systems involve a large number of degrees of freedom $N_h$, which makes their numerical solution costly to compute. In a parametrised PDE context, this is particularly prohibitive since one aims at recovering solutions in real-time for a number of different physical or geometrical configurations, giving rise to a parameter-dependent high-fidelity solution $\bfmu_h(\boldsymbol{\mu})$. ROMs seek to mitigate this issue by building and subsequently solving a cheaper reduced system with a much lower number of degrees of freedom $N \ll N_h$, giving rise to the reduced order solutions $\bfmu_N(\boldsymbol{\mu})$. High-fidelity solutions can then be approximated from the reduced order ones via some transform $\bfmu_h(\boldsymbol{\mu}) \approx \phi \left( \bfmu_N(\boldsymbol{\mu}) \right)$, where $\phi$ is a linear or nonlinear mapping dependent on the choice of ROM employed. Such methods can be constructed in an intrusive or non-intrusive fashion, depending on whether knowledge of the high-fidelity system is required for the computation of $\bfmu_N(\boldsymbol{\mu})$.

§.§ Projection-based intrusive methods

Intrusive methods require knowledge of the high-fidelity system. For ease of illustration, we consider the case of a linear high-fidelity system:

_h(μ) _h(μ) = _h(μ),
where $\boldsymbol{\mu} \in \mathbb{R}^{N_\mu}$ is the parameter, $\bfmA_h \in \mathbb{R}^{N_h \times N_h}$ is the stiffness matrix, $\bfmu_h \in \mathbb{R}^{N_h \times  N_{V_F}}$ is the high-fidelity solution, and $\bfmf_h \in \mathbb{R}^{N_h \times N_{V_F}}$ is the forcing term[For ease of notation, in the following we drop the dependence on the parameters $\boldsymbol{\mu}$.].

A ROM aims at replacing <ref> with a cheaper reduced system. Linear MOR techniques express an approximation of the high-fidelity solution $\tilde{\bfmu}_h$ as a linear expansion over some chosen basis functions $\{\boldsymbol{\psi}_i\}_{i=1}^N$ with $\boldsymbol{\psi}_i \in \mathbb{R}^{{N_h}}$. Introducing the matrix $\bfmV = \left[ \boldsymbol{\psi}_1 | \cdots | \boldsymbol{\psi}_N \right]$ one can concisely write $\bfmu_h \approx \bfmV \bfmu_N$. POD is a popular choice for constructing a basis, extracting the principal components from a series of high-fidelity solutions for different parameter realisations, known as snapshots. POD provides the best rank $N$ subspace for approximating the snapshots in a least-squares sense [69]. Once the set of basis functions has been selected, it remains to calculate the reduced coefficients $\bfmu_N$. Intrusive projection-based MOR approaches exploit the knowledge of the high-fidelity system by imposing the projected residual of the reduced order solution onto $\bfmV$ is zero i.e. one solves for $\bfmu_N$ the system

^T (_h - _h _N) = 0.
An efficient application in the many-query context requires <ref> to exhibit an affine parametric dependence in order to circumvent $\mathbb{R}^{N_h}$-dependent online operations. This is usually not fulfilled when dealing with complex and/or nonlinear problems.

§.§ Non-intrusive methods

When the high-fidelity system is not available, e.g. black-box or commercial solvers, intrusive MOR cannot be pursued and the reduced coefficients have to be recovered in a non-intrusive manner. For example, POD with interpolation (PODI) recovers the reduced coefficients from the projection of the snapshots onto the POD basis [15, 23]. A similar regression-based  machine learning approach, POD-NN, is also possible by training a network to predict the reduced coefficients, with the projected coefficients used as the training dataset [47, 7].

Although non-intrusive, PODI, POD-NN and other such methods still represent linear approaches since they approximate the high-fidelity coefficients via a linear scheme $\bfmu_h \approx \bfmV \bfmu_N$. This can be problematic, since many important physics phenomena exhibit a slow decay in the Kolmogorov $n$-width, making linear methods inefficient for such cases [36, 73]. As a result, nonlinear approaches have arisen which seek to avoid this limitation by instead prescribing the more general form of $\bfmu_h(\boldsymbol{\mu}) \approx \psi \left( \bfmu_N(\boldsymbol{\mu}) \right)$, where $\psi$ is a nonlinear mapping. A huge variety of nonlinear compression schemes have been explored for standard MOR: using local reduced bases to construct a piecewise linear scheme [3], considering nonlinear compression via kernel POD [24, 54], shifted POD [87] or registration methods [95].

Machine learning approaches are also possible, with the advantage of typically being cheap at inference. In this direction, Gaussian process regression (GPR) approaches have been heavily investigated [37, 38, 103, 52, 39]. In contrast to standard MOR, these approaches learn a probability distribution and are also advantageous in terms of uncertainty quantification [18]. However, GPR models are expensive with cubic scaling, and direct application of GPR to large datasets is infeasible, instead requiring techniques such as local approximate GPR [34, 35].

On the other hand, autoencoders have also arisen as a hugely popular means of nonlinear compression and do not suffer from such poor scaling. The method was originally introduced as a nonlinear generalisation of principal component analysis [57], making the approach a generalisation of POD-based methods. Recently, it has been used in many MOR applications with great success [71, 45, 62, 88, 33, 32], which makes it a hugely promising direction that we pursue in this work.

§ MULTIFIDELITY REDUCED ORDER MODELS

By removing any dependence on the generation procedure for the training data, non-intrusive ROMs can be used as black-box approaches, allowing e.g. to augment a small experimental dataset with simulated data [31, 30, 59]. In fact, even without embedding any physics into the ROM, it is still possible to learn the dynamics of the system. Nonetheless, standard approaches typically require all training data to be generated or measured on the same grid, and with the same accuracy. Whilst often assumed for algorithmic convenience, these requirements are highly undesirable for two key reasons: firstly, generating highly-fidelity data for training is expensive and secondly, experimental data is often multifidelity data.

§.§ Generating high-fidelity data for training is expensive
Several problems require fine discretisations to obtain accurate results and resolve multiscale systems.
However, ROMs need to sample the parameter space well in order to make reasonable predictions. Thus, for standard MOR this involves the expensive computation of a large number of high-fidelity solutions. By allowing a ROM to train on data of different fidelities, one has the potential to greatly reduce the computational burden by computing only a small number of expensive (high-fidelity) simulations, and augmenting the data with a number of cheap (low-fidelity) simulations. This reflection has lead to the development of a large number of multifidelity ROMs [39, 52, 70, 82, 19, 28, 43], aiming to alleviate the computational burden of generating training data.

Low-fidelity simulated data can be generated in a number of ways, either by reducing the accuracy of the computational method e.g. using $\mathbb{P}1$ instead of $\mathbb{P}2$ finite elements, considering a coarsened discretisation, or by simplifying the physics. The central task for multifidelity MOR is to determine how to combine data arising from a number of different fidelities. A common assumption is that all fidelities are prespecified and separate models are built for each of the $F$ different fidelities. A common approach has been to adopt a multifidelity GPR based approach for MOR [82, 39, 52, 70, 19], commonly known as cokriging [53, 2], expressing each of the $F$ solution fidelities as a combination of $F$ independent Gaussian processes. Similarly, multigrid approaches have also been heavily explored for parametric PDEs [96, 41, 5, 68, 42, 43], relying on a hierarchy of models, and iteratively computing corrections through a sequence of $F$ coarser discretisations.

§.§ Experimental data is often multifidelity data
Other than synthetic data, ROMs are of utmost importance also when dealing with experimental data, e.g. coming from sensors which may move or break. This can lead to a high number of resolutions $F$, particularly for contexts such as climate modelling where data is dependent on weather stations, the number and locations of which can vary greatly [26, 25]. As a concrete example, the average number of stations used in the Global Precipitation Climatology Center [8] full-data product over Africa has varied from under 250 to over 3250 stations between 1901 and 2013 [25]. Even with a coarse yearly temporal resolution, this would still give a total of $F=112$ different resolutions.

This is particularly problematic when one considers multifidelity ROMs train a new model for each fidelity. Furthermore, in many cases it is not possible to upscale or downscale the model predictions to unseen fidelities (e.g. finer discretisation), or to leverage training data from other previously unseen fidelities either.

Thus, in the experimental context, a single model capable of training and being evaluated on complexity arbitrary discretisations is required as depicted in <ref>, i.e. there should be no dependence in the model architecture on the discretisations of the training samples. The ability to evaluate and train with data on any discretisation is known as resolution invariance [55].

None of the MOR approaches discussed thus far are resolution-invariant, underscoring the need for the development of a new, resolution-invariant ROM.

GFN-ROM as a resolution-invariant ROM, capable of handling data from any discretisation, both in training and in testing modes.

§ RELATED RESOLUTION-INVARIANT METHODS

There exist a number of resolution-invariant approaches relevant for PDEs problems. Of particular interest are PINNs, NOs, and GNNs.

§.§ Physics-informed neural networks

Standard ROM approaches usually do not require any explicit positional information, simply relying on the snapshot matrix, thereby making resolution-invariant extensions not straightforward. However, already when considering the approximation of a single snapshot, it is possible to map directly from positional information to a PDE solution i.e. $\bfmx \mapsto \phi(\bfmx) \approx u(\bfmx)$, where $\phi$ is a function to fit.

This is precisely the idea of PINNs [81, 83, 20], whereby a neural network $\phi$ is trained via a physics-informed loss i.e. the loss is taken to be the residual of the network approximation for the governing PDE, and no training data is required.  Whilst PINNs can perform well on several problems [84, 80, 85, 92, 50, 4, 22, 51], they suffer from a number of issues such as “spectral bias”, where they are incapable of learning functions with high frequencies [100], convergence discrepancies among their loss components [100] and stiff gradient flow dynamics [99]. Furthermore, their extension towards parametrised problems is in practice unsuitable due to the complexity of the training phase, rendering them not useful in the ROM context [65, 89]. For example, the POD-PINN approach trains a neural network to predict the reduced coefficients from the parameters by minimising the residual of the reduced problem, making it an intrusive approach requiring an efficient means of evaluating the residual of the reduced problem [17].

§.§ Neural operators

The inability of PINNs to learn several instances for parametric PDEs is due to the fact that such methods are not designed to directly model the solution operator. Unlike PINNs and standard machine learning methods which learn a map between evaluations of two functions, neural operators learn maps between the function spaces themselves [55]. By attempting to model explicitly the solution operator itself, neural operators remove the burdensome dependence of typical machine learning methods on the discretisation of the data, meaning they are inherently capable of achieving super- or sub-resolution [64, 86, 94, 55] and also typically benefit from desirable properties such as discretisation invariance [55]. Neural operators have shown impressive performance [86, 64, 40, 16, 12], even learning different physics at the same time [58, 74]. However, many state-of-the-art neural operators such as Fourier neural operators [64] and convolutional neural operators [86] require structured data as input. Additionally, such methods are designed for learning maps between two spatially dependent functions, and not from a global parameter to a PDE solution, as in the simplified yet more common MOR setting. With the exception of the general neural operator transformer [40], existing neural operators therefore do not take global parameters as input.

§.§ Graph neural networks

PDEs are often posed on complex domains, giving rise to unstructured meshes when computing numerical solutions. Similarly, experimental data is also rarely present on structured grids. As a result, training data for MOR is best considered as graph-based data, necessitating methods that can handle such unstructured meshes. GNNs have arisen in recent years as powerful tools to learn on graphs, and have already celebrated remarkable successes in modelling complex physics problems [60, 44, 9]. These methods aim to enforce some geometric structure to a problem by embedding geometric priors into the model. Much work (reviewed in [101, 102, 66]) has been dedicated to the the generalisation of pooling, unpooling and convolutional layers for graph-based data. GNNs are also not without limitations, with known issues for many graph convolutions such as reduced expressive power [21] and oversmoothing [48, 63, 91]. In the MOR context, one seeks a map from a global parameter to a PDE solution, which can be represented as a graph. In GNN parlance, this represents a graph unpooling operation. However, currently there exist very few unpooling methods. In the major GNN libraries in Python, namely TF-GNN [27], PyG [29] and DGL [98], there currently exists only one single graph unpooling layer, the $k$-nn interpolation layer, which is non-trainable and therefore clearly unsuitable for MOR.

§ GRAPH FEEDFORWARD NETWORK (GFN)

As we discussed before, machine learning approaches for MOR usually exploit neural network layers which are not suitable for dealing with multifidelity data. Specifically, in the GNN context there is a lack of suitable unpooling layers for mapping from an input vector to an output graph, allowing for multi-resolution approaches. Given that feedforward networks have been heavily used in MOR, regardless of their potential inconsistencies, in this work we present GFNs, a generalisation of feedforward networks capable of leveraging data on any discretisation.

Prior to presenting the method, we fix some notation and conventions for ease of reading, also providing their graphical illustrations in <ref>.

* A mesh $\mathcal{M}$ is an ordered collection of nodes (co-ordinates), without any duplicates[We represent $\mathcal{M}$ as a set, whose ordering is arbitrary and not important here, provided it is fixed for a given set.].

* $i_{\inmeshlg{}}$ denotes an index associated to a node in $\mathcal{M}$, defined by the ordering of the mesh $\mathcal{M}$. Thus, the position of a node of index $i_{\inmeshlg{}}$ is given as $\mathcal{M}\left[ i_{\inmeshlg{}} \right]$.

* $i_{\inmeshlg{o}} \leftarrow j_{\inmeshlg{n}}$ means that the node of index $i_{\inmeshlg{o}}$ in the mesh $\mathcal{M}_o$ is the nearest neighbour of the node of index $j_{\inmeshlg{n}}$ in the mesh $\mathcal{M}_n$, i.e., $i_{\inmeshlg{o}} = \operatorname{argmin}_{k_{\inmesh{o}}} \left|{\mathcal{M}_{n}} \left[ j_{\inmeshlg{n}} \right] - {\mathcal{M}_{o}}\left[ k_{\inmeshlg{o}} \right] \right|$.

* $i_{\inmeshlg{o}} \rightarrow j_{\inmeshlg{n}}$ means that the node of index $j_{\inmeshlg{n}}$ in the mesh $\mathcal{M}_n$ is the nearest neighbour[Note that we could write $i_{\inmeshlg{o}} \rightarrow j_{\inmeshlg{n}}$ as $j_{\inmeshlg{n}} \leftarrow i_{\inmeshlg{o}}$  and likewise $i_{\inmeshlg{o}} \leftarrow j_{\inmeshlg{n}}$ as $j_{\inmeshlg{n}} \rightarrow i_{\inmeshlg{o}}$. However, for readability we adopt the convention that nodes pertaining to $\mathcal{M}_o$ are given on the left side of the expression.] of the node of index $i_{\inmeshlg{o}}$ in the mesh $\mathcal{M}_o$, i.e., $j_{\inmeshlg{n}} = \operatorname{argmin}_{k_{\inmesh{n}}} \left|{\mathcal{M}_{n}} \left[ k_{\inmeshlg{n}} \right] - {\mathcal{M}_{o}}\left[ i_{\inmeshlg{o}} \right] \right|$.

* $i_{\inmeshlg{o}} \eitherarrow j_{\inmeshlg{n}}$ means that either $i_{\inmeshlg{o}} \leftarrow j_{\inmeshlg{n}}$ or $i_{\inmeshlg{o}} \rightarrow j_{\inmeshlg{n}}$.

* $i_{\inmeshlg{o}} \longleftrightarrow j_{\inmeshlg{n}}$ means that both $i_{\inmeshlg{o}} \leftarrow j_{\inmeshlg{n}}$ and $i_{\inmeshlg{o}} \rightarrow j_{\inmeshlg{n}}$.

Illustration of the arrow notations used in this work.

§.§ Feedforward networks for multifidelity data

Numerous state-of-the-art ROMs such as GCA-ROM [77] and DL-ROM [33] employ autoencoder-based approaches for MOR.
In its simplest possible form, a standard feedforward network for an autoencoder task on an unstructured mesh $\mathcal{M}_o$ with $|\mathcal{M}_o|=N_o$ nodes consists of a single-layer encoder and single-layer decoder defined as
enc(_o)_i    =    σ( ∑_j_o=1^N_o W^e_ij_o u_j_o + b^e_i ),     ∀i=1,⋯,L,

dec()_i_o    =    ∑_j=1^L W^d_i_oj z_j + b^d_i_o,     ∀i_o=1,⋯,N_o,
where $\bfmu_{\inmeshlg{o}} \in \mathbb{R}^{N_o}$ is an input sample on the mesh $\mathcal{M}_o$, $\bfmz \in \mathbb{R}^L$ is a latent vector, $\sigma$ is an activation function, $\bfmW^e \in \mathbb{R}^{L \times N_o}$ and $\bfmb^e \in \mathbb{R}^L$ are the encoder weights and biases, respectively, and $\bfmW^d \in \mathbb{R}^{N_o \times L}$ and $\bfmb^d \in \mathbb{R}^{N_o}$ are the decoder weights and biases, respectively. The reconstruction loss to optimise is simply given as $|| \operatorname{dec}\left(\operatorname{enc}(\bfmu_{\inmeshlg{o}}) \right) - \bfmu_{\inmeshlg{o}}||_2^2$, thus defining the machine learning task. We note that the encoder and decoder weights $W^e_{ji_{\inmesh{o}}}$, $W^d_{i_{\inmesh{o}}j}$ and the decoder bias $b^d_{i_{\inmesh{o}}}$ are all associated to a node $i_{\inmeshlg{o}}$. Consequently, we can regard these weights and biases as belonging to a node in the mesh $\mathcal{M}_o$, meaning that each has a spatial location given by the node's position. An illustration of this is shown for an autoencoder in <ref>.

Single-layer graph feedforward autoencoder architecture. As shown via shading, the columns of the encoder weight matrix and the rows of the decoder weight matrix are associated to the individual nodes in the mesh $\mathcal{M}_o$.

Now, we suppose that we have a trained autoencoder defined above performing well on the unstructured mesh $\mathcal{M}_o$, and we wish to use this autoencoder on a new similar mesh $\mathcal{M}_n$ with $|\mathcal{M}_n|=N_n$ nodes. This is not possible for standard feedforward network, even if it is intuitively clear that one should be able to leverage the information from the previous model. As a result, we introduce the notion of a GFN and use the nearest neighbour information between the meshes $\mathcal{M}_o$ and $\mathcal{M}_n$ to generalise the autoencoder, with provable guarantees and desirable properties. Specifically, GFNs create a means of transforming weights and biases of the encoder and decoder $(\bfmW^e, ~\bfmb^e, ~\bfmW^d \text{ and } \bfmW^d)$ for the training mesh $\mathcal{M}_o$ to new weights and biases $( \tilde{\bfmW}^e, ~\tilde{\bfmb}^e, ~\tilde{\bfmW}^d \text{ and } \tilde{\bfmW}^d )$ for the new mesh $\mathcal{M}_n$, allowing for the evaluation of a feedforward network on multifidelity data.
To achieve this transform from $\mathcal{M}_o$ to $\mathcal{M}_n$, which we denote coherently as $\mathcal{M}_o \rightarrow \mathcal{M}_n$, while retaining good performance, one can leverage the similarity of proximate nodes as illustrated in <ref>.

GFN approach to transform weights between meshes. Light blue shades denote the original mesh with its associated $\bfmW^e$ and $\bfmW^d$. The arrows depict how the new weight matrices (in darker blue shades) $\tilde{\bfmW}^e$ and $\tilde{\bfmW}^d$ can be created from the original weights and used for prediction on the new mesh.

Mathematically, the GFN transforms for $\mathcal{M}_o \rightarrow \mathcal{M}_n$ are defined as

\begin{equation}
\begin{aligned}
\label{eq:gfn}
\tilde{W}^e_{ij_{\inmesh{n}}} &= \sum_{\forall k_{\inmesh{o}} \text{ s.t } k_{\inmesh{o}} \eitherarrow j_{\inmesh{n}}} \frac{W^e_{ik_{\inmesh{o}}}}{\lvert \{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} ~\eitherarrow~ h_{\inmeshlg{n}} \}\rvert}, \\
\tilde{b}^e_{i} &= b^e_{i}, \\
\tilde{W}^d_{i_{\inmesh{n}}j} &= \underset{\forall k_{\inmesh{o}} \text{ s.t } k_{\inmesh{o}} \eitherarrow i_{\inmesh{n}}}{\operatorname{mean}} {W}^d_{k_{\inmesh{o}}j}, \\
\tilde{b}^d_{i_{\inmesh{n}}} &=  \underset{\forall k_{\inmesh{o}} \text{ s.t } k_{\inmesh{o}} \eitherarrow i_{\inmesh{n}}}{\operatorname{mean}} {b}^d_{k_{\inmesh{o}}}.
\end{aligned}
\end{equation}

Note that trivially $\tilde{b}^e_i = b^e_i$ since the encoder bias is not associated with the graph, but rather with the latent vector. These transforms require nearest neighbour computations for all nodes in $\mathcal{M}_o$ and $\mathcal{M}_n$, leading to algorithm with time complexity of $O(N_o \log N_n + N_n \log N_o)$, using the $k$-d tree method for nearest neighbour computations [11]. The predictions on the new mesh $\mathcal{M}_n$ then read analogously to before as

enc(_n)_i    =    σ( ∑_j_n=1^N_n W̃^e_ij_n u_j_n + b^e_i ),     ∀i=1,⋯,L,

dec()_i_n    =    ∑_j=1^L W̃^d_i_nj z_j + b̃^d_i_n,     ∀i_n=1,⋯,N_n,
where $\bfmu_{\inmeshlg{n}} \in \mathbb{R}^{N_n}$ is now a new input sample on the mesh $\mathcal{M}_n$ and the new weights and biases $\tilde{\bfmW}^e \in \mathbb{R}^{L \times N_n}$, $\tilde{\bfmb}^e \in \mathbb{R}^L$, $\tilde{\bfmW}^d \in \mathbb{R}^{N_n \times L}$ and $\tilde{\bfmb}^d \in \mathbb{R}^{N_n}$ are given in <ref>. Note that the transformed weights are now associated to each of the nodes on the new mesh $\mathcal{M}_n$.

To keep track of meshes and transforms, we denote the presented GFN approach of transforming weights from a mesh $\mathcal{M}_o$ to a new mesh $\mathcal{M}_n$ as

^e, ^d, ^d = GFN^ℳ_o →ℳ_n(^e, ^d, ^d).
Similarly, we also write $\operatorname{enc}^{\mathcal{M}_o \to \mathcal{M}_n}$ and $\operatorname{dec}^{\mathcal{M}_o \to \mathcal{M}_n}$ to track the applied transforms.

For autoencoders with multiple hidden layers, the GFN approach trivially generalises and only needs to be applied to the first encoder layer and the last decoder layer, since these are the only layers directly interacting with the changing meshes (for further details see [app:mult_layers]Appendix <ref>app:mult_layers).

This approach of transforming via GFN thus defines a means of allowing an autoencoder, which has been trained on single fidelity data using standard feedforward networks only, to be evaluated on multifidelity data.

§.§ Algorithmic details

In this section, we highlight two notable instances of the possible transforms $\mathcal{M}_o \to \mathcal{M}_n$ which lead to computational savings: the expansive and the agglomerative cases. Demonstrative examples of agglomerative, expansive, neither agglomerative nor expansive and both agglomerative and expansive cases are shown in <ref>.

Examples of the different types of possible transforms.

§.§.§ Expansive case

The nearest neighbour of each node $i_{\inmeshlg{o}}$ in $\mathcal{M}_n$ has $i_{\inmeshlg{o}}$ as its nearest neighbour in $\mathcal{M}_o$, i.e., the nearest neighbour of the nearest neighbour of each node $i_{\inmeshlg{o}}$ is itself, following our notation:

∀i_o,   i_o →j_n i_o ←j_n.

In such case, the transform can be simplified (see the proof in [subsec:exp_proof]Appendix <ref>subsec:exp_proof) as:

\begin{equation}
\label{eq:gfn_exp}
\begin{aligned}
\tilde{W}^e_{ij_{\inmesh{n}}} &= \frac{1}{|\{h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow h_{\inmeshlg{n}} \}|} {W}^e_{ik_{\inmesh{o}}}, \qquad \text{ where}\quad k_{\inmeshlg{o}} \leftarrow j_{\inmeshlg{n}}, \\
\tilde{W}^d_{i_{\inmesh{n}}j} &= {W}^d_{k_{\inmesh{o}}j}, \qquad \text{ where}\quad k_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{n}}, \\
\tilde{b}^d_{i_{\inmesh{n}}} &= {b}^d_{k_{\inmesh{o}}},  \qquad \text{ where}\quad k_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{n}}.
\end{aligned}
\end{equation}

We note that the absence of rightward arrows means we are only required to find nearest neighbours of nodes in the mesh $\mathcal{M}_n$ and not those in $\mathcal{M}_o$, reducing the algorithmic time complexity to $O(N_n \log N_o)$. The implementation of the expansive transformation with this complexity is shown in [alg:gfn_expand]Algorithm <ref>alg:gfn_expand.

Simplified GFN transform in the expansive case, as given in <ref>.

Expand$\bfmW^e \in \mathbb{R}^{L \times N_o}, \bfmW^d \in \mathbb{R}^{N_o \times L}, \bfmb^d \in \mathbb{R}^{N_o}, \mathcal{M}_o, \mathcal{M}_n$
Initialise $\tilde{\bfmW}^e \in \mathbb{R}^{L \times N_n}, \tilde{\bfmW}^d \in \mathbb{R}^{N_n \times L}, \tilde{\bfmb}^d \in \mathbb{R}^{N_n}$
$\operatorname{counts} \gets \boldsymbol{0} \in \mathbb{R}^{N_o}$
$i_{\inmeshlg{n}}=1 \text{ to } N_n$
$j_{\inmeshlg{o}} \gets \operatorname{nearest-neighbour}^{\inmeshlg{o}}(\mathcal{M}_n \left[ i_{\inmeshlg{n}} \right])$
$\operatorname{counts}[j_{\inmeshlg{o}}] \gets \operatorname{counts}[j_{\inmeshlg{o}}]+1$
$i_{\inmeshlg{n}}=1 \text{ to } N_n$
$j_{\inmeshlg{o}} \gets \operatorname{nearest-neighbour}^{\inmeshlg{o}}(\mathcal{M}_n \left[ i_{\inmeshlg{n}} \right])$
$\tilde{\bfmW}^e \left[:,~i_{\inmeshlg{n}}\right] \gets \bfmW^e \left[:,~j_{\inmeshlg{o}}\right] \mathbin{/} \operatorname{counts}[j_{\inmeshlg{o}}]$
$\tilde{\bfmW}^d \left[i_{\inmeshlg{n}},~:\right] \gets {\bfmW}^d \left[j_{\inmeshlg{o}},~:\right]$
$\tilde{\bfmb}^d \left[ i_{\inmeshlg{n}}\right] \gets{\bfmb}^d \left[ j_{\inmeshlg{o}}\right]$
$\tilde{\bfmW}^e, \tilde{\bfmW}^d, \tilde{\bfmb}^d$

§.§.§ Agglomerative case

The nearest neighbour of each node $j_{\inmeshlg{n}}$ in $\mathcal{M}_o$ has $j_{\inmeshlg{n}}$ as its nearest neighbour in $\mathcal{M}_n$, i.e., the nearest neighbour of the nearest neighbour of each node $j_{\inmeshlg{n}}$ is itself, following our notation:

∀j_n,   i_o ←j_n i_o →j_n.
In such case, the transform can be simplified (see the proof in [subsec:agg_proof]Appendix <ref>subsec:agg_proof) as:

\begin{equation}
\label{eq:gfn_agg}
\begin{aligned}
\tilde{W}^e_{ij_{\inmesh{n}}} &=& \sum_{\forall k_{\inmesh{o}} \text{ s.t. }  k_{\inmesh{o}} \rightarrow j_{\inmesh{n}}} W^e_{ik_{\inmesh{o}}}, \\
\tilde{W}^d_{i_{\inmesh{n}}j} &=& \underset{\forall k_{\inmesh{o}} \text{ s.t. }  k_{\inmesh{o}} \rightarrow i_{\inmesh{n}}}{\operatorname{mean}} W^d_{k_{\inmesh{o}}j}, \\
\tilde{b}^d_{i_{\inmesh{n}}} &=&  \underset{\forall k_{\inmesh{o}} \text{ s.t. }  k_{\inmesh{o}} \rightarrow i_{\inmesh{n}}}{\operatorname{mean}} b^d_{k_{\inmesh{o}}}.
\end{aligned}
\end{equation}

We note that the absence of leftward arrows means we are only required to find nearest neighbours of nodes in the mesh $\mathcal{M}_o$ and not those in $\mathcal{M}_n$, reducing the algorithmic time complexity to $O(N_o \log N_n)$. The implementation of the expansive transformation with this complexity is shown in [alg:gfn_agglomerate]Algorithm <ref>alg:gfn_agglomerate.

Simplified GFN transform in the agglomerative case, as given in <ref>.

Agglomerate$\bfmW^e \in \mathbb{R}^{L \times N_o}, \bfmW^d \in \mathbb{R}^{N_o \times L}, \bfmb^d \in \mathbb{R}^{N_o}, \mathcal{M}_o, \mathcal{M}_n$
Initialise $\tilde{\bfmW}^e \in \mathbb{R}^{L \times N_n}, \tilde{\bfmW}^d \in \mathbb{R}^{N_n \times L}, \tilde{\bfmb}^d \in \mathbb{R}^{N_n}$
$\operatorname{counts} \gets \boldsymbol{0} \in \mathbb{R}^{N_n}$
$i_{\inmeshlg{o}}=1 \text{ to } N_o$
$j_{\inmeshlg{n}} \gets \operatorname{nearest-neighbour}^{\inmeshlg{n}}(\mathcal{M}_o \left[ i_{\inmeshlg{o}} \right])$
$\operatorname{counts}[j_{\inmeshlg{n}}] \gets \operatorname{counts}[j_{\inmeshlg{n}}]+1$
$\tilde{\bfmW}^e \left[:,~j_{\inmeshlg{n}}\right] \gets \bfmW^e \left[ :,~i_{\inmeshlg{o}}\right]$
$\tilde{\bfmW}^d \left[j_{\inmeshlg{n}},~:\right] \gets \left(\left( \operatorname{counts}[j_{\inmesh{n}}] - 1\right) \tilde{\boldsymbol{W}}^d \left[j_{\inmeshlg{n}},~:\right] + {\boldsymbol{W}}^d \left[i_{\inmeshlg{o}},~:\right]\right) \mathbin{/} \operatorname{counts}[j_{\inmeshlg{n}}]$
$\tilde{\bfmb}^d \left[ j_{\inmeshlg{n}}\right] \gets \left(\left( \operatorname{counts}[j_{\inmeshlg{n}}] - 1\right) \tilde{\boldsymbol{b}}^d \left[ j_{\inmeshlg{n}}\right] + {\boldsymbol{b}}^d \left[ i_{\inmeshlg{o}}\right]\right) \mathbin{/} \operatorname{counts}[j_{\inmeshlg{n}}]$
$\tilde{\bfmW}^e, \tilde{\bfmW}^d, \tilde{\bfmb}^d$

§.§.§ Computational advantages

For hierarchical meshes, it follows that if $\mathcal{M}_o \subseteq \mathcal{M}_n$ then the transform $\mathcal{M}_o \rightarrow \mathcal{M}_n$ is expansive, whereas if $\mathcal{M}_n \subseteq \mathcal{M}_o$ then the transform $\mathcal{M}_o \rightarrow \mathcal{M}_n$ is agglomerative. Moreover, if the transform is both expansive and agglomerative then it must hold that $\forall i_{\inmeshlg{o}}$, $i_{\inmeshlg{o}} \longleftrightarrow j_{\inmeshlg{n}}$ and $\forall j_{\inmeshlg{n}}$, $i_{\inmeshlg{o}} \longleftrightarrow j_{\inmeshlg{n}}$, i.e. there is a one-to-one correspondence between nodes in the original mesh $\mathcal{M}_o$ and in the new mesh $\mathcal{M}_n$, and the GFN transforms reduce to nearest neighbour interpolation of the weights.

However, for general meshes a transform can be both, neither or one of expansive and agglomerative. Despite this, notably any general GFN transform can be expressed as a succession of an expansive update and an agglomerative update, (proof in [app:gfn_decomp]Appendix <ref>app:gfn_decomp) i.e. <ref> is equivalent to defining the new mesh

ℳ_ = ℳ_o ∪{ℳ_n [i_n]  s.t.  j_o ←i_n  but not  j_o →i_n },
and computing the transform in two steps, as

\begin{equation}
\begin{rcases}
\hat{W}^e_{ij_{\inmesh{}}} &= \frac{1}{|\{h_{\inmesh{}} \text{ s.t. } k_{\inmesh{o}} \leftarrow h_{\inmesh{}} \}|} {W}^e_{ik_{\inmesh{o}}}, \qquad \text{ where}\quad k_{\inmeshlg{o}} \leftarrow j_{\inmeshlg{}}, \\
\hat{W}^d_{i_{\inmesh{}}j} &= {W}^d_{k_{\inmesh{o}}j}, \qquad \text{ where}\quad k_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{}}, \\
{{\hat{b}}^d}_{i_{\inmesh{}}} &= {b}^d_{k_{\inmesh{o}}},  \qquad \text{ where}\quad k_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{}},
\end{rcases} \quad \text{Expansion step}
\end{equation}

followed by

\begin{equation}
\begin{rcases}
\tilde{W}^e_{ij_{\inmesh{n}}} &= \sum_{\forall k_{\inmesh{}} \text{ s.t. }  k_{\inmesh{}} \rightarrow j_{\inmesh{n}}} {\hat{W}}^e_{ik_{\inmesh{}}}, \\
\tilde{W}^d_{i_{\inmesh{n}}j} &= \underset{\forall k_{\inmesh{}} \text{ s.t. }  k_{\inmesh{}} \rightarrow i_{\inmesh{n}}}{\operatorname{mean}} {\hat{W}}^d_{k_{\inmesh{}}j}, \\
\tilde{b}^d_{i_{\inmesh{n}}} &=  \underset{\forall k_{\inmesh{}} \text{ s.t. }  k_{\inmesh{}} \rightarrow i_{\inmesh{n}}}{\operatorname{mean}} {\hat{b}}^d_{k_{\inmesh{}}}.
\end{rcases} \quad \text{Agglomerative step}
\end{equation}

The general GFN algorithm can thus be easily expressed as in [alg:gfn_general]Algorithm <ref>alg:gfn_general, with a demonstrative illustration of this interpretation in <ref>.

General GFN transform leveraging an expansion step followed by an agglomerative update, as given in <ref>.

GFN$\bfmW^e \in \mathbb{R}^{L \times N_o}, \bfmW^d \in \mathbb{R}^{N_n \times L}, \bfmb^d \in \mathbb{R}^{N_o}, \mathcal{M}_o, \mathcal{M}_n$
$\mathcal{M} = \mathcal{M}_o \cup \{\mathcal{M}_n \left[i_{\inmeshlg{n}}\right] \text{ s.t. } j_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{n}} \text{ but not } j_{\inmeshlg{o}} \rightarrow i_{\inmeshlg{n}} \}$
$\hat{\bfmW}^e, \hat{\bfmW}^d, \hat{\bfmb}^d \gets \Call{Expand}{\bfmW^e, \bfmW^d, \bfmb^d, \mathcal{M}_o, \mathcal{M}}$
$\tilde{\bfmW}^e, \tilde{\bfmW}^d, \tilde{\bfmb}^d \gets \Call{Agglomerate}{\hat{\bfmW}^e, \hat{\bfmW}^d, \hat{\bfmb}^d, \mathcal{M}, \mathcal{M}_n}$
$\tilde{\bfmW}^e, \tilde{\bfmW}^d, \tilde{\bfmb}^d$

The GFN transform can be thought about as applying an expansive update followed by an agglomerative update.

§ GRAPH FEEDFORWARD NETWORK REDUCED ORDER MODEL (GFN-ROM)

Having defined GFNs as a means of evaluating trained autoencoders on multifidelity data, we leverage them to propose a reduced order model capable of being evaluated on arbitrary meshes.
As depicted in <ref>, the proposed GFN-ROM architecture consists of two components: an autoencoder and a mapper. The autoencoder is based on GFNs, while the mapper can be any neural network model which maps the fixed-size vector parameters directly to the fixed-size latent representation. Notice that, within this parametric context, the mapper does not need to be defined using GFNs since it does not deal with graph-structured data.

GFN-ROM is inspired by recent state-of-the-art autoencoder-based models such as DL-ROM [33], which uses standard convolutions, and GCA-ROM [77], which uses graph convolutions. In fact, GFN-ROM can be regarded as an extension of such methods. Like both models, we train the GFN-ROM model based on a weighted sum of two losses: the reconstruction loss and the mapper loss.
Given a GFN-ROM model with weights defined on a mesh $\mathcal{M}_o$, the per-example loss functions are defined for the $t$th training sample $\left( \boldsymbol{\mu}_t, ~\bfmu_{\inmeshlgtwo{n}{t}}(\boldsymbol{\mu}_t) \right)$, defined on the $t$-th mesh $\mathcal{M}_n^t$ and for the $t$-th parameter value $\boldsymbol{\mu}_t$, as
ℒ_recon[μ_t,_nt(μ_t)]    =    1/|ℳ_n^t| || dec^ℳ_o →ℳ_n^t(enc^ℳ_o →ℳ_n^t(_nt(μ_t)) ) - _nt(μ_t) ||_2^2,

ℒ_map[μ_t, _nt(μ_t)]    =    1/L || enc^ℳ_o →ℳ_n^t(_nt(μ_t)) - map(μ_t)||_2^2.
The overall loss is given as a weighted mean of these losses over a set $\{ \left( \boldsymbol{\mu}_t, ~ \bfmu_{\inmeshlgtwo{n}{t}}(\boldsymbol{\mu}_t) \right) \}_{t=1}^{T}$ of $T$ training examples
ℒ = 1/T ∑_t=1^T |ℳ_n^t|/∑_s=1^T |ℳ_n^s| ( ℒ_recon[μ_t,_nt(μ_t)] + ωℒ_map[μ_t,_nt(μ_t)] ),
where we have introduced the mapper weight hyperparameter $\omega$.
We note that at inference, the encoder is no longer required and, for a parameter $\boldsymbol{\mu}_t$ and mesh $\mathcal{M}_n^t$, a GFN-ROM model with weights defined on a mesh $\mathcal{M}_o$ predicts $\operatorname{dec}^{\mathcal{M}_o \to \mathcal{M}_n^t}(\operatorname{map}(\boldsymbol{\mu}_t))$.

GFN-ROM architecture, consisting of an autoencoder and a mapper. During inference, predictions can be computed by decoding the mapper output.

§.§ Main features of the architecture

The proposed GFN-ROM architecture benefits from a number of desirable properties that we discuss in the following sections, which make it an attractive tool for applications with multifidelity data.

§.§.§ Self-consistency
The GFN transforms are self-consistent approaches, meaning that a series of two transforms from $\mathcal{M}_o$ to $\mathcal{M}_n$ and then back to $\mathcal{M}_o$ recovers exactly the original weights $\bfmW^e, \bfmW^d, \bfmb^d$ provided that the transform $\mathcal{M}_o$ to $\mathcal{M}_n$ is expansive, i.e.
$$\bfmW^e, \bfmW^d, \bfmb^d = \operatorname{GFN}^{\mathcal{M}_n \to \mathcal{M}_o}( \operatorname{GFN}^{\mathcal{M}_o \to \mathcal{M}_n}(\bfmW^e, \bfmW^d, \bfmb^d) ).$$

The self-consistency, a proof of which is given in [app:self_consistency]Appendix <ref>app:self_consistency, is a marked advantage in comparison to any other interpolation methods such as Kriging, linear regression, spline, and $k$NN (with the exception of the case $k=1$). This means that no information is lost during the expansive transform and the original weights can be directly recovered, unlike for other interpolation-based approaches.

§.§.§ Bound on super- and sub-resolution error of GFN-ROM

GFN-ROM is a resolution-invariant architecture, meaning one is often interested in evaluating the model on a finer discretisation (super-resolution) or coarser discretisation (sub-resolution) compared to a reference resolution. A particularly attractive feature of the architecture is that one can bound the performance on any testing fidelity given certain information. Specifically, given a bound $\delta$ on the difference between the features of nearest neighbours, and a bound $\tau$ on the performance of GFN-ROM on the known mesh $\mathcal{M}_o$, i.e.[We use the notation $u(x_{i_{\inmesh{o}}})$, equivalent to $u_{i_{\inmesh{o}}}$, in order to emphasise the dependence of $u$ on spatial position.]
∀i_o,   i_o j_n | u(x_i_o ) - u(x_j_n) | ≤δ,

∀i_o,  | u(x_i_o) - dec^ℳ_o →ℳ_o(map(μ))_i_o | ≤τ,
one can show that the performance of GFN-ROM on the new mesh $\mathcal{M}_n$ can be bounded as:
∀i_n,  | u(x_i_n) - dec^ℳ_o →ℳ_n(map(μ))_i_n | ≤τ+ δ.
A proof is given in [subsec:error_gfnrom]Appendix <ref>subsec:error_gfnrom. We note that there is no dependence at all on the weights and biases, only on the mesh similarities and the GFN-ROM performance on the reference resolution.

§.§.§ Bound on super- and sub-resolution error of the GFN-ROM mapper

Given the same bound $\delta$ as in <ref>, and a bound $\alpha$ on the performance of the mapper on the known mesh $\mathcal{M}_o$, i.e.∀i,  | map(μ)_i - enc^ℳ_o →ℳ_o(_o)_i | ≤α,
one can show that the performance of the mapper on the new mesh $\mathcal{M}_n$ can be bounded as:
∀i,  | map(μ)_i - enc^ℳ_o →ℳ_n(_n)_i | ≤α+ δC^P+1 ||^e||_∞.
A proof is given in [subsec:error_mapper]Appendix <ref>subsec:error_mapper.

§.§.§ Bound on super- and sub-resolution error of the GFN-ROM autoencoder

Given the same bound $\delta$ as in <ref>, and a bound $\beta$ on the performance of the autoencoder on the known mesh $\mathcal{M}_o$, i.e.∀i_o, | u(x_i_o) - dec^ℳ_o →ℳ_o(enc^ℳ_o →ℳ_o(_o))_i_o | ≤β,
one can show that the performance of the autoencoder on the new mesh $\mathcal{M}_n$ can be bounded as:

\begin{equation}
\begin{aligned}
\forall i_{\inmeshlg{n}},\quad \left| u(x_{i_{\inmesh{n}}}) - \operatorname{dec}^{\mathcal{M}_o \to \mathcal{M}_n}(\operatorname{enc}^{\mathcal{M}_o \to \mathcal{M}_n}(\bfmu_{\inmeshlg{n}}))_{i_{\inmesh{n}}} \right| \leq \beta + \delta +  \delta C^{Q+1} ||\bfmW^{d}||_\infty ||\bfmW^{e}||_\infty.
\end{aligned}
\end{equation}

A proof is given in [subsec:error_autoenc]Appendix <ref>subsec:error_autoenc. The dependence on the norm of the weights in the bound suggests that regularisation is important to allow for good generalisation for the autoencoder.

§.§ Link between GFN-ROM and novel architectures

Following the discussion in [subsec:nos]Section <ref>subsec:nos and [subsec:gnns]Section <ref>subsec:gnns, we discuss here how GFN-ROM is connected with state-of-the-art resolution invariant approaches.

§.§.§ GNNs and GFN-ROM

Concerning graph neural networks, we can interpret a GFN encoder layer and a GFN decoder layer as two different types of GNN: namely, a graph pooling layer and a graph unpooling layer, respectively. We can consider the input to the autoencoder, which is some PDE solution $u(x, \boldsymbol{\mu})$ evaluated on some discretized domain, as a fully connected graph. The GFN encoder layer is then tasked with computing a vector representation summarising the input graph, i.e. a graph pooling task. The GFN decoder layer has precisely the opposite task of reconstructing a graph given a summarised representation i.e. graph unpooling. At inference, the overall GFN-ROM architecture can thus be thought of as a graph unpooling model, where the PDE solution is identified by the parameters $\boldsymbol{\mu}$.

§.§.§ NOs and GFN-ROM

A parameterised PDE with a parameter $a \in \mathcal{A}$ and solution $u \in \mathcal{U}$, where $\mathcal{A}$ and $\mathcal{U}$ are Banach spaces, defines a mapping $\mathcal{G}^\dag: \mathcal{A} \to \mathcal{U}$, i.e. $u(x)=\mathcal{G}^\dag[a](x)$. A neural operator aims to directly approximate this mapping $\mathcal{G}^\dag$ by means of an operator $\mathcal{G}_\theta \approx \mathcal{G}^\dag$ with trainable parameters $\theta$.
A single-layer neural operator can be written as
G_θ= 𝒬 ∘σ(W+𝒦 + b) ∘𝒫,
where $\mathcal{P}$, $\mathcal{Q}$ are the local lifting and projection mappings, $W$ is a local linear operator (matrix), $\mathcal{K}$ is an integral kernel operator, $b$ is a bias function and $\sigma$ is a pointwise activation function. We refer to Kovachki et al. for further details [55]. Considering a simplified neural operator with $W=0$ and no lifting or projection mappings, we obtain

G_θ[v](y) = σ( ∫_D κ(x,y) v(x)  dx + b(y) ),
where $\kappa$ is a kernel function and $D$ is the domain of integration. Usually, one does not have complete access to the function $v(x)$, but only to its evaluations at certain points, meaning this integral has to be numerically computed. Supposing we are given function evaluations on a mesh $\mathcal{M}$, the numerical approximation of the integral reads as
G_θ[v](y) ≈σ(∑_i_ h_i_ κ(x_i_,y) v(x_i_) + b(y) ),
where $h_{i_{\inmesh{}}}$ is a factor depending on the integration method.

GFN-ROM is closely related to NOs. In fact, in the expansive case, GFN-ROM can be recast as a particular instance of NOs, since we can directly obtain both the encoder and decoder. To make the link with GFN-ROM, we start with the mesh $\mathcal{M}_o$, the associated weights and biases ($\bfmW^e$, $\bfmb^e$, $\bfmW^d$, $\bfmb^d$), and choose a particular kernel $\kappa$.

Starting from <ref>, discarding the $y$-dependence and taking $D=\Omega$, $b=b^e_i$, and $\kappa(x) = W^e_{ik_{\inmesh{o}}}/|B_{k_{\inmesh{o}}}|$ where $k_{\inmeshlg{o}} \leftarrow x$, and $B_{k_{\inmesh{o}}} = \{x \in \Omega\ |\ k_{\inmeshlg{o}} \leftarrow x \}$, one can write[We abuse notation to write $k_{\inmeshlg{o}} \leftarrow x$ since $x$ is a position and not an index, but the meaning is analogous to before, i.e. $k_{\inmeshlg{o}} = \operatorname{argmin}_{h_{\inmesh{o}}} \left|{\mathcal{M}_{o}}\left[ h_{\inmeshlg{o}} \right]  - x \right|$.]

\begin{equation}
\begin{aligned}
G_\theta[u] &= \sigma \left( \int_D \kappa(x) u(x) \, \mathrm{d}x + b \right),\\
&= \sigma \left( \int_\Omega \frac{W^e_{ik_{\inmesh{o}}}}{|B_{k_{\inmesh{o}}}|} u(x) \, \mathrm{d}x  + b^e_i \right), \qquad k_{\inmeshlg{o}} \leftarrow x,\\
&= \sigma \left( \sum_{k_{\inmesh{o}}=1}^{N_o} W^e_{ik_{\inmesh{o}}} \frac{1}{|B_{k_{\inmesh{o}}}|}\int_{B_{k_{\inmesh{o}}}} u(x) \, \mathrm{d}x + b^e_i \right).
\end{aligned}
\end{equation}

Approximating the integral with observations of $u(x)$ on a mesh $\mathcal{M}_n$, with $k_{\inmeshlg{o}} \leftarrow j_{\inmeshlg{n}}$, gives the formula for the GFN-ROM encoder in the expansive case

\begin{equation}
\begin{aligned}
G_\theta[u] &= \sigma \left( \sum_{k_{\inmesh{o}}=1}^{N_o} W^e_{ik_{\inmesh{o}}} ~\underset{\forall j_{\inmesh{n}} \text{ s.t } k_{\inmesh{o}} \leftarrow j_{\inmesh{n}}}{\operatorname{mean}} u_{j_{\inmesh{n}}} + b^e_i \right), \\
&= \sigma \left( \sum_{j_{\inmesh{n}}=1}^{N_n} \frac{W^e_{ik_{\inmesh{o}}}}{\lvert \{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow h_{\inmeshlg{n}} \}\rvert}  u_{j_{\inmesh{n}}} + b^e_i \right),  \qquad \text{ where}\quad k_{\inmeshlg{o}} \leftarrow j_{\inmeshlg{n}},
\\
&= \operatorname{enc}^{\mathcal{M}_o \to \mathcal{M}_n}(\bfmu_{\inmeshlg{n}})_i.
\end{aligned}
\end{equation}

Similarly, the decoder can be obtained from <ref> without any $x$-dependence.
Taking $\tilde{\kappa} = |D| \kappa$, setting $\sigma$ to the identity function, $\tilde{\kappa}(y)=\bfmW^d_{k_{\inmesh{o}},:}$ where $k_{\inmeshlg{o}} \leftarrow y$, and $b(y)=b^d_{k_{\inmesh{o}}}$ where $k_{\inmeshlg{o}} \leftarrow y$, one obtains

\begin{equation}
\begin{aligned}
G_\theta[\bfmz](y) &= \sigma \left( \int_D \kappa(y) \bfmz \, \mathrm{d}x + b(y) \right), \\
&= \sigma \left( \tilde{\kappa}(y) \bfmz + b(y) \right),\\
&= \sum_{j=1}^L {W}^d_{k_{\inmesh{o}}j} z_j + {b}^d_{k_{\inmesh{o}}}, \qquad k_{\inmeshlg{o}} \leftarrow y.
\end{aligned}
\end{equation}

Considering $y=\mathcal{M}_n[i_{\inmeshlg{n}}]$, we recover the GFN-ROM decoder in the expansive case as

\begin{equation}
\begin{aligned}
G_\theta[\bfmz](y) &= \sum_{j=1}^L {W}^d_{k_{\inmesh{o}}j} z_j + {b}^d_{k_{\inmesh{o}}},\\
&= \operatorname{dec}^{\mathcal{M}_o \to \mathcal{M}_n}(\bfmz)_{i_{\inmesh{n}}}.
\end{aligned}
\end{equation}

Like the decoder, the GFN-ROM model at inference in the expansive case can be thought of as a neural operator, but with the mapper acting as the lifting function. Thus, in such settings, GFN-ROM encoder and decoder are very closely linked to resolution-invariant neural operators with a choice of piecewise constant kernel function.

§.§.§ ROM-related properties

In addition to the above proofs, the GFN-ROM architecture benefits from a number of typical properties important for reduced order models. Namely, it is purely data-driven and non-intrusive, not requiring any knowledge of the generation process for the training data. It is furthermore a nonlinear method, leveraging an autoencoder for nonlinear compression. The method is suitable for unstructured meshes, which is vital in order to deal with general PDEs, where the presence of complex domains typically lead to unstructured meshes. Most importantly, the method is a completely multifidelity approach and is capable of inference on arbitrary meshes.

§ ADAPTIVE MULTIFIDELITY TRAINING

The GFN method proposes a means of transferring weights from one graph $\mathcal{M}_o$ to a new graph $\mathcal{M}_n$ as $\tilde{\bfmW}^e, \tilde{\bfmW}^d, \tilde{\bfmb}^d = \operatorname{GFN}^{\mathcal{M}_o \to \mathcal{M}_n}(\bfmW^e, \bfmW^d, \bfmb^d)$, allowing one to evaluate a feedforward-based approach on new meshes. Whilst this is already a major advantage, it does not prescribe how to effectively train in a multifidelity setting.

In this section, we discuss means of training in multifidelity applications. We aim at learning from a series of $T$ meshes $\mathcal{M}_n^1,\cdots,\mathcal{M}_n^T$ by training the GFN-ROM model with weights and biases corresponding to a fixed or an adaptive mesh $\mathcal{M}_o$. We describe below the two approaches, the difference between which is illustrated in <ref>.

§.§ Fixed mesh
The GFN transforms are entirely differentiable, meaning that for weights and biases $\bfmW^e$, $\bfmW^d$ and $\bfmb^d$ associated on a fixed mesh $\mathcal{M}_o$, it is possible to compute gradients $\nabla_{\boldsymbol{W}^e} \mathcal{L}$, $\nabla_{\boldsymbol{W}^d} \mathcal{L}$ and $\nabla_{\boldsymbol{b}^d} \mathcal{L}$ when training on new meshes. As a result, it is in fact possible to use the standard gradient-based optimisation methods to learn the optimal weights $\bfmW^e$, $\bfmW^d$ and $\bfmb^d$, which are always associated to the mesh $\mathcal{M}_o$. If one considers the mesh $\mathcal{M}_o$ as being a suitable sampling of the domain, this approach therefore suffices as a means of training the model in a multifidelity setting. However, the choice of mesh $\mathcal{M}_o$ is extremely important for the performance of the model. The fixed nature of $\mathcal{M}_o$ is also undesirable since, like interpolation approaches onto fixed grids, it does not allow the model to adapt the mesh $\mathcal{M}_o$ to better match the training data it sees.

§.§ Adaptive mesh

GFN-ROM differs from interpolation onto a fixed grid since it is possible to transform weights and biases from a coarse mesh into weights and biases for a finer mesh. Thus, one can always overwrite the old weights and biases corresponding to the old mesh with new weights and biases of the new mesh, and then optimise the new weights. As a result, the approach leads to a number of advantages in allowing the architecture to adapt its number of parameters during training. That is, instead of having a single fixed mesh $\mathcal{M}_o$ associated to the fixed weights and biases for the model, we replace it by a master mesh $\mathcal{M}^t$ which can change during each training sample $t=1,\cdots,T$. Importantly, the sequence of master meshes during training must be hierarchically increasing i.e. $\mathcal{M}^0 \subseteq \cdots \subseteq \mathcal{M}^T$, since we do not wish to lose information at any point.

A simple scheme for adaptive training could be setting $\mathcal{M}^t = \mathcal{M}^{t-1} \cup \mathcal{M}_n^t$. However, this approach is suboptimal since it necessarily adds every node seen during training to the master mesh. This means that the master mesh can become prohibitively large, and result in a worse time complexity of $O((N_o+N_n)\log N_o + (N_o+N_n)\log N_n)$ for the transform $\mathcal{M}_o \to \mathcal{M}_n$.

Instead, a better approach is possible without incurring the computational issues of the previous algorithm. As discussed in <ref> and summarised in [alg:gfn_general]Algorithm <ref>alg:gfn_general, the GFN transforms can be decomposed into an expansive and agglomerative part via the choice of a suitable intermediate mesh given in <ref>. This intermediate mesh is already computed during the GFN transforms, meaning it does not represent an extra operation. Furthermore, the intermediate mesh is larger than the original mesh. Therefore, it is possible to use such mesh to develop a computationally efficient adaptive approach i.e., $\mathcal{M}^t = \mathcal{M}^{t-1} \cup  \{\mathcal{M}_n^t \left[i_{\inmeshlgtwo{n}{t}}\right] \text{ s.t. } j_{\inmeshlgtwo{}{t-1}} \leftarrow i_{\inmeshlgtwo{n}{t}} \text{ but not } j_{\inmeshlgtwo{}{t-1}} \rightarrow i_{\inmeshlgtwo{n}{t}} \}$.

Training with  fixed or adaptive mesh methods when given a starting mesh $\mathcal{M}_o$ and new mesh $\mathcal{M}_n$. The fixed mesh approach updates the original mesh weights, whilst the adaptive mesh approach updates weights on a new mesh $\mathcal{M}$, created from information from both meshes.

For a transform $\mathcal{M}_o \to \mathcal{M}_n$, we retain the time complexity of the non-adaptive algorithm. Furthermore, this approach only adds nodes to the master mesh when the master mesh undersamples[Undersampling here refers purely to nearest-neighbour information i.e., areas where a node in the master mesh has multiple nearest neighbours are considered undersampled.]. We note that because the weights and biases are changing shape along with the master mesh during training for adaptive methods, one cannot use Adam, RMSProp or any momentum-based optimisation methods for updates[In practice, we often know all of the training fidelities we shall see during training. In such cases, the final mesh $\mathcal{M}^T$ can be precomputed and the model trained as for a fixed mesh.].

The training loop for GFN-ROM is given in [alg:gfn_train]Algorithm <ref>alg:gfn_train, illustrating the difference between the adaptive and fixed mesh modes.

Training algorithm for GFN-ROM, showing both adaptive and fixed mesh approaches.

$T$ training samples of parameter and solution pairs i.e. $\{ \left( \boldsymbol{\mu}_t,~ \bfmu_{\inmeshlgtwo{n}{t}}(\boldsymbol{\mu}_t) \right) \}_{t=1}^{T}$
Initialise a mesh $\mathcal{M}_o$ with associated GFN weights ${\bfmW}^e, {\bfmW}^d, {\bfmb}^d$
Initialise $\theta$: all parameters needed for GFN-ROM excluding ${\bfmW}^e, {\bfmW}^d, {\bfmb}^d$
$t\gets 1, \dots, T$
$\mathcal{M} = \mathcal{M}_o \cup \{\mathcal{M}^t_n \left[i_{\inmeshlgtwo{n}{t}}\right] \text{ s.t. } j_{\inmeshlg{o}} \leftarrow i_{\inmeshlgtwo{n}{t}} \text{ but not } j_{\inmeshlg{o}} \rightarrow i_{\inmeshlgtwo{n}{t}} \}$
adaptive adaptive mode
${\bfmW}^e, {\bfmW}^d, {\bfmb}^d \gets \Call{Expand}{\bfmW^e, \bfmW^d, \bfmb^d, \mathcal{M}_o, \mathcal{M}}$
$\mathcal{M}_o \gets \mathcal{M}$
$\tilde{\bfmW}^e, \tilde{\bfmW}^d, \tilde{\bfmb}^d \gets \Call{Agglomerate}{{\bfmW}^e, {\bfmW}^d, {\bfmb}^d, \mathcal{M}_o, \mathcal{M}_n}$
fixed mode
$\hat{\bfmW}^e, \hat{\bfmW}^d, \hat{\bfmb}^d \gets \Call{Expand}{\bfmW^e, \bfmW^d, \bfmb^d, \mathcal{M}_o, \mathcal{M}}$
$\tilde{\bfmW}^e, \tilde{\bfmW}^d, \tilde{\bfmb}^d \gets \Call{Agglomerate}{\hat{\bfmW}^e, \hat{\bfmW}^d, \hat{\bfmb}^d, \mathcal{M}, \mathcal{M}_n}$
Compute loss $\mathcal{L}(\boldsymbol{\mu}_t, \bfmu_{\inmeshlgtwo{n}{t}}; \tilde{\bfmW}^e, \tilde{\bfmW}^d, \tilde{\bfmb}^d, \theta)$
Update parameters ${\bfmW}^e, {\bfmW}^d, {\bfmb}^d, \theta$ using gradients $\nabla_{\boldsymbol{W}^e} \mathcal{L},~ \nabla_{\boldsymbol{W}^d} \mathcal{L},~ \nabla_{\boldsymbol{b}^d} \mathcal{L},~ \nabla_{\theta} \mathcal{L}$

§ NUMERICAL SIMULATIONS

Having introduced the GFN-ROM architecture, providing guarantees on super- and sub-resolution errors, we examine its performance on three challenging benchmarks for MOR. We compare the results of the methodology to intrusive and non-intrusive ROMs: namely, POD-Galerkin (POD-G) and POD with projection (representing the best possible linear method), POD-NN [47, 76] and GCA-ROM [77]. For fairness, we strive to make the MOR baselines for the different models as comparable as possible by adopting similar architectural and hyperparameter choices between methods, as shown in [app:hyperparams]Appendix <ref>app:hyperparams. Importantly, we do not fine tune the architecture or seek the best possible initialisation for each specific problem. As such, we stress that these results do not represent the best possible performance, but illustrate that GFN-ROM is able to directly handle several challenging benchmarks with no specific fine tuning. As a general rule, we have chosen a latent/reduced dimension $N,L= \lfloor 1.5\times N_\mu \rfloor$, and the small-data regime with a train/test split of 30/70.

Each of the problems examined here are selected to showcase a different physical or computational complexity. The first two benchmarks have been previously investigated with GCA-ROM [77], which we use as a baseline for nonlinear reduction – namely the Graetz problem and an advection-dominated problem. Furthermore, we examine a Stokes flow benchmark with a 7-dimensional parameter space to illustrate the model's capability of dealing with high-dimensional physical and geometric parameter spaces, and the low data regime with sparse samplings of the parameter space. Training data is generated via finite element methods using the RBniCS package [90], which is built on top of the FEniCS package [1, 67]. For each problem, we consider four meshes in order to test the multifidelity performance of our architecture[The meshes are generated from the largest mesh by using the moving front subsampling algorithm [61].].

The results for single fidelity experiments across our baselines are shown in <ref>.

POD-G     POD-Proj     POD-NN     GCA-ROM     4cGFN-ROM

Large     Large     Large     Large     Large     Medium     Small     Tiny
[HTML]E5E3E3Graetz     [HTML]E5E3E33.44     [HTML]E5E3E33.44     [HTML]E5E3E33.54     [HTML]E5E3E30.74     [HTML]E5E3E31.02     [HTML]E5E3E30.88     [HTML]E5E3E30.98     [HTML]E5E3E31.28

Advection     27.09     25.66     30.58     4.73     4.73     4.48     7.22     12.35

[HTML]E5E3E3Stokes     [HTML]E5E3E32.45     [HTML]E5E3E32.45     [HTML]E5E3E32.94     [HTML]E5E3E34.63     [HTML]E5E3E34.24     [HTML]E5E3E34.12     [HTML]E5E3E34.44     [HTML]E5E3E35.55

Mean relative errors (%) evaluated on large mesh.

Since GFN-ROM is capable of being trained on multifidelity data, we investigate its generalisation ability across different resolutions. In <ref>, we report the change in performance of GFN-ROM when one substitutes half of the training data with cheaper, lower-fidelity data.

Large &     Large &     Large &     Medium &     Medium &     Small &

Medium      Small     Tiny     Small     Tiny     Tiny
[HTML]E5E3E3Graetz     [HTML]E5E3E30.96 (+0.06)     [HTML]E5E3E30.98 (+0.04)     [HTML]E5E3E31.40 (-0.37)     [HTML]E5E3E34.44 (-3.56)     [HTML]E5E3E31.03 (-0.15)     [HTML]E5E3E31.09 (-0.11)

Advection     5.02 (-0.30)     5.35 (-0.62)     5.63 (-0.91)     5.22 (-0.73)     6.27 (-1.79)     8.77 (-1.55)

[HTML]E5E3E3Stokes     [HTML]E5E3E33.63 (+0.61)     [HTML]E5E3E34.94 (-0.70)     [HTML]E5E3E34.44 (-0.19)     [HTML]E5E3E34.51 (-0.39)     [HTML]E5E3E35.56 (-1.45)     [HTML]E5E3E35.65 (-1.21)

Mean relative errors (%) evaluated on large mesh, and the change when training only with the finer of the two in <ref>. A positive change indicates reduced error and thus better performance.

GFN-ROM is a more lightweight and flexible method compared to existing state-of-the-art ROMs. We show the method's computational efficiency with respect to GCA-ROM in <ref>. GFN-ROM offers large savings in terms of training time (up to 63x reduction) and in terms of the number of trainable parameters (up to 25x reduction). All simulations have been performed on a workstation equipped with an Nvidia Quadro RTX 4000 GPU.

4cGFN-ROM     GCA-ROM

Large     Medium     Small     Tiny     Large
[HTML]E5E3E3     [HTML]E5E3E3Training time (s)     [HTML]E5E3E3119      [HTML]E5E3E346     [HTML]E5E3E331     [HTML]E5E3E326     [HTML]E5E3E3600

[HTML]E5E3E3    [HTML]E5E3E3Trainable parameters     [HTML]E5E3E32898761     [HTML]E5E3E3911004     [HTML]E5E3E3311910     [HTML]E5E3E3115821     [HTML]E5E3E32898795

[HTML]E5E3E3-3*Graetz     [HTML]E5E3E3Mesh nodes     [HTML]E5E3E37205     [HTML]E5E3E32248 (31%)     [HTML]E5E3E3754 (10%)     [HTML]E5E3E3265 (4%)     [HTML]E5E3E37205
[HTML]FFFFFF    [HTML]FFFFFFTraining time (s)     [HTML]FFFFFF124     [HTML]FFFFFF53     [HTML]FFFFFF34     [HTML]FFFFFF23     [HTML]FFFFFF408.9

[HTML]FFFFFF    [HTML]FFFFFFTrainable parameters     [HTML]FFFFFF3538757     [HTML]FFFFFF1110702     [HTML]FFFFFF387298     [HTML]FFFFFF140282     [HTML]FFFFFF3538791

[HTML]FFFFFF-3*Advection    [HTML]FFFFFFMesh nodes     [HTML]FFFFFF8801     [HTML]FFFFFF2746 (31%)     [HTML]FFFFFF942 (11%)     [HTML]FFFFFF326 (4%)     [HTML]FFFFFF8801
[HTML]E5E3E3     [HTML]E5E3E3Training time (s)     [HTML]E5E3E3248     [HTML]E5E3E390     [HTML]E5E3E345     [HTML]E5E3E331     [HTML]E5E3E31949

[HTML]E5E3E3    [HTML]E5E3E3Trainable parameters     [HTML]E5E3E32827589     [HTML]E5E3E3905596     [HTML]E5E3E3316126      [HTML]E5E3E3118032     [HTML]E5E3E32827623

[HTML]E5E3E3-3*Stokes    [HTML]E5E3E3Mesh nodes     [HTML]E5E3E37019     [HTML]E5E3E32226 (32%)     [HTML]E5E3E3756 (11%)     [HTML]E5E3E3262 (4%)     [HTML]E5E3E37019

Computational efficiency of GFN-ROM in comparison to GCA-ROM, and its scaling w.r.t. the mesh size.

§.§ Graetz problem

The Graetz problem represents a common benchmark in MOR, combining a forced heat diffusion with horizontal heat advection in a parametrised geometry. Specifically, we study a steady-state Graetz problem on the parametrised geometry $\Omega(\mu_1) = \Omega^1 \cup \Omega^2(\mu_1) \in \mathbb{R}^2$, with $\Omega^1 = [0,1]\times[0,1]$ and $\Omega^2(\mu_1)=[1,1+\mu_1]\times[0,1]$  illustrated in <ref>, and discretised by four meshes as summarised in Table <ref>.

(A) The parameterised domain for the Graetz problem. (B) Mean relative errors for GFN-ROM trained on the large mesh for the Graetz problem. Squares indicate the parameter realisations seen during training.

Homogeneous Dirichlet conditions are imposed on $\Gamma_D = \Gamma_1 \cup \Gamma_5 \cup \Gamma_6$, whilst a constant Dirichlet boundary condition is imposed on $\Gamma_G =\Gamma_2(\mu_1) \cup \Gamma_4(\mu_1)$ as $u(\boldsymbol{\mu})|_{\Gamma_G}=1$. The remaining boundary $\Gamma_3(\mu_1)$ is given a homogeneous von Neumann boundary condition. The weak formulation of the Graetz problem reads as follows: given $\boldsymbol{\mu} \in \mathbb{P}$, find $u(\boldsymbol{\mu}) \in \mathbb{U}(\mu_1)$
s.t.μ_2 ∫_Ω(μ_1) ∇u ·∇v  dΩ+ ∫_Ω(μ_1) y(1-y) ∂u/∂x v  dΩ= 0,   ∀v ∈𝕍,
where $\mathbb{U}(\mu_1) = \{v \in H^1(\Omega(\mu_1)): v |_{\Gamma_D}=0, ~ v |_{\Gamma_G}=1 \}$ and $\mathbb{V}(\mu_1) = \{v \in H^1(\Omega(\mu_1)): v |_{\Gamma_D \cup \Gamma_G} =0\}$. Note that due to the parametrised domain, these function spaces are dependent on the geometric parameter.
This problem is parametrised by $\boldsymbol{\mu} = (\mu_1, \mu_2)$, where $\mu_1$ is a geometrical parameter governing the length of $\Omega^2$, and $\mu_2$ is a physical parameter, namely the diffusivity coefficient. We consider a dataset with 200 snapshots computed on a uniform grid of $\mathbb{P}$, with 10 and 20 equispaced values for $\mu_1 \in [1,3]$ and $\mu_2 \in [0.01, 0.1]$, respectively.

Training and testing the performance on GFN-ROM on the largest mesh, we achieve a mean relative error of 1.02% over the full dataset. As illustrated in <ref>, we can see GFN-ROM is well able to capture the overall solution. Note that this performance is very good, particularly considering only 60 parameter realisations are seen for training, and the model learns a latent representation of size $L=3$. Interestingly, for the Graetz problem we see that errors are much higher for small values of $\mu_2$ i.e. in the low diffusivity region, where the advection dominates. Moreover, few training samples are representing this regime, making it more difficult for GFN-ROM to learn how to model the region.

Nonetheless, the overall performance on GFN-ROM is extremely good on this benchmark, achieving over a 3x reduction in mean relative error compared to linear methods (see <ref>), and comparable performance to GCA-ROM.

Training on large mesh.

Training on medium mesh.

Training on small mesh.

Training on tiny mesh.

GFN-ROM solution fields for the Graetz problem when trained on (A) large mesh, (B) medium mesh, (C) small mesh, and (D) tiny mesh. The mesh points are shown by white pluses.

Moreover, we can further improve the performance of GFN-ROM by leveraging its resolution invariance. In fact, when training on a coarser representation, e.g. the medium mesh, the mean relative error is reduced to 0.88%, illustrating the new opportunities arising from GFN-ROM's ability to perform super-resolution. Even when training on the coarsest (tiny) discretisation, the mean relative error increases only slightly by 0.26% to 1.28% despite this representing only 4% of the training data. In <ref>, we show the predictions of the four different GFM-ROM models trained on the four different meshes for a parameter realisation $\boldsymbol{\mu}=(1.44,0.06)$ not seen during training. The model predictions are qualitatively indistinguishable, showing that GFN-ROM is capable of extracting enough information even on cheaper computational data without incurring any penalty on performance. Even further generalisations are possible thanks to the  possibility of training on mixed-fidelity data. GFN-ROM shows excellent generalisation properties and performance is minorly impacted when substituting high-fidelity data for lower-fidelity data, as seen in <ref>. For example, the performance even improves from 1.02% when training solely on the large mesh to 0.96% when substituting half of the large mesh dataset with cheaper data from the medium mesh (representing a 35% reduction in training data).

§.§ Advection-dominated problem

Advection-dominated problems, due to the slow decay of the Kolmogorov $n$-width [36, 73], represent an important benchmarking case in the ROM context to compare linear and nonlinear reduction strategies. We consider a fixed domain $\Omega=[0,1]\times[0,1]\in\mathbb{R}^2$, illustrated in <ref> with homogeneous Dirichlet boundary conditions. We generate four meshes summarised in Table <ref>.

(A) The parameterised domain for the advection problem. (B) Mean relative errors for GFN-ROM trained on the large mesh for the advection problem. Squares indicate the parameter realisations seen during training.

The weak formulation of the advection problem reads as follows:
given $\boldsymbol{\mu} \in \mathbb{P}$, find $u(\boldsymbol{\mu}) \in \mathbb{V}$ s.t.D(μ_1) ∫_Ω∇u ·∇v  dΩ+ ∫_Ω(β ·∇u) v  dΩ= ∫_Ωv  dΩ,   ∀v ∈𝕍= H^1_0(Ω).
A total of 100 snapshots are generated for this problem on a uniform grid of $\mathbb{P}$, with 10 and 10 equispaced values for $\mu_1 \in [0,6]$ and $\mu_2 \in [-1, 1]$, respectively. The first parameter $\mu_1$ determines the diffusion coefficient via $\operatorname{D}(\mu_1)=10^{-\mu_1}$ and the second parameter $\mu_2$ determines the strength and sign of the transport $\boldsymbol{\beta}(\mu_2)=(\mu_2,\mu_2)$.

As expected, POD-based techniques perform poorly due to their linear nature on this benchmark, with POD-G only achieving a mean relative error of 27.09% over the full dataset. GFN-ROM's nonlinear nature allows it to achieve a far superior mean relative error of only 4.73% when training on the large mesh, which matches GCA-ROM's performance. We show the relative errors over the parameter space in <ref>. We stress that this performance is possible even in the low data regime with only 30 training examples and a small latent representation of size $L=3$. Contrarily to the Graetz problem, it can be seen that poorer performances occur for lower Péclet numbers ($\operatorname{Pe}\approx \lvert \frac{\mu_2}{10^{-\mu_1}} \rvert$) compared to higher Péclet numbers i.e. the model performs better when transport dominates. This is likely due to the fact that the model training data is skewed towards examples where the advection in the opposite direction dominates, making it less adept at predicting solutions where diffusion has a stronger impact and the $\mu_2$ is positive.

On this benchmark, GFN-ROM's resolution invariant nature allows for training on cheaper meshes, and outperforming  GCA-ROM, achieving a mean relative error of 4.48% when training on the medium mesh (69% reduction in training data). In <ref>, we show the predictions of the GFN-ROM models trained on the four meshes for a parameter $\boldsymbol{\mu}=(3.33,-0.78)$ not seen during training. Qualitatively, we can see that all models but tiny are capable of learning well-enough the solution dynamics. The tradeoff to consider is of course given by the performance of large models and the computational efficiency of the smaller ones. Similarly to the Graetz benchmark, we also observed good generalisation for the multifidelity cases, as shown in <ref>.

Training on large mesh.

Training on medium mesh.

Training on small mesh.

Training on tiny mesh.

GFN-ROM solution fields for the advection problem when trained on (A) large mesh, (B) medium mesh, (C) small mesh, and (D) tiny mesh. The mesh points are shown by white pluses.

§.§ Stokes problem

We now consider a more complex benchmark in computational fluid dynamics, modelling a steady-state Stokes flow problem with high-dimensional physical/geometric parameter space and sparse sampling of it. This setting makes the problem particularly challenging to learn, and an important test case for the ability of a ROM to deal with the low-data regime. For simplicity, we consider the case of a viscosity of fixed value $1$, but we allow for a variable forcing term $\bfmf = (0, \mu_6)$, governed by a physical parameter $\mu_6$. We also consider a total of 6 geometrical parameters $\boldsymbol{\mu}_{\text{geo}} = (\mu_0, \mu_1, \mu_2, \mu_3, \mu_4, \mu_5)$, to deform the domain depicted in <ref>. We generate four meshes summarised in <ref>.

[left] at (-1, 3) $(0, \mu_0+ \mu_2 + \mu_3)$;
[right] at (1, 3) $(\mu_1, \mu_0+ \mu_2 + \mu_3)$;
[right] at (3, 1) $(\mu_1 + \mu_4, \mu_0+ \mu_2 + \mu_4\tan(\mu_5))$;
[right] at (3, -1) $(\mu_1 + \mu_4, \mu_2 + \mu_4\tan(\mu_5))$;
[left] at (-1, 1) $(0, \mu_0+\mu_2)$;
[left] at (-1, -1) $(0, \mu_2)$;
[left] at (-1, -3) $(0, 0)$;
[right] at (1, -3) $(\mu_1, 0)$;
[green!55!blue] (0,-1) rectangle (3,1);
[green!55!blue] (-1,-3) rectangle (1,3);
[white] at (1,0) $\Omega(\boldsymbol{\mu}_{\text{geo}})$;
(-1,-3) – (-1,3) node[pos=0.5,left] $\Gamma_{\text{w}}$;
(-1,3) – (1,3) node[white,pos=0.5,below] $\Gamma_{\text{N}}$;
(-1,-3) – (1,-3) node[white,pos=0.5,above] $\Gamma_{\text{N}}$;
(1,-1) – (3,-1) node[pos=0.5,below] $\Gamma_{\text{w}}$;
(1,1) – (3,1) node[pos=0.5,above] $\Gamma_{\text{w}}$;
(1,-3) – (1,-1);
(1,1) – (1,3);
(3,-1) – (3,1) node[pos=0.5,right] $\Gamma_{\text{in}}$;

figureThe domain for the Stokes problem.

As boundary conditions, the inlet $\Gamma_{\text{in}}$ is given a Dirichlet boundary condition of $\bfmu_{in} = [4(y-1)(2-y), 0]$, the walls $\Gamma_{\text{w}}$ are given zero Dirichlet boundary conditions and the outlet $\Gamma_{\text{N}}$ is left as a homogeneous Neumann boundary condition. A weak formulation for this problem reads: for a given parameter $\boldsymbol{\mu} \in\mathbb{P}$, find $\boldsymbol{u}(\boldsymbol{\mu}) \in\mathbb{U}$, $p \in\mathbb{M}$ s.t.

∫_Ω ∇u : ∇v   dΩ- ∫_Ω p ∇·v   dΩ= ∫_Ω f ·v   dΩ,   ∀v ∈𝕍,

∫_Ω q ∇·u   dΩ= 0,   ∀q ∈𝕄,

where the parameter dependent function spaces for the velocity are defined as $\mathbb{V}(\boldsymbol{\mu}_{\text{geo}}) = \{\bfmv \in (H^1(\Omega(\boldsymbol{\mu}_{\text{geo}}))^2: \bfmv|_{\Gamma_{w}\cup\Gamma_{in}}=\boldsymbol{0} \}$ and $\mathbb{U}(\boldsymbol{\mu}_{\text{geo}}) = \{\bfmv \in (H^1(\Omega(\boldsymbol{\mu}_{\text{geo}})))^2: \bfmv|_{\Gamma_{w}}=\boldsymbol{0}, \bfmv|_{\Gamma_{in}} = \bfmu_{in}\}$, and the function space for the pressure is defined as $\mathbb{M}(\boldsymbol{\mu}_{\text{geo}}) = L^2(\Omega(\boldsymbol{\mu}_{\text{geo}}))$. The problem uses a mixed finite element discretisation with the velocity and pressure as solution variables, meaning that the inf-sup condition is necessary for the well posedness of this problem and as a result the supremiser operator $T^{\mu}: \mathbb{M} \rightarrow \mathbb{V}$ is used (see [6] for more details).

We consider a sparse sampling with solutions computed only for the following parameter choices: $\mu_0, \mu_1, \mu_2, \mu_3, \mu_4 \in \{0.5,1.5 \}, \mu_5 \in \{-\frac{\pi}{6}, \frac{\pi}{6} \}, \mu_6 \in \{-10,-8,-6,-4,-2,0,2,4,6,8,10\}$, leading to a total of 704 high-fidelity solutions, and the field of interest is the magnitude of the the velocity $|\bfmu|$.

Histogram of the mean relative error for GFN-ROM trained on the large mesh for the Stokes problem.

Training on the large mesh, GFN-ROM achieves a mean relative error of 4.24%, outperforming GCA-ROM. We stress that this performance is possible whilst remaining in the low data regime, undersampling the parameter space. Despite this, the model is well able to capture the solution. As shown in <ref>, model performance is found to be very similar on average for the possible geometrical parameters, while, as expected, the results are more sensitive to the physical parameter $\mu_6$, with slightly worse approximation for the case $\mu_6=0$ where the forcing term vanishes.

GFN-ROM shows even further performance improvements by training on cheaper computational meshes, achieving a mean relative error of 4.12% for the medium mesh (representing a 68% reduction in the number of nodes). Even when training on the coarsest discretisation with just 4% of the data, performance only increases to 5.55% and is still able to capture the behaviour with respect to the physical and geometric parameters not seen during training, respectively in <ref> and <ref>.

Super-resolution via GFN-ROM when training on the tiny mesh and evaluating on large mesh for variations in the physical parameter $\mu_6$ showcasing downward, none and upward forcings.

Super-resolution via GFN-ROM when training on the tiny mesh and evaluating on large mesh for variations in the geometrical parameters.

Once again, the complex parametric variability and the scarce amount of information (training on 211 parameter realisations for the 7-dimensional parameter space) do not compromise the generalisation capability of the methodology.

We also see excellent results for multifidelity training. In general, performances are only minorly affected when substituting half of the data for cheaper data, underlining the possibility of training on cheaper data whilst retaining excellent performance. It is even possible to achieve a better performance of 3.63% by replacing 50% of data on the large mesh with cheaper data on the medium mesh (representing a 34% reduction in training data), as detailed in <ref>.

§ CONCLUSION

In this work, we presented a novel  graph-based resolution-invariant ROM, as a powerful tool to investigate multifidelity applications in the MOR context. We showed that it is possible to apply feedforward networks for graphical data by means of the graph feedforward network, allowing for the extension of many single-fidelity architectures to the multifidelity setting. Furthermore, we showed that this extension comes with provable guarantees on the performance in terms of error bounds.

Beyond such theoretical guarantees, we have also demonstrated that our architecture performs well in practice on three challenging benchmarks involving parametrised domains, advection-dominated problems and high-dimensional parameter spaces. Across all benchmarks, GFN-ROM achieved satisfactory performances, especially considering we worked within a low-data regime. Indeed, like GCA-ROM, our architecture has been trained on only 30% of the dataset, and has been exploited to learn a small latent representation.
GFN-ROM achieved comparable or better performance than GCA-ROM whilst being more flexible, interpretable and lightweight, and applicable for multifidelity data.

We showed that GFN-ROM can perform super-resolution, and demonstrated that within the proposed architecture, it is possible to substitute high-fidelity data for lower-fidelity data without incurring performance penalties. In practice, we showed that often performance can even be improved when training on cheaper computational data. Our work thus highlighted the importance of multifidelity ROMs, showing that computational
savings can be made by training on smaller meshes without deterioration in performance.

Future work will investigate further possibilities with graph feedforward networks, including tackling time-dependent problems, problems with spatially dependent parameters, enforcing more smoothness in GFN-ROM predictions, forcing greater dependence on local connections and automating finding the best training mesh (adapting the number of model parameters during training). We also plan to investigate multimodal modelling by using information from the Fourier domain, and the generalisation of convolutional neural operators to unstructured grids.

§ ACKNOWLEDGMENT

FP acknowledges the “GO for IT” program within the CRUI fund for the project “Reduced order method for nonlinear PDEs enhanced by machine learning”. This work has been conducted within the research activities of the consortium iNEST (Interconnected North-East Innovation Ecosystem), Piano Nazionale di Ripresa e Resilienza (PNRR) – Missione 4 Componente 2, Investimento 1.5 – D.D. 1058 23/06/2022, ECS00000043, supported by the European Union's NextGenerationEU program.

[1]
M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells.
The FEniCS project version 1.5.
Archive of Numerical Software, 3(100), 2015.

[2]
M. A. Alvarez, L. Rosasco, N. D. Lawrence, et al.
Kernels for vector-valued functions: A review.
Foundations and Trends® in Machine Learning,
4(3):195–266, 2012.

[3]
D. Amsallem, M. J. Zahr, and C. Farhat.
Nonlinear model order reduction based on local reduced-order bases.
International Journal for Numerical Methods in Engineering,
92(10):891–916, 2012.

[4]
A. Arzani and S. T. Dawson.
Data-driven cardiovascular flow modelling: examples and
Journal of the Royal Society Interface, 18(175):20200802, 2021.

[5]
J. Ballani, D. Kressner, and M. D. Peters.
Multilevel tensor approximation of PDEs with random data.
Stochastics and Partial Differential Equations: Analysis and
Computations, 5(3):400–427, 2017.

[6]
F. Ballarin, A. Manzoni, A. Quarteroni, and G. Rozza.
Supremizer stabilization of POD–Galerkin approximation of
parametrized steady incompressible Navier–Stokes equations.
International Journal for Numerical Methods in Engineering,
102(5):1136–1161, 2015.

[7]
J. L. Barnett, C. Farhat, and Y. Maday.
Neural-network-augmented projection-based model order reduction for
mitigating the Kolmogorov barrier to reducibility of CFD models.
arXiv preprint arXiv:2212.08939, 2022.

[8]
A. Becker, P. Finger, A. Meyer-Christoffer, B. Rudolf, K. Schamm, U. Schneider,
and M. Ziese.
A description of the global land-surface precipitation data products
of the global precipitation climatology centre with sample applications
including centennial (trend) analysis from 1901–present.
Earth System Science Data, 5(1):71–99, 2013.

[9]
F. D. A. Belbute-Peres, T. Economon, and Z. Kolter.
Combining differentiable PDE solvers and graph neural networks for
fluid flow prediction.
In International Conference on Machine Learning, pages
2402–2411. PMLR, 2020.

[10]
P. Benner, M. Ohlberger, A. Cohen, and K. Willcox.
Model reduction and approximation: Theory and algorithms.
SIAM, 2017.

[11]
J. L. Bentley.
Multidimensional binary search trees used for associative searching.
Communications of the ACM, 18(9):509–517, 1975.

[12]
B. Bonev, T. Kurth, C. Hundt, J. Pathak, M. Baust, K. Kashinath, and
A. Anandkumar.
Spherical Fourier neural operators: Learning stable dynamics on
the sphere.
In International Conference on Machine Learning, pages
2806–2823. PMLR, 2023.

[13]
N. Boullé and A. Townsend.
A Mathematical Guide to Operator Learning, 2023.

[14]
S. L. Brunton, J. L. Proctor, and J. N. Kutz.
Discovering governing equations from data by sparse identification of
nonlinear dynamical systems.
Proceedings of the national academy of sciences,
113(15):3932–3937, 2016.

[15]
T. Bui-Thanh, M. Damodaran, and K. Willcox.
Proper orthogonal decomposition extensions for parametric
applications in compressible aerodynamics.
In 21st AIAA Applied Aerodynamics Conference, page 4213, 2003.

[16]
S. Cao.
Choose a transformer: Fourier or Galerkin.
Advances in Neural Information Processing Systems,
34:24924–24940, 2021.

[17]
W. Chen, Q. Wang, J. S. Hesthaven, and C. Zhang.
Physics-informed machine learning for reduced-order modeling of
nonlinear problems.
Journal of Computational Physics, 446:110666, 2021.

[18]
L. Cicci, S. Fresca, M. Guo, A. Manzoni, and P. Zunino.
Uncertainty quantification for nonlinear solid mechanics using
reduced order models with Gaussian process regression.
arXiv preprint arXiv:2302.08216, 2023.

[19]
P. Conti, M. Guo, A. Manzoni, and J. S. Hesthaven.
Multi-fidelity surrogate modeling using long short-term memory
Computer Methods in Applied Mechanics and Engineering,
404:115811, 2023.

[20]
S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, and F. Piccialli.
Scientific machine learning through physics–informed neural
networks: Where we are and what’s next.
Journal of Scientific Computing, 92(3):88, 2022.

[21]
T. Danel, P. Spurek, J. Tabor, M. Śmieja, Ł. Struski, A. Słowik, and
Ł. Maziarka.
Spatial graph convolutional networks.
In Neural Information Processing: 27th International Conference,
ICONIP 2020, Bangkok, Thailand, November 18–22, 2020, Proceedings, Part V,
pages 668–675. Springer, 2020.

[22]
T. De Ryck, A. D. Jagtap, and S. Mishra.
Error estimates for physics informed neural networks approximating
the navier-stokes equations.
arXiv preprint arXiv:2203.09346, 2022.

[23]
N. Demo, M. Tezzele, and G. Rozza.
A non-intrusive approach for the reconstruction of POD modal
coefficients through active subspaces.
Comptes Rendus Mécanique, 347(11):873–881, 2019.

[24]
P. Díez, A. Muixí, S. Zlotnik, and A. García-González.
Nonlinear dimensionality reduction for parametric problems: A
kernel proper orthogonal decomposition.
International Journal for Numerical Methods in Engineering,
122(24):7306–7327, 2021.

[25]
T. Dinku.
Challenges with availability and quality of climate data in Africa.
In Extreme hydrology and climate variability, pages 71–80.
Elsevier, 2019.

[26]
T. Dinku, M. C. Thomson, R. Cousin, J. del Corral, P. Ceccato, J. Hansen, and
S. J. Connor.
Enhancing national climate services (ENACTS) for development in
Climate and Development, 10(7):664–672, 2018.

[27]
O. Ferludin, A. Eigenwillig, M. Blais, D. Zelle, J. Pfeifer,
A. Sanchez-Gonzalez, W. L. S. Li, S. Abu-El-Haija, P. Battaglia,
N. Bulut, J. Halcrow, F. M. G. de Almeida, P. Gonnet, L. Jiang, P. Kothari,
S. Lattanzi, A. Linhares, B. Mayer, V. Mirrokni, J. Palowitch, M. Paradkar,
J. She, A. Tsitsulin, K. Villela, L. Wang, D. Wong, and B. Perozzi.
TF-GNN: graph neural networks in tensorflow.
CoRR, abs/2207.03522, 2023.

[28]
G. M. Fernández-Godino.
Review of multi-fidelity models.
Advances in Computational Science and Engineering,
1(4):351–400, 2023.

[29]
M. Fey and J. E. Lenssen.
Fast graph representation learning with PyTorch Geometric.
In ICLR Workshop on Representation Learning on Graphs and
Manifolds, 2019.

[30]
K. Fidkowski.
Quantifying uncertainties in radiation hydrodynamics models.
Ann Arbor, 1001:48109, 2014.

[31]
A. I. Forrester.
Black-box calibration for complex-system simulation.
Philosophical Transactions of The Royal Society A: Mathematical,
Physical and Engineering Sciences, 368(1924):3567–3579, 2010.

[32]
N. R. Franco, S. Fresca, F. Tombari, and A. Manzoni.
Deep learning-based surrogate models for parametrized PDEs:
Handling geometric variability through graph neural networks.
Chaos: An Interdisciplinary Journal of Nonlinear Science,
33(12):123121, 2023.

[33]
S. Fresca, L. Dede’, and A. Manzoni.
A comprehensive deep learning-based approach to reduced order
modeling of nonlinear time-dependent parametrized PDEs.
Journal of Scientific Computing, 87:1–36, 2021.

[34]
R. B. Gramacy.
laGP: large-scale spatial modeling via local approximate Gaussian
processes in R.
Journal of Statistical Software, 72:1–46, 2016.

[35]
R. B. Gramacy.
Surrogates: Gaussian process modeling, design, and
optimization for the applied sciences.
CRC press, 2020.

[36]
C. Greif and K. Urban.
Decay of the Kolmogorov n-width for wave problems.
Applied Mathematics Letters, 96:216–222, 2019.

[37]
M. Guo and J. S. Hesthaven.
Reduced order modeling for nonlinear structural analysis using
Gaussian process regression.
Computer methods in applied mechanics and engineering,
341:807–826, 2018.

[38]
M. Guo and J. S. Hesthaven.
Data-driven reduced order modeling for time-dependent problems.
Computer methods in applied mechanics and engineering,
345:75–99, 2019.

[39]
M. Guo, A. Manzoni, M. Amendt, P. Conti, and J. S. Hesthaven.
Multi-fidelity regression using artificial neural networks:
Efficient approximation of parameter-dependent output quantities.
Computer Methods in Applied Mechanics and Engineering,
389:114378, 2022.

[40]
Z. Hao, Z. Wang, H. Su, C. Ying, Y. Dong, S. Liu, Z. Cheng, J. Song, and
J. Zhu.
GNOT: A general neural operator transformer for operator learning.
In International Conference on Machine Learning, pages
12556–12569. PMLR, 2023.

[41]
H. Harbrecht, M. Peters, and M. Siebenmorgen.
Multilevel accelerated quadrature for PDEs with log-normally
distributed diffusion coefficient.
SIAM/ASA Journal on Uncertainty Quantification, 4(1):520–551,

[42]
C. Heiß, I. Gühring, and M. Eigel.
A neural multilevel method for high-dimensional parametric PDEs.
In The Symbiosis of Deep Learning and Differential Equations,

[43]
C. Heiß, I. Gühring, and M. Eigel.
Multilevel CNNs for parametric PDEs.
arXiv preprint arXiv:2304.00388, 2023.

[44]
Q. Hernández, A. Badías, F. Chinesta, and E. Cueto.
Thermodynamics-informed graph neural networks.
arXiv preprint arXiv:2203.01874, 2022.

[45]
Q. Hernandez, A. Badias, D. Gonzalez, F. Chinesta, and E. Cueto.
Deep learning of thermodynamics-aware reduced-order models from data.
Computer Methods in Applied Mechanics and Engineering,
379:113763, 2021.

[46]
J. S. Hesthaven, G. Rozza, and B. Stamm.
Certified Reduced Basis Methods for Parametrized Partial
Differential Equations.
SpringerBriefs in Mathematics. Springer International Publishing,

[47]
J. S. Hesthaven and S. Ubbiali.
Non-intrusive reduced order modeling of nonlinear problems using
neural networks.
Journal of Computational Physics, 363:55–78, 2018.

[48]
W. Huang, Y. Rong, T. Xu, F. Sun, and J. Huang.
Tackling over-smoothing for general graph convolutional networks.
arXiv preprint arXiv:2008.09864, 2020.

[49]
W. Jiang and J. Luo.
Graph neural network for traffic forecasting: A survey.
Expert Systems with Applications, page 117921, 2022.

[50]
X. Jin, S. Cai, H. Li, and G. E. Karniadakis.
Nsfnets (navier-stokes flow nets): Physics-informed neural networks
for the incompressible navier-stokes equations.
Journal of Computational Physics, 426:109951, 2021.

[51]
A. Kashefi and T. Mukerji.
Physics-informed pointnet: A deep learning solver for steady-state
incompressible flows and thermal fields on multiple sets of irregular
arXiv preprint arXiv:2202.05476, 2022.

[52]
M. Kast, M. Guo, and J. S. Hesthaven.
A non-intrusive multifidelity method for the reduced order modeling
of nonlinear problems.
Computer Methods in Applied Mechanics and Engineering,
364:112947, 2020.

[53]
M. C. Kennedy and A. O'Hagan.
Predicting the output from a complex computer code when fast
approximations are available.
Biometrika, 87(1):1–13, 2000.

[54]
M. Khamlich, F. Pichi, and G. Rozza.
Optimal Transport-inspired Deep Learning Framework for
Slow-Decaying Problems: Exploiting Sinkhorn Loss and Wasserstein
Kernel, 2023.

[55]
N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and
A. Anandkumar.
Neural operator: Learning maps between function spaces with
applications to pdes.
Journal of Machine Learning Research, 24(89):1–97, 2023.

[56]
B. Kramer, B. Peherstorfer, and K. E. Willcox.
Learning nonlinear reduced models from data with operator inference.
Annual Review of Fluid Mechanics, 56, 2024.

[57]
M. A. Kramer.
Nonlinear principal component analysis using autoassociative neural
AIChE Journal, 37(2):233–243, 1991.

[58]
T. Kurth, S. Subramanian, P. Harrington, J. Pathak, M. Mardani, D. Hall,
A. Miele, K. Kashinath, and A. Anandkumar.
FourCastNet: Accelerating global high-resolution weather
forecasting using adaptive Fourier neural operators.
In Proceedings of the Platform for Advanced Scientific Computing
Conference, pages 1–11, 2023.

[59]
Y. Kuya, K. Takeda, X. Zhang, and A. I. Forrester.
Multifidelity surrogate modeling of experimental and computational
aerodynamic data sets.
AIAA journal, 49(2):289–298, 2011.

[60]
R. Lam, A. Sanchez-Gonzalez, M. Willson, P. Wirnsberger, M. Fortunato, F. Alet,
S. Ravuri, T. Ewalds, Z. Eaton-Rosen, W. Hu, et al.
Learning skillful medium-range global weather forecasting.
Science, page eadi2336, 2023.

[61]
A. P. Lawrence, M. E. Nielsen, and B. Fornberg.
Node subsampling for multilevel meshfree elliptic PDE solvers.
arXiv preprint arXiv:2303.09080, 2023.

[62]
K. Lee and K. T. Carlberg.
Model reduction of dynamical systems on nonlinear manifolds using
deep convolutional autoencoders.
Journal of Computational Physics, 404:108973, 2020.

[63]
Q. Li, Z. Han, and X.-M. Wu.
Deeper insights into graph convolutional networks for semi-supervised
In Proceedings of the AAAI Conference on Artificial
Intelligence, volume 32, 2018.

[64]
Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and
A. Anandkumar.
Fourier neural operator for parametric partial differential
arXiv preprint arXiv:2010.08895, 2020.

[65]
Z. Li, H. Zheng, N. Kovachki, D. Jin, H. Chen, B. Liu, K. Azizzadenesheli, and
A. Anandkumar.
Physics-informed neural operator for learning partial differential
arXiv preprint arXiv:2111.03794, 2021.

[66]
C. Liu, Y. Zhan, J. Wu, C. Li, B. Du, W. Hu, T. Liu, and D. Tao.
Graph pooling for graph neural networks: Progress, challenges, and
arXiv preprint arXiv:2204.07321, 2022.

[67]
A. Logg, K.-A. Mardal, and G. Wells, editors.
Automated Solution of Differential Equations by the
Finite Element Method: The FEniCS Book, volume 84 of Lecture
Notes in Computational Science and Engineering.
Springer, Berlin, Heidelberg, 2012.

[68]
K. O. Lye, S. Mishra, and R. Molinaro.
A multi-level procedure for enhancing accuracy of machine learning
European Journal of Applied Mathematics, 32(3):436–469, 2021.

[69]
A. Manzoni, F. Negri, and A. Quarteroni.
Dimensionality reduction of parameter-dependent problems through
proper orthogonal decomposition.
Annals of Mathematical Sciences and Applications,
1(2):341–377, 2016.

[70]
X. Meng and G. E. Karniadakis.
A composite neural network that learns from multi-fidelity data:
Application to function approximation and inverse PDE problems.
Journal of Computational Physics, 401:109020, 2020.

[71]
B. Moya, A. Badias, D. Gonzalez, F. Chinesta, and E. Cueto.
Physics perception in sloshing scenes with guaranteed thermodynamic
IEEE Transactions on Pattern Analysis and Machine Intelligence,
45(2):2136–2150, 2022.

[72]
S. K. Mukkavilli, D. S. Civitarese, J. Schmude, J. Jakubik, A. Jones,
N. Nguyen, C. Phillips, S. Roy, S. Singh, C. Watson, et al.
AI foundation models for weather and climate: Applications,
design, and implementation.
arXiv preprint arXiv:2309.10808, 2023.

[73]
M. Ohlberger and S. Rave.
Reduced basis methods: Success, limitations and future challenges.
arXiv preprint arXiv:1511.02021, 2015.

[74]
J. Pathak, S. Subramanian, P. Harrington, S. Raja, A. Chattopadhyay,
M. Mardani, T. Kurth, D. Hall, Z. Li, K. Azizzadenesheli, et al.
FourCastNet: A global data-driven high-resolution weather model
using adaptive Fourier neural operators.
arXiv preprint arXiv:2202.11214, 2022.

[75]
B. Peherstorfer and K. Willcox.
Data-driven operator inference for nonintrusive projection-based
model reduction.
Computer Methods in Applied Mechanics and Engineering,
306:196–215, 2016.

[76]
F. Pichi, F. Ballarin, G. Rozza, and J. S. Hesthaven.
An artificial neural network approach to bifurcating phenomena in
computational fluid dynamics.
Computers & Fluids, 254:105813, 2023.

[77]
F. Pichi, B. Moya, and J. S. Hesthaven.
A graph convolutional autoencoder approach to model order reduction
for parametrized PDEs.
arXiv preprint arXiv:2305.08573, 2023.

[78]
E. Qian, I.-G. Farcas, and K. Willcox.
Reduced operator inference for nonlinear partial differential
SIAM Journal on Scientific Computing, 44(4):A1934–A1959, 2022.

[79]
A. Quarteroni, A. Manzoni, and F. Negri.
Reduced basis methods for partial differential equations: An
introduction, volume 92.
Springer, 2015.

[80]
M. Raissi, H. Babaee, and P. Givi.
Deep learning of turbulent scalar mixing.
Physical Review Fluids, 4(12):124501, 2019.

[81]
M. Raissi and G. E. Karniadakis.
Hidden physics models: Machine learning of nonlinear partial
differential equations.
Journal of Computational Physics, 357:125–141, 2018.

[82]
M. Raissi, P. Perdikaris, and G. E. Karniadakis.
Inferring solutions of differential equations using noisy
multi-fidelity data.
Journal of Computational Physics, 335:736–746, 2017.

[83]
M. Raissi, P. Perdikaris, and G. E. Karniadakis.
Physics-informed neural networks: A deep learning framework for
solving forward and inverse problems involving nonlinear partial differential
Journal of Computational Physics, 378:686–707, 2019.

[84]
M. Raissi, Z. Wang, M. S. Triantafyllou, and G. E. Karniadakis.
Deep learning of vortex-induced vibrations.
Journal of Fluid Mechanics, 861:119–137, 2019.

[85]
M. Raissi, A. Yazdani, and G. E. Karniadakis.
Hidden fluid mechanics: Learning velocity and pressure fields from
flow visualizations.
Science, 367(6481):1026–1030, 2020.

[86]
B. Raonić, R. Molinaro, T. Rohner, S. Mishra, and E. de Bezenac.
Convolutional neural operators.
arXiv preprint arXiv:2302.01178, 2023.

[87]
J. Reiss, P. Schulze, J. Sesterhenn, and V. Mehrmann.
The shifted proper orthogonal decomposition: A mode decomposition
for multiple transport phenomena.
SIAM Journal on Scientific Computing, 40(3):A1322–A1344, 2018.

[88]
F. Romor, G. Stabile, and G. Rozza.
Non-linear manifold reduced-order models with convolutional
autoencoders and reduced over-collocation method.
Journal of Scientific Computing, 94(3):74, 2023.

[89]
S. G. Rosofsky, H. Al Majed, and E. Huerta.
Applications of physics informed neural operators.
Machine Learning: Science and Technology, 4(2):025022, 2023.

[90]
G. Rozza, F. Ballarin, L. Scandurra, and F. Pichi.
Real Time Reduced Order Computational Mechanics:
Parametric PDEs Worked Out Problems, volume 5 of SISSA Springer
Springer Nature Switzerland, Cham, 2024.

[91]
T. K. Rusch, M. M. Bronstein, and S. Mishra.
A survey on oversmoothing in graph neural networks.
arXiv preprint arXiv:2303.10993, 2023.

[92]
L. Sun, H. Gao, S. Pan, and J.-X. Wang.
Surrogate modeling for fluid flows based on physics-constrained deep
learning without simulation data.
Computer Methods in Applied Mechanics and Engineering,
361:112732, 2020.

[93]
M. Sun, S. Zhao, C. Gilvary, O. Elemento, J. Zhou, and F. Wang.
Graph convolutional networks for computational drug development and
Briefings in Bioinformatics, 21(3):919–935, 2020.

[94]
Y. Sun, C. Moya, G. Lin, and M. Yue.
DeepGraphONet: A deep graph operator network to learn and
zero-shot transfer the dynamic response of networked systems.
arXiv preprint arXiv:2209.10622, 2022.

[95]
T. Taddei.
A registration method for model order reduction: data compression and
geometry reduction.
SIAM Journal on Scientific Computing, 42(2):A997–A1027, 2020.

[96]
A. L. Teckentrup, P. Jantsch, C. G. Webster, and M. Gunzburger.
A multilevel stochastic collocation method for partial differential
equations with random input data.
SIAM/ASA Journal on Uncertainty Quantification,
3(1):1046–1074, 2015.

[97]
J. Wang, S. Zhang, Y. Xiao, and R. Song.
A review on graph neural network methods in financial applications.
arXiv preprint arXiv:2111.15367, 2021.

[98]
M. Wang, D. Zheng, Z. Ye, Q. Gan, M. Li, X. Song, J. Zhou, C. Ma, L. Yu,
Y. Gai, T. Xiao, T. He, G. Karypis, J. Li, and Z. Zhang.
Deep graph library: A graph-centric, highly-performant package for
graph neural networks.
arXiv preprint arXiv:1909.01315, 2019.

[99]
S. Wang, Y. Teng, and P. Perdikaris.
Understanding and mitigating gradient flow pathologies in
physics-informed neural networks.
SIAM Journal on Scientific Computing, 43(5):A3055–A3081, 2021.

[100]
S. Wang, X. Yu, and P. Perdikaris.
When and why pinns fail to train: A neural tangent kernel
Journal of Computational Physics, 449:110768, 2022.

[101]
S. Zhang, H. Tong, J. Xu, and R. Maciejewski.
Graph convolutional networks: Algorithms, applications and open
In Computational Data and Social Networks: 7th International
Conference, CSoNet 2018, Shanghai, China, December 18–20, 2018, Proceedings
7, pages 79–91. Springer, 2018.

[102]
S. Zhang, H. Tong, J. Xu, and R. Maciejewski.
Graph convolutional networks: A comprehensive review.
Computational Social Networks, 6(1):1–23, 2019.

[103]
Z. Zhang, M. Guo, and J. S. Hesthaven.
Model order reduction for large-scale structures with local
Computer Methods in Applied Mechanics and Engineering,
353:491–515, 2019.

§ PROOFS

§.§ Extension to multiple layers

For ease of notation we define each of the outputs of the $\operatorname{GFN}^{\mathcal{M}_o \to \mathcal{M}_n}$ function as
GFN_W^e^ℳ_o →ℳ_n(^e),  GFN_W^d^ℳ_o →ℳ_n(^d),  GFN_b^d^ℳ_o →ℳ_n(^d) = GFN^ℳ_o →ℳ_n(^e, ^d, ^d).
The single-layer GFN autoencoder given in <ref> and <ref> can then be rewritten as
enc^ℳ_o →ℳ_n(_n)_i    =    σ( ∑_j_n=1^N_n GFN_W^e^ℳ_o →ℳ_n(^e)_ij_n u_j_n + b^e_i ),     ∀i=1,⋯,L,

dec^ℳ_o →ℳ_n()_i_n    =    ∑_j=1^L GFN_W^d^ℳ_o →ℳ_n(^d)_i_nj z_j + b^d_i_n,     ∀i_n=1,⋯,N_n.
This can be extended to $Q+1$ hidden layers as
enc_0^ℳ_o →ℳ_n()_i    =    σ( ∑_j_n=1^N_n GFN_W^e^ℳ_o →ℳ_n(^e)_ij_n x_j_n + b^e_i ) ,

enc_P⋯1()_i    =    σ( ∑_l_P=1^L_P W^(P)_i,l_P σ( ⋯σ( ∑_l_1=1^L_1 W^(1)_l_2, l_1 x_l_1 + b^(1)_l_2 ) + ⋯) + b^(P)_i ),

dec_Q⋯P+1()_i    =    σ( ∑_l_Q=1^L_Q W^(Q)_i, l_Q σ( ⋯σ( ∑_l_P+1=1^L_P+1 W^(P+1)_l_P+2, l_P+1 x_l_P+1 + b^(P+1)_l_P+2 ) + ⋯) + b^(Q)_i ),

dec_Q+1^ℳ_o →ℳ_n()_i_n    =    ∑_j=1^L_Q+1 GFN_W^d^ℳ_o →ℳ_n(^d)_i_nj x_j + GFN_b^d^ℳ_o →ℳ_n(^d)_i_n,

where $0 \leq P\leq Q$, $\mathcal{M}_n$ is a new mesh, $L_{i}$ represents the output sizes for the $(i+1)$-th network layer with $L_{P+1}=L$ and $L_{Q+2}=N_n$, $\operatorname{enc}^{\mathcal{M}_o \to \mathcal{M}_n}(\bfmx)_i = \operatorname{enc}_{P\cdots 1}(\operatorname{enc}^{\mathcal{M}_o \to \mathcal{M}_n}_{0}(\bfmx))_i, \forall i=1,\cdots,L$ is the encoded latent representation of an input $\bfmx \in \mathbb{R}^{N_n}$ on mesh $\mathcal{M}_n$ and $\operatorname{dec}^{\mathcal{M}_o \to \mathcal{M}_n}(\bfmz)_{i_{\inmesh{n}}}=\operatorname{dec}^{\mathcal{M}_o \to \mathcal{M}_n}_{Q+1}(\operatorname{dec}_{Q\cdots P+1}(\bfmz))_{i_{\inmesh{n}}}, \forall i_{\inmeshlg{n}}=1,\cdots,N_n$ is the reconstructed solution on $\mathcal{M}_n$ from a given latent representation $\bfmz \in \mathbb{R}^{L}$.

§.§ Simplified expansive equations

We consider the GFN transforms $\mathcal{M}_o \to \mathcal{M}_n$ and show that they can be simplified under the expansive condition given in <ref>.

First of all, we show that the condition implies
∀i_n,   { k_o  s.t.  k_o i_n } = { k_o  s.t.  k_o ←i_n },
which is a set containing a unique element $k_{\inmeshlg{o}}$ defined as $k_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{n}}$.

We have to show that $\{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \eitherarrow i_{\inmeshlg{n}} \} = \{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{n}} \} \cup \{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \rightarrow i_{\inmeshlg{n}} \} = \{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{n}} \}$. If $\{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \rightarrow i_{\inmeshlg{n}} \} = \emptyset$, this is trivially the case. Therefore, in the following we consider only the case where it is not.

Trivially, $|\{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{n}} \}|=1$ since a node $i_{\inmeshlg{n}}$ can only have one nearest neighbour. In fact, additionally we can show it is also the case that $|\{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \rightarrow i_{\inmeshlg{n}} \}|=1$, which can be proven by contradiction. Suppose that $|\{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \rightarrow i_{\inmeshlg{n}} \}|=R \neq 1$. This means that there exist multiple distinct nodes $\left\{k^r_{\inmeshlg{o}}\right\}_{r=1}^R$, each satisfying $k^r_{\inmeshlg{o}} \rightarrow i_{\inmeshlg{n}}$. Via the expansive condition in <ref>, it follows that each of these nodes must also satisfy $k^r_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{n}}$. This is not possible since $i_{\inmeshlg{n}}$ can only have a single nearest neighbour, meaning therefore $|\{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \rightarrow i_{\inmeshlg{n}} \}|=1$.

We note that the single node in the set $|\{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \rightarrow i_{\inmeshlg{n}} \}|$ must also satisfy $k_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{n}}$ due to <ref>. Therefore, it is in fact the exact same node as the node in the set $|\{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{n}} \}|$, thus concluding the proof.

Secondly, we show that for a given $k_{\inmeshlg{o}}$
|{ h_n  s.t.  k_o h_n }|= |{ h_n  s.t.  k_o ←h_n }|.

Due to <ref>, it is the case that $\{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \rightarrow h_{\inmeshlg{n}} \} \subseteq \{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow h_{\inmeshlg{n}} \}$, concluding the proof.

The updates can therefore be simplified to
W̃^e_ij_n    =    1/|{h_n  s.t.  k_o ←h_n }| W^e_ik_o,      where  k_o ←j_n,

W̃^d_i_nj    =    W^d_k_oj,      where  k_o ←i_n,

b̃^d_i_n    =    b^d_k_o,       where  k_o ←i_n,
which therefore concludes the proof.

§.§ Simplified agglomerative equations

We consider the GFN transforms $\mathcal{M}_o \to \mathcal{M}_n$ and show that they can be simplified under the agglomerative condition given in <ref>.

First of all, we show that the condition implies
∀j_n,   { k_o  s.t.  k_o j_n } = { k_o  s.t.  k_o →j_n }.

Due to <ref>, it is the case that $\{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow j_{\inmeshlg{n}} \} \subseteq \{ k_{\inmeshlg{o}} \text{ s.t. } k_{\inmeshlg{o}} \rightarrow j_{\inmeshlg{n}} \}$, concluding the proof.

Secondly, we show that for a given $k_{\inmeshlg{o}}$
|{ h_n  s.t.  k_o h_n }|= 1.

We want to show that
$$|\{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \eitherarrow h_{\inmeshlg{n}} \}| = |\{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \rightarrow h_{\inmeshlg{n}} \} \cup \{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow h_{\inmeshlg{n}} \}| = 1 .$$
If $\{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow h_{\inmeshlg{n}} \} = \emptyset$, this is trivially the case since $k_{\inmeshlg{o}}$ only has one nearest neighbour i.e. $\lvert \{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \rightarrow h_{\inmeshlg{n}} \}\rvert = 1$. Therefore, in the following we consider only the case where it is not.

With this, we can show it is also the case that $|\{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow h_{\inmeshlg{n}} \}|=1$, which can be proven by contradiction. Suppose that $|\{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow h_{\inmeshlg{n}} \}|=R \neq 1$. This means that there exist multiple distinct nodes $\left\{h^r_{\inmeshlg{n}}\right\}_{r=1}^R$, each satisfying $k_{\inmeshlg{o}} \leftarrow h^r_{\inmeshlg{n}}$. Via the agglomerative condition in <ref>, it follows that each of these nodes must also satisfy $k_{\inmeshlg{o}} \rightarrow h^r_{\inmeshlg{n}}$. This is not possible since $k_{\inmeshlg{o}}$ can only have a single nearest neighbour, meaning therefore $|\{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \rightarrow h_{\inmeshlg{n}} \}|=1$.

We note that the single node in the set $|\{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow h_{\inmeshlg{n}} \}|$ must also satisfy $k_{\inmeshlg{o}} \rightarrow h_{\inmeshlg{n}}$ due to <ref>. Therefore, it is in fact the exact same node as the node in the set $|\{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} \rightarrow h_{\inmeshlg{n}} \}|$, thus concluding the proof.

The updates can therefore be simplified to
W̃^e_ij_n    =    ∑_∀k_o  s.t.   k_o →j_n W^e_ik_o,

W̃^d_i_nj    =    ∀k_o  s.t.   k_o →i_nmean W^d_k_oj,

b̃^d_i_n    =     ∀k_o  s.t.   k_o →i_nmean b^d_k_o,
which therefore concludes the proof.

§.§ Self-consistency

We consider a sequence of two transforms, firstly an expansive transform $\mathcal{M}_o \to \mathcal{M}_n$, and then secondly a transform $\mathcal{M}_n \to \mathcal{M}_o$. Note that if $\mathcal{M}_o \to \mathcal{M}_n$ is expansive, then it must follow by definition that $\mathcal{M}_n \to \mathcal{M}_o$ is agglomerative.
The first transform $\mathcal{M}_o \to \mathcal{M}_n$ is expansive. Following <ref>, this gives the first update as
Ŵ^e_ij_n    =    1/|{h_n  s.t.  k_o ←h_n }| W^e_ik_o,      where  k_o ←j_n,

Ŵ^d_i_nj    =    W^d_k_oj,      where  k_o ←i_n,

b̂^d_i_n    =    b^d_k_o,       where  k_o ←i_n.
The second transform $\mathcal{M}_o \leftarrow \mathcal{M}_n$ is agglomerative (note that the transform is not $\mathcal{M}_o \to \mathcal{M}_n$). Following <ref>, this gives the second update as
W̃^e_ij_o    =    ∑_∀k_n  s.t.  j_o ←k_n Ŵ^e_ik_n,

W̃^d_i_oj    =    ∀k_n  s.t.  i_o ←k_nmean Ŵ^d_k_nj,

b̃^d_i_o    =     ∀k_n  s.t.  i_o ←k_nmean b^d_k_n.
Subbing everything in gives
W̃^e_ij_o    =    ∑_∀k_n  s.t.  j_o ←k_n 1/|{h_n  s.t.  l_o ←h_n }| W^e_il_o,      where  l_o ←k_n,

=    ∑_∀k_n  s.t.  j_o ←k_n 1/|{h_n  s.t.  j_o ←h_n }| W^e_ij_o,

=    W^e_ij_o,

W̃^d_i_oj    =    ∀k_n  s.t.  i_o ←k_nmean W^d_l_oj,      where  l_o ←k_n,

=    ∀k_n  s.t.  i_o ←k_nmean W^d_i_oj,

=    W^d_i_oj,

b̃^d_i_o    =     ∀k_n  s.t.  i_o ←k_nmean b^d_l_o,       where  l_o ←k_n,

=    ∀k_n  s.t.  i_o ←k_nmean b^d_i_o,

=    b^d_i_o.
This concludes the proof.

§.§ General GFN transform as a composition of an expansive and agglomerative transform

Splitting up the summations, the general GFN transforms given in <ref> can be rewritten as

W̃^e_ij_n    = [ ∑_∀k_o  s.t  k_o →j_n
+ ∑_∀k_o  s.t  k_o ←j_n
but not k_o →j_n ] W^e_ik_o/|{ h_n  s.t.  k_o    h_n }|,

W̃^d_i_nj    = 1/|{ k_o  s.t  k_o i_n }| [ ∑_∀k_o  s.t  k_o →i_n + ∑_∀k_o  s.t  k_o ←i_n
but not k_o →i_n ] W^d_k_oj,

b̃^d_i_n    = 1/|{ k_o  s.t  k_o i_n }| [ ∑_∀k_o  s.t  k_o →i_n + ∑_∀k_o  s.t  k_o ←i_n
but not k_o →i_n ] b^d_k_o.

By defining an auxiliary mesh
ℳ_ = ℳ_o ∪{ℳ_n [i_n]  s.t.  j_o ←i_n  but not  j_o →i_n },
we can convert the summations to
∑_∀k_o  s.t  k_o →i_n ⟶∑_∀l_  s.t  l_ →i_n
and l_ ∈ℳ_o,   k_o ←l_,

∑_∀k_o  s.t  k_o ←i_n
but not k_o →i_n ⟶∑_∀l_  s.t  l_ →i_n
and l_ ∈ℳ_n,   k_o ←l_.
Subbing this in, this gives an updated expression for the general GFN transforms of

W̃^e_ij_n    = [ ∑_∀l_  s.t  l_ →j_n
and l_ ∈ℳ_o
+ ∑_∀l_  s.t  l_ →j_n
and l_ ∈ℳ_n ] W^e_ik_o/|{ h_n  s.t.  k_o    h_n }|,   k_o ←l_,

W̃^d_i_nj    = 1/|{ g_o  s.t  g_o i_n }| [ ∑_∀l_  s.t  l_ →i_n
and l_ ∈ℳ_o
+ ∑_∀l_  s.t  l_ →i_n
and l_ ∈ℳ_n ] W^d_k_oj,   k_o ←l_,

b̃^d_i_n    = 1/|{ g_o  s.t  g_o i_n }| [ ∑_∀l_  s.t  l_ →i_n
and l_ ∈ℳ_o
+ ∑_∀l_m  s.t  l_ →i_n
and l_ ∈ℳ_n ] b^d_i_n,   k_o ←l_.

Since it is the case by construction of $\mathcal{M}$ that $|\{ h_{\inmeshlg{n}} \text{ s.t. } k_{\inmeshlg{o}} ~\eitherarrow~ h_{\inmeshlg{n}} \}| = |\{ h_{\inmeshlg{}} \text{ s.t. } k_{\inmeshlg{o}} \leftarrow h_{\inmeshlg{}} \}|$ and $|\{ g_{\inmeshlg{o}} \text{ s.t } g_{\inmeshlg{o}} \eitherarrow i_{\inmeshlg{n}} \}| = |\{ l_{\inmeshlg{}} \text{ s.t } l_{\inmeshlg{}} \rightarrow i_{\inmeshlg{n}} \}|$, we can further rewrite these expressions as

W̃^e_ij_n    = ∑_∀l_  s.t  l_ →j_n W^e_ik_o/|{ h_  s.t.  k_o ←h_ }|,   k_o ←l_,

W̃^d_i_nj    = ∀l_  s.t  l_ →i_nmean W^d_k_oj,   k_o ←l_,

b̃^d_i_n    = ∀l_  s.t  l_ →i_nmean b^d_k_o,   k_o ←l_.

We can finally compute the transform in two steps, as

\begin{equation*}
\begin{rcases}
\hat{W}^e_{ij_{\inmesh{}}} &= \frac{1}{|\{h_{\inmesh{}} \text{ s.t. } k_{\inmesh{o}} \leftarrow h_{\inmesh{}} \}|} {W}^e_{ik_{\inmesh{o}}}, \qquad \text{ where}\quad k_{\inmeshlg{o}} \leftarrow j_{\inmeshlg{}}, \\
\hat{W}^d_{i_{\inmesh{}}j} &= {W}^d_{k_{\inmesh{o}}j}, \qquad \text{ where}\quad k_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{}}, \\
{{\hat{b}}^d}_{i_{\inmesh{}}} &= {b}^d_{k_{\inmesh{o}}},  \qquad \text{ where}\quad k_{\inmeshlg{o}} \leftarrow i_{\inmeshlg{}},
\end{rcases} \quad \text{Expansion step}
\end{equation*}

followed by

\begin{equation*}
\begin{rcases}
\tilde{W}^e_{ij_{\inmesh{n}}} &= \sum_{\forall k_{\inmesh{}} \text{ s.t. }  k_{\inmesh{}} \rightarrow j_{\inmesh{n}}} {\hat{W}}^e_{ik_{\inmesh{}}}, \\
\tilde{W}^d_{i_{\inmesh{n}}j} &= \underset{\forall k_{\inmesh{}} \text{ s.t. }  k_{\inmesh{}} \rightarrow i_{\inmesh{n}}}{\operatorname{mean}} {\hat{W}}^d_{k_{\inmesh{}}j}, \\
\tilde{b}^d_{i_{\inmesh{n}}} &=  \underset{\forall k_{\inmesh{}} \text{ s.t. }  k_{\inmesh{}} \rightarrow i_{\inmesh{n}}}{\operatorname{mean}} {\hat{b}}^d_{k_{\inmesh{}}},
\end{rcases} \quad \text{Agglomerative step}
\end{equation*}

which concludes the proof.

§.§ Super- and sub-resolution error bound for GFN-ROM

We seek a bound on the following quantity
| u(x_i_n) - dec^ℳ_o →ℳ_n(map(μ))_i_n |, ∀i_n.
We can rewrite the expression and bound it as
| u(x_i_n) - dec^ℳ_o →ℳ_n(map(μ))_i_n |

=    | u(x_i_n) - ∑_j=1^L W̃^d_i_nj map(μ)_j + b̃^d_i_n |,

=    | u(x_i_n) - ∑_j=1^L ∀k_o  s.t  k_o i_nmean W^d_k_oj map(μ)_j + ∀k_o  s.t  k_o i_nmean b^d_k_o |,

=    | u(x_i_n) - ∀k_o  s.t  k_o i_nmean ( ∑_j=1^L W^d_k_oj map(μ)_j + b^d_k_o ) |

=    | u(x_i_n) - ∀k_o  s.t  k_o i_nmean dec^ℳ_o →ℳ_o(map(μ))_k_o |,

=    | ∀k_o  s.t  k_o i_nmean ( u(x_i_n) - u(x_k_o) + u(x_k_o) - dec^ℳ_o →ℳ_o(map(μ))_k_o ) |,

≤    ∀k_o  s.t  k_o i_nmean | u(x_i_n) - u(x_k_o) + u(x_k_o) - dec^ℳ_o →ℳ_o(map(μ))_k_o |,

≤    ∀k_o  s.t  k_o i_nmean ( | u(x_i_n) - u(x_k_o) | + | u(x_k_o) - dec^ℳ_o →ℳ_o(map(μ))_k_o | ),

≤    ∀k_o  s.t  k_o i_nmean ( δ+ τ),

=    δ+ τ.
This concludes the proof. For the multiple layer case, the result is also the same.

§.§ Super- and sub-resolution error bound for GFN-ROM mapper

We seek a bound on
| map(μ)_i - enc^ℳ_o →ℳ_n(_n)_i |, ∀i.
We can rewrite the expression as
| map(μ)_i - enc^ℳ_o →ℳ_o(_o)_i + enc^ℳ_o →ℳ_o(_o)_i - enc^ℳ_o →ℳ_n(_n)_i |

≤    | map(μ)_i - enc^ℳ_o →ℳ_o(_o)_i | + | enc^ℳ_o →ℳ_o(_o)_i - enc^ℳ_o →ℳ_n(_n)_i |,

≤    α+ | enc^ℳ_o →ℳ_o(_o)_i - enc^ℳ_o →ℳ_n(_n)_i |,

=    α+ | σ( ∑_j_o=1^N_o W_ij_o^e u(x_j_o) + b^e_i ) - σ(∑_k_n=1^N_n W̃_k_ni^e u(x_k_n) + b̃^e_i) |.
We assume a Lipschitz condition with Lipschitz constant $C$ for the activation function $\sigma$, which allows us to continue bounding as
≤    α+ C | ∑_j_o=1^N_o W_ij_o^e u(x_j_o) + b^e_i - (∑_k_n=1^N_n W̃_ik_n^e u(x_k_n) + b̃^e_i) |,

=    α+ C | ∑_j_o=1^N_o W_ij_o^e u(x_j_o) - ∑_k_n=1^N_n W̃_ik_n^e u(x_k_n) |,

=    α+ C | ∑_j_o=1^N_o W_ij_o^e u(x_j_o) - ∑_k_n=1^N_n   ∑_∀m_o  s.t  m_o k_n W^e_im_o/|{ h_n  s.t.  m_o    h_n }| u(x_k_n) |,

=    α+ C | ∑_j_o=1^N_o W_ij_o^e ( u(x_j_o) - ∑_∀k_n  s.t  j_o k_n u(x_k_n)/|{ h_n  s.t.  j_o    h_n }| ) |,

=    α+ C | ∑_j_o=1^N_o W_ij_o^e/|{ h_n  s.t.  j_o    h_n }| ∑_∀k_n  s.t  j_o k_n (  u(x_j_o) - u(x_k_n) ) |,

=    α+ C | ∑_j_o=1^N_o W_ij_o^e/|{ h_n  s.t.  j_o    h_n }| ∑_∀k_n  s.t  j_o k_n δ|,

≤    α+ δC ∑_j_o=1^N_o |W_ij_o^e|.
To obtain a bound independent of $i$, one can bound as
≤α+ δC ||^e||_∞.
This concludes the proof. For the multiple layer case, the result is
≤    α+ δC^P+1 ∑_l_P=1^L_P |W^(P)_i,l_P|  ∑_l_P-1=1^L_P-1 |W^(P-1)_l_P,l_P-1| ⋯∑_l_1=1^L_1 |W^(1)_l_2,l_1|   ∑_j_o=1^N_o |W_l_1,j_o^e|,

≤    α+ δC^P+1 ||^e||_∞∏_p=1^P ||^(p)||_∞.

§.§ Super- and sub-resolution error bound for GFN-ROM autoencoder

We seek a bound on
| u(x_i_n) - dec^ℳ_o →ℳ_n(enc^ℳ_o →ℳ_n(_n))_i_n |, ∀i_n.
To begin bounding this, we can rewrite the expression as
=    | u(x_i_n) - ∀k_o  s.t  k_o i_nmean u(x_k_o) .

+ ∀k_o  s.t  k_o i_nmean u(x_k_o) - ∀k_o  s.t  k_o i_nmean dec^ℳ_o →ℳ_o(enc^ℳ_o →ℳ_o(_o))_k_o

. + ∀k_o  s.t  k_o i_nmean dec^ℳ_o →ℳ_o(enc^ℳ_o →ℳ_o(_o))_k_o - dec^ℳ_o →ℳ_n(enc^ℳ_o →ℳ_n(_n))_i_n |,

≤    | u(x_i_n) - ∀k_o  s.t  k_o i_nmean u(x_k_o) |

+ | ∀k_o  s.t  k_o i_nmean u(x_k_o) - ∀k_o  s.t  k_o i_nmean dec^ℳ_o →ℳ_o(enc^ℳ_o →ℳ_o(_o))_k_o |

+ | ∀k_o  s.t  k_o i_nmean dec^ℳ_o →ℳ_o(enc^ℳ_o →ℳ_o(_o))_k_o - dec^ℳ_o →ℳ_n(enc^ℳ_o →ℳ_n(_n))_i_n |,

≤    δ+ β+ | ∀k_o  s.t  k_o i_nmean dec^ℳ_o →ℳ_o(enc^ℳ_o →ℳ_o(_o))_k_o - dec^ℳ_o →ℳ_n(enc^ℳ_o →ℳ_n(_n))_i_n |.
We consider the last term in more detail. We have
| ∀k_o  s.t  k_o i_nmean dec^ℳ_o →ℳ_o(enc^ℳ_o →ℳ_o(_o))_k_o - dec^ℳ_o →ℳ_n(enc^ℳ_o →ℳ_n(_n))_i_n |

=    | ∀k_o  s.t  k_o i_nmean ( ∑_j=1^L W^d_k_oj enc^ℳ_o →ℳ_o(_o)_j + b^d_k_o ) .

. - ∑_m=1^L W̃^d_i_nm enc^ℳ_o →ℳ_n(_n)_m - b̃^d_i_n |,

=    | ∀k_o  s.t  k_o i_nmean ( ∑_j=1^L W^d_k_oj enc^ℳ_o →ℳ_o(_o)_j + b^d_k_o ) .

. - ∑_m=1^L ∀k_o  s.t  k_o i_nmean ( W^d_k_om ) enc^ℳ_o →ℳ_n(_n)_m - ∀k_o  s.t  k_o i_nmean ( b^d_k_o ) |,

=    | ∀k_o  s.t  k_o i_nmean ∑_j=1^L W^d_k_oj ( enc^ℳ_o →ℳ_o(_o)_j -enc^ℳ_o →ℳ_n(_n)_j) |,

≤    ∀k_o  s.t  k_o i_nmean ∑_j=1^L | W^d_k_oj | | enc^ℳ_o →ℳ_o(_o)_j -enc^ℳ_o →ℳ_n(_n)_j |.
We assume a Lipschitz condition with Lipschitz constant $C$ for the activation function $\sigma$, allowing us to bound
| enc^ℳ_o →ℳ_o(_o)_i -enc^ℳ_o →ℳ_n(_n)_i |

=    | σ( ∑_j_o=1^N_o W_ij_o^e u(x_j_o) + b^e_i ) - σ(∑_k_n=1^N_n W̃_k_ni^e u(x_k_n) + b̃^e_i) |,

≤    C | ∑_j_o=1^N_o W_ij_o^e u(x_j_o) + b^e_i - (∑_k_n=1^N_n W̃_ik_n^e u(x_k_n) + b̃^e_i) |,

=    C | ∑_j_o=1^N_o W_ij_o^e u(x_j_o) - ∑_k_n=1^N_n W̃_ik_n^e u(x_k_n) |,

=    C | ∑_j_o=1^N_o W_ij_o^e u(x_j_o) - ∑_k_n=1^N_n   ∑_∀m_o  s.t  m_o k_n W^e_im_o/|{ h_n  s.t.  m_o    h_n }| u(x_k_n) |,

=    C | ∑_j_o=1^N_o W_ij_o^e ( u(x_j_o) - ∑_∀k_n  s.t  j_o k_n u(x_k_n)/|{ h_n  s.t.  j_o    h_n }| ) |,

=    C | ∑_j_o=1^N_o W_ij_o^e/|{ h_n  s.t.  j_o    h_n }| ∑_∀k_n  s.t  j_o k_n (  u(x_j_o) - u(x_k_n) ) |,

=    C | ∑_j_o=1^N_o W_ij_o^e/|{ h_n  s.t.  j_o    h_n }| ∑_∀k_n  s.t  j_o k_n δ|,

≤    δC ∑_j_o=1^N_o |W_ij_o^e|.
Substituting these expressions into the original bound, we obtain
≤    δ+ β+ δC ∀k_o  s.t  k_o i_nmean ∑_j=1^L | W^d_k_oj | ∑_m_o=1^N_o |W_jm_o^e|.
To obtain a bound independent of $i$, one can bound as
≤    δ+ β+ δC ||^d||_∞||^e||_∞.
This concludes the proof. For the multiple layer case, the result is
≤    β+ δ+ δC^Q+1 ∀k_o  s.t  k_o i_nmean ∑_j=1^L_Q+1 |W^d_k_oj| ∑_l_Q=1^L_Q |W^(Q)_j,l_Q| ⋯∑_l_1=1^L_1 |W^(1)_l_2,l_1|   ∑_j_o=1^N_o |W_j_ol_1^e|,

=    β+ δ+ δC^Q+1 ||^d||_∞||^e||_∞∏_p=1^Q ||^(p)||_∞.

§ TRAINING DETAILS

§.§ Network architecture and hyperparameter choices

Network architecture and hyperparameter choices for the MOR methods.

POD-NN     GCA-ROM     GFN-ROM
[HTML]E5E3E3 Bottleneck size ($N$)     [HTML]E5E3E3 $\lfloor 1.5\times N_\mu \rfloor$     [HTML]E5E3E3 $\lfloor 1.5\times N_\mu \rfloor$     [HTML]E5E3E3 $\lfloor 1.5\times N_\mu \rfloor$

Mapper weight ($\omega$)     -     10     10

[HTML]E5E3E3 Optimiser     [HTML]E5E3E3 Adam     [HTML]E5E3E3 Adam     [HTML]E5E3E3 Adam$^a$

Learning rate     $10^{-3}$     $10^{-3}$     $10^{-3}$

[HTML]E5E3E3 L2 regularisation     [HTML]E5E3E3 $10^{-5}$     [HTML]E5E3E3 $10^{-5}$     [HTML]E5E3E3 $10^{-5}$

Training epochs     5000     5000     5000

[HTML]E5E3E3 Autoencoder sizes     [HTML]E5E3E3 -     [HTML]E5E3E3 [$|\mathcal{M}|$, 200, $N$]     [HTML]E5E3E3 [$|\mathcal{M}|$, 200, $N$]

Train/test split     30/70     30/70     30/70

[HTML]E5E3E3Autoencoder activation     [HTML]E5E3E3 -     [HTML]E5E3E3 ELU$^b$     [HTML]E5E3E3 Tanh$^b$

Mapper sizes     [$N_\mu$, 50, 50, 50, 50, $N$]     [$N_\mu$, 50, 50, 50, 50, $N$]     [$N_\mu$, 50, 50, 50, 50, $N$]

[HTML]E5E3E3 Mapper activation     [HTML]E5E3E3 Tanh     [HTML]E5E3E3 Tanh$^c$     [HTML]E5E3E3 Tanh

* POD-G is also undertaken with $N= \lfloor 1.5\times N_\mu \rfloor$ degrees of freedom and a train/test split of 30/70.

a We employ the precomputed adaptive method so that we can use Adam.

b Activation not applied to the last decoder layer.

c Activation not applied to the last mapper layer.