# A Comparison of Neural Network Architectures for Data-Driven Reduced-Order Modeling

The popularity of deep convolutional autoencoders (CAEs) has engendered effective reduced-order models (ROMs) for the simulation of large-scale dynamical systems. However, it is not known whether deep CAEs provide superior performance in all ROM scenarios. To elucidate this, the effect of autoencoder architecture on its associated ROM is studied through the comparison of deep CAEs against two alternatives: a simple fully connected autoencoder, and a novel graph convolutional autoencoder. Through benchmark experiments, it is shown that the superior autoencoder architecture for a given ROM application is highly dependent on the size of the latent space and the structure of the snapshot data, with the proposed architecture demonstrating benefits on data with irregular connectivity when the latent space is sufficiently large.

## Authors

• 4 publications
• 11 publications
• 15 publications
• 20 publications
07/13/2020

### Drum Beats and Where To Find Them: Sampling Drum Patterns from a Latent Space

This paper presents a large dataset of drum patterns and compares two di...
06/04/2021

### Neural Network Surrogate Models for Absorptivity and Emissivity Spectra of Multiple Elements

Simulations of high energy density physics are expensive in terms of com...
06/21/2019

### Connectivity-Optimized Representation Learning via Persistent Homology

We study the problem of learning representations with controllable conne...
05/05/2021

### Novelty Detection and Analysis of Traffic Scenario Infrastructures in the Latent Space of a Vision Transformer-Based Triplet Autoencoder

Detecting unknown and untested scenarios is crucial for scenario-based t...
03/21/2019

### Deep Learning with Anatomical Priors: Imitating Enhanced Autoencoders in Latent Space for Improved Pelvic Bone Segmentation in MRI

We propose a 2D Encoder-Decoder based deep learning architecture for sem...
10/22/2021

### Deep Convolutional Autoencoders as Generic Feature Extractors in Seismological Applications

The idea of using a deep autoencoder to encode seismic waveform features...
01/25/2021

### Deep Learning-Based Autoencoder for Data-Driven Modeling of an RF Photoinjector

We adopt a data-driven approach to model the longitudinal phase-space di...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

High-fidelity computational models are indispensable tools for the prediction and analysis of physical systems which are governed by parameterized partial differential equations (PDEs). As more and more industries are relying on simulated results to inform their daily operations, the significant amount of computational resources demanded by such models is becoming increasingly prohibitive. Indeed, actions which increase model fidelity such as refining the spatio-temporal resolution can also lead to an explosion of dimensionality, making use of the full-order model (FOM) infeasible in real-time or many-query scenarios. To remedy this, emphasis has been placed on reduced-order models (ROMs) which approximate the high-fidelity, full-order models at any desired configuration of parameters. Indeed, an appropriately built ROM can enable important applications such as uncertainty quantification or predictive control, which respectively require access to the FOM solution at thousands of points or near-instantaneous access to approximate solutions.

The governing assumption inherent in the design of ROMs is that the FOM solution lies in a submanifold of the ambient space which is intrinsically low-dimensional. This is reasonable to posit when the FOM solution has been generated from the semi-discretization of an evolving quantity over a dense set of nodal points, as the values of this quantity in each dimension are strongly determined by the original (presumably lower-dimensional) dynamics. When this is the case, efficient approximation of the solution submanifold can lead to massive decreases in simulation time with relatively small approximation errors, since restriction to this space is lossless in principle. On the other hand, the problem of computing this submanifold is generally highly nontrivial, which has led to a wide variety of ROMs appearing in practice.

To enable this speedup in approximation of the FOM solution, ROMs are divided into offline and online stages. The offline stage is where the majority of expensive computation takes place, and may include operations such as collecting snapshots of the high-fidelity model, generating a reduced basis, or training a neural network. On the other hand, an effective ROM has an online stage that is designed to be as fast as possible, involving only a small amount of function evaluations or the solution of small systems of equations. Classically, most ROMs have relied on some notion of reduced basis to accomplish this task, which is designed to capture all of the important information about the high-fidelity system. For example, the well studied method of proper orthogonal decomposition (POD) described in Section 2 uses a number of solution snapshots to generate a reduced basis of size whose span has minimal reconstruction error. The advantage of this is that the high-fidelity problem can then be projected into the reduced space, so that only its low-dimensional representation must be manipulated at solving time. On the other hand, POD like many reduced-basis methods is a fundamentally linear technique, and therefore struggles in the presence of advection-dominated phenomena where the solution does not decay over time.

### 1.1. Related Work

Desire to break the accuracy barrier of linear methods along with the rapid advancement of neural network technology has reignited interest in machine learning methods for model order reduction. In particular, it was first noticed in Milano and Koumoutsakos

[milano2002]

that the reduced basis built with a linear multi-layer perceptron (MLP) model is equivalent to POD, and that nonlinear MLP offers improved prediction and reconstruction abilities at additional computational cost. Later, this idea was extended in works such as

[kashima2016, hartman2017] which established that a fully connected autoencoder architecture can be trained to provide a nonlinear analogue of projection which often realizes significant performance gains over linear methods such as POD.

Following the mainstream success of large-scale parameter sharing network architectures such as convolutional neural networks, the focus of nonlinear ROMs began to shift in this direction. Since many high-fidelity PDE systems such as large-scale fluid simulations involve a huge number of degrees of freedom, the ability to train an effective network at a fraction of the memory cost for a fully connected network is very appealing and led to several important breakthroughs. The work of Lee and Carlberg in

[lee2020] demonstrated that deep convolutional autoencoders (CAEs) can overcome the Kolmogorov width barrier for advection-dominated systems which limits linear ROM methods, leading to a variety of similar ROMs based on this architecture seen in e.g. [fukami2020, eivazi2020, fresca2021, maulik2021]. Moreover, works such as [eivazi2020, fresca2021] have experimented with entirely data-driven ROMs based on deep CAEs, and seen success using either fully connected or recurrent long short term networks to simulate the reduced low-dimensional dynamics. At present, there are a large variety of nonlinear ROMs based on fully connected and convolutional autoencoder networks, and this remains an active area of research.

Despite the now widespread use of autoencoder networks for reduced-order modeling, there have not yet been studies which compare the performances of fully connected and convolutional networks in a ROM setting. In fact, many high-fidelity simulations take place on domains which are unstructured and irregular, so it is somewhat surprising that CAEs still perform well for this type of problem (see e.g. [fresca2021, Remark 5]). On the other hand, there is clearly a conceptual issue with applying the standard convolutional neural network (c.f. Section 2) to input data with variable connectivity; since the convolutional layers in this architecture accept only rectangular data, the input must be reshaped in a way which can change its neighborhood relationships, meaning the convolution operation may no longer be localized in space. This is a known issue which has been thoroughly considered in the literature on graph convolutional neural networks (see e.g. [wu2020] and references therein), which have been successful in a wide range of classification tasks.

### 1.2. Contributions

The goal of this work is to augment the ROM literature with a comparison between autoencoder-ROM architectures based on fully connected and convolutional neural networks. In particular, we compare the most popular deep convolutional autoencoder architecture from e.g. [lee2020, fresca2021]

to two alternatives: a basic fully connected autoencoder with batch normalization

[ioffe2015], and a novel autoencoder based on the GCNII graph convolution operation of [chen2020]. Through benchmark experiments on regular and irregular data, we study the effectiveness of each architecture in compressing/decompressing high-fidelity PDE data, and evaluate its performance as part of a data-driven ROM. Results show that the superior architecture is highly dependent on the size of the reduced latent space as well as the modeling task involved. Moreover, some of the assumed advantages of the CAE (such as reduced memory consumption and easier training) are not satisfied in every case, and for some applications it is advantageous to use one of the other available architectures. It is our hope that the experiments provided here will enable modelers to make a more informed choice regarding which network architecture they employ for their autoencoder-based ROM.

## 2. Preliminaries

Consider a full-order model which is a dynamical system of parameterized ODEs,

 (1) ˙x(t,μ)=f(t,x(t,μ),μ),x(0,μ)=x0(μ).

Here is time (),

is the vector of model parameters,

is the parameterized state function with initial value , and is the state velocity. In practice, models (1) arise naturally from the semi-discretization of a parameterized system of PDE, where the discretization of the domain creates an -length vector (or array) representing a quantity defined at the nodal points. Clearly, a semi-discretization of size has the potential to create a large amount of redundant dimensionality, as a unique dimension is created for every nodal point while the solution remains characterized by the same number of parameters.

To effectively handle the high-dimensional ODE (1), the overarching goal of reduced-order models is to predict approximate solutions without solving the governing equations in the high-dimensional space . In principle, this is possible as long as the assignment is unique for each , meaning that the space of solutions spans a low-dimensional manifold in the ambient (c.f. Figure 1). Therefore, the ROM should provide fast access to the solution manifold

 S={x(t,μ)|t∈[0,T],μ∈D}⊂RN,

which contains integral curves of the high-dimensional FOM. Note that the uniqueness assumption on implies that the coordinates form a global chart for and hence is intrinsically -dimensional (at most), which is typically much less than the ambient dimension .

Traditionally, many ROMs are generated through some form of Proper Orthogonal Decomposition (POD). To apply POD, some number of solutions to (1) are first collected and stored in a snapshot matrix where is the number of training samples, , represents the component of the solution , and nonuniqueness is allowed in the parameters

. A singular value decomposition of the snapshot matrix

is then performed, and the first columns of (corresponding to the largest eigenvalues of ) are chosen to serve as a reduced basis for the solution space. Geometrically, this procedure amounts to choosing as basis an -plane of POD modes in which minimizes the projection error over the snapshot data. A ROM can then be generated from this by projecting (1) to this -plane. More precisely, writing the approximate solution for some and using the FOM along with the fact that leads to the system

 ˙^x(t,μ)=U⊤nf(t,Un^x,μ),^x(0,μ)=^x0(μ)=U⊤nx0(μ),

which requires the solution of only equations.

POD based techniques are known to be highly effective in many cases, especially when the problem dynamics are diffusion-dominated. On the other hand, advection-dominated problems exhibit a slowly decaying Kolmogorov n-width

 dn(S)=infMnsupx∈Sinf^x∈Mn∥x−^x∥,

meaning a very large number of POD modes is required for accurate approximation of their solutions. This reveals a fundamental issue with linear approximation spaces which motivated recent study into neural network based ROMs. By now, it is well known that the barrier posed by can be surmounted using deep convolutional autoencoders, which are now reviewed.

### 2.1. Convolutional Autoencoder Networks

Recall that a fully connected neural network (FCNN) is a function of the input vector such that at layer ,

 (2) xℓ\coloneqqfℓ(xℓ−1)=σℓ(Wℓxℓ−1+bℓ),

where

is a nonlinear activation function acting row-wise on its input,

is the weight matrix of the layer and is its bias. FCNNs are extremely expressive due to their large amount of parameters, but for the same reason they are typically difficult to train and incur a high memory cost. On the other hand, a convolutional neural network (CNN) can be considered a special case of FCNN where weights are shared and determined by a translation-invariant kernel function. More precisely, a convolutional neural network is a function of a multidimensional input such that at layer ,

 xℓ,i=σℓ(Cin∑j=1x(ℓ−1),j∗Wjℓ,i+bℓ,i),

where denotes the spatially-independent “channel dimension” of , is the convolution kernel at layer , and is a bias object which is constant for each fixed . Note that the kernel has dimension equal to that of the input plus one, since it depends on the shape of as well as the number of channels at layers . Here plays the role of the discrete convolution operator, defined for two-dimensional, one-channel arrays as

 (x∗W)αβ=∑γ,δx(sα+γ)(sβ+δ)w(L−1−γ)(M−1−δ),

where

is the stride length,

, and all Greek indices are compactly supported and begin at zero. Note from this definition of the discrete convolution that the ranges of in the output depend on the constants

as well as the shape of the input. To determine this more precisely, consider projecting the tensor

along some coordinate direction and computing the length of this projection in terms of the lengths of the corresponding projections of . If the projection of along this direction has length and the projection of has length , then it is straightforward to show [dumoulin2016, Relationship 5] that the output projection has length . Therefore, convolution will naturally decrease the dimensionality of the input data when

, creating a downsampling effect which is highly useful for techniques such as image compression. Note that this downsampling can be further controlled by padding the input

with an array of zeros, in which case the size of the output changes based on other calculable rules (see e.g. [dumoulin2016]).

Since the discrete convolution (2) decreases the dimensionality of the input data in a controlled way, it is reasonable to make use of it when constructing convolutional autoencoder (CAE) networks. In particular, the main idea behind CAE networks is to stack convolutional layers in such a way that the input is successively downsampled (or encoded) to a low-dimensional representation, and then upsampled (or decoded) back to its original form. Combined with an appropriate notion of reconstruction error to be minimized, this encourages the neural network to learn an approximation to the identity mapping which is useful for tasks such as compressed sensing and PDE approximation provided the decoder can be decoupled from the encoder . In the context of ROM, the notion of data downsampling serves as a nonlinear form of projection analogous to the truncated eigenvalue decomposition in POD. In principle, as more characteristic features are extracted by the convolutional layers in the network, less dimensionality is required to encode the original solution . Therefore, repeated convolution enables the most basic features of to be represented with a latent space of small dimension, which can be exploited by the network for accurate reconstruction of the solution manifold. Moreover, the assumption of uniqueness in the solution mapping guarantees that this idea is loss-less whenever the latent dimension is at least as large as .

On the other hand, if convolution is to be useful for ROM, it is necessary to have a way to upsample (or decode) the latent representation for . Although the convolution operation is clearly not invertible, it can be transposed in a way which is useful for this purpose. In particular, note that discrete convolution can be considered a matrix operation by lexicographically “unrolling” the inputs into a vector and similarly unrolling the kernel into a sparse matrix which acts on that vector. The operator then has a well defined transpose, which effectively broadcasts the (unrolled) low-dimensional input to a high-dimensional output which can be rolled back into its original shape. In practice, this operation is accomplished by simply exchanging the forward and backward passes of the relevant convolutional layer, as it can be shown (see [dumoulin2016, Section 4]) that this is operationally equivalent to computing the transposed convolution. Doing this several times in succession produces the required upsampling effect, and traditionally makes up the decoder part of the CAE structure.

### 2.2. Autoencoder-based ROMs

As mentioned before, using a well trained CAE in place of the linear reconstruction when designing a ROM can provide a huge benefit to accuracy. In particular, examples from [lee2020, fresca2021] show that deep CAEs can approach the theoretical minimum error with only a few latent variables, and enable the same (or better) accuracy as POD with a much smaller latent space. Moreover, autoencoder-based ROMs do not require any more information than standard POD. Indeed, developing a ROM using the decoder mapping means seeking which satisfies

 ~x(t,μ)=g(^x(t,μ)),

where , is a latent representation of the state , and is a nonlinear mapping from the latent space to the ambient . In particular, if can be constructed such that its image is the solution submanifold , then the comparatively small latent representation is sufficient for complete characterization of the state . As is typically trained using only the snapshots in , when this provides the potential for an accurate nonlinear mapping using the same information necessary for POD.

Using the approximation in the FOM (1), it follows that where denotes the derivative operator (i.e. Jacobian matrix) of . Therefore, a reduced order model involving only can be generated by assuming that the approximate state obeys the same dynamics as the original state . In particular, if and has full rank at every , it follows quickly that

 (3) ˙^x(t,μ)=g′(^x)+f(t,g(^x),μ),~x(0,μ)=~x0(μ)=g(^x(0,μ)),

where denotes the Moore-Penrose pseudoinverse of . This is a well posed low-dimensional dynamical system that is feasible to solve when and its derivatives are known. Note that the ROM (3) takes place on the tangent bundle to the solution manifold , so that is a low-dimensional integral curve in this space. Moreover, it is clear that choosing the mapping immediately recovers the standard POD ROM from before since is a linear mapping.

Solving the ROM (3) is often accomplished using Newton or quasi-Newton methods based on minimizing the residual (e.g. [lee2020, fu2018]), though recent work in [eivazi2020, fresca2021] has demonstrated that a purely data-driven approach can also be effective for end-to-end reduced-order modeling of PDEs without appealing to conventional numerical solvers. More precisely, instead of relying on a low-dimensional system of differential equations to compute the latent representation , it is reasonable to consider computing this quantity with another neural network which is fully connected and relatively simple. In this case, the autoencoder can be trained alongside the mapping through the joint loss

 (4) L(x,t,μ)=∥x−g∘h(x)∥2+∥^x(t,μ)−h(x)∥2.

The first term in is simply the reconstruction error coming from the autoencoder, while the second term encourages the latent representation to agree with the encoded state . This approach has the advantage of requiring only function evaluations in its online stage, at the potential cost of model generalizability. On the other hand, experiments in [fresca2021] as well as Section 4 of this work show that data-driven ROMs still produce competitive results in practice while maintaining a lower online cost.

It is important to keep in mind that the existence of the decoder mapping is independent of the numerical technique used to generate it, and therefore the use of the convolutional propagation rule (2) is a design decision rather than a necessity. In fact, at present there are an enormous number of different neural network architectures available for accomplishing both predictive and generative tasks, each with their own unique properties. Therefore, it is reasonable to question whether or not deep CAEs based on (2) provide the best choice of network architecture for ROM applications. The goal of the sequel is to introduce an alternative to the conventional CAE, along with some experiments which compare the performance of this CAE-ROM to other autoencoder-based ROMs involving fully connected and graph convolutional neural networks.

## 3. A Graph Convolutional Autoencoder for ROM

The convolutional propagation rule (2) comes with inherent advantages and disadvantages that become readily apparent in a ROM setting. One very positive quality of this rule is that if the convolutional filter is of size , then its associated discrete convolution is precisely -localized in space. This locality ensures that only the features of nearby nodes are convolved together, increasing the interpretability and expressiveness of the neural network. On the other hand, it is clear that (2) does not translate immediately to irregular domains, where each node has a variable number of neighbors in its -hop neighborhood. Because of this, bothersome tricks are needed to apply the CAE based on (2) to PDE data defined on irregular or unstructured domains. In particular, one popular approach (c.f. [fresca2021]) is to simply pad and reshape the solution. More precisely, each feature of can be flattened into a vector, zero-padded until it has length for some , and then reshaped to a square of size . This square is then fed into the CAE, generating an output of the same shape which can be flattened and truncated to produce the required feature of . While this procedure works surprisingly well in many cases, it has the potential to create a very high overhead cost due to artificial padding of the solution. Moreover, it fails to preserve the beneficial locality property of convolution, since nodes which were originally far apart may end up close together after reshaping. Remedying these issues requires a new autoencoder architecture based on graph convolutional neural networks, which will now be discussed.

### 3.1. Graph Convolutional Neural Networks

While the discrete convolution in (2) is in some sense a direct translation of the continuous version, it is only meaningful for regular, structured data. However, PDE simulations based on high-fidelity numerical techniques such as finite element methods frequently involve domains which are irregular and unstructured. In this case, developing an autoencoder-based ROM requires a notion of convolution which is appropriate for data with arbitrary connectivity. On the other hand, it is not immediately obvious how to define a useful notion of convolution when the domain in question is unstructured. Because discrete procedures typically obey a “no free lunch” scenario with respect to their continuous counterparts, it is an active area of research to determine which properties of the continuous convolution are most beneficial when preserved in the discrete setting. Since unstructured domains can be effectively modeled as mathematical graphs, this has motivated a large amount of activity in the area of graph convolutional neural networks (GCNNs). The next goal is to recall the basics of graph convolution and propose an autoencoder architecture which makes use of this technology.

To that end, consider an undirected graph with (weighted or unweighted) adjacency matrix . Recall that is symmetric with entries defined by

 aij={wij∃eij∈Econnectingviandvj0otherwise,

where are the weight values associated to each edge. Given this information, there are a wide variety of approaches to defining a meaningful convolution between two signals defined at the nodes of (see e.g. [wu2020]). One highly successful approach is to make use of the fact that convolution is just ordinary multiplication in Fourier space. In particular, recall the combinatorial Laplacian which can be computed as where is the diagonal degree matrix with entries . It is straightforward to check that is symmetric and positive semi-definite (see e.g. [bollobas2013]), therefore has a complete set of real eigenvalues

with associated orthonormal eigenvectors

known as the Fourier modes of . Since any signal can be considered an -length array of -vectors denoted

, this yields the discrete Fourier transform

where is the vectorized component function of , and each entry . The advantage of this notion is a straightforward extension of convolution to graphical signals through multiplication in frequency space. In particular, it follows that for

 fi∗gj=U(U⊤fi⊙U⊤gj),

where denotes the element-wise Hadamard product and is the inverse Fourier transform on . Note the usual convention that matrix products occur before the Hadamard product.

This expression provides a mathematically precise (though indirect) definition of spatial graph convolution through the graph Fourier transform, but carries the significant drawback of requiring matrix multiplication with the Fourier basis which is generically non-sparse. To remedy this, Defferrard et al. [defferrard2016] proposed considering matrix-valued filters which can be recursively computed from the Laplacian . In particular, consider a convolution kernel where is an order polynomial in the eigenvalues of the normalized Laplacian and is a vector of learnable coefficient parameters. Such kernels are easily computed if the eigenvalues are known, but do not obviously eliminate the need for multiplication by the Fourier basis . On the other hand, if can be computed recursively from , then for any signal the convolution can be computed in operations, making the spectral approach more feasible for practical neural network applications. As an added benefit of the work in [defferrard2016], note that order- polynomial kernels are also precisely -localized in space. In particular, the component of the kernel localized at vertex is given by

 (gθ(L)δi)j=(gθ(L))ij=K∑k=0θk(Lk)ij,

where is the nodal basis vector and whenever lies outside the -hop neighborhood of (c.f. [hammond2011, Lemma 5.2]). This shows that the frequency-based convolution maintains the desirable property of spatial locality when the filter is a polynomial in .

Further simplification to this idea was introduced by Kipf and Welling [kipf2016] who chose , and with a learnable parameter to obtain the first-order expression . By renormalizing the initial graph to include self-loops (i.e. adding identity to the adjacency and degree matrices), they proposed the graph convolutional network (GCN) which obeys the propagation rule

 (5) xℓ+1=σ(~PxℓWℓ),

where is a signal on the graph with channels, is a (potentially nonsquare) weight matrix containing the learnable parameters and

 ~P=~D−1/2~A~D−1/2=(D+I)−1/2(A+I)(D+I)−1/2.

Networks of this form are simple and fast to compute, but were later shown [wu2019] to be only as expressive as a fixed polynomial filter on the spectral domain of the self-looped graph with order equal to the depth of the network. While the propagation rule (5) displayed good performance on small-scale graphs, it is now known that simple GCNs are notorious for overfitting and oversmoothing the data, as their atomic operation is related to heat diffusion on the domain. This severely limits their expressive power, as GCNNs based on (5) are depth-limited in practice and cannot effectively reconstruct complex signals when used in a ROM setting.

On the other hand, many other graph convolutional operations are available at present, some of which were developed specifically to remedy this aspect of the original GCN. One such operation was introduced in [chen2020] and directly extends approach of [kipf2016]. In particular, the authors of [chen2020] define the “GCNII” or GCN2 network, whose propagation rule is

 (6) xℓ+1=σ[((1−αℓ)~Pxℓ+αℓx0)((1−βℓ)I+βℓWℓ)],

where

are hyperparameters and

is the component-wise ReLU function

. The terms in this operation are motivated by the ideas of residual connection and identity mapping popularized in the well known ResNet architecture (c.f.

[he2016]), and the authors prove that -layer networks of this form can express an arbitrary -order polynomial in the self-looped graph Laplacian. While the hyperparameters are typically constant in practice, it is noted in [chen2020] that improved performance is attained when decreases with increasing network depth. Therefore, the experiments in Section 4 use the recommended values of and where is a constant hyperparameter.

With this, it is now possible to use the propagation rule (6) as the basis for a convolutional autoencoder. Although spectral graph convolution does not produce any natural downsampling analogous to the CNN case, it is still reasonable to expect that a succession of convolutional layers will learn to extract important features of the data it accepts as input. Therefore, it should be possible to recover the decoder mapping necessary for ROM from this approach as well. Moreover, in contrast to the discrete convolution (2) the GCN2 operation is invariant to label permutation and remains -localized on unstructured data. Therefore, it is reasonable to suspect that a CAE-ROM based on (6) will exhibit superior performance when the PDE to be simulated takes place on an unstructured domain. The remainder of the work will investigate whether or not this is the case, beginning with a description of the GCNN autoencoder that will be used.

###### Remark 3.1.

Note that there are a wide variety of graph convolutional operations which can be used as the basis for a convolutional autoencoder on unstructured data. During the course of this work, many such operations were invesigated including the ones proposed in [bresson2017, gilmer2017, qi2017, ma2019, wang2019]. Through experimentation, it was determined that the graph convolution proposed in [chen2020] is the most effective for the present purpose.

### 3.2. Graph Convolutional Autoencoder Architecture

Recall that the standard CAE is structured so that each layer downsamples the input in space at the same rate that it extrudes the input in the channel dimension. In particular, it is common to use convolution kernels with a stride of two and “same” zero padding (c.f. [dumoulin2016, Section 4.4]), so that the spatial dimension of decreases by a factor of four at each layer. While this occurs, the channel dimension is simultaneously increased by a factor of four, so that the total dimensionality of remains constant until the last layer(s) of the encoder. The final part of the encoder is often fully connected, taking the simplified output of the convolutional network and mapping it directly to the encoded representation in some small latent space . To recover from its encoding, this entire block of layers is simply “transposed” to form the decoder mapping , so that the autoencoder comprises the mapping as desired. The learnable parameters of the mappings are occasionally shared (resulting in a “tied autoencoder”), though most frequently they are not.

To study whether autoencoder-based ROMs benefit from the core ideas behind GCNNs, we propose the ROM architecture displayed in Figure 2. Similar to [fresca2021], this ROM is entirely data-driven and is comprised of two neural networks: a specially designed CAE based on the operation (6) whose decoder is the required mapping , and a small-scale FCNN which simulates the low-dimensional mapping . In this way, the component networks are trained directly from snapshot data during the offline stage, so that the online stage requires only function evaluations. Recall that the input is not downsampled by the convolutional procedure (6), so the feature dimension is similarly held constant. Moreover, use of the GCN2 operation ensures that each convolutional layer is skip-connected to the first, i.e. each convolutional layer retains a fraction of the pre-convolutional input. Therefore, the proposed autoencoder is defended against the troublesome oversmoothing known to occur in deep GCNs. To evaluate the performance of our GCNN-ROM, it is compared directly to two analogous autoencoder-based data-driven ROMs which use the same prediction network

and loss function (

4) but different autoencoder architectures. The first is the standard deep CAE based on (2) and designed according to the common conventions of CAEs on regular data, and the second is a simple feed-forward autoencoder which uses batch-normalized fully connected layers throughout. The results of a comparison between these ROMs are given in the next Section.

## 4. Numerical Examples

This section details a performance comparison between ROMs using autoencoder architectures based on fully connected (FCNN), convolutional (CNN), and graph convolutional (GCNN) networks. Through examples of PDEs on structured and unstructured domains, ROM performance is evaluated on two distinct tasks. The first task measures how well each architecture can approximate PDE solutions at new parameter values (referred to as the prediction problem), and the second measures how well they can compress and decompress PDE data (referred to as the compression problem).

### 4.1. Data Preprocessing and Evaluation

Begin with a snapshot array each row of which is a solution at a particular parameter configuration found in the (consistently ordered) parameter matrix . We first choose integers such that and separate and row-wise into training and validation sets, taking care that no value of appears in both . After separation, the training data is normalized “column-wise” to the interval . Momentarily dropping the subscript , this means computing the transformation

 Ms↦Ms−MminMmax−MminSs↦Ss−SminSmax−Smin,

where (resp. ) are arrays containing the element-wise maxima and minima of (resp. over the sampling dimension (e.g. for every relevant ), and the division is performed element-wise. Once the training data has been normalized, the validation data is normalized in a similar fashion based on the arrays computed from the training set.

Once the relevant neural networks have been trained, performance is measured through the empirical relative reconstruction of the FOM solution over the validation set. In particular, denoting the reconstructed solution by , the empirical relative reconstruction error is computed as

 Rℓi(x,~x)=∑s∥xs−~xs∥i∑s∥xs∥i,

where . Note that the reconstructed solution for the compression problem while

for the prediction problem. In addition to the reconstruction error, it is helpful to report the memory required to store the trained models as well as the average wall-clock time per training epoch. While the memory cost is a property of the network, note that the wall-clock time is hardware and implementation dependent with significant fluctuation throughout the course of training. On the other hand, it provides a useful comparison (along with the plots of the loss) to estimate the difference in computational cost across autoencoder architectures. Pseudocode for the training procedure is given in Algorithm

1

. Note that an early stopping criterion is employed, so that training is stopped when the validation loss does not decrease over a predefined number of consecutive epochs. All experiments are performed using Python 3.8.10 running PyTorch 1.8.1

[torch] on an early 2015 model MacBook Pro with 2.7 GHz Dual-Core Intel Core i5 processor and 8GB of RAM. The implementation of GCN2 convolution is modified from source code found in the PyTorch Geometric package [fey2019].

### 4.2. Parameterized 1-D Inviscid Burger’s Equation

The first experiment considers the parameterized one-dimensional inviscid Burgers’ equation, which is a common model for fluids whose motion can develop discontinuities. In particular, let be a vector of parameters and consider the initial value problem,

 (7) wt+12(w2)x=0.02eμ2x,w(a,t,μ)=μ1,w(x,0,μ)=1,

where subscripts denote partial differentiation and represents the conserved quantity of the fluid for . As a first test, it is interesting to examine how the CNN, GCNN, and FCNN autoencoder-based ROMs perform on this relatively simple hyperbolic equation. First, note that this problem is easily converted to the form (1) by discretizing the interval , in which case the solution at each nodal point is , and is a vector for each of length equal to thee number of nodes . With this, the equations (7) can be discretized and solved for various parameter values using a simple upwind forward Euler scheme. In particular, given stepsizes , let denote the solution at time step and spatial node . Then, simple forward differencing gives

 wk+1i−wkiΔt+12(w2)ki−(w2)ki−1Δx =0.02eμ2x,

which expresses the value at in terms of its value at . For the experiment shown, the discretization parameters chosen are for and for . The parameters used for training are where and , while the validation parameters are the midpoints where and . A precise description of the autoencoder architectures employed can be found in Tables 1, 2, and 3. The initial learning rates for the ADAM optimizer are chosen as 0.0025 for in the GCNN case, 0.001 in the CNN case, and 0.0005 in the FCNN case. The early stopping criterion for this example is 200 epochs and the mini-batch size is 20.

As seen in Table 4, the results for both the prediction and compression problems are significantly different across architectures and latent space sizes. Indeed, the CNN autoencoder is consistently the most accurate on the prediction problem while also being the most time consuming and memory consumptive, and its performance does not improve with increased latent dimension. On the other hand, the performance of the GCNN-ROM increases quite a bit when the latent dimension is increased from 3 to 10, but remains significantly less accurate than the ROM based on CNN or FCNN. Interestingly, on the prediction problem the performance of all networks appears to be bottlenecked by the fully connected network which simulates the low-dimensional dynamics, as their accuracy does not increase when the latent dimension increases from to . This suggests that a ROM based on Newton or quasi-Newton methods may produce better results in this case.

Conversely, the compression data shows that all network architectures benefit from an increased latent space dimension when the goal is to encode and decode, with the GCNN autoencoder benefiting the most and the CNN autoencoder benefiting the least. Again, the GCNN architecture struggles when the latent space dimension is set to the theoretical minimum , however its performance exceeds CNN and becomes competitive with that of FCNN once the latent dimension is increased to 32, at half of the memory cost.

###### Remark 4.1.

Instead of increasing the latent dimension , it is reasonable to consider increasing the depth of the autoencoder and prediction networks. However, we find that increasing the depth of the GCNN autoencoder beyond does not significantly improve its performance on either the prediction or compression problems. On the other hand, the GCNN ROM benefits modestly from an increase in depth of the prediction network from 4 to 8 layers, although the CNN and FCNN based ROMs do not. Therefore, all prediction experiments use 8 fully connected layers in the GCNN case and 4 in the CNN and FCNN cases.

Figure 3 illustrates the results at various values of for the parameter . Clearly, the architecture based on CNN is best equipped for solving the prediction problem in this case. On the other hand, when it is evident that all architectures can compress and reconstruct solution data well, although the FCNN and GCNN yield the best results. Figure 4 shows the corresponding plots of the losses incurred during network training. Note that all architectures have comparable training and validation losses for the compression problem, while all but the GCNN overfit the data slightly on the prediction problem.

### 4.3. Parameterized 2-D Heat Equation

The next experiment involves the parameterized two-dimensional heat equation on a uniform grid. In particular, let , , and consider the problem

 (8) ut−Δu=0,u(0,y,t)=−0.5,u(1,y,t)=μ1cos(μ2y),u(x,y,0)=0,

where is the temperature profile. Multiplying the first equation by a test function and integrating by parts yields the standard weak equation

 ∫Ωutv+∫Ω∇u⋅∇v=0,forallv∈H10(Ω),

which can be discretized and solved using standard finite element techniques. In particular, choosing a finite element basis so that the discretized solution has coefficient vector yields the linear system

 Mut=−Su,

where are the symmetric mass and stiffness matrices

 (9) mij=∫Ωφiφj,sij=∫Ω∇φi⋅∇φj.

For the present experiment, is semi-discretized using bilinear -elements over a regular grid to yield of dimension , and the IFISS software package [elman2007] is used to assemble and solve the boundary value problem (8) using a backward Euler method with stepsize . The system is integrated for a time period of 4 seconds at various parameters . Since no classical solution to (8) exists at , solution snapshots are stored starting from . Some solutions along this heat flow contained in the validation set are displayed in Figure 6.

The parameters appearing in the boundary conditions are again chosen to lie in a lattice. In particular, the 9 values used for training are for and the 4 validation values are the lattice midpoints . A sample from the validation set along with its ROM reconstructions for the case is shown in Figure 6. The precise architectures used are displayed in Tables 1, 5, and 6. The initial learning rates are for the GCNN case, for the CNN case, and for the FCNN case. The early stopping criterion for network training is 100 epochs, and the mini-batch size is .

The results for this experiment are displayed in Table 7, which clearly illustrates the drawbacks of using the standard CNN architecture along with solution reshaping. Not only is the computation time increased due to the introduction of almost 2000 fake nodes, but the accuracy of the CNN ROM is significantly lower than the FCNN-ROM despite its comparable or larger memory cost. On the other hand, the architecture based on GCNN is both the least memory consumptive and the most accurate on each problem when . Again, results show that the GCNN-ROM is most effective on the compression problem, demonstrating significantly improved accuracy when the latent dimension is not too small. Conversely, if a latent dimension equal to the theoretical minimum is desired, then the ROM based on deep CNN is still superior, potentially due to its multiple levels of filtration with relatively large kernels. An example ROM reconstruction along with pointwise error plots is displayed in Figure 7. Interestingly, Figure 5 shows that all autoencoder networks exhibit significant overfitting on this dataset, suggesting that none of the models has sufficient generalizability to produce optimal results on new data.

### 4.4. Unsteady 2-D Navier-Stokes Equations

The final experiment considers the unsteady Navier-Stokes equations on a rectangular domain with a slightly off-center cylindrical hole, which is a common benchmark model for computational fluid dynamics (c.f [schafer1996]). The domain is pictured in Figure 8, and the governing system is

 ut−νΔu+(u⋅∇)u+∇p=0,∇⋅u=0,u|t=0=0,

subject to the boundary conditions

where are the velocity and pressure functions and is the viscosity parameter. Recall the Reynolds number

 Re=|u|meanLν,

where is the average speed along the flow and is the characteristic length of the flow configuration. The goal of this experiment is to see how well the ROM-autoencoder architectures can compress and predict solutions to these equations for a given parameter pair .

To generate the training and validation data in this case, three values of are selected corresponding to for . This yields solutions with corresponding . The flow is simulated starting from rest over the time interval , using Netgen/NGSolve with Taylor-Hood finite elements for spatial discretization and implicit-explicit Euler with for time stepping. Solution snapshots are stored once every time steps, yielding a total of snapshots per example. Note that domain is represented by an unstructured irregular mesh whose element density is higher near the cylindrical hole, yielding two solution features (the components of velocity) with dimension equal to the number of nodes. For the present experiments, the data corresponding to are used to train the ROM networks and the data correseponding to are used for validation. The architectures used in this case are given in Tables 1, 8, and 9, and the early stopping criterion is 50 epochs. The initial learning rate used to train each network is .

The results of this experiment are given in Table 10, making it clear that the trick of reshaping the nodal features in order to apply CNN carries an associated accuracy cost. Indeed, the performance of the CNN-ROM on the prediction problem is inferior to both the GCNN and FCNN-ROMs (provided the latent space is not too small) despite its much higher computational cost. Moreover, on this highly irregular domain the CNN-based architecture does not even offer benefits in the low-dimensional case, where the architecture based on FCNN is roughly 5 times faster while maintaining superior accuracy. As before, the GCNN-ROM struggles when is small but gains in performance as the latent dimension increases. In fact, results indicate that the GCNN is superior to the FCNN on the compression problem when while also being faster to train and an order of magnitude less memory consumptive. Although is far from the “best case scenario” of , it is still much smaller than the original dimension of and may be useful in a real-time or many-query setting. Interestingly, Figure 9 shows that the CNN-ROM does not overfit on the prediction or compression problems, indicating that the standard CNN may not be expressive enough to handle complicated irregular data that has been aggressively reshaped. Similarly, the GCNN-ROM does not overfit during the compression examples, although its performance scales with the dimension of the latent space. Figures 10, 11, and 12 provide a visual comparison of the ROM reconstructions and their pointwise error at . Finally, note that for this problem it may be infeasible in the case to generate an accurate prediction ROM using any method, since it is still an open question whether the mapping is unique in this case.

## 5. Conclusion

Neural network-based strategies for the reduced-order modeling of PDEs have been discussed, and a novel autoencoder architecture based on graph convolution has been proposed for this purpose. Through interesting benchmark problems, the performance of three fundamentally different autoencoder-based ROM architectures has been compared, demonstrating that each architecture has advantages and disadvantages depending on the task involved and the connectivity pattern of the input data. In general, results indicate that the proposed graph convolutional autoencoder is superior for compression and decompression tasks provided the dimension of the latent space is sufficiently large. Moreover, it seems that the standard CAE is highly performing at low values of but should be avoided when the model domain is irregular. Finally, the fully connected autoencoder exhibits strong overall performance at a much higher memory cost than the alternatives. Future work will examine how the proposed GCNN autoencoder can be combined with traditional PDE-based numerical solvers to generate solutions at unknown parameter configurations without the need for a fully connected prediction network, potentially eliminating the noted accuracy bottleneck arising from this part of the ROM.

## Acknowledgments

This work is partially supported by U.S. Department of Energy Scientific Discovery through Advanced Computing under grants DE-SC0020270 and DE-SC0020418.