1 Introduction↩︎

Over the past few decades, numerical simulations of physical systems have provided invaluable insight, significantly enhancing our understanding of complex phenomena. These systems are typically modeled using Partial Differential Equations (PDEs) for which closed-form solutions are rarely attainable. Instead, they must be solved through high-fidelity numerical frameworks such as Finite Difference (FD), Finite Element Methods (FEM), and Finite Volume (FV) techniques.

The complexity of these systems does not translate directly into a more difficult implementation of those numerical schemes, but rather arises from the challenge of managing diverse temporal and spatial scales. These scales, often differing by several orders of magnitude, do not evolve separately; they interact dynamically due to the couplings arising from the nonlinear terms, leading to phenomena such as the energy cascade observed in viscous dissipative systems like turbulence. To achieve accurate solutions, a simulation must resolve the entire spectrum of scales, from the largest structures down to the smallest dissipative eddies. This requirement results in discretized systems with a massive number of Degrees of Freedom (DoF).

Despite the continuous advancement of High-Performance Computing (HPC) facilities over the last thirty years [3], [4], certain realistic solutions remain virtually unattainable within reasonable time frames. A prominent example is the Direct Numerical Simulation (DNS) of an entire aircraft under cruise conditions, which remains computationally prohibitive even with modern HPC hardware. Even in scenarios where DNS is technically feasible, its application is often impractical for multi-query contexts such as structural parameter optimization or for real-time prediction models embedded in control systems.

In this context Reduced Order Models (ROMs) comes into play. Their objective is to speed up the numerical solutions, hence drastically decreasing the computational burden, at the cost of losing some degree of accuracy [5]. Within this framework, classical turbulence modeling approaches such as Large Eddy Simulation (LES) [6] and Reynolds-Averaged Navier-Stokes (RANS) [7] can be interpreted as specialized forms of model reduction. In these approaches, only the most energetic large-scale structures of the flow are explicitly resolved, while the influence of the unresolved small scales is accounted for through various closure models. This simplification enables the use of coarser computational grids, significantly reducing the required degrees of freedom and, consequently, the time-to-solution.

Dissipative systems, which arise from nonlinear PDEs, exhibit a notable property: their long-time dynamics reside on a simpler, finite-dimensional structure known as the attractor [8]. The attractor is an invariant subset with respect to the semigroup described by the solution map of the dissipative systems, thus encapsulating all essential information regarding the non-trivial dynamics of the system [9]. Given that this attractor is finite-dimensional, the long-term dynamics following transient behaviors can be effectively represented by a system of Ordinary Differential Equations (ODEs) [10]. This finite-dimensional aspect can be articulated through the concept of an inertial manifold, defined as a smooth finite-dimensional manifold \(\mathcal{M}\) that remains invariant under the semigroup associated with the solution operator. Additionally, the inertial manifold contains the global attractor and attracts all trajectories at an exponential rate [9]–[11]. For certain dissipative dynamical systems, such as the Navier-Stokes equations, reaction-diffusion equations, and damped hyperbolic equations, the existence of an inertial manifold has been established [12], [13]. The presence of these inertial manifolds thus lends rigorous theoretical support to the Manifold Hypothesis, which posits that the set of all solutions of a parametric PDE lies on, or can be well-approximated by, a low-dimensional smooth manifold \(\mathcal{M}\) embedded within the high-dimensional state space [9], [14]. Within this framework, Reduced Order Models (ROMs) seek to identify a mapping from the high-dimensional physical space, where the manifold is embedded, to a low-dimensional representation [5] commonly referred to in ROM terminology as the latent space [2], [15], [16].

Once a proper reduced representation is established, a latent dynamical model is required to evolve the low-dimensional state over time [17]. The low dimensionality of the latent space is what makes the time integration of the reduced system computationally cheaper, allowing for a much faster solution of the physical problem. Traditionally, Reduced Basis (RB) methods approximate this solution manifold as a linear combination of basis functions, typically derived via Proper Orthogonal Decomposition (POD), Dynamic Mode Decomposition (DMD) [18], or Greedy sampling algorithms.

However, linear RB methods face significant limitations, most notably the slow decay of the Kolmogorov \(n\)-width for advection-dominated problems [19]. In such cases, the linear subspace spanned by the reduced basis fails to efficiently represent the nonlinear manifold, requiring an excessively large number of basis functions to maintain acceptable accuracy and thereby diminishing the computational advantage of the reduction.

The advent of Machine Learning has significantly transformed the landscape of Model Order Reduction (MOR), providing powerful tools to circumvent the restrictions inherent in traditional Reduced Basis methods [20]. By utilizing specialized Neural Network architectures, one can establish nonlinear mappings that connect the high-fidelity full-order space with a low-dimensional latent representation. These nonlinear mappings facilitate superior information compression compared to linear frameworks, thereby addressing the challenges posed by the slow decay of the Kolmogorov \(n\)-width in advection-dominated regimes [20]. Currently, several architectures are prominently used for this purpose, most notably Graph Neural Networks (GNNs) [21], [22], fully connected networks (MLPs) [2], Convolutional Neural Networks (CNNs) [15], [23]–[25], and Variational Autoencoders (VAEs) [26]–[28].

Regarding the temporal evolution of the system, a wide range of neural architectures have been explored to approximate the latent dynamics. These range from traditional Recurrent Neural Networks (RNNs), such as Long Short-Term Memory (LSTM) and Gated Recurrent Units (GRU) [29]–[31], to more recent frameworks like Neural Ordinary Differential Equations (Neural ODEs) [32], [33], Koopman-based operators [34], [35] and Transformers [26], [27]. For a comprehensive overview of these methodologies and additional references, the reader is referred to [17].

The governing equations of the latent dynamics can be rigorously derived under the assumption that the solution set of the parametric PDE constitutes a differentiable manifold. This derivation leverages the properties of manifold parametrization via diffeomorphic mapping which, when restricted to the manifold, are differentiable, invertible, and have a differentiable inverse too. In this framework, if \(\mathbf{u} \in \mathcal{M}\) represents the high-fidelity state on the manifold \(\mathcal{M}\), and \(\boldsymbol{z} \in U_\alpha\) is its latent representation, the latent velocity \(\dot{\boldsymbol{z} \in U_\alpha}\) is formally linked to the physical velocity \(\mathbf{f}(\mathbf{u})\) through the Jacobian of the mapping.

While this framework is mathematically solid, a significant challenge arises in practice: in standard autoencoder architectures, the decoder \(D\) does not strictly serve as the exact inverse of the encoder \(E\) (i.e., \(E \circ D \neq I\)). This discrepancy introduces a fundamental structural inconsistency that can invalidate the mathematically derived relations governing the latent evolution. It is important to highlight that this issue, intrinsic to the approximation of non-linear mapping through neural networks, is absent in linear RB-POD methods, where the projection and reconstruction operators are duals by construction.

Although recent literature has proposed fully connected AEs designed to enforce the invertibility required for consistent manifold parametrization [1], [2], this kind of architecture faces significant practical hurdles. Specifically, they suffer from poor scalability in terms of the number of trainable parameters relative to the input vector dimension and exhibit low sensitivity to local topological features of the physical domain. These limitations motivate the transition towards convolutional architectures [36].

In this work, we extend the symmetry properties developed by [1], [2] to hybrid 1D convolutional autoencoders. The term hybrid refers to the structure of the network: while the encoder consists of traditional convolutional layers, similar to those previously utilized in convolutional AEs, the decoder is constructed using a fully connected layer with matrices parameterized by a user-defined number of tunable parameters. These hybrid convolutional AEs, combined with a simple and straightforward latent dynamical model, are employed to develop ROMs for 1D dynamical systems, which serve as classic academic test cases.

The remainder of this work is organized as follows. Section 2 establishes the theoretical foundation through a formal derivation of the latent governing equations. Building upon these principles, Section 3 details the construction of symmetric Convolutional Autoencoders designed to satisfy the representation consistency requirement. The integration of these architectures into a coherent Reduced Order Model is addressed in Section 4, which outlines the interaction between the AE and latent dynamical models alongside the decoupled training protocols. The practical utility of the proposed framework is evaluated in Section 5, which describes the numerical experiments and physical test cases selected for validation. Section 6 presents a comprehensive performance analysis, comparing the symmetric ROM architecture against standard model order reduction techniques. Finally, Section 7 summarizes the study’s findings and discusses current limitations while proposing directions for future extensions.

2 Methodology↩︎

2.1 Latent dynamics↩︎

In this section we derive the latent governing equation as a projection of the FOM equation system, by exploiting the basics of differential geometry. To do so, we adopt the mathematical framework and notation introduced by [17], which is specifically tailored for applications in nonlinear model order reduction.

Consider the system of ODEs resulting from the semi-discretization of a parametric PDE, where the parameters are collected in a vector \(\boldsymbol{\mu} \in \mathcal{P}\): \[\label{eq:32FOM32equation} \left\{\begin{array}{l} \dot{\mathbf{u}}_h(t ; \boldsymbol{\mu})=\mathbf{f}_h\left(t, \mathbf{u}_h(t ; \boldsymbol{\mu}) ; \boldsymbol{\mu}\right), \quad t \in\left(t_0, T\right], \\ \mathbf{u}_h\left(t_0 ; \boldsymbol{\mu}\right)=\mathbf{u}_{0, h}(\boldsymbol{\mu}). \end{array}\right.\tag{1}\] Here, the subscript \(h\) denotes a high-fidelity discretization method (e.g., FEM, FVM, or FD). The set of time-dependent, parametric solutions \(\mathbf{u}_h \in \mathbb{R}^{N_h}\) forms a solution manifold: \[S_h = \left\{ \mathbf{u}_h(t; \boldsymbol{\mu}) : t \in [t_0, T], \, \boldsymbol{\mu} \in \mathcal{P} \right\} \subset \mathbb{R}^{N_h}\] with \(N_h\) denoting the number of degrees of freedom of the discretized dynamical system. Assuming \(S_h\) to be a differentiable \(m\)-dimensional manifold with \(m = \dim(\mathcal{P}) + 1\) there exists a family of injective parameterizations \[\label{eq:32exact32decoder} \mathbf{x}_\alpha : U_\alpha \subset \mathbb{R}^m \to S_h \subset \mathbb{R}^{N_h},\tag{2}\] where each \(\mathbf{x}_\alpha\) maps a local coordinate \(\boldsymbol{z} \in U_\alpha\) to a point on the manifold \(S_h\). To provide a wider view, the map \(\mathbf{x}_\alpha\) could be approximated by the decoder \(D\) of an AE, while \(\mathbf{x}_\alpha ^{-1}\) is approximated by the encoder \(E\).

Let \(\beta : I \to \mathbf{x}_\alpha (U_\alpha)\) denote the trajectory traced on the manifold by the high-fidelity solution \(\mathbf{u}_h(t; \boldsymbol{\mu}^*)\) for a fixed parameter instance \(\boldsymbol{\mu}^*\). This trajectory is defined on an open interval \(I \subseteq (t_0, T]\) such that \(\beta(t) = \mathbf{u}_h(t; \boldsymbol{\mu}^*) \in \mathbf{x}_\alpha (U_\alpha)\) for all \(t \in I\). The tangent vector to \(\beta\) is the temporal derivative \(\dot{\beta}(t) = \dot{\mathbf{u}}_h(t; \boldsymbol{\mu}^*)\). If \(\gamma: I \to U_\alpha\) represents the preimage of \(\beta\) under the map \(\mathbf{x}_\alpha\), then \(\dot{\boldsymbol{z} }(t; \boldsymbol{\mu}^*)\) is the tangent vector to the curve \(\gamma\) at time \(t\). The structural relationship between \(\gamma\), \(\beta\), and \(\mathbf{x}_\alpha\) is given by the composition \(\beta = \mathbf{x}_\alpha \circ \gamma\). We denote by \(d\mathbf{x}_\alpha: T_{\boldsymbol{z}} U_\alpha \to T_{\mathbf{u}_h} \mathcal{M}\) the differential of the manifold parametrization. By the definition of the differential and the application of the chain rule, we obtain: \[\dot{\mathbf{u}}_h(t; \boldsymbol{\mu}^*) = \dot{\beta}(t) = d\mathbf{x}_\alpha (\gamma(t))[\dot{\gamma}(t)] =d\mathbf{x}_\alpha (\mathbf{z}(t; \boldsymbol{\mu}^*))[\dot{\mathbf{z}}(t; \boldsymbol{\mu}^*)].\] While this relation was derived for a specific trajectory associated with \(\boldsymbol{\mu}^*\) , the result generalizes to any parameter \(\boldsymbol{\mu}\in\mathcal{P}\). Applying the change of variables to the FOM in Equation 1 yields an implicit expression for the latent time derivative: \[d\mathbf{x}_\alpha (\mathbf{z}(t; \boldsymbol{\mu}))[\dot{\mathbf{z}}(t; \boldsymbol{\mu})] = \mathbf{f}_h\left(t, \mathbf{x}_\alpha (\boldsymbol{z}(t ; \boldsymbol{\mu})) ; \boldsymbol{\mu}\right).\] Deriving an explicit evolution law requires pre-multiplying both sides by the inverse of the differential of the manifold parametrization \((d \mathbf{x}_\alpha)^{-1}\). Diffeomorphism properties ensure that this inverse exists and is equivalent to the differential of the inverse of the manifold parametrization [37], namely \((d \mathbf{x}_\alpha)^{-1} = d(\mathbf{x}_\alpha^{-1})\). Projecting the high-fidelity governing equations into the latent space via such identity defines the evolution of the reduced state variables: \[\label{eq:32latent32equation} \dot{\boldsymbol{z}}(\boldsymbol{z},t ; \boldsymbol{\mu})=d \mathbf{x} _\alpha ^{-1} (\boldsymbol{z}(t ; \boldsymbol{\mu})) [\mathbf{f}_h\left(t, \mathbf{x}_\alpha (\boldsymbol{z}(t ; \boldsymbol{\mu})) ; \boldsymbol{\mu}\right)], \quad t \in I.\tag{3}\] The system described by Equation 3 evolves autonomously once the map \(\mathbf{x}_\alpha\) is established. Despite the reduced dimensionality of the latent state \(\boldsymbol{z}\), the integration process still necessitates the evaluation of high-dimensional vectors, including the reconstructed state \(\mathbf{x}_\alpha(\boldsymbol{z})\) and the full-order operator \(\mathbf{f}_h\). In practice, this requirement can make the direct integration of the latent equations more computationally demanding than the original high-fidelity system [25].

The ROM circumvents this bottleneck by approximating the right-hand side of the latent equations with a suitable NN (later named Latent Flow Network). While \(\mathbf{x}_\alpha\) remains a linear map and reduces to the standard POD basis when the manifold \(S_h\) is an affine space, non-linear manifolds require more sophisticated approximations. For such tasks, autoencoders have proven to be an ideal candidate for capturing complex, non-linear manifold structures.

3 AutoEncoders↩︎

AutoEncoders constitute a broad class of neural network architectures that learn, in an unsupervised manner, to compress data into a lower-dimensional representation. The encoder, denoted by \(E\), maps the input data to a latent space by discovering an appropriate reduced representation. The complementary component, the decoder \(D\), approximates the inverse of this mapping and reconstructs the input from its compressed latent code, ideally recovering the original high-dimensional data. Common variants of this framework include variational Autoencoders [38], denoising Autoencoders [39], and contractive Autoencoders [40]. Architecturally, these models are distinguished primarily by the types of layers they employ: convolutional AutoEncoders are built from convolutional layers, whereas fully connected AutoEncoders are composed of dense (fully connected) layers.

Owing to their high expressive power, AEs can be employed to tackle a broad range of tasks. In many applications, however, it is desirable to bias the AE toward solutions of a particular type, even when alternative representations achieve the same reconstruction error. One simple example is the reduction of an incompressible fluid velocity field, for which it is important not only that the reconstructed velocity looks similar to the original field, but also that the divergence free condition is preserved. Such preferences can be weakly enforced by modifying the loss function in an ad hoc manner [41] or strongly by carefully constraining the architecture of the AE [42]. Otto et al. [1] were the first to introduce AEs that are intrinsically representation-consistent, in the sense that any latent point remains invariant under the composition of decoder and encoder. By imposing a specific structural constraint on the AE, they effectively realized a nonlinear projector onto the slow manifold along which the reduced dynamics evolve. A subsequent theoretical analysis of these architectures was provided by Brivio et al. [2], who derived upper and lower bounds on the reconstruction error in terms of the Lipschitz constants of the activation functions and POD error bounds at each encoder/decoder layer stage. Both studies [1], [2] were restricted to fully connected AE architectures.

Let us consider an AE composed by an encoder \(E\) and a decoder \(D\). We assume that the encoder and decoder comprise the same number of convolutional layers, and we denote this common depth by \(n_d\). In order to preserve symmetry (\(E \circ D = id\)) the following conditions are needed [2]:

• Bi-lipschitz activation function.

• Orthogonal / Bi-orthogonal linear layer pairs.

Let layer \(1\) denote the first layer of the encoder and layer \(n_d\) denote the last layer. The opposite indexing applies to the decoder, as schematically shown in Figure 1. Since the decoder is a mirrored “copy” of the encoder, a relation must hold between the last encoder layer and the last decoder layer (according to the enumeration given in Figure 1), the first encoder layer and the first decoder layer, and, more generally, between the \(i\)-th encoder layer and the \(i\)-th decoder layer. If the layers are fully connected and represented by the weight matrices \[E_i \in \mathbb{R}^{n_{i+1} \times n_i}, \quad D_{i} \in \mathbb{R}^{n_i \times n_{i+1}},\] with \(n_{i+1} \le n_i\), then bi-orthogonality means \[E_i D_{i} = \mathbb{I} \in \mathbb{R}^{n_{i+1} \times n_{i+1}}. \label{eq:32bi-ortho32condition}\tag{4}\]

Several strategies have been proposed to enforce the condition defined in Equation 4 . For instance, Brivio and Franco [2] implemented a matrix decomposition approach where the encoder and decoder weight matrices are expressed as a product of five distinct matrices belonging to specific classes, such as diagonal, orthogonal, and general matrices. Alternatively, Otto et al. [1] pursued a different methodology by defining the decoder weight matrices as the generalized inverse of the corresponding encoder matrices. It is important to note that while the Moore-Penrose pseudoinverse is unique, the generalized inverse is not, providing additional degrees of freedom in the architectural design.

If \(E_i\) is an orthogonal matrix, the orthogonality condition implies that \(D_{i} = E_i^T\), which automatically satisfies Equation 4 .

While the theoretical significance of symmetric fully connected autoencoders is considerable, their practical utility is often constrained by inherent limitations. In contrast, convolutional neural networks and autoencoders have emerged as the standard for large-scale applications. With the goal of creating symmetric AEs that can make a meaningful impact in practical applications, we are taking a step forward by developing 1D convolutional symmetric AEs.

3.1 Convolutional Symmetric AutoEncoders↩︎

In this section, we extend the notion of symmetric AEs to convolutional architectures. For the sake of simplicity, we will assume that there are no biases involved. The inclusion of bias is straightforward, as explained in previous works on symmetric AEs [1], [2].

The convolution of a function against a kernel is a linear operator. Let \(\boldsymbol{k} \in \mathbb{R}^{2n_k + 1}\) denote the (discrete) kernel: \[\boldsymbol{k} = [k_{-n_k}, k_{-n_k +1}, ..., k_0, ... ,k_{n_k -1}, k_{n_k}].\] In the discrete form, the convolution of a vector \(\boldsymbol{x}\) with a kernel \(\boldsymbol{k}\) of size \(2n_k + 1\) can be expressed, away from the head and tail of the vector \(\boldsymbol{x}\), as \[z_i = (\boldsymbol{x} \star \boldsymbol{k})_i = \sum _{j =- n_k} ^{n_k} x_{i + j} k_{j}.\] Assuming periodic (circular) boundary conditions and no stride, the matrix \(\mathbf{C}\) representing the convolution is the so-called circulant matrix. For the sake of clarity, we start with the simple case of one channel input-output layer. Writing the convolution in matrix notation \[\boldsymbol{z} = \mathbf{C} \boldsymbol{x}\] where \(\mathbf{C} \in \mathbb{R}^{n_i \times n_i}\) is a circulant matrix defined as \[c_{i,j} = \begin{cases} k_{j-i}, \quad |j-i|\le n_k \\ 0, \qquad \,|j-i|> n_k \\ \end{cases}.\] This formulation can be extended to the case where a stride greater than one is applied, while retaining, for simplicity, the assumption of periodic padding. Let \(n_{s}\) denote the stride length. \[c_{i,j} = \begin{cases} k_{j-i}, \quad |j-i|\le n_k \quad \text{and} \quad \lfloor (i-1) / n_s \rfloor = 0 \\ 0, \qquad \, \text{elsewhere}\\ \end{cases}.\] Given a vector \(\boldsymbol{x} \in \mathbb{R}^{n_i }\). Convolving it with \(\boldsymbol{k}\) using a stride of \(n_s\) produces an output vector \(\boldsymbol{z} \in \mathbb{R}^{\tfrac{n_i}{n_s}}\).

If we replace the convolution operation with a matrix–vector multiplication, the result is a vector of dimension \(n_i\), where only \(n_i / n_s\) entries are nonzero. \[\boldsymbol{z}' = \mathbf{C} \boldsymbol{x} \in \mathbb{R}^{n_i }.\] Specifically, \[z'_i =0 , \quad \text{if} \quad \lfloor (i-1) / n_s \rfloor \ne 0.\] We can build a suitable permutation matrix to push on the bottom the zero rows of the matrix \(\mathbf{C}\): \[\mathbf{P} \begin{pmatrix} \mathbf{U} \\ \mathbf{0} \\ \end{pmatrix} \boldsymbol{x}= \boldsymbol{z}',\] with \(\mathbf{U} \in \mathbb{R}^{n_i/n_s \times n_i}\). We can consequently write: \[\label{eq:32single32LayOrigin} \begin{pmatrix} \mathbf{U} \\ \mathbf{0} \\ \end{pmatrix} \boldsymbol{x}= \mathbf{P}^T \boldsymbol{z}',\tag{5}\] where \(\mathbf{P}^T \boldsymbol{z}'\) has the last \(n_i - n_i/n_s\) rows equal to zero. Since \(\mathbf{U}\) is a rectangular matrix, it admits the generalized right inverse \(\mathbf{U}^+\), which is a matrix so that \[\mathbf{U} \mathbf{U} ^ + = \mathbb{I} \in \mathbb{R}^{n_{i+1} \times n_{i+1}}.\] Now we introduce the decoder weight matrix defined as \(D = \left(\mathbf{U}^+ , \,\, \mathbf{0}^T\right) \mathbf{P}^T\), which gives as output the vector \(\boldsymbol{x}'\) \[\boldsymbol{x}'= D\boldsymbol{z}' = \left(\mathbf{U}^+,\,\, \mathbf{0}\right) \mathbf{P}^T \boldsymbol{z}'.\] Notice that in general \(\boldsymbol{x}' \ne \boldsymbol{x}\). This peculiar choice of the decoder guarantees that the non-zero elements of the vector \(\boldsymbol{z}'\) are mapped to themselves through the matrix \(E D\), which, expanded, reads as \[E D = \mathbf{P} \begin{pmatrix} \mathbf{U} \\ \mathbf{0} \\ \end{pmatrix} \left(\mathbf{U}^+ , \,\, \mathbf{0}^T\right) \mathbf{P}^T = \mathbf{P} \begin{pmatrix} \mathbf{U}\mathbf{U}^+ & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \\ \end{pmatrix} \mathbf{P}^T = \mathbf{P} \begin{pmatrix} \mathbb{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \\ \end{pmatrix} \mathbf{P}^T.\] Indeed, clustering in \(\bar{\boldsymbol{z}} =( \mathbf{P}^T \boldsymbol{z}')_{1:n_i/n_s} = (\mathbf{P}_{:, 1:n_i/n_s})^T \boldsymbol{z}'\) all the non-zero elements of \(\boldsymbol{z}'\), it is clear that \((\mathbf{P}_{:, 1:n_i/n_s})^T E_i D_{i}\mathbf{P}_{:, 1:n_i/n_s}\) maps \(\bar{\boldsymbol{z}}\) in \(\bar{\boldsymbol{z}}\) and that \((\mathbf{P}_{:, 1:n_i/n_s})^T E_i D_{i}\mathbf{P}_{:, 1:n_i/n_s} = \mathbb{I} \in \mathbb{R}^{n_{i+1} \times n_{i+1}}\).

The extension to convolutional layers with multiple input and output channels is achieved by employing a specific mathematical reformulation that simplifies the formal derivation. To ease the notation, the layer index will be dropped. The size of the \(j\)-th input channel vector \(\boldsymbol{x}_j\) is be denoted by \(n\) and \(l\) is the number of non zero elements in the \(i\)-th output channel vector \(\boldsymbol{z}_i\). Indeed, let us start by assuming \(n_{in}^{ch}, n_{out}^{ch}>1\), and \(n_{out}^{ch} \ge n_{in}^{ch}\), where \(n_{out}^{ch}\) and \(n_{in}^{ch}\) denote the number of output and input channels for a given convolutional layer. The \(i\)-th output channel vector is the sum of the convolved input channels: \[\boldsymbol{z}_i' = \sum _ {j=1} ^{n^{ch}_{in}} \mathbf{C} _{i, j} \boldsymbol{x} _j,\] \(i = 1, ..., n_{out}^{ch}\), and \(\boldsymbol{z}_i ' \in \mathbb{R} ^ {n}\), \(\boldsymbol{x}_j \in \mathbb{R} ^ {n}\). \(\mathbf{C} _{i, j}\) is a single convolution matrix, which represents a convolution operation. Since \(j = 1, ..., n_{in}^{ch}\), a single pass through a multi channels layer involves \(n_{in}^{ch} \times n_{out}^{ch}\) convolutions. Following the same notation of the single channel, we can write \[\boldsymbol{z}_i' = \sum _ {j=1} ^{n^{ch}_{in}} \mathbf{P} \begin{pmatrix} \mathbf{U}_{i, j} \\ \mathbf{0} \\ \end{pmatrix} \boldsymbol{x} _j \quad \mathrm{and} \; \mathrm{so} \quad \mathbf{P} ^ T \boldsymbol{z}_i' = \sum _ {j=1} ^{n^{ch}_{in}} \begin{pmatrix} \mathbf{U}_{i, j} \\ \mathbf{0} \\ \end{pmatrix} \boldsymbol{x} _j,\] with \(\mathbf{U}_{i, j} \in \mathbb{R} ^ {l \times n }\). Assume it is possible to write \[\boldsymbol{x}_j' = \sum _ {k=1} ^{n^{ch}_{out}} \left(\mathbf{U}_{j, k}^+ , \,\, \mathbf{0}^T\right) \mathbf{P}^T \boldsymbol{z}_k',\] where \(\mathbf{U}_{i, j} ^ + \in \mathbb{R} ^ {n \times l}\) then \[\begin{align} \hat{\boldsymbol{z}}_i &= \sum _ {j=1} ^{n^{ch}_{in}} \mathbf{P} \begin{pmatrix} \mathbf{U}_{i, j} \\ \mathbf{0} \\ \end{pmatrix} \boldsymbol{x}' _j \\ & = \sum _ {j=1} ^{n^{ch}_{in}} \mathbf{P} \begin{pmatrix} \mathbf{U}_{i, j} \\ \mathbf{0} \\ \end{pmatrix} \sum _ {k=1} ^{n^{ch}_{out}} \left(\mathbf{U}_{j, k}^+ , \,\, \mathbf{0}^T\right) \mathbf{P}^T \boldsymbol{z}_k' \\ & = \sum _ {j=1} ^{n^{ch}_{in}} \sum _ {k=1} ^{n^{ch}_{out}} \mathbf{P} \begin{pmatrix} \mathbf{U}_{i, j} \\ \mathbf{0} \\ \end{pmatrix} \left(\mathbf{U}_{j, k}^+ , \,\, \mathbf{0}^T\right) \mathbf{P}^T \boldsymbol{z}_k' \\ & = \sum _ {j=1} ^{n^{ch}_{in}} \sum _ {k=1} ^{n^{ch}_{out}} \mathbf{P} \begin{pmatrix} \mathbf{U}_{i, j}\mathbf{U}_{j, k}^+ & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \\ \end{pmatrix} \mathbf{P}^T \boldsymbol{z}_k'. \end{align}\] To have the identity, we have to impose the identity \(\hat{\boldsymbol{z}}_i = \boldsymbol{z}_i'\), leading to: \[\sum _ {j=1} ^{n^{ch}_{in}} \mathbf{U}_{i, j}\mathbf{U}_{j, i}^+ = \begin{cases} \mathbb{I} \quad \mathrm{for} \quad k= i\\ \mathbf{0} \quad \mathrm{for} \quad k \ne i \end{cases}.\] We can get rid of the double sum in all previous relations if we consider the matrices \(\mathbf{U}_{j, i}\) and \(\mathbf{U}_{j, i} ^ +\) as blocks of a bigger matrix \(\overline{\overline{\mathbf{U}}} \in \mathbb{R}^{(l \,n^{ch}_{out})\times (n\, n^{ch}_{in})}\) and \(\overline{\overline{\mathbf{U}}}^+\in \mathbb{R}^{(n\, n^{ch}_{in}) \times (l \, n^{ch}_{out})}\) which are respectively defined as \[\overline{\overline{\mathbf{U}}} = \begin{pmatrix} \mathbf{U}_{1,1} & \mathbf{U}_{1,2} & \cdots & \mathbf{U}_{1,n^{ch}_{in}} \\ \mathbf{U}_{2,1} & \mathbf{U}_{2,2} & \cdots & \mathbf{U}_{2,n^{ch}_{in}} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{U}_{n^{ch}_{out},1} & \mathbf{U}_{n^{ch}_{out},2} & \cdots & \mathbf{U}_{n^{ch}_{out},n^{ch}_{in}} \end{pmatrix} \quad \mathrm{and} \quad \overline{\overline{\mathbf{U}}}^+ = \begin{pmatrix} \mathbf{U^+}_{1,1} & \mathbf{U^+}_{1,2} & \cdots & \mathbf{U^+}_{1,n^{ch}_{out}} \\ \mathbf{U^+}_{2,1} & \mathbf{U^+}_{2,2} & \cdots & \mathbf{U^+}_{2,n^{ch}_{out}} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{U^+}_{n^{ch}_{in},1} & \mathbf{U^+}_{n^{ch}_{in},2} & \cdots & \mathbf{U^+}_{n^{ch}_{in},n^{ch}_{out}} \end{pmatrix}.\] To maintain the consistency it is also necessary to concatenate all the input vectors together, as well as the output vectors. Capital letters denote the vectors obtained by concatenation of the corresponding lowercase letters. For instance: \[\boldsymbol{X} = [\boldsymbol{x}_1^T, \boldsymbol{x}_2^T, ..., \boldsymbol{x}_{n_{in}^{ch}}^T]^ T \in \mathbb{R} ^{n\, n^{ch}_{in}},\] and \[\boldsymbol{Z}' = [(\mathbf{P}^T\boldsymbol{z'}_1^T)_{1:l}, (\mathbf{P}^T\boldsymbol{z'}_2^T)_{1:l}, ..., (\mathbf{P}^T\boldsymbol{z'}_{n_{out}^{ch}})_{1:l}^T]^ T \in \mathbb{R} ^{l \, n^{ch}_{out}}.\]

Under this formulation, the case of one input and one output channel is a special case of the multiple input-output channels. The condition for the existence of the generalized right inverse (\(\overline{\overline{\mathbf{U}}}^+\)) ties together the number of input-output channels and the dimension of the input-output vectors. Specifically, \(n^{ch}_{out} \,\, l \le n\,\, n^{ch}_{in}\).

3.2 Building the generalised right inverse↩︎

Several strategies can be employed to satisfy the bi-orthogonality condition defined in Equation 4 . A preliminary approach, following the work of [1], involves parameterizing the generalized inverse as: \[\overline{\overline{\mathbf{U^+}}} = \overline{\overline{\mathbf{D}}}(\overline{\overline{\mathbf{U}}} \,\overline{\overline{\mathbf{D}}})^{-1},\] where \(\overline{\overline{\mathbf{D}}}\) is a trainable weight matrix. While a judicious construction of \(\overline{\overline{\mathbf{D}}}\) allows for a low-rank parameterization of the decoder, numerical experiments indicated that this formulation suffers from poor training convergence. Specifically, the loss function frequently stagnated, or the model reached a pathological state where the weights failed to represent a physically meaningful manifold.

To overcome these limitations, we propose an alternative formulation based on orthogonal projectors. In this framework, the decoding stage is not restricted to a linear mapping but is instead decomposed into a linear component and a non-linear term: \[\label{eq:32projecion32layer} \boldsymbol{Z}'_o = D(\boldsymbol{Z}'_i) = \underbrace{\overline{\overline{\mathbf{U}}} ^{P} \boldsymbol{Z}'_i }_{\text{Linear}}+ \underbrace{(I - \overline{\overline{\mathbf{U}}} ^{P} \overline{\overline{\mathbf{U}}})\boldsymbol{f}(\boldsymbol{Z}'_i)}_{\text{Non Linear}}\tag{6}\] where \(\overline{\overline{\mathbf{U}}} ^{P}\) denotes the Moore-Penrose pseudoinverse. Here, the non-linear contribution is governed by a function \(\boldsymbol{f}\), which in our implementation consists of classical deconvolutional layers. This architecture offers several theoretical advantages:

• If \(\boldsymbol{f}\) is the null function, the decoder becomes parameter-free, reducing to a Symmetric Orthogonal Autoencoder (SOAE) as discussed in [2].

• When \(\boldsymbol{f}(\boldsymbol{z}'_i) \neq \boldsymbol{0}\), it provides the decoder with additional degrees of freedom to capture complex features of the physical solution without violating the bi-orthogonality constraint.

• The non-linear term is designed to lie within the null space of the encoder. Consequently, the identity \(E \circ D = id\) is preserved regardless of the complexity of \(\boldsymbol{f}\). Keeping in mind that \(\overline{\overline{\mathbf{U}}} \overline{\overline{\mathbf{U}}} ^{P} = \mathbb{I}\), this is formally demonstrated by: \[\begin{gather} \boldsymbol{Z} = E(\boldsymbol{X}) = \overline{\overline{\mathbf{U}}} \boldsymbol{X} \\ \boldsymbol{X}' = D(\boldsymbol{Z} ) \end{gather}\] \[\boldsymbol{Z}' = E(\boldsymbol{X}') = E( D(\boldsymbol{Z} )) = \underbrace{\overline{\overline{\mathbf{U}}} \overline{\overline{\mathbf{U}}} ^{P} }_{\mathbf{I}} \boldsymbol{Z} + \underbrace{\overline{\overline{\mathbf{U}}} (I - \overline{\overline{\mathbf{U}}} ^{P} \overline{\overline{\mathbf{U}}})}_{\boldsymbol{0}} \boldsymbol{f}(\boldsymbol{Z}) = \boldsymbol{Z}.\] This ensures that any point decoded from the latent space is perfectly preserved upon re-encoding, providing a consistent and stable coordinate system for the reduced-order model.

3.3 Activation functions↩︎

As previously discussed, the selection of appropriate activation functions is critical for maintaining a globally representation consistent AE. To ensure that the composition of the encoder \(E\) and decoder \(D\) yields the identity mapping across multiple layers, each activation function must be invertible. Let \(\sigma\) denote a Bi-Lipschitz, invertible activation function that operates element-wise. The global consistency of the network is preserved through the following chain of compositions: \[E \circ D = E_1 \circ \sigma \circ \cdots \circ \sigma \circ E_n \circ D_n \circ \sigma^{-1} \circ \cdots \circ \sigma^{-1} D_1 = E_1 \circ D_1 = id.\] It is clear from the expanded chain that the overall composition results into the identity. The invertible activation functions already existing in the literature are the LeakyReLU [2] and the Hyperbolic activation [1], [2]. We enlarge this collection by proposing a new invertible and Bi-Lipschitz activation function resembling the hyperbolic tangent that we name Rational activation function: \[\label{eq:32rational32AF} y = \sigma_\alpha(x) = \alpha \left( x + \frac{x}{1+x^2}\right),\tag{7}\] where \(\alpha\) accounts as hyperparameter. The inverse is analytic and can be written as:\[x = \sigma ^{-1}_\alpha(y) = \left( - \frac{q}{2} + \sqrt{\Delta} \right)^{1/3} + \left( - \frac{q}{2} - \sqrt{\Delta} \right)^{1/3},\] with \[\begin{align} t&=\frac{y}{\alpha}, & p &=2 -\frac{t^2}{3}, \\ q&=-\frac{2t^3+9t}{27}, & \Delta &= \left(\frac{q}{2}\right)^2 + \left(\frac{p}{3}\right)^2 .\\ \end{align}\]

The upper and lower Lipschitz constants of the rational activation function are \[L = 2\alpha \quad \mathrm{and} \quad \eta = \alpha \frac{7}{8}\] such that \[\eta\|x - y \| \le \|\sigma_\alpha(x) - \sigma_\alpha(y)\| \le L \| x - y \|, \qquad \forall x, y \in \mathbb{R}.\] We conduct all the numerical experiments with \(\alpha = 0.5\). A plot of the rational activation function and its inverse, for the given value of the hyperparameter \(\alpha\) is given in Figure 2.

4 Data-driven approach↩︎

Having introduced the enhanced version of convolutional autoencoders, which are distinguished by their advantageous parameterization properties, we can now employ them for manifold learning tasks to obtain a reduced representation of high-dimensional data. Additionally, we can leverage the universal approximation capabilities of feedforward neural networks to model the function that defines the latent flow described in Equation 3 .

To construct the reduced representation of the dynamical systems governed by Equation 1 , we integrate the aforementioned architectures within a decoupled two-stage training framework. In the first stage, the autoencoder is trained on the available snapshot data to minimize reconstruction error. Once the AE has converged, its weights are frozen, and the encoder is utilized to project the high-fidelity trajectories into the latent space. In the second stage, these latent trajectories serve as the ground truth for training the LF-net, ensuring that the learned flow accurately captures the time evolution of the reduced system.

4.1 AutoEncoders↩︎

The maps \(\mathbf{x}_\alpha\) are not known a priori, but an approximation of them can be learned by means of a data-driven approach. Indeed, let us assume to have access to a certain number of realizations of the solution of Equation 1 for a given set of parameters \(\boldsymbol{P} \subset \mathcal{P}\). We call this set of snapshots \[\mathcal{D} =\{ \mathbf{u}_h(t; \boldsymbol{\mu}) | t \in \Theta , \, \mu \in \boldsymbol{P} \},\] where \(\Theta\) is a set of discrete time instants. By exploiting the manifold hypothesis, we assume that all the points in the dataset live in a low dimensional manifold, which we may want to learn by using a suitable neural network architecture such as an AE. Specifically, the decoder \(D: \mathbb{R}^l \rightarrow \mathbb{R}^{N_h}\) tends to approximate the map \(\mathbf{x}_\alpha\) introduced in Equation 2 , while the encoder \(E: \mathbb{R}^{N_h} \rightarrow \mathbb{R}^{l}\) is an approximation of its inverse. The AutoEncoder is defined as the composition: \[AE = D \circ E : \mathbb{R}^{N_h} \rightarrow \mathbb{R}^{N_h}.\] Training the AE on the dataset \(\mathcal{D}\) consists in finding the maps \(E\) and \(D\) which minimize the reconstruction error. Namely, \[\label{eq:32minimization32formulation} \min_{E, D} \text{MSE} (\boldsymbol{u}, \,\, AE(\boldsymbol{u})), \,\, \boldsymbol{u} \in \mathcal{D}.\tag{8}\] The combined action of decoder and encoder provide an approximation of the solution manifold: \[\widetilde{S}_h ^n\approx \{AE(u) | \, u \in S_h\}.\] Since \(\mathbf{u}_h \simeq \widetilde{\mathbf{u}} \in \widetilde{S}_h\) it is possible to assume that \(\mathbf{u}_h\) satisfies Equation 1 as well, leading to: \[\label{eq:32approx32FOM32eq46} \left\{\begin{array}{l} \dot{\widetilde{\mathbf{u}}}(\boldsymbol{z}(t ; \boldsymbol{\mu}))=\mathbf{f}_h\left(t, \widetilde{\mathbf{u}}(\boldsymbol{z}(t ; \boldsymbol{\mu})) ; \boldsymbol{\mu}\right), \quad t \in\left(t_0, T\right], \\ \widetilde{\mathbf{u}}\left(\boldsymbol{z}(t_0 ; \boldsymbol{\mu})\right)=\widetilde{\mathbf{u}}_{0}(\boldsymbol{\mu}). \end{array}\right.\tag{9}\] After some careful manipulations, the latent governing equations produced by the AEs maps read: \[\label{eq:32latent32equation32AE} \dot{\boldsymbol{z}}(\boldsymbol{z}, t ; \boldsymbol{\mu})) = (dEdD)^{-1} dE \mathbf{f}_h\left(t, \widetilde{\mathbf{u}}(\boldsymbol{z}(t ; \boldsymbol{\mu})) ; \boldsymbol{\mu}\right).\tag{10}\] Before proceeding any further, several considerations are warranted. First, a striking similarity can be observed between Equation 10 , derived from true manifold parametrizations, and Equation 3 . Although these two equations differ only by the additional term \((dE \, dD)^{-1}\) present in Equation 10 , this term vanishes if the AE belongs to the family of symmetric AEs, thereby reducing Equation 10 to Equation 3 . Despite the fact that Equation 10 is globally low-dimensional and allows for the autonomous evolution of the dynamics, the temporal integration of such a system would offer no computational advantage, as evaluating the RHS of Equation 3 requires intermediate computation of high-dimensional vectors. Any type of nonlinearity appearing in the flow \(\mathbf{f}_h\!\left(t, \widetilde{\mathbf{u}}(\boldsymbol{z}(t;\boldsymbol{\mu}))\right)\) or in the model order reduction technique itself, in fact, leads to such a caveat. To circumvent this issue, an appropriate neural network architecture can be introduced to approximate the latent flow, thereby providing the efficiency required for a high-performance ROM. It is worth stressing that the proposed approach is therefore particularly well suited for non-linear parametric PDEs. First, non-linear PDEs often suffer from a slow decay of the Kolmogorov \(n\)-width (e.g., turbulence, wave propagation), thus highlighting the need for non-linear model order reduction techniques such as AEs. Secondly, the need for a dedicated neural network to advance the latent dynamics in time is intrinsically related to the non-linearity of the FOM equations themselves.

4.2 Latent Flow Network (LF-net)↩︎

The latent flow is approximated by a Feed-Forward Neural Network, denoted as \(h_\theta\). According to the dynamics established in Equation 3 , the flow typically depends on the current latent state \(\boldsymbol{z}\), the parameter vector \(\boldsymbol{\mu}\), and, in non-autonomous cases, the explicit time \(t\). To account for these dependencies, these variables are concatenated into a single vector, which is given as input to the network, as illustrated in Figure 3.

Specifically, the latent state is obtained through the encoding mapping \(\boldsymbol{z} = E(\mathbf{u}_h(t; \boldsymbol{\mu}))\). The ground-truth latent velocity \(\dot{\boldsymbol{z}}\) is derived by computing the differential of the encoder with respect to the high-fidelity state \(\mathbf{u}_h\), evaluated at \(\mathbf{u}_h(t; \boldsymbol{\mu})\), and applying to the high-fidelity dynamics \(\mathbf{f}_h\) as follows: \[\dot{\boldsymbol{z}}(\boldsymbol{z}, t; \boldsymbol{\mu}) =dE\left(\mathbf{u}_h(t; \boldsymbol{\mu})\right) \cdot \mathbf{f}_h\left(t, \mathbf{u}_h(t ; \boldsymbol{\mu}) ; \boldsymbol{\mu}\right).\] The LF-net is trained by minimizing a Mean Squared Error (MSE) loss function, which penalizes the discrepancy between the network’s prediction and the projected high-fidelity dynamics: \[\min_{h_\theta} \text{MSE} \{h_\theta(E(\boldsymbol{u}_h), \boldsymbol{\mu}, t) , dE \mathbf{f}_h(t, \boldsymbol{u}_h; \boldsymbol{\mu}) \} \,\, \boldsymbol{u}_h \in \mathcal{D}.\]

4.3 ROM Deployment↩︎

Following the offline training of both the AE and the LF-net, the ROM can be deployed for efficient online predictions. The deployment workflow consists of two main steps: initial condition projection and latent trajectory evolution. First, the high-fidelity initial condition \(\mathbf{u}_{0, h}(\boldsymbol{\mu})\) is projected into the latent space using the learned encoder \(E\): \[\boldsymbol{z}(0 ; \boldsymbol{\mu}) = E(\mathbf{u}_{0, h}(\boldsymbol{\mu})).\] Subsequently, the latent state is evolved in time by integrating the learned latent flow \(h_\theta\). This is achieved using a suitable numerical integration scheme (e.g., explicit or implicit Euler, Heun, or higher-order Runge-Kutta methods): \[\label{eq:32ROM32solution} \widetilde{\boldsymbol{z}}(t; \boldsymbol{\mu}) = \boldsymbol{z}(0 ; \boldsymbol{\mu}) + \int_{0}^{t} h_\theta(\widetilde{\boldsymbol{z}}(\tau), \boldsymbol{\mu}, \tau) \, d\tau.\tag{11}\] In general, integrating \(h_\theta\) does not equate to solve Equation 10 . The primary reason for this is that \(h_\theta\) is designed to approximate the projection of the tangent vector within the latent tangent space learned by the encoder. However, in the case of a general autoencoder, this projection does not adequately represent the latent space governing equation, as it neglects the term \((dEdD)^{-1}\), accounting for the representation-inconsistent parameterization that has been learned. Conversely, if the encoder and decoder are chosen from families of functions such that \(E \circ D = id\) (i.e., a representation-consistent formulation), then Equation 11 is the solution of Equation 10 .

5 Numerical experiments↩︎

To evaluate the proposed Convolutional Symmetric Autoencoders within a Model Order Reduction framework, the pipeline detailed in Section 4 is applied to three canonical one-dimensional dynamical systems, whose spatio-temporal solutions are depicted in Figure 4. Each system was selected to represent a distinct challenge, thereby demonstrating the versatility of the architecture across different levels of complexity. The assessment focuses on both reconstruction accuracy, computational cost and latent space stability. We define latent stability as the capacity of the reduced model to produce solutions that are physically consistent. In deterministic, non-chaotic systems, this means that the latent trajectories remain close to the ground truth. However, in chaotic systems, this condition becomes impractical because, by nature, trajectories that initially start close will eventually diverge. In these cases, stability is characterized by the ability of the latent trajectories to remain within the invariant set, which is represented by the projection of the high-dimensional attractor into the latent space. The analysis begins with Linear Advection, a benchmark for advection-dominated problems where traditional linear reduction techniques, such as POD, typically underperform. This is followed by the parametric Viscous Burgers’ equation, where the viscosity serves as the primary parameter. Finally, the model is tested against the Kuramoto-Sivashinsky equation, representing a chaotic system.

5.1 Linear Advection↩︎

The one-dimensional Linear Advection equation describes the transport of a passive scalar field within a spatial domain of size \(L\) over a temporal horizon \(T\). The governing equations, assuming periodic boundary conditions, are defined as: \[\label{eq:32Advection32equation} \begin{cases} \frac{\partial u}{\partial t} + v \frac{\partial u}{\partial x} = 0, \quad x \in [0, L], t \in [0, T]\\ u(x, 0) = u_0(x) \\ u(0, t) = u(L, t) \end{cases}\tag{12}\] where \(v\) denotes the constant transport velocity. Despite its mathematical simplicity and linearity, this equation remains a significant challenge for classical Model Order Reduction techniques. In advection-dominated regimes, the Kolmogorov \(n\)-width exhibits a slow decay rate, meaning that linear subspace methods, such as POD, require an excessively large number of modes to achieve accurate compression. Recently it was shown that the decay of the Kolmogorov \(n\)-width is also largely affected by the regularity and smoothness of initial and boundary conditions [19].

Equation 12 is solved numerically using a second-order central finite difference scheme as spatial discretization and a \(4^{\mathrm{th}}\)-order five stages Runge-Kutta method for the time integration. For this benchmark, the physical parameters, specifically the advection velocity \(v\), the simulation time \(T\), and the domain length \(L\), were set to unity. The initial condition \(u_0(x)\) was defined as a Gaussian centered in the middle of the domain with a normalized standard deviation of \(\sigma/L = 0.05\). The spatial interval \([0, L]\) was partitioned into \(N_h + 1\) equidistant points with a grid spacing of \(h = L/N_h\). Let \(\boldsymbol{u}_h \in \mathbb{R}^{N_h}\) represent the discrete solution vector, where \(u_i \approx u(x_i, t)\) for \(x_i = ih\) and \(i = [0, 1, ..., N_h-1]\). The resulting semi-discrete system is given by: \[\label{eq:32semi95discrete95advection} \frac{d \boldsymbol{u}}{d t} + \mathbf{A} \boldsymbol{u} = 0.\tag{13}\] The dataset used for training and testing consists of snapshots derived from these numerical solutions. To ensure the generalization of the model, \(80\%\) of the data was allocated for training, while the remaining \(20\%\) was reserved for testing and validation.

5.2 Viscous Burgers Equation↩︎

The Viscous Burgers’ equation is a fundamental one-dimensional non-linear parametric PDE that serves as a benchmark for modeling advection-diffusion phenomena. In this study, we consider the equation in its conservative form under periodic boundary conditions: \[\begin{cases} \frac{\partial u}{\partial t} + \frac{\partial }{\partial x} \left( \frac{u ^ 2}{2} \right) = \nu \frac{\partial ^2 u}{\partial x ^2} \quad x \in [-L, L],\, t \in [0, T],\, \nu \in (0.001, 0.1],\\ u(x, 0) = u_0(x) \\ \text{Periodic B.C.s} \end{cases}\] where the viscosity \(\nu \in (0.001, 0.1]\) acts as the system parameter. The non-linear convective flux is primarily responsible for the development of steep gradients and shock-like discontinuities. Such localized features are notoriously difficult for global linear bases to capture. This limitation often manifests itself as spurious oscillations around the discontinuity [16], [24].

The governing equation is discretized using a FD framework with \(L=1\) and \(T=1\). Spatial derivatives are approximated using a second-order central difference scheme on a uniform grid with \(N_h + 1\) discretization points, resulting in a uniform mesh size \(h = 2L/N_h\). For temporal integration, an explicit Euler scheme is employed.

To construct the ROM, the parameter interval for \(\nu\) is sampled at \(15\) logarithmically spaced points. This sampling strategy ensures that the model is exposed to a wide range of physical regimes, from diffusion-dominated to advection-dominated flows. Out of the \(15\) generated trajectories, \(12\) are used for the training phase, while the remaining 3 are reserved for testing the model.

5.3 Kuramoto-Sivashinsky Equation↩︎

The final test case considers the Kuramoto-Sivashinsky (KS) equation, which introduces a significantly higher level of complexity due to its chaotic dynamics. Originally developed to model the instabilities of laminar flame fronts, the KS equation is one of the simplest PDEs known to exhibit high-dimensional chaos [43]. The governing equation, under periodic boundary conditions, is given by: \[\label{eq:ks95equation} \begin{cases} \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \frac{\partial^2 u}{\partial x^2} + \nu \frac{\partial^4 u}{\partial x^4} = 0, & x \in [0, L], \;t \in [0, T] \\ u(x, 0) = u_0(x) \\ u(0, t) = u(L, t) . \end{cases}\tag{14}\] The qualitative behavior of the solution, ranging from a bursting regime to fully developed chaos, is strictly dependent on the viscosity \(\nu\) and the domain length \(L\) [44]. In this study, we focus on the chaotic regime, which is readily achieved by selecting a sufficiently large domain \(L\) [45]. Indeed the number of modes that growth exponentially increases linearly with \(L\). In the chaotic regime, the nonlinear term plays a stabilizing role by facilitating an energy cascade, effectively redistributing energy from unstable low-frequency modes to damped high-frequency modes [46].

The system is solved numerically using a pseudo-spectral Fourier method for spatial discretization, coupled with an Exponential Time Differencing Fourth-order Runge-Kutta (ETDRK4) scheme for time marching, as detailed in [47]. A more comprehensive discussion on the properties of chaotic systems and their reduced-order modeling is provided in Appendix 8.

6 Results↩︎

6.1 LA equation↩︎

The study of the Linear Advection problem is motivated by two primary factors. First, for non-linear reduction techniques, it facilitates the generation of a two-dimensional latent space, which lends itself to straightforward visualization and interpretability. Second, it provides a rigorous benchmark for evaluating the relative performance of three distinct methodologies: POD, classical Convolutional Autoencoders, and their symmetric counterparts.

To construct the POD-based ROM, the numerical solutions are organized into a snapshot matrix \(\boldsymbol{X} = [\boldsymbol{u}_1, \boldsymbol{u}_2, \dots, \boldsymbol{u}_{N_T}] \in \mathbb{R}^{N_h \times N_T}\), where \(\boldsymbol{u}_i\) denotes the discrete solution to Equation 13 at time \(t_i\). By computing the Singular Value Decomposition (SVD) of \(\boldsymbol{X}\), the POD modes \(\boldsymbol{\varphi}_i\) (the left singular vectors) are extracted to span a reduced-order subspace. The governing equation is then projected onto this subspace to derive the evolution equations for the POD coefficients, defined as \(c_i = \langle \boldsymbol{\varphi}_i, \boldsymbol{u} \rangle\): \[\label{eq:pod95evolution} \frac{d \mathbf{c}}{d t} + \boldsymbol{\Phi}^T \mathbf{A} \boldsymbol{\Phi} \mathbf{c} = 0\tag{15}\] where \(\boldsymbol{\Phi} = [\boldsymbol{\varphi}_1, \boldsymbol{\varphi}_2, \dots, \boldsymbol{\varphi}_l]\) represents the basis truncated at dimension \(l\). According to the Eckart–Young theorem, the truncated SVD provides the optimal \(l\)-rank approximation of \(\boldsymbol{X}\) in the Frobenius norm. Furthermore, the total reconstruction error is bounded by the sum of the squares of the discarded singular values \(\sigma_i\) from \(i = l+1\) to \(r\), where \(r = \min\{N_h, N_T\}\) is the rank of the snapshot matrix [48].

To estimate the truncated dimension \(l\) required for an accurate POD-ROM in the unit domain problem described in Section 5.1, we analyze the cumulative energy carried by the singular values of \(\boldsymbol{X}\). As illustrated in Figure 5, the cumulative sum of the squared singular values suggests that approximately \(l \simeq 20\) modes are sufficient to capture the dominant (\(99.82\%\)) information within the dataset.

To quantify the performance of the ROMs, we define the relative \(L_2\) reconstruction error as: \[\label{eq:32epsilon32u} \epsilon_{\boldsymbol{u}}(t) = \frac{\| \boldsymbol{u}(t) - \tilde{\boldsymbol{u}}(t) \|_{L_2}}{\| \boldsymbol{u}(0) \|_{L_2}}\tag{16}\] where \(\tilde{\boldsymbol{u}}(t)\) is the solution predicted by the ROM. For the POD-ROM, the reconstruction is given by \(\tilde{\boldsymbol{u}}(t) = \boldsymbol{\Phi} \mathbf{c}(t)\). As illustrated in Figure [fig: Linear Advection POD error], a POD basis of \(l=20\) is required to maintain the relative error below \(1\%\) throughout the simulation. Reducing the number of modes to \(l < 10\) results in significantly higher errors, underscoring the limitations of linear reduction techniques even for relatively simple advection problems.

The limitations of POD motivate the use of Autoencoders for non-linear manifold learning. Given that the solution to Equation 12 is periodic and governed by a single parameter (velocity), the FOM solution traces a closed curve embedded in \(\mathbb{R}^{N_h}\). While this manifold is locally one-dimensional, its periodic nature prevents a global, injective parametrization into \(\mathbb{R}\) [49]. Consequently, a minimum latent dimension of \(l=2\) is required for the AE to capture the topological feature of the solution.

To evaluate the robustness of the ROMs, based on classical and symmetric CAEs, we generated \(10\) instances of the classic ROMs and \(10\) instances of the symmetric ROMs to assess the reproducibility of the results. With a little abuse of notation we denote with classic and symmetric ROMs, the ROMs that employ classic and symmetric CAEs, respectively. For the symmetric variant, we utilize a parameter-free decoder that serves as the pseudoinverse of the encoder (i.e. \(\boldsymbol{f} = 0\) in Equation 6 ), motivated by the simplicity of the problem at hand.

To facilitate a direct comparison between models trained independently, which may produce topologically equivalent but geometrically distinct latent representations, we map all latent trajectories onto a unit circle. For a high-dimensional trajectory \(\boldsymbol{u}_h(t, \boldsymbol{\mu})\), first we project it onto the latent space (\(\boldsymbol{z}(t, \boldsymbol{\mu}) = E(\boldsymbol{u}_h(t, \boldsymbol{\mu}))\)), and then define the trajectory center of mass \(\bar{\boldsymbol{z}}(\boldsymbol{\mu})\) as: \[\bar{\boldsymbol{z}}(\boldsymbol{\mu}) = \frac{\int_0^T \boldsymbol{z}(t, \boldsymbol{\mu}) \| \dot{\boldsymbol{z}}(t, \boldsymbol{\mu}) \| dt}{\int_0^T \| \dot{\boldsymbol{z}}(t, \boldsymbol{\mu}) \| dt}.\] Next, introduce polar coordinates relative to the center: \[\theta(\boldsymbol{x})=\text{atan2}{\left(\frac{x_2}{x_1}\right)}, \qquad \boldsymbol{x}\in \mathbb{R}^2\] and define the radial profile of the trajectory as the function \[\rho(\alpha)=\| \boldsymbol{z}(t, \boldsymbol{\mu}) - \bar{\boldsymbol{z}}(\boldsymbol{\mu}) \|, \quad \text{such that} \quad \theta(\boldsymbol{z}(t, \boldsymbol{\mu}) - \bar{\boldsymbol{z}}(\boldsymbol{\mu}))=\alpha.\] In other words, \(\rho(\alpha)\) gives the distance from the origin of the true latent trajectory along direction \(\alpha\). It is now possible to map each point \(\boldsymbol{z}(t, \boldsymbol{\mu})\) into a unique point of the unit circle \(\hat{\boldsymbol{z}}\), by means of a normalization function \(\mathcal{N}\): \[\hat{\boldsymbol{z}}(t, \boldsymbol{\mu}) = \frac{\boldsymbol{z}(t, \boldsymbol{\mu}) - \bar{\boldsymbol{z}}(\boldsymbol{\mu})}{\rho(\theta(\boldsymbol{z}(t, \boldsymbol{\mu}) - \bar{\boldsymbol{z}}(\boldsymbol{\mu})))} = \mathcal{N}(\boldsymbol{z}(t, \boldsymbol{\mu}); E, \boldsymbol{u}_h).\] This normalization function \(\mathcal{N}\) is solely determined by the encoder \(E\) and the reference solution \(\boldsymbol{u}_h\), as the shifting constant \(\bar{\boldsymbol{z}}(\boldsymbol{\mu})\) and the scaling function \(\rho\) are defined by the true latent trajectory \(\boldsymbol{z}(t, \boldsymbol{\mu})\). More generally, for any point \(\boldsymbol{x}\in \mathbb{R}^2\) we define \[\boldsymbol{\hat{x}} = \mathcal{N}(\boldsymbol{x}; E, \boldsymbol{u}_h) = \frac{\boldsymbol{x} - \bar{\boldsymbol{z}}(\boldsymbol{\mu})}{\rho(\theta(\boldsymbol{x} - \bar{\boldsymbol{z}}(\boldsymbol{\mu})))}.\] Notably, the normalization function maps the predicted ROM trajectories onto the unit circle only when they accurately follow the ground truth. Any deviation from the true latent trajectory is thus reflected as a departure from the unit circle of the normalized ROMs prediction.

After the training phase, the ROMs were evolved for \(10\) periods (\(t \in [0, 10T]\)) to test the stability of the latent flow.

The resulting normalized trajectories are visualized in Figure 7. The \(10\) trajectories generated by the symmetric ROMs tend to stay closer to the unit circle, where they should belong. This is not the case for the classic one because, soon after the second period, some of the trajectories depart from the unit circle and start wandering around. A clear evolution of the error defined by Equation 16 and the error in latent space, defined as \[\epsilon _ {\boldsymbol{z}}(t) = \| \hat{\boldsymbol{z}}(t) - \mathcal{N}(E(\boldsymbol{u}_h(t))) \|\] is illustrated in Figure 8. The symmetric ROMs maintain an error profile that is significantly lower than that of the classical models. We attribute this superior performance to three primary factors:

1. Having a pair of encoder-decoder models that are consistent provides a latent space that is a true coordinate system. Any latent point which is decoded from the latent space to the physical space and then encoded it again, it results in the same input point.

2. The constraint \(E \circ D = id\) acts as a powerful regularizer during training. This leads to a smoother latent manifold that is inherently more amenable to interpolation and temporal integration by the latent flow network.

3. The symmetric constraint promotes a lower Lipschitz constant for the latent flow [17]. This reduces the sensitivity of the dynamics to small numerical perturbations, ensuring long-term stability in periodic regimes.

6.2 Viscous Burger Equation↩︎

In the second test case, we explore the application of symmetric CAEs in the context of non-linear parametric PDE. Given that the only two parameters involved are time and viscosity, we can choose a two-dimensional latent space size (\(l=2\)). Multiple models were trained, but we focus on the two best-performing models.

Figure 9 illustrates the reduced trajectories in the latent space. The lack of smoothness in the latent space produced by the classical CAE is particularly evident when examining viscosity values from the test trajectories. The solid lines in Figure 9 represent the projection of the FOM solution \(\boldsymbol{z} = E(\boldsymbol{u}_h)\) through the encoder. Among these, three lines stand out as notably thicker: these correspond to the test set. When employing a standard convolutional AE for compression, the test trajectories exhibit significant wiggling, indicating irregularities in regions distant from the training points. In contrast, when projecting the curve using the Symmetric CAE, these irregularities are substantially reduced, further demonstrating the regularizing effect inherently introduced by this architecture.

This issue also impacts the chain of prediction capabilities within the ROM. Specifically, the LFNet relies on the knowledge of the tangent vector (the derivative taken along the trajectory) at a particular point in the latent space. When the derivative of a given noisy curve is computed, it tends to amplify the effects of the noise, resulting in even more irregular derivatives (tangent vectors). These perturbations cause the predicted trajectory by classic ROM (computed via Equation 11 ) and depicted by the dashed lines in the left plot of Figure 9, to deviate from the reference trajectory due to the accumulation of large errors. In this framework, the overlap of the solid and dashed lines for the symmetric ROM can be explained by better accuracy in approximating the smoother time derivative, leading to slower error propagation.

The improved regularity of the latent trajectories translates into superior physical reconstruction. Figures 10 and 11 compare the FOM solutions with the ROM predictions and their corresponding punctual errors for viscosity values in both the training and testing sets. For the chosen training viscosity value (Figure 10), the classical architecture notably fails to capture the precise position of the shock front, the most fundamental feature of the flow in this regime. Furthermore, as the viscosity increases and the importance of the non linear term decreases accordingly, the diffusive nature of the PDE results in a more smeared shock profile. In these higher-viscosity test cases (Figure 11), the classical ROM manifests significantly larger errors for \(t \ge 0.2\), precisely when the diffusion process begins to dominate. The symmetric ROM, however, maintains high fidelity across both regimes.

In order to make a general comparison, we introduce the absolute error: \[\label{eq:32absolute32rec32error} \epsilon_{\boldsymbol{u}} ^{\text{abs}} (t; \nu) = \|\boldsymbol{u}_h(t; \nu) - D(\boldsymbol{\hat{z}}(t; \nu)) \|_{L_2}.\tag{17}\] The performance of the classical and symmetric ROMs is evaluated via the mean absolute error, defined as: \[\overline{\epsilon_{\boldsymbol{u}}^{\text{abs}}}(\nu) = \frac{1}{T} \int_0^T \epsilon_{\boldsymbol{u}}^{\text{abs}}(t; \nu) , dt.\] As illustrated in Figure 12, a distinct trend characterizes the symmetric architecture: the mean error consistently decreases as the viscosity parameter increases across both training and testing ranges. By contrast, the classical ROM lacks this monotonic behavior, displaying a less predictable error distribution across the parametric space. Detailed analysis of the temporal evolution reveals the main source of these errors. In low-viscosity regimes, the classical ROM fails to resolve precise shock locations, resulting in peak error values at the beginning of the simulation. This behavior is clearly visible in the bottom row of Figure 10. Conversely, in viscosity-dominated scenarios, the classical model exhibits its largest inaccuracies toward the end of the temporal domain, as shown in Figure 11. While the symmetric ROM follows similar qualitative patterns, it maintains significantly lower error magnitudes on average. This superior performance, summarized in Figure 12, underscores the greater reliability preserved by the symmetric CAE-ROM throughout the investigated parametric range.

6.3 Kuramoto-Sivashinsky Equation↩︎

The application of ROMs to chaotic systems, such as the Kuramoto-Sivashinsky equation, represents a significant challenge due to the inherent unpredictability of the dynamics. In chaotic regimes, trajectories exhibit a sensitive dependence on initial conditions, where small perturbations grow exponentially. Consequently, instantaneous field comparisons, standard in non-chaotic cases, lose their physical significance after a short integration interval. This unpredictability adds complication in understanding weather the reduced model is able to actually capture the true essence of the system being modeled (invariant properties of the chaotic attractor). We should then be equipped with tools that reveal the essence of the dynamics and extract quantities that are truly representative of the chaotic system at hand. These tools can be recast into two macro groups: stability analysis (Covariant Lyapunov vectors and Lyapunov exponents) and statistical analysis (evolution of the mean variables, Probability Density Functions (PDFs), energy spectrum, etc.). The latter are extensively employed to analyze turbulent fluid fields, which can be considered one of the prime examples of chaos [50].

In this study, we fix the KS parameters to \(L = 20\pi\) and \(\nu = 1\) to ensure a fully developed chaotic regime [23], characterized by a leading Lyapunov Exponent (LE) \(\lambda = 0.085\), which is associated to a characteristic Lyapunov Time (LT) \(\tau_ \lambda = 1/\lambda\). The space domain is discretized using \(N_h = 512\) nodes, and the equation is solved using a time step \(\Delta t = 0.25\), for a total time of \(T= 2.5 \times 10^4\), corresponding to \(2125\) LTs. We chose a latent space size \(l=24\), based on the study made by Özalp and Magri [23]. To evaluate the performance of the classic and symmetric ROMs, we move beyond trajectory-based errors and assess the models’ ability to capture the flow’s statistical and spectral characteristics.

The classical CAEs exhibited a more pronounced lack of robustness when subjected to the chaotic regime of the KS equation (compared, for instance, to the Linear Advection test case). Of the three instances trained, two suffered from instability, with their solutions effectively “exploding” after only a few LTs. In these instances, the norm of the latent state solution \(\boldsymbol{z}\)—calculated via Equation 11 —increases rapidly, escaping the latent attractor (either the global attractor or the inertial manifold) showing an unphysical behavior and drifting into regions of the latent space that were never encountered during the training phase. Decoding these regions far from the inertial manifold, leads to an unphysical high-dimensional solution that depends strongly on the weights learned and it is unpredictable in principle. Even the single instance that achieved partial stability remained highly sensitive, manifesting instabilities for several specific initial conditions.

In contrast, the models based on the symmetric AE architecture proved to be remarkably robust. All seven generated instances maintained stable behavior throughout the integration period, representing a significant improvement over the classical approach. Therefore, the consistent stability of the symmetric models confirms that the structural constraints effectively regularize the latent flow, keeping the dynamics confined to the attractor.

The prediction horizon is defined as the temporal interval beyond which the ROM trajectory noticeably diverges from the ground-truth high-fidelity solution. In the existing literature, surrogate models for the KS equation—particularly those utilizing autoregressive approaches on FOM snapshots (e.g., RNNs, LSTMs, and Reservoir Computing)—report prediction horizons ranging from \(1\) to \(5\) Lyapunov Times [51]–[53]. When Reservoir Computing is applied specifically to latent dynamics, this horizon is often reported around \(1.25\) LT [52].

Starting from a selected initial condition, we integrated the reduced-order models—utilizing both the classical CAE (Figure 13) and the symmetric architecture (Figure 14)—for a total duration of \(12\) LT. Both models display the characteristic features of the KS solutions, with strides that merge and split over time. The prediction horizon stands at about \(1\) and \(2\) LT for classic and symmetric ROMs respectively, which is not too far from \(1.25\) found by Özalp and Magri [23].

The decision to employ a Latent Flow Network as the temporal propagator, rather than traditional Recurrent Neural Networks or Reservoir Computing, provides distinct advantages in terms of both training efficiency and model robustness. Primarily, the LFNet avoids the need for Backpropagation Through Time, which is a requirement for training LSTM or GRU-like architectures. While increasing the unroll length in RNNs can enhance predictive performance, it simultaneously increases memory overhead and computational cost, often rendering long-range dependency training prohibitive [54]. Furthermore, while Reservoir Computing methods circumvent traditional problems of BPTT [51], their performance is notoriously sensitive to the calibration of hyperparameters, such as the spectral radius and Tikhonov regularization parameter [55]. By shifting the objective to the supervised learning of the latent state’s time derivative, the LFNet is a “plug-and-play” surrogate. This approach allows the network to internalize the underlying dynamics during the training phase, effectively eliminating the exhaustive parameter tuning associated with RC while maintaining equal or even superior predictive accuracy.

Another useful analysis concerns the Lyapunov exponents, , which measure the sensitivity of a solution to initial conditions by quantifying the average rate of growth or decay of a small perturbation along a given direction. For a rigorous theoretical treatment and numerical algorithms to compute LEs, we refer the reader to [56]–[58]. An high level overview about the methods to compute the LEs from the FOM equation is provided in Appendix 9. To compute the LEs of the ROM, we leverage the automatic differentiation capabilities of PyTorch [59]. Given that the latent dynamics of the ROM are governed by \[\label{eq:32ROM32eq32LE} \dot{\mathbf{z}}(t; \boldsymbol{\mu}) = h_\theta(\mathbf{z}, \boldsymbol{\mu}, t)\tag{18}\] the perturbation evolves according to: \[\label{eq:32ROM32variational32equation} \dot{\mathbf{v}}(t) = d_ {\boldsymbol{z}} h_\theta(\mathbf{z}, \boldsymbol{\mu}, t) \mathbf{v}\tag{19}\] where \(d_ {\boldsymbol{z}} h_\theta(\mathbf{z}, \boldsymbol{\mu}, t)\) can be easily computed using backpropagation. Equation 18 and 19 are solved together using a RK4 scheme, with time stepping for equation 18 which is half of the time step used to solve equation 19 . Furthermore, Appendix 10 provides a mathematical proof demonstrating that, under specific conditions, the LEs of the ROM and FOM are equivalent. As established in Appendix 8.1, the LS serves as the “fingerprint” of a chaotic system, encapsulating critical information regarding its dimensionality, entropy, and long-term predictability. Figure 15 illustrates the comparison of the LEs computed by the FOM with those derived from the ROMs. The models accurately capture the positive LEs. However, as noted in [51], [60], ROMs are unable to represent the multiplicity of the zero LEs. This behavior is suggested in [60] but has not been elaborated upon. As discussed in Appendix 11, the zero LEs are linked to the number of continuous symmetries present in the dynamical system. In particular, the Kuramoto-Sivashinsky equation exhibits three continuous symmetries: time-translation invariance, space-translation invariance, and Galilean invariance. Neural networks do not maintain symmetries such as space translation or Galilean invariance [61]. Consequently, the only symmetry retained in the latent dynamical system is the time-translation invariance associated with the Lyapunov vector \(\dot{\mathbf{z}}(\mathbf{z},t ; \boldsymbol{\mu})=d\boldsymbol{x}_\alpha ^{-1}(\mathbf{f}_h\left(t, \mathbf{x}_\alpha (\mathbf{z}(t ; \boldsymbol{\mu})) ; \boldsymbol{\mu}\right))\). Given that the ROMs compute \(l\) LEs while omitting the two zero LEs, it is evident that the ROM produced two spurious LEs. By manually eliminating them, you can observe the overlapping between the true and predicted spectra in Figure 15.

We computed the Kaplan-Yorke (KY) dimension taking into account the two missing LEs. The KY dimension and the value of the first LE is reported in table 1, along with the relative errors.

To assess the long-term ergodic properties of the reduced-order models, we computed the Probability Density Function of the kinetic energy over an extensive integration period of \(2083\) LT. The kinetic energy is defined as: \[k(t) = \frac{1}{2} \int _ 0 ^ L u^2(x, t) dx.\] The statistical characterization of the classical ROM was notably hampered by its inherent numerical instability; as previously discussed, certain initial conditions triggered divergent trajectories. These unstable instances resulted in unphysical state reconstructions, which produced an erroneous PDF with a pronounced bias toward high-energy states. As illustrated in Figure 16, while both the classical and symmetric ROMs capture the general distribution of the high-fidelity kinetic energy, the classical model exhibits spurious high-energy configurations that are absent in the FOM data.

To conclude the spectral characterization of the Kuramoto-Sivashinsky system, we compute the averaged kinetic energy spectrum utilizing Welch’s method. This analysis allows for a precise evaluation of how the reduced-order models distribute energy across different spatial scales. The symmetric ROM demonstrates a superior capability in capturing the energy content across the entire frequency spectrum, whereas the classical CAE-based ROM exhibits significant discrepancies. These inaccuracies are particularly evident in the low-frequency regime, which corresponds to the large-scale coherent structures.

7 Conclusions↩︎

In this work, we begin by mathematically deriving the equation of the reduced model, leveraging some basic knowledge from differential geometry. Next, we adopt tools from the machine learning domain that, without any external assistance—i.e., in an unsupervised manner—are capable of learning an approximation of the manifold’s parametrization by examining a subset of points belonging to that manifold. Unfortunately, the standard versions of these tools do not fully satisfy all the properties that parametrization should possess, particularly the concept of representation consistency. Recent literature has introduced augmented neural network types, specifically autoencoders with fully connected layers, that address the representation consistency property. This paper extends this property to hybrid convolutional autoencoders. In this context, “hybrid” refers to an encoder composed of convolutional layers, while the decoder is constructed with fully connected neural networks designed to be parametrized by a limited number of weights—either matching the encoder or, potentially, none. This approach allows the encoding process to retain the advantages of CNN layers, such as a number of trainable parameters that is independent of the input vector size, shared weights, and a limited receptive field, while permitting the decoder to leverage the full expressiveness of feedforward neural networks for the reconstruction task. Although there are convolutional neural networks capable of processing not only vectors but also 2D arrays, our implementation focuses solely on the 1D version. Addressing the extension to 2D cases presents challenges that will be explored in future work. Consequently, this study is limited to 1D PDEs.

Here, the symmetric AEs, are integrated inside a ROM framework, to see if benefits in terms of performance can be achieved once you get closer to "ideal" maps. The ROM comprises both an AE and a network, which we called Latent Flow Network, able to predict the latent flow at any point and for any given value of the parameters. Indeed, three different cases were investigated, namely the Linear Advection (LA), the Viscous Burger (VB) equation and the Kuramoto-Sivashinsky (KS) equation. Each of them, carries its own complexity, like advection dominating problem (LA), parameters dependency (VB) or chaos (KS). Throughout these test cases, the symmetric CAE performed equal or better compared to the classical version of the CAEs. In particular, it can be seen that introducing this representation consistency constraint, and in particular the formulation based on projector, acts as a powerful regularizer, producing a latent space which is smoother and easier to be learnt, by any latent dynamic model. Specifically, for LA and VB, ROM based on symmetric AE appeared to be more accurate. Moreover, employing symmetric AEs, results in more stable ROMs (like models that are easier to train or models that doesn’t explode). This was particularly clear, when a reduced model was created for the chaotic KS equation. Indeed, it was not so easy to deliver a working model with the classical CAE, while, with the introduction of symmetric AEs we were able to provide viable reduced models for each trained model. For the chaotic KS problem, particular attention is given to the spectra, specifically, to the Lyapunov spectrum. Here we proved first that the FOM and the ROM, possess the same non zeros LEs, and then we explained why, with NN it is difficult to obtain the exact multiplicity of the zero LEs.

8 Chaotic systems↩︎

8.1 Lyapunov Exponents (LEs) and related quantities↩︎

Let \(\boldsymbol{u}(t, \boldsymbol{\mu})\) be the solution of the autonomous dynamical system \(\dot{\boldsymbol{u}}(t, \boldsymbol{\mu}) = \boldsymbol{f}(\boldsymbol{u}(t, \boldsymbol{\mu}))\), \(\boldsymbol{u}(t, \boldsymbol{\mu}) = \Phi (\boldsymbol{u}_0;t, \boldsymbol{\mu})\) with \(\Phi: \mathcal{M} \times \mathbb{R} \times \Omega \rightarrow \mathcal{M}\) the solution map, \(\mathcal{M}\) the solution manifold and \(\Omega \subseteq \mathbb{R} ^ p\), the parameters space.

Suppose that we introduce a small perturbation (of magnitude \(\varepsilon_0\)), along a direction \(\boldsymbol{t} \in T_{\boldsymbol{u}} \mathcal{M}\), that we call \(\boldsymbol{v}_0 = \varepsilon_0 \boldsymbol{t}\). We can study how the perturbation evolves along the trajectory \(\boldsymbol{u}(t; \boldsymbol{\mu})\) by linearizing the dynamical’s system governing equation around the reference solution \(\boldsymbol{u}(t; \boldsymbol{\mu})\). By expanding in Taylor series the flow \(\boldsymbol{f}(\boldsymbol{u}(t, \boldsymbol{\mu}))\) the perturbed dynamical system can be written as \[\cancel{\dot{\boldsymbol{u}}(t, \boldsymbol{\mu})} + \dot{\boldsymbol{v}}(t) = \cancel{\boldsymbol{f}(\boldsymbol{u}(t, \boldsymbol{\mu}))} + d_ {\boldsymbol{u}} \boldsymbol{f}(\boldsymbol{u}(t, \boldsymbol{\mu})) \boldsymbol{v} + o(\varepsilon \boldsymbol{1}),\] and so \[\label{eq:32variational32equation} \dot{\boldsymbol{v}}(t) \approx d_ {\boldsymbol{u}} \boldsymbol{f}(\boldsymbol{u}(t, \boldsymbol{\mu})).\tag{20}\] Equation 20 is called variational equation and it is a linear differential equation in the variable \(\boldsymbol{v}\) with time-dependent coefficient matrix \(d_ {\boldsymbol{u}} \boldsymbol{f}(\boldsymbol{u}(t, \boldsymbol{\mu}))\). The solution to equation 20 can be written using the solution map \(\Phi\) as: \[\boldsymbol{v}(t) = d_ {\boldsymbol{u}} \Phi(\boldsymbol{u}_0; t, \boldsymbol{\mu}) \boldsymbol{v}_0.\]

By tracking the evolution of one perturbation, you will get with probability one, the largest LE [56]. This is the most important LE, since its sign is indicative of the nature of the system under investigation. The Multiplicative Ergodic Theorem (MET) guarantees the existence of the LEs and of a base \(\{\boldsymbol{e}_1, \boldsymbol{e}_2, ..., \boldsymbol{e}_n\}\) of \(T_{\boldsymbol{u}_0}\mathcal{M}\), such that the LEs associated to each element of the base [57]: \[\lambda(d_{\boldsymbol{u}}\Phi,\boldsymbol{e}_1) \ge \lambda(d_{\boldsymbol{u}}\Phi,\boldsymbol{e}_2) \ge ... \ge \lambda(d_{\boldsymbol{u}}\Phi,\boldsymbol{e}_n)\] The collection of all the LEs form the Lyapunov Spectrum. The Lyapnunov Spectrum provides useful information about the inertial manifold and about the global attractor embedded within it.

While definition 2 and 3 looks similar they are quite different mathematical objects. Indeed, while the global attractor could be a fractal set, the inertial manifold is required to be smooth, or Lipschitz. Another important difference concerns the rate at which trajectories that do not initially lie on these manifolds approach them. For the global attractor, the distance between a trajectory and the attractor is known to vanish asymptotically as time tends to infinity. In contrast, trajectories approach the inertial manifold at an exponential rate in time.. Last but not least, the global attractor is a subset of the inertial manifold.

For some dissipative systems (KS equation, Ginsburg-Landau equation, Cahn-Hilliard equations), it can be proved the existence of a inertial manifold [9]. The Lyapunov exponents are tightly related with these manifolds. For instance, the Kaplan-York (KY) dimension gives an upper bound to the fractal dimension of the global attractor. \[\label{eq:32KY32dimension} D_{KY} = k + \frac{\sum _ {i=1} ^k \lambda_i}{|\lambda _{k+1}|}\tag{21}\]

Where \(k\) is the index such that \(\sum _ {i=1} ^k \lambda_i \ge 0\) and \(\sum _ {i=1} ^{k + 1} \lambda_i < 0\).

A chaoticity measure of the dynamical system is given by the Kolmogorov-Sinai entropy, which can be interpreted as the average amount of new information generated by the system per unit of time, which means how unpredictable is the evolution of the system itself [62]. The Pessin’s theorem indeed states that the Kolmogorov-Sinai entropy can be expressed as the sum of the positive LEs [56].

9 FOM LEs of KS equation↩︎

The variational equation for the KS system represented by Equation 14 can be written as: \[\label{eq:32variational32equation32KS} \frac{\partial v}{\partial t} + \frac{\partial u v}{\partial x} + \frac{\partial^2 v}{\partial x ^2} + \nu \frac{\partial^4 u}{\partial v ^4}= 0.\tag{22}\] Both equation 14 and equation 22 are solved together, using a spectral scheme coupled with ETDRK4 time stepping scheme. Specifically, in order to evaluate the non linear term in between the time discretization nodes, as ETDRK4 requires, we integrated equation 14 by using half of the time steps used to solve 22 . Then, the LEs for the FOM are computed using the algorithm described by [56]. A validation of this solver was performed by comparing the computed spectra with the reference one reported in [23], for the same value of the domain length \(L\). The validation is not reported here.

10 Equivalence of LE↩︎

To show the equivalence of the LEs found by the FOM and ROM system’s governing equations we consider the same framework introduced in section 8.1. Suppose that \[T_{\boldsymbol{u}} \mathcal{M} = \text{span}(\{\boldsymbol{\xi}_1(\boldsymbol{u}),\, ..., \, \boldsymbol{\xi}_l(\boldsymbol{u})\})\] where \(\boldsymbol{\xi}_i(\boldsymbol{u})\) is the \(i\)-th covariant Lyapunov vector associated with the \(i\)-th Lyapunov exponent. Let us denote \(E:\mathcal{M} \rightarrow \Sigma \subseteq \mathbb{R}^l\) a (full rank) manifold parameterization, and \(\varphi = E \circ \Phi \circ E^{-1}\), the solution operator in the latent space. It is possible to apply Definition 1 to the solution map defined by \(\varphi\) in the latent space: \[\label{eq:32reduced32lambda} \lambda ^r (d_{\boldsymbol{u}}\varphi,\boldsymbol{v})=\limsup_{t\to\infty}\frac{1}{t}\ln\|d_{\boldsymbol{z}}\varphi \boldsymbol{v}\|.\tag{23}\]

Let \(\boldsymbol{w} \in T_{\boldsymbol{u}} \mathcal{M}\), and hence \(\boldsymbol{v} = d_{\boldsymbol{u}} E \boldsymbol{w}\). We want to prove that \(\lambda ^r (d_{\boldsymbol{u}}\varphi,\boldsymbol{v}) = \lambda (d_{\boldsymbol{u}}\Phi,\boldsymbol{w})\). We can then expand the terms inside the norm in Equation (23 ) as: \[\label{eq:32norm} \| dE d\Phi d E^{-1} dE \boldsymbol{w} \|,\tag{24}\] since \(\boldsymbol{w} \in T_{\boldsymbol{u}} \mathcal{M}\), \(\boldsymbol{w} = d E^{-1} dE \boldsymbol{w}\). Let \(\sigma_{min}(\boldsymbol{u})\) and \(\sigma_{max}(\boldsymbol{u})\) be maximum and minimum modules of the singular values of the linear map \(dE\) at the point \(\boldsymbol{u} \in \mathcal{M}\). Let us introduce the definitions \[\overline{\sigma}_{min} = \inf_{\boldsymbol{u}\in \mathcal{M}} \left( \sigma_{min}(\boldsymbol{u}) \right) \quad \mathrm{and} \quad \overline{\sigma}_{max} = \sup_{\boldsymbol{u}\in \mathcal{M}} \left( \sigma_{max}(\boldsymbol{u}) \right).\] Being \(E\) full rank, it is guaranteed that \(\overline{\sigma}_{min}>0\) and a linear map between finite spaces is always bounded. The argument of the logarithm in equation (23 ) can then be bounded as \[\overline{\sigma }_{min} \|d\Phi \boldsymbol{w} \| \le \|dE d\Phi \boldsymbol{w} \| \leq \overline{\sigma }_ {max} \|d\Phi \boldsymbol{w} \|.\] By applying the logarithm to the above expression and exploiting the monotonicity of the logarithm, the same inequalities hold: \[\ln (\overline{\sigma} _ {min} \|d\Phi \boldsymbol{w} )\| \le \ln( \|dEd\Phi \boldsymbol{w} \|) \leq \ln(\overline{\sigma} _ {max} \|d\Phi \boldsymbol{w} \|).\] Since the logarithm of a product is the sum of the logarithm of each factor, by the squeezing theorem, we can prove that \(\lambda ^r (d_{\boldsymbol{u}}\varphi,dE\boldsymbol{w}) = \lambda (d_{\boldsymbol{u}}\Phi,\boldsymbol{w})\): \[\limsup_{t\to\infty}\frac{1}{t} \ln (\overline{\sigma} _ {min} \|d\Phi \boldsymbol{w} \| ) \le \limsup_{t\to\infty}\frac{1}{t} \ln( \|d\varphi \boldsymbol{v} \|) \\ \leq \limsup_{t\to\infty}\frac{1}{t} \ln(\overline{\sigma} _ {max} \|d\Phi \boldsymbol{w} \|)\] \[\lambda (d_{\boldsymbol{u}}\Phi,\boldsymbol{w}) \le\lambda ^r (d_{\boldsymbol{u}}\varphi,dE\boldsymbol{w}) \leq \lambda (d_{\boldsymbol{u}}\Phi,\boldsymbol{w})\]

11 Symmetries and Lyapunov Exponents in Infinite-Dimensional Systems↩︎

Let us consider an abstract evolution equation on a Banach (or Hilbert) space \(\mathcal{H}\): \[\frac{\partial u}{\partial t} = F(u), \qquad u(t) \in \mathcal{H},\] where \(F : \mathcal{H} \to \mathcal{H}\) is a sufficiently smooth (e.g.,\(C^1\)) operator. We denote by \(\Phi(t,u_0)\) the associated solution map, so that \(u(t)=\Phi(t,u_0)\) solves the equation with initial condition \(u(0)=u_0\).

Let \(\{g_\alpha\}_{\alpha \in \mathbb{R}}\) be a one-parameter group of transformations \[g_\alpha : \mathcal{H} \to \mathcal{H},\] such that: \[\begin{align} g_0 &= \mathrm{Id}, \\ g_{\alpha+\beta} &= g_\alpha \circ g_\beta, \end{align}\] and the map \(\alpha \mapsto g_\alpha u\) is \(C^1\) for all \(u \in \mathcal{H}\).

We say that the evolution equation is equivariant under this group if \[F(g_\alpha u) = g_\alpha F(u), \qquad \forall u \in \mathcal{H}, \;\alpha \in \mathbb{R}.\]

Equivalently, the flow satisfies \[\Phi(t, g_\alpha u_0) = g_\alpha \Phi(t, u_0), \qquad \forall t \in \mathbb{R}.\]

Since the symmetry is continuous, we define the infinitesimal generator \[v(u) := \left.\frac{d}{d\alpha} g_\alpha u \right|_{\alpha=0}.\]

Under the above assumptions, \(v : \mathcal{H} \to \mathcal{H}\) is a well-defined vector field tangent to the group orbit \(\{ g_\alpha u : \alpha \in \mathbb{R} \}\). Differentiating the equivariance relation \(\Phi(t, g_\alpha u_0) = g_\alpha \Phi(t, u_0)\) with respect to \(\alpha\) at \(\alpha=0\), we obtain \[d_u \Phi(t,u_0)\, v(u_0) = v(\Phi(t,u_0)).\]

Thus, the vector \[v(t) := v(u(t)) = v(\Phi(t,u_0))\] is a solution of the linearized equation \[\dot{v} = dF(u(t))\, v.\]

In particular, the subspace spanned by the symmetry generator is invariant under the linearized dynamics.

The Lyapunov exponent associated with this direction is defined as \[\lambda_v(u_0) = \limsup_{t \to \infty} \frac{1}{t} \log \| D_u \Phi(t,u_0)\, v(u_0) \|.\]

Using the identity above, this becomes \[\lambda_v(u_0) = \limsup_{t \to \infty} \frac{1}{t} \log \| v(\Phi(t,u_0)) \|.\]

Therefore, the Lyapunov exponent associated with the symmetry direction is determined entirely by the growth of the norm of the generator along the trajectory.

In particular, if there exist constants \(C > 0\) and \(p \ge 0\) such that \(\| v(\Phi(t,u_0)) \| \le C (1+t)^p\), then, as shown in[63], \(\lambda_v(u_0) = 0\).

12 NNs hyper parameters↩︎

Hyperparameter tuning for the autoencoder was performed exclusively on the viscous Burgers’ equation, with the resulting architecture subsequently applied to both the linear advection and Kuramoto–Sivashinsky (KS) problems. Given the relative simplicity of the linear advection case, we set the function \(f\) in Equation 6 to zero yielding a weight-free decoder. For the more complex viscous Burgers’ and KS cases, we defined \(f\) as a ConvTranspose1d layer from PyTorch

References↩︎