Kinetic equations are crucial for modeling non-equilibrium phenomena, but their computational complexity is a challenge. This paper presents a data-driven approach using reduced order models (ROM) to efficiently model non-equilibrium flows in kinetic equations by comparing two ROM approaches: Proper Orthogonal Decomposition (POD) and autoencoder neural networks (AE). While AE initially demonstrate higher accuracy, POD's precision improves as more modes are considered. Notably, our work recognizes that the classical POD-MOR approach, although capable of accurately representing the non-linear solution manifold of the kinetic equation, may not provide a parsimonious model of the data due to the inherently non-linear nature of the data manifold. We demonstrate how AEs are used in finding the intrinsic dimension of a system and to allow correlating the intrinsic quantities with macroscopic quantities that have a physical interpretation.

Keywords: Model order reduction, data-driven methods, kinetic equations, neural autoencoder networks, proper orthogonal decomposition, Sod shock tube, Boltzmann-BGK

§ INTRODUCTION

Kinetic equations are widely used in science and engineering [20, 25, 26, 33]. They allow the modeling of deviations from an equilibrium model which is given by an underlying macroscopic equation like the Euler equations, providing detailed insight into fundamental physical processes [35]. However, kinetic equations are often characterized by a large dimensional phase space, making them computationally expensive to solve and sometimes even unfeasible for realistic applications [35].

Investing in solving kinetic equations is only beneficial if large deviations from equilibrium are present [35]. Striking a balance between a fast but inaccurate equilibrium solver and a slow but accurate non-equilibrium solver remains an open challenge. Our work aims to address this challenge by providing a proof-of-concept for a data-driven solution to efficient modeling of flows in different non-equilibrium regimes.

In the field of non-equilibrium gas flows, several standard methods to discretize the high dimensional phase space exist. Particle-based Monte Carlo methods are only tractable in the free flight regime and out of scope in the transition regime of moderate non-equilibrium unless special techniques are used [5]. The straightforward Discrete Velocity Method (DVM) uses a pointwise discretization of the velocity space, potentially leading to a large number of equations [30]. Specially tailored moment models are based on the expansion of the particle distribution function and lead to a set of extended fluid dynamical equations [35]. However, it is by no means clear a-priori how many equations are sufficient and which variables are optimal [19, 34].

To tackle the computational complexity of kinetic equations, recently, reduced order models (ROM) have been introduced, enabling reductions in computational complexity by orders of magnitudes [1, 7, 8, 9]. Two different approaches have been followed in the literature. The classical offline-online decomposition as used by [1] involves a two-stage procedure. In the offline stage, the full order model (FOM) is assessed to create a database, which is then utilized to generate a data-dependent basis through proper orthogonal decomposition (POD). This basis allows an efficient description of the FOM on a low-dimensional linear subspace during the online phase. On the other hand, the online adaptive basis method called dynamic low-rank approximation [13] constructs the low dimensional linear basis during the online phase itself, eliminating the need to evaluate the expensive FOM. It has been successfully applied to kinetic equations in the works by Einkemmer et al. [7, 8, 9]. However, the additional complexity of updating the basis during the evolution makes it less online efficient than the classical offline-online approach shown in [15] for a shallow water moment model.

In this work, we adopt the same offline strategies as in [1]. Specifically, we sample data for a classical test case called Sod shock tube using a discrete velocity model as our FOM and compare the compression of the linear reduced subspace created by POD with a non-linear description provided by neural autoencoder networks. Neural networks, based on the universal approximation theorem [29], allow for the approximation of a wide range of function classes and appear promising in identifying the intrinsic dimension of a system. However, the non-linear relation between macroscopic model equations and the discrete velocity model hinders the determination of these dimensions using linear reduction methods like the POD.

This paper aims to utilize these data-driven model reduction techniques to reduce the number of describing variables and equations and determine how many and which variables are useful in specific test cases. For the non-vanishing Knudsen number, we expect to need more non-equilibrium variables with corresponding balance laws, while in the limit of vanishing Knudsen number, we expect to recover the Euler equations, given by conservation laws for mass, momentum, and energy. To the knowledge of the authors, this is the first paper aiming to bridge the gap between equilibrium and non-equilibrium flows using neural networks in this way.

The long-term objective of this line of work is to enable dynamically adapting the model by varying the number of variables during the online phase, paving the way for more efficient and accurate model adaptive simulations of kinetic equations.

The organization of the paper is as follows: In Section 2, we introduce the 1D model problem and the reference data used for model reduction. Section 3 describes the two model reduction techniques used in this study: Proper Orthogonal Decomposition (POD) and Autoencoder Networks. The results are presented in Section 4, and the paper concludes with a summary in Section 5.

§ THE BOLTZMANN-BGK MODEL AND DATA

This paper considers a proof-of-concept of using reduced models for the solution approximation of the 1D Boltzmann-BGK equation [2] for monoatomic, ideal gases

\begin{equation}
\label{Eq:BGK}
\partial_t f + c \partial_x f = \frac{1}{\tau} (f_M - f),
\end{equation}

which is a potentially high-dimensional equation for the unknown probability density function $f(t,x,c)$, where $t \in \mathbb{R}^+$ is the time, $x \in \mathbb{R}$ is the spatial variable, and $c \in \mathbb{R}$ the microscopic particle velocity. For simplicity we consider the one-dimensional case in this paper, but the results can be extended to the multi-dimensional case.

Computing solutions and generating data of the Boltzmann-BGK model is essential for industrial and scientific applications, but often so computationally prohibitive that a large number of test cases is not feasible. To reduce time and cost during the data generating process, experiments or numerical simulations can be replaced by reduced-order models (ROMs).

For standard continuum flows the widely-used Euler equations can be applied, but more rarefied regimes require different extended fluid dynamical models. Rarefaction levels are distinguished with the help of the Knudsen number $\Kn$ defined by the ratio of the mean free path length of the particles $\lambda$ over a reference length $l$ as

\begin{equation}
\Kn = \frac{\lambda}{l}.
\end{equation}

The right-hand side of the BGK collision operator (<ref>) models the relaxation with relaxation time $\tau \in \mathbb{R}^+$ towards the equilibrium Maxwellian distribution $f_M(t,x,c)$ given by

\begin{equation}
f_M(t,x,c) = \frac{\rho(t,x)}{(2\pi R T(t,x))^{\frac{3}{2}}}\exp \left(-\frac{(c - u(t,x))^2}{2 R T(t,x)}\right),
\end{equation}

where $\rho(t,x)$, $v(t,x)$ and $T(t,x)$ are density, bulk velocity, and temperature of the flow, respectively. $R$ is the universal gas constant.
In this work, we consider the relaxation time $\tau$ a parameter and set it equal to the Knudsen number, $\tau = \Kn$, however, the relaxation time can also be changed, e.g., to depend on the gas density and temperature in addition.

For practical computations, we consider macroscopic moments of the distribution function, which are given by multiplying the distribution function with the co-called collision invariants $(1,c,\frac{1}{2} c^2)$ and integrating in velocity space

\begin{eqnarray}
\rho(t,x) &=& \int f(t,x,c) \,\mathrm{d}c,\label{Eq:Moments1}\\
\quad\rho(t,x) u(t,x) &=& \int c f(t,x,c) \,\mathrm{d}c,\label{Eq:Moments2}\\
\quad E(t,x) &=& \int \frac{1}{2}c^2 f(t,x,c)  \,\mathrm{d}c,\label{Eq:Moments3}
\end{eqnarray}

where $E$ denotes the total energy. The temperature $T(t,x)$ and the pressure $p(t,x)$ can be obtained by

\begin{equation}
T(t,x) = \frac{2E(t,x)}{3\rho(t,x)} - \frac{u(t,x)^2}{3} \quad\textrm{and}\quad p(t,x) = \rho(t,x) T(t,x).
\end{equation}

<Ref> illustrates the relation between the macroscopic moments and the distribution function $f(t,x,c)$ at a certain position in time and space. The density $\rho(t,x)$ is the integral of the distribution function, which is centered around the macroscopic velocity $u(t,x)$, and the mean deviation is related to the temperature $T(t,x)$.

Illustration of the macroscopic moments corresponding to an example distribution function.

The Boltzmann-BGK equation (<ref>) is in equilibrium when $f = f_M$. Multiplying the equilibrium solution with the collision invariants and integrating in velocity space, one finds the Euler equations of classical gas dynamics

\begin{align}
\partial_t&\rho + \partial_x(\rho u) = 0, \label{Eq:Conservation1} \\
\partial_t&(\rho u) + \partial_x(\rho u^2 + p) = 0, \label{Eq:Conservation2}\\
\partial_t&E + \partial_x(u(E+p)) = 0, \label{Eq:Conservation3}
\end{align}

which are conservation laws for mass, momentum, and energy, respectively.

For distribution functions further away from equilibrium, for example due to a larger relaxation time $\tau$ and a significantly large Knudsen number $\Kn$, the Euler equations do not give accurate results. In this case, additional equations can be used, which are derived by the so-called method of moments [14, 35]. This effectively leads to an extended set of equations, called moment model. It is possible to preserve important properties like hyperbolicity with moment models [11, 17]. The additional equations (for example for the heat flux and higher-order moments) add complexity, but allow for more accurate solutions [18, 34]. However, it is often unclear a-priori, how many equations are needed for an efficiently accurate and computationally feasible solution. In this work, we aim to give a proof-of-concept for a data-based identification of the necessary number of variables, called the intrinsic physical dimension.

§.§ Sod shock tube test case and reference data

Sod's shock tube is a well-established test case in the field of rarefied gases [19]. It uses discontinuous initial conditions based on equilibrium values

\begin{equation}
\label{eqn:shock_tube_problem}
\begin{cases} (\rho_L, u_L, p_L) = (1,0,1) &\mbox{if } x < 0.5, \\ (\rho_R, u_R, p_R) = (0.125,0,0.1) & \mbox{if } x > 0.5, \end{cases}
\end{equation}

corresponding to a jump in density and pressure at $x=0.5$ due to a diaphragm at that position, which is removed at time $t=0$.

The problem setup at $t=0$ is shown in <ref>, which is split into two regions left and right of the diaphragm.

Problem setup for the 1D Sod shock tube. A diaphragm at the center is initially separating the domain in two regions, where initial conditions for density $\rho$, macroscopic velocity $u$, and pressure $p$ are indicated.

For the generation of reference data we employ a discrete velocity method (DVM) [27], which uses a pointwise microscopic velocity space discretization

\begin{equation}
\partial_t f_{k}(t,x) = -(c_k) \partial_x f_{k}(t,x) + \frac{1}{\tau}\left({M_f}_{k}(t,x) - f_{k}(t,x)\right),
\label{Eq:Discrete BGK}
\end{equation}

where a uniform grid in velocity space is considered with $c_k = k\Delta c$ to discretize the distribution function $f_k(t,x) = f(t,x,c_k)$, for some $k \in \mathbb{Z}$. After a subsequent discretization in space, the DVM <ref> leads to a coupled ODE system in time than can be solved with standard methods.

For the numerical reference data, we use $N_x = 200$ spatial points $x_k \in [0, 1]$, $N_c = 40$ discrete velocities $c_j \in [-10, 10]$ and $N_t = 25$ time steps $t^n \in [0, 0.12]$ summarized in <ref>. It is possible to choose another range for the discrete velocity points, but in typical applications the range of the bulk velocity is not known, such that one has to include a safety margin. We therefore chose the domain $[-10, 10]$. The goal of the model order reduction is now to reduce the complexity of the computation using lower dimensional models. For that matter, it is not relevant what the actual error of the numerical reference data is or if is fully converged. It is fair to say that a full reference solution might easily take into account more spatial points, time steps, and discrete velocities, which makes it even more necessary to reduce the complexity.

For the model order reduction later, we consider two different Knudsen numbers for Sod's shock-tube test case: $\textrm{Kn} = 0.00001$ for a small Knudsen number in the hydrodynamic regime and $\textrm{Kn} = 0.01$ for a relatively large Knudsen number in the rarefied regime.

To understand the behavior of the reference solutions for non-vanishing Knudsen numbers, we first describe the solution for vanishing Knudsen number in equilibrium, i.e., $\Kn = 0$, which can be obtained using the method of characteristics and the Rankine Hugoniot jump conditions connecting the states before and after the shocks [24].

Starting from the initial condition in <ref>, the solution evolves for $t> 0$ and five regions are formed that are depicted in <ref> [32]. A rarefaction wave is moving to the left between $x_1$ and $x_2$. The contact discontinuity is located at $x_3$, where the macroscopic velocity $u$ and the pressure $p$ are continuous in contrast to the density $\rho$ and the energy $E$. $x_4$ is the position of the shock wave.

In non-equilibrium, i.e., for solutions evolving with Knudsen numbers $\Kn > 0$, the solution does not have discontinuities due to the finite relaxation time $\tau$. <Ref> shows the reference solutions $f(t,x,c)$ at $t_0=0$, $t_1=0.06$ and $t_3=0.12$ for the two levels of rarefaction considered in this paper: $\Kn = 0.00001$ and $\Kn =0.01$. Increasing the Knudsen number leads to a smoother transition from region 1 to region 5 with a less pronounced shock front.

Sod's shock tube ($t=0$) containing two regions with corresponding initial conditions.

Sod's shock tube ($t>0$) containing the rarefaction wave, the contact discontinuity and the shock wave.

Discrete reference solutions $f(t,x,c)$ for $\Kn=0.00001$ (top row) and $\Kn=0.01$ (bottom row), for times $t_0=0$ (left), $t_1=0.06$ (center) and $t_2=0.12$ (right).

Macroscopic quantities $\rho$, $\rho u$, $E$ ($t=0.12$) for $\Kn=0.00001$ ($\hy$) and $\Kn=0.01$ ($\rare$).

Sod shock tube and reference data. Initial conditions (a); equilibrium solution (b); Reference solutions in rarefied and hydrodynamic regime (c); Macroscopic quantities at $t=0.12s$ (d).

Problem setup for the Boltzmann-BGK model in Sod's shock tube.

Variable       Number of nodes \(i\)      Domain extension    Step size (uniform)

\(x\) 		   	200        [0, 1]   	    0.005

\(c\)             40      		    [-10,10]   	    \(\approx\) 0.51282051

\(t\)   	   	25            	[0,0.12]   	      0.005

§ METHODS

In this section we present two common methods used for reducing the dimensionality of high dimensional data: (1) the proper orthogonal decomposition (POD) and (2) neural autoencoder networks (AE). The methods will be used to parameterise the high dimensional data stemming from the DVM simulation, using a linear mapping in case of the POD and a non-linear mapping in case of AE. Although the classical POD-MOR approach shows that linear mappings are sufficient to describe the non-linear solution manifold of the BGK equation to a good accuracy [1], it is in general not sufficient to determine a parsimonious model of the full model data, since the data manifold can be non-linear. Here, neural autoencoder networks can be used as they are capable to find the intrinsic dimension of a system.

§.§ Proper orthogonal decomposition

The proper orthogonal decomposition [31] (POD) approximates the data with help of dyadic pairs:

\begin{equation}
\label{eq:POD-approx}
f(x,t,c_i) \approx \sum_{k=1}^r \hat{f}_k(x,t)\psi_k(c_i) \qquad \text{for } r\ll N_c\,.
\end{equation}

The pairs $\{(\hat{f}_k(x,t),\psi_k(c_i))\}_{k=1,\dots,r}$ are the structures in the data that contain the most energy and they are chosen to minimize the gap between the data and the reconstruction <ref>.
In the following, $\psi_k(c_i)$ are termed POD-modes and $\hat{f}_k(x,t)$ the corresponding reduced variables.

For notation, we define $f^{(i)}(x,t) = f(x,t,c_i)$ and the vectors $f(x,t) = (f^{(1)}(x,t), \dots, f^{(\Nc)}(x,t))$ and $\psi_k = (\psi_k(c_1),\dots, \psi_k(c_{\Nc}))$.
The proper orthogonal decomposition computes the solution of the minimization problem:

\begin{equation}
\min_{\psi_k}\norm{f(x,t)-\sum_{k=1}^r\langle f(x,t),\psi_k\rangle\psi_k}_2^2\, \quad\text{such that}\quad\langle \psi_k,\psi_l\rangle =\delta_{kl}.
\end{equation}

Technically one can solve this optimization problem using a singular value decomposition (SVD) of the so-called snapshot matrix [22]:

\begin{equation}
\label{eq:snapshotmatrix}
\mathbf{F}=\begin{bmatrix}
f^{(1)}(x_1, t_1) &   \cdots & f^{(\Nc)}(x_1, t_1) \\
\vdots &  \ddots & \vdots \\
f^{(1)}(x_{\Nx}, t_1) &   \cdots & f^{(\Nc)}(x_{\Nx}, t_1) \\
f^{(1)}(x_1, t_2) &   \cdots & f^{(\Nc)}(x_1, t_2) \\
\vdots & \ddots &  \vdots \\
f^{(1)}(x_{\Nx}, t_{\Nt}) &   \cdots & f^{(\Nc)}(x_m, t_{\Nt}) \\
\end{bmatrix}\in\mathbb{R}^{(\Nx\Nt)\times\Nc}\,.
\end{equation}

The snapshot matrix collects the time and space discrete distribution function in its columns. Each column holds the time-spatial values for a different discrete velocity. Performing an SVD factorizes $\matr{F}$ as

\begin{equation}
\matr{F}=\matr{\Phi\Sigma \Psi^T},
\end{equation}

with diagonal matrix $\matr{\Sigma}=\diag(\sigma_1,\dots,\sigma_m)$, $m=\min(\Nx\Nt,\Nc)$, containing the singular values $\sigma_1\ge \sigma_2\ge\dots\ge \sigma_m\ge0$ and $\matr{\Phi}\in\mathbb{R}^{\Nx\Nt\times m}, \matr{\Psi}\in\mathbb{R}^{\Nc\times m}$ are orthogonal matrices containing the left and right singular vectors, respectively.

The first $r$ columns of the truncated $\matr{\Psi}_r=[\psi_1,\dots,\psi_r]\in\mathbb{R}^{\Nc\times r}$ contain the POD modes in <ref>. Together with $\matr{\Sigma}_r=\diag(\sigma_1,\dots,\sigma_r)$ and the $r$ leading left singular vectors $\matr{\Phi}_r\in\mathbb{R}^{(\Nx\Nt)\times r}$ they yield the rank $r$-term approximation $\matr{F}_r$ of the snapshot matrix given by

\begin{equation}
\matr{F}_r:= \matr{\Phi}_r\matr{\Sigma}_r \matr{\Psi}_r^T.
\end{equation}

According to the Eckart-Young-Mirsky theorem [6, 28] $\matr{F}_r$ is the best rank $r$ approximation and the resulting error in the Frobenius norm is rigorously computed from the trailing singular values

\begin{equation}
\norm{\matr{F}-\matr{F}_r}_\mathrm{F}^2=\sum_{k=r+1}^m\sigma_k^2\,.
\end{equation}

A common choice for $r$ is to truncate after a certain energy percentage is reached in the reduced system compared to the full system:

\begin{equation}\label{e:ePOD}
\end{equation}

In this paper, POD is used to compare with the autoencoder from the next section.

§.§ Autoencoders

Neural networks, particularly autoencoder networks, have become widely used tools for dimension reduction [36]. A comprehensive introduction to autoencoder networks can be found in [12]. Here we only summarize briefly the common idea of autoencoder networks and give the specific details of our implementation.

An autoencoder aims to reproduce the input data while compressing it through an information bottleneck. It consists of two main components:

* The encoder, denoted as $\enc$, maps the input data ${f}$ to points ${\hat{f}}$ in a lower-dimensional latent space: $\enc \colon \mathbb{R}^{M}\to\mathbb{R}^r, {f}\mapsto {\hat{f}} = \enc({f})$, $r\ll M$.

* The decoder, denoted as $\dec$, reconstructs the input space from the latent representation: $\dec \colon\mathbb{R}^{r}\to\mathbb{R}^{M}, {\hat{f}} \mapsto \dec({\hat{f}}) = {\tilde{f}}$, $r\ll M$.

Note that the dimension of the latent space is denoted by $r$, to match the rank of the POD approximation.

The autoencoder is defined as the composition of both parts: ${\tilde{f}} = \dec \left( \enc \left( {f} \right) \right)$.
For our purpose, we identify the discrete velocity space as the input dimension $M=\Nc$. Thus, the autoencoder maps each time-spatial value of the distribution function $f(x,t)\in \mathbb{R}^{\Nc}$ onto a smaller latent space $\hat f(x,t)\in \mathbb{R}^{r}$, which parameterizes the necessary physical information of the system.

The goal of the optimization procedure is to determine $\dec$ and $\enc$ such that the reconstruction error over a set of training/testing data contained in $\matr{F}$ is minimized. The reconstruction error is defined as:

\begin{equation}
\label{eq:loss}
\mathcal{L} = \frac{1}{\Nc}\norm{f(x,t)- \dec\circ\enc\circ f(x,t)}_2^2   \,.
\end{equation}

The reconstruction error is the sum of the two-norm of the discrete velocities vector of the difference between the input data $f$ and the reconstructed data $\tilde{f}$ that has been squeezed through the informational bottleneck. The assumption is, that if the original data can be represented well while the information went through a smaller latent space, there exists a physical law in the latent space that describes the system sufficiently. The intrinsic latent dimension $r=\intrDim$ which is sufficient to describe the data is then called the intrinsic physical dimension similar to the intrinsic dimension defined in [23]. Such a reduced model is then termed parsimonious because it explains the data with a minimum number of variables.

In the training procedure, the functions $\enc$ and $\dec$ are determined by trainable parameters of the network, referred to as weights and biases. The networks are constructed using a composition of layers $\enc=L_1\circ L_2 \circ \dots \circ L_N$. Typically, each layer $L_n\colon \mathbb{R}^i \to \mathbb{R}^o$ in the network consists of an affine linear mapping $\vect{x}\mapsto h_n(\matr{W}_n {x}+{b}_n)$, where $\matr{W}_n\in\mathbb{R}^{o,i}$ represent the weights, $\vect{b}_n \in\mathbb{R}^o$ denote the biases, and $h_n$ are predefined non-linear functions. The configuration of the input and output dimensions $i$ and $o$ for each layer, the choice of activation function, and the number of layers collectively determine the architecture of the network. The choice of these so-called hyper-parameters is often difficult and a matter of trial and error.

In our studies we have exploited fully connected neural autoencoder networks (FCNN) and convolutional autoencoder networks (CNN). However, in this manuscript we restrict ourselves to the results of the fully connected network, since it gave structurally the best results.
We have studied a variety of different activation functions, hidden layers, batch sizes and depths of the network. The best results concerning the validation error and acceptable training time where obtained by the network
defined in <ref>. A comprehensive study of the parameter optimization is attached to the manuscript in <ref>.

Parameter     hydrodynamic     rarefied

Layer Sizes     [$\Nc$, 30, $r$, 30, $\Nc$]    [$\Nc$, 40, $r$, 40, $\Nc$]

Activation function     ELU     ReLU

Loss function     MSE <ref>     MSE <ref>

Optimizer     Adam     Adam

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

Epochs     3000     3000

Batch size     4     4

Hyper-parameters of Fully Connected Autoencoder Network (FCNN) for the hydrodynamic and rarefied regime.

Before the training we initialize the weights of the network using the standard initialization implemented in . Thus the weights are randomly uniform distributed between $m^{-1/2}$ and $m^{1/2}$ with $m$ being the number of input nodes in the layer.
Our network is trained by splitting the data consisting of $\Nx\times\Nt$ samples in a testing and training set with a 80/20 split over 3000 epochs using a batch size of 4. In each epoch the network is updated using the Adam optimizer with a learning rate of $10^{-5}$. More information about hyperparameters and training of the network can be found in the appendix <ref>.

§ RESULTS

We reconstruct the full order model (FOM) solution with the help of POD and an autoencoder, the FCNN, for which the selection of hyperparameters and the training are described in the previous section. Note that we apply both model reduction techniques to reconstruct both the rarefied reference data and the hydrodynamic reference data. The goal is to later determine the intrinsic dimension of the data for both cases.
We therefore compare the two dimension reduction techniques by means of different measures. The intrinsic variables obtained from POD and the FCNN will be referred to as \(\idhy\) and \(\idrare\), where the former describes the intrinsic variables when reducing the hydrodynamic data and the latter when reducing the rarefied data.

For the purpose of comparing the results we define the $L_2$-error

\begin{equation}
\mathcal{E}_\text{rel}=\frac{\norm{\matr{F}-\tilde{\matr{F}}}_2}{\norm{\matr{F}}_2}\,,
\end{equation}

where $\matr{F}$ the reference data is given in <ref> and the reconstructed data $\tilde{\matr{F}}$ is either $\matr{F}_r$ in case of the POD or the FCNN predictions with $r$ latent variables for every $(x,t,c)$ in the data set.

§.§ Singular value decay of reference data

As a first step, we perform a POD with the hydrodynamic data and with the rarefied data. The obtained singular values \(\sigma\), as well as the cumulative energy (cusum-e) defined in <ref>, are shown in <ref>. As expected, more modes are necessary in the rarefied regime compared to the hydrodynamic regime. With a total of \(\intrPOD=4\) intrinsic variables, a cumulative energy of over \(99\%\) can be achieved for the hydrodynamic regime.
The cumulative energy of the singular values of the rarefied regime only reaches above \(99\%\) with \(\intrPOD=6\) singular values.

For the POD we define the intrinsic dimension $\intrPOD$ as the smallest truncation rank $r$ of the reduced system, at which the cumulative energy $E_\text{cum}$ defined in <ref> reaches $99\%$. Although this choice is arbitrary it is a common practice in classical MOR to truncate <ref>, whenever $99\%$ of the cumulative energy is reached.

In <ref>, we further see that the rate at which the singular values drop is approximately exponential in both regimes, which has been also observed by [1]. Consequently, a rapid decay of the Kolmogorov N-width is indicated.
Note that the singular value decay is similar for both domains but not exactly the same, thus leading to an expected increase in the number of intrinsic variables in the rarefied regime necessary to achieve similar \(L_2\)-errors.

Hydrodynamic regime.

Rarefied regime.

Singular value decay \(\sigma\) and cumulative energy increase for the number of singular variables \(k\) in the hydrodynamic regime (a) and in the rarefied regime (b). A black cross marker corresponds to over 99% cumulative energy.

It is important to note that the parameter $\intrPOD$ is not expected to precisely match the actual intrinsic dimension $\intrDim$ of the solution manifold. The intrinsic dimension represents the minimum number of variables required to accurately describe the system's exact solution manifold. This discrepancy arises because the solution manifold is fundamentally nonlinear, making it challenging to adequately capture with a parsimonious linear approximation.

For the FCNN, the intrinsic dimension is defined as the smallest number of intrinsic variables that minimizes the error. In well-trained models, the FCNN's intrinsic dimension should ideally align with $\intrDim$.

From a fluid mechanics perspective, the hydrodynamic regime theoretically requires only $\intrDim=3$ intrinsic variables. This is because near-equilibrium flows in this regime can be effectively characterized by three macroscopic quantities: density $\rho$, macroscopic velocity $u$, and total energy $E$, as outlined in <ref>-<ref>, see also [1, 19].

Conversely, the rarefied regime demands a larger intrinsic dimension, denoted as $\intrDim$. This is due to the need for more than only the equilibrium Maxwellian distribution function to describe the microscopic velocities adequately. Therefore, we initially set $\intrDim=3$ intrinsic variables ($\idhy$) for the FCNN in the hydrodynamic case and choose $\intrDim=5$ intrinsic variables ($\idrare$) for the rarefied regime. This choice aligns with extended fluid dynamic models as described in [19, 35].

We note that each FCNN with different latent space dimension $r$ needs to be trained separately. This is different from the POD, where the decomposition is only performed once. Thereafter the approximation quality is given by the truncation rank $r$.

§.§ Variation of the number of intrinsic variables

The variation of the number of intrinsic variables \(r\) in <ref> sheds light on the performance of the autoencoder with different bottleneck layer sizes. In the case of the POD $r$ is the truncation rank of the decomposition <ref> and the latent space dimension in case of the FCNN. To this end, $r$ is varied for both the POD and the FCNN over $r \in \{1,2,3,4,8,16,32\}$ for the hydrodynamic case and over $r \in \{1,2,4,5,8,16,32\}$ for the rarefied case. We note that the loss of information when applying POD goes exponentially to zero with increasing $r$, which is not surprising when consulting the Eckard-Young Theorem  [6].
Note that the FCNN is retrained for each different $r$. By changing $r$, i.e. widening the bottleneck layer, a gain or loss of capacity occurs that can be connected to stability during training.

The $L_2$-Error over the variation of the latent space dimension/truncation rank \(r\) using FCNN/POD for the hydrodynamic regime (left) and the rarefied regime (right).

Both for the hydrodynamic and the rarefied regime, POD initially yields a larger error than the FCNN for small number of intrinsic variables $r$. Not surprisingly, the POD accuracy increases with the number of singular values taken into account until the error reaches machine precision. The FCNN error decreases as well and then reaches a plateau, with a typical remaining error due to the network architecture and training. For the previously identified values $\intrDim=3$ in the hydrodynamic case and $\intrDim=5$ in the rarefied case, the FCNN results in a more accurate approximation than the POD.

We note that when testing the FCNN against POD and fixing \(r\) the FCNN is limited by the estimation error of the training and performs under its abilities. However, POD uses five to six times more parameters than the FCNN while the deterministic character enables POD to achieve any possible accuracy, which was not observed with the neural network.

In the following, we consider the intrinsic variables with constant values $\intrDim=3$ in the hydrodynamic case and $\intrDim=5$ in the rarefied case.

This leads to the number of trainable parameters of the POD and the FCNN shown in <ref>. We can see that the FCNN achieves a relatively small error with a small number of parameters in comparison with the POD for this choice of the number of intrinsic variables.

Amount of parameters used to reconstruct \(f\), the number of intrinsic variables \(p\) and the corresponding $L_2$-Error for POD and FCNN, both for hydrodynamic (\(\hy\)) and rarefied (\(\rare\)) regimes.

Algorithm     2cParameters \(\frepar\)     2cInt. variables \(p\)    2c$L_2$-error
\(\hy\)   \(\rare\)   \(\hy\)   \(\rare\)   \(\hy\)   \(\rare\)

POD         15129     25225     3     5     0.0205     0.0087

FCNN 	    2683     3725     3     5     0.0008     0.0009

Next, a qualitative analysis with the actual reconstructions is presented. From computations of the \(L_2\)-error over time \(t\), which are not shown due to conciseness, it became clear that the time step that contributes most to the error is the last time step in case of POD, while the FCNN distributes the error more evenly over all time steps.

§.§ Reconstruction quality

The reconstructed solutions compared to the full order model (FOM) at the final time step $t=0.12s$ are  given in <ref> . Because of the small overall errors indicated in <ref> both the POD and the FCNN reproduce the FOM solution without any visible differences at first sight.

Comparison of the FOM solutions \(f\) (left column) with reconstructions \(\tilde{f}\) obtained from POD (middle column) and FCNN (right column) at end time \(t=0.12s\) for \(x\in [0.375,0.75]\) in the hydrodynamic regime (top row) and in the rarefied regime (bottom row).

Comparison of FOM macroscopic quantities \(\rho\), \(\rho u\) and \(E\) with POD and the FCNN for the hydrodynamic regime (top row) and the rarefied regime (bottom row) at time \(t=0.12s\).

As shown in <ref>, the more subtle information loss from the model reduction can unfold in actual differences in the macroscopic quantities \(\rho\), \(\rho u\) and \(E\).
Overall, <ref> shows that the errors are larger for the hydrodynamic regime (top row), most notably for the momentum and energy of the POD model close to the contact discontinuity. However, the position of the shock is well approximated. In contrast, the FCNN model yields a very good agreement in the hydrodynamic case.
For the rarefied regime, both models approximate the FOM solution very well. The lack of sharp shock structures in the full model and the increased intrinsic dimension $\intrDim=5$ combined seem to notably influence the accuracy.

§.§ Conservation properties

The physical consistency of the reduced \(\tilde{f}\), in terms of conservation of mass, momentum, and energy, is a critical criterion for its validity. Hence, conservation properties are analyzed in the following. We note that conservation of mass, momentum, and energy is not directly built-in by means of a specifically tailored loss function. Even though we can expect to recover some conservation properties as they are implicitly built into the numerical reference data. The investigation of different loss functions to improve upon this is left for future work.

We investigate the conservation properties by means of the derivative of the cell-averaged conserved quantities mass, momentum, and total energy, defined exemplary as

\begin{equation}
\frac{\mathrm{d}}{\mathrm{d}t}\int \rho(x,t)\, \mathrm{d}x\Delta t  =\overline{\dot{\rho}}\mathrm{,}
\end{equation}

for the mass. Note that a derivative $\overline{\dot{\rho}}= 0$ denotes conservation of mass, for example. We expect the conservation of mass and total energy, while the momentum increases with a constant rate, due to the boundary conditions of the test case, featuring a larger pressure on the left-hand side of the domain.

<Ref> shows the evolution of the derivatives of mass, momentum, and total energy as a function of time for the hydrodynamic regime (top row) and the rarefied regime (bottom row).

Comparison of the conservative properties of reconstructions obtained from POD and the FCNNs with the conservative properties of the FOM solution using the mass, momentum, and energy derivative. A value of $0$ indicates conservation. FCNN approximately conserves mass and energy, while the momentum increases with the correct rate of the test case.

Indeed, <ref> indicates that conservation of mass and energy are achieved with reasonable accuracy for the FCNN, while the error is larger for the POD reconstruction for both regimes. Also the increase in momentum of the full order model (FOM) is accurately described by the FCNN with a larger error for the POD method for both regimes. Overall, the errors are slightly smaller for the rarefied case compared to the hydrodynamic case, which might be due to the higher capacity of the neural network and more modes for the POD using five intrinsic variables in the rarefied regime.

§.§ Physical interpretability

An important question in the context of model reduction with neural networks is the interpretability of the results, because of the usual black-box nature of neural networks [10].
Especially when benchmarking neural networks for model order reduction with POD, evaluating the interpretability of the intrinsic variables is important, since POD by construction achieves a so-called physically interpretable decomposition of the input data [3] as outlined previously.

Following the assumption that the hydrodynamic case can be fully described in terms of three macroscopic quantities and that the rarefied case is reasonably describable in a similar way with an extended set of five variables, we test the intrinsic variables \(\idhy\) and \(\idrare\) for similarities and investigate if they match any macroscopic quantity. Two related macroscopic quantities, namely the temperature \(T\) and macroscopic velocity \(u\), are added to the three macroscopic variables $\rho$, $\rho u$, and $E$. In <ref> and <ref>, these are plotted first over the whole domain of \(x\) and \(t\) and for the end time \(t=0.12s\) for both regimes. Similarly, we plot both the FCNN's and the POD's first 3 intrinsic variables \(\idhy\) of the hydrodynamic case and 5 intrinsic variables \(\idrare\) of the rarefied test case

\begin{align*}
\idhy &= [h_0(x,t),h_1(x,t),h_2(x,t)],\\
\idrare &= [r_0(x,t),r_1(x,t),r_2(x,t),r_3(x,t),r_4(x,t)],
\end{align*}

depicted in <ref> and <ref> respectively.

Macroscopic quantities of the hydrodynamic case obtained by the FOM. Density \(\rho\), momentum \(\rho u\), total energy \(E\), temperature \(T\), and velocity \(u\) over time \(t\) and space \(x\) in the top row and at \(t=0.12\) in the bottom row.

Intrinsic variables \(h_0(x,t)\), \(h_1(x,t)\) and \(h_2(x,t)\) of hydrodynamic case obtained by the FCNN. Top row depicts the whole \((x,t)\) domain, bottom row is for \(t=0.12\).

First three intrinsic variables \(h_0(x,t)\), \(h_1(x,t)\) and \(h_2(x,t)\) of hydrodynamic case obtained by the POD. Top row depicts the whole \((x,t)\) domain, bottom row is for \(t=0.12\).

Macroscopic quantities of the rarefied case obtained by the FOM. Density \(\rho\), momentum \(\rho u\), total energy \(E\), temperature \(T\), and velocity \(u\) over time \(t\) and space \(x\) in the top row and at \(t=0.12\) in bottom row.

Intrinsic variables \(r_0(x,t)\), \(r_1(x,t)\), \(r_2(x,t)\), \(r_3(x,t)\) and \(r_4(x,t)\) of rarefied case obtained by the FCNN. Top row depicts the whole \((x,t)\) domain, bottom row is for \(t=0.12\).

First five intrinsic variables \(r_0(x,t)\), \(r_1(x,t)\), \(r_2(x,t)\), \(r_3(x,t)\) and \(r_4(x,t)\) of rarefied case obtained by the POD. Top row depicts the whole \((x,t)\) domain, bottom row is for \(t=0.12\).

Strikingly, most intrinsic variables of the FCNN in <ref> and <ref> and the POD in <ref> and <ref> appear to be a combination of the five intrinsic variables shown in <ref> and <ref>, respectively.
In particular consider FCNN in the hydrodynamic case by comparing <ref> and <ref>. The rarefaction wave, shock wave, and contact discontinuity, which can be identified in \(h_0\) reflect a combination of those found in the density \(\rho\) and the total energy \(E\). Furthermore, \(h_1\) seems to model the negative momentum \(\rho u\), with different boundary values. The temperature \(T\) appears linked to \(h_2\), where the same fluctuation appears.
Similar results hold true for the POD variables in <ref>, where especially the first two variables closely resemble the density $\rho$ and the momentum $\rho u$. Interestingly, this does not relate to very good conservation properties of the POD in comparison with FCNN as shown in the previous section and <ref>.
Considering the FCNN in the rarefied case, we compare <ref> and <ref>. Here \(r_3\) clearly reflects the shape of the density \(\rho\). Moreover, the peak of the velocity \(u\) can be observed in \(r_0\). For the other intrinsic variables of \(r_1\), \(r_2\) and \(r_4\) a clear discernability of macroscopic quantities is difficult to observe and might require linear or nonlinear combinations. Additionally, it is possible that those intrinsic variables resemble non-equilibrium variables not present among the macroscopic variables. Considering the POD results in <ref>, we again see a very good agreement of the first intrinsic variables with density and momentum.

For more physical insight into the relation between macroscopic variables and intrinsic variables, the Pearson correlation between all variable combinations is computed in [12] for the FCNN. Note that we expect similar results for the POD based on the previous results, but do not present them here for conciseness. The Pearson correlation coefficient $r_{X,Y} = r_{Y,X}$ is a measure of linear correlation between two sets of data, here represented by a macroscopic variable $X \in \{\rho, \rho u, E, T, u\}$ and an intrinsic variable $Y \in \{r_0, r_1, r_2,\}$ for the hydrodynamic case and $Y \in \{r_0, r_1, r_2, r_3, r_4\}$ for the rarefied case. It is commonly defined as

\begin{equation}
r_{X,Y} = \frac{\sum \left(x_i - \Bar{x} \right)\left(y_i - \Bar{y} \right)}{\sum \sqrt{\left(x_i - \Bar{x} \right)}\sqrt{\sum \left(y_i - \Bar{y} \right)}}.
\end{equation}

Note that $r_{X,Y}\in [-1,1]$, with $r_{X,Y} = 0$ meaning that there is no correlation between both data sets, $r_{X,Y} = 1$ indicating a perfect correlation, and $r_{X,Y} = -1$ indicating a perfect anti-correlation.

The Pearson correlation coefficients for the hydrodynamic test case are presented in <ref>. As predicted by the previous analysis, there appears to be an almost perfect correlation of the first intrinsic variable $r_0$ with the density $\rho$. This means that the FCNN succeeds at identifying the density $\rho$ precisely as an internal variable. Note that $r_0$ also correlates almost perfectly with the energy $E$, as the energy depends linearly on $\rho$. Additionally, there is a relatively strong linear correlation between $r_2$ and both the momentum $\rho u$ as well as the temperature $T$. The intrinsic variable $r_2$ on the other hand seems to be correlated with all variables.

The Pearson correlation coefficients for the rarefied test case are presented in <ref>. In agreement with the previous results, there is no clear correlation of most of the intrinsic variables. An exception is $r_0$, which is anti-correlated with $u$ and $r_3$, which correlates almost perfectly with the energy $E$. Note that only linear correlations are tested here. For further analysis on the disentanglement of the intrinsic variables, it might be more suitable to consider nonlinear correlations beyond the Pearson correlation coefficients. One option that could be explored in future work would be to train the functional relation between the intrinsic variables and macroscopic variables, for example by means of symbolic regression [4].

Pearson correlation between macroscopic quantities and intrinsic variables for the hydrodynamic case.

Pearson correlation between macroscopic quantities and intrinsic variables for the rarefied case.

§ CONCLUSION

This paper marks the first comparison of velocity space model reduction techniques for rarefied flows using proper orthogonal decomposition (POD) and a fully connected neural network (FCNN).

As physically expected, the rarefied regime needs more modes than the hydrodynamic regime. Choosing three and five intrinsic variables for the hydrodynamic and rarefied case, respectively, leads to less than one percent error.
The FCNN is initially more accurate but has a remaining error even for larger number of intrinsic variables, while POD achieves subsequently higher accuracy  with increasing the number of parameters.
The resulting errors of the macroscopic variables are small, especially in the smoother rarefied case.

Even though not strictly enforced, the FCNN approximately exhibits the correct conservation of mass, momentum, and energy, while POD has a slightly larger error.

The correlation of intrinsic variables and macroscopic variables is investigated by means of the evolution of the reconstructed values and the pairwise correlation for the FCNN. The density is directly included in the latent space while the relation with other macroscopic variables seems to be more complex.

In addition to optimizing the neural network's performance, our research endeavors to enhance its predictive capabilities by incorporating fundamental physical properties such as conservation and interpretability. This objective includes the exploration of different loss functions and the regularization techniques.
Furthermore, future work can utilize the model reduction approach in real-world scenarios, including simulations and parameter predictions.
Lastly, the presented concepts could be applied in related fields, such as shallow flows [16] or fusion plasmas [21].

§ AUTHOR CONTRIBUTION STATEMENT (CREDIT)

All authors contributed to this publication.
The authors contributions differ in the following points:

Julian Koellermeier:     initial idea, methodology, setup of test cases, reference data, writing

Philipp Krah:     initial idea, methodology, implementation, visualization, writing

Julius Reiss:     initial idea, supervision

Zachary Schellin:     implementation, training, visualization

§ CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.

§ CODE AND DATA AVAILABILITY

All data and scripts to reproduce the results can be obtained from the authors on reasonable request.

§ ACKNOWLEDGEMENT

The authors would like to acknowledge the financial support of the CogniGron research center and the Ubbo Emmius Funds (University of Groningen).
Centre de Calcul Intensif d’Aix-Marseille is acknowledged for granting access to its high performance computing resources.
P. Krah acknowledges partial funding from the Agence Nationale de la Recherche (ANR), project CM2E, grant ANR-20-CE46-0010-01.

[1]
F. Bernard, A. Iollo, and S. Riffaud.
Reduced-order model for the BGK equation based on POD and optimal
Journal of Computational Physics, 373:545–570, 2018.

[2]
P. L. Bhatnagar, E. P. Gross, and M. Krook.
A model for collision processes in gases. 1. Small amplitude
processes in charged and neutral one-component systems.
Physical Review, 94:511–525, 1954.

[3]
S. L. Brunton and J. N. Kutz.
Data driven science and engineering.
Cambridge University Press, 2019.

[4]
M. Cranmer.
Interpretable Machine Learning for Science with PySR and
ArXiv, 2023.

[5]
K. Debrabant, G. Samaey, and P. Zieliński.
A micro-macro acceleration method for the monte carlo simulation of
stochastic differential equations.
SIAM Journal on Numerical Analysis, 55(6):2745–2786, 2017.

[6]
C. Eckart and G. Young.
The approximation of one matrix by another of lower rank.
Psychometrika, 1(3):211–218, 1936.

[7]
L. Einkemmer.
A low-rank algorithm for weakly compressible flow.
SIAM Journal on Scientific Computing, 41(5):A2795–A2814, 2019.

[8]
L. Einkemmer, J. Hu, and Y. Wang.
An asymptotic-preserving dynamical low-rank method for the
multi-scale multi-dimensional linear transport equation.
Journal of Computational Physics, 439:110353, Aug. 2021.

[9]
L. Einkemmer, J. Hu, and L. Ying.
An efficient dynamical low-rank algorithm for the boltzmann-bgk
equation close to the compressible viscous flow regime.
SIAM Journal on Scientific Computing, 43(5):B1057–B1080, 2021.

[10]
F. Fan, J. Xiong, M. Li, and G. Wang.
On interpretability of artificial neural networks: A survey.
IEEE Transactions on Radiation and Plasma Medical Sciences,
5:741–760, 2020.

[11]
Y. Fan, J. Koellermeier, J. Li, R. Li, and M. Torrilhon.
Model reduction of kinetic equations by operator projection.
Journal of Statistical Physics, 162(2):457–486, 2016.

[12]
I. Goodfellow, Y. Bengio, and A. Courville.
Deep Learning.
MIT Press, 2016.

[13]
O. Koch and C. Lubich.
Dynamical low-rank approximation.
SIAM Journal on Matrix Analysis and Applications,
29(2):434–454, 2007.

[14]
J. Koellermeier.
Derivation and numerical solution of hyperbolic moment equations
for rarefied gas flows.
PhD thesis, RWTH Aachen University, 2017.

[15]
J. Koellermeier, P. Krah, and J. Kusch.
Macro-micro decomposition for consistent and conservative model order
reduction of hyperbolic shallow water moment equations: A study using
pod-galerkin and dynamical low rank approximation, 2023.

[16]
J. Koellermeier and M. Rominger.
Analysis and numerical simulation of hyperbolic shallow water moment
Communications in Computational Physics, 28(3):1038–1084,

[17]
J. Koellermeier, R. P. Schaerer, and M. Torrilhon.
A framework for hyperbolic approximation of kinetic equations using
quadrature-based projection methods.
Kinetic and Related Models, 7(3):531–549, 2014.

[18]
J. Koellermeier and M. Torrilhon.
Numerical solution of hyperbolic moment models for the boltzmann
European Journal of Mechanics - B/Fluids, 64:41–46, July 2017.

[19]
J. Koellermeier and M. Torrilhon.
Numerical study of partially conservative moment equations in kinetic
Communications in Computational Physics, 21(4):981–1011, 2017.

[20]
J. Koellermeier and M. Torrilhon.
Two-dimensional simulation of rarefied gas flows using
quadrature-based moment equations.
Multiscale Modeling and Simulation, 16(2):1059–1084, 2018.

[21]
P. Krah, X.-Y. Yin, J. Bergmann, J.-C. Nave, and K. Schneider.
A characteristic mapping method for vlasov–poisson with extreme
resolution properties.
in preparation, 2023.

[22]
K. Kunisch and S. Volkwein.
Control of the burgers equation by a reduced-order approach using
proper orthogonal decomposition.
Journal of optimization theory and applications, 102:345–371,

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

[24]
R. J. LeVeque.
Finite Volume Methods for Hyperbolic Problems.
Cambridge Texts in Applied Mathematics. Cambridge University Press,

[25]
V. Maes, W. Dekeyser, J. Koellermeier, M. Baelmans, and G. Samaey.
Hilbert expansion based fluid models for kinetic equations
describing neutral particles in the plasma edge of a fusion device.
Physics of Plasmas, 30(6):063907, 06 2023.

[26]
R. G. McClarren and C. D. Hauck.
Robust and accurate filtered spherical harmonics expansions for
radiative transfer.
Journal of Computational Physics, 229(16):5597–5614, 2010.

[27]
L. Mieussens.
Discrete velocity model and implicit scheme for the BGK equation of
rarefied gas dynamics.
Mathematical Models and Methods in Applied Sciences,
10(08):1121–1149, 2000.

[28]
L. Mirsky.
Symmetric gauge functions and unitarily invariant norms.
The quarterly journal of mathematics, 11(1):50–59, 1960.

[29]
A. Pinkus.
Approximation theory of the mlp model in neural networks.
Acta Numerica, 8:143–195, 1999.

[30]
S. B. &. C. Prigent.
Local discrete velocity grids for multi-species rarefied flow
Communications in Computational Physics, 28(4):1274–1304, June

[31]
L. Sirovich.
Turbulence and the dynamics of coherent structures. part i-iii.
Quarterly of Applied Mathematics, 45(3):561–571, 1987.

[32]
G. A. Sod.
Review. A Survey of Several Finite Difference Methods for Systems of
Nonlinear Hyperbolic Conservation Laws.
Journal of Computational Physics, 27(1):1–31, Apr. 1978.

[33]
H. Struchtrup and M. Torrilhon.
Higher-order effects in rarefied channel flows.
Physical Review E, 78(4):46301, 2008.

[34]
M. Torrilhon.
Convergence study of moment approximations for boundary value
problems of the Boltzmann-BGK equation.
Communications in Computational Physics, 18(3):529–557, 2015.

[35]
M. Torrilhon.
Modeling nonequilibrium gas flow based on moment equations.
Annual Review of Fluid Mechanics, 48(1):429–458, 2016.

[36]
L. Van Der Maaten, E. Postma, J. Van den Herik, et al.
Dimensionality reduction: a comparative.
J Mach Learn Res, 10(66-71), 2009.

§ APPENDIX: HYPERPARAMETERS FOR THE FULLY CONNECTED AUTOENCODER

In this section, we describe the tests that have been conducted to optimize our hyperparameters. The hyperparameters include: the number of layers, i.e., depth; the number of nodes per hidden layer, i.e., width, batch size, and non-linear activation functions; the number of epochs for training; and the learning rate. Experiments are evaluated through: the validation error which estimates the model's ability to generalize; the training error which estimates the optimization to training data;
and the $\L2$-error, which gives an estimate of how well the model performs on the whole dataset and hence is applied as a comparative metric against POD.

To start with a working model, an estimate over the initial hyperparameters is done, which are summarized in <Ref>. These include a mini-batch size of 16, the width of the bottleneck layer is 3 and 5  for the hydrodynamic case \(\hy\) and the rarefied case \(\rare\), respectively, and a learning rate of 0.0001. LeakyReLU is applied as an activation function for the output, input, and any hidden layer besides the output of the bottleneck layer, and is referred to as activation hidden. The hyperbolic tangent is applied as an activation function for the output of the last hidden layer in the encoder which outputs the code, referred to as the activation code. Moreover, 2000 initial number of epochs are used. This might appear exaggerated but is justified by the small amount of input data and the small size of the network which yields fast training.

The initial selection of batch size, bottleneck size, number of epochs, learning rate, and applied activation functions.

Mini-batch size       Intr. dim.       Epochs    Learn. rate     Activ. hidden/code
16 		   	3/5         2000   	    0.0001     LeakyReLU/Tanh

Five designs for finding an optimal number of layers, i.e., depth, are explored. These are as follows:

* 10 layers with layer widths: 40, 40, 20, 10, 5, 3/5, 5, 10, 20, 40, 40.

* 8 layers with layer widths: 40, 40, 20 , 10, 3/5, 10, 20, 40, 40.

* 6 layers with layer widths: 40, 40, 20 , 3/5, 20, 40, 40.

* 4 layers with layer widths: 40, 40, 3/5, 40, 40.

* 2 layers with layer widths: 40, 3/5, 40.

The model's depth is determined in a primary step because it determines the model's representational capacity and therefore can initiate over- and underfitting at an early stage in the hyperparameter search. The results of the experimentation are shown in <ref> and <ref> for both rarefaction levels.

Results for the variation of the depth. Given are minimum values of training and validation error as well as the \(\L2\) error. The minima were reached around the last 50 epochs of the training. The \(\L2\) error is evaluated with the model at the last epoch.

Depth     2cMinimum training error     2cMinimum validation error     2c\(\L2\) error
\(\hy\)   \(\rare\)   \(\hy\)   \(\rare\)   \(\hy\)   \(\rare\)

10    1.53e-7     5.96e-7     2.22e-7     5.19e-7     0.0048     0.0091

8     1.17e-7     2.05e-7     1.58e-7     2.32e-7     0.0041     0.0054

6     9.76e-8     1.40e-7     1.49e-7     1.72e-7     0.0038     0.0045

4     6.29e-8     1.52e-7     7.74e-8     1.61e-7     0.0031     0.0048

2     1.29e-6     3.29e-6     1.37e-6     3.42e-6     0.0136     0.0217

For the hydrodynamic case \(\hy\), the lowest validation error of 7.74e-8 and an \(\L2\) error of 0.0031 is reached with 4 layers and constitutes the best-performing design. Additionally, as seen in <ref>(left), a design that exceeds 4 layers results in slight overfitting from the 500th epoch. Less than 4 layers do not reach the validation error and \(\L2\) error of the other designs, yielding the conclusion, that the capacity is too low. Overfitting occurs with 4 layers only after the 1000th epoch and is of smaller magnitude compared to the other three models that show overfitting.

For rarefied case \(\rare\), the lowest validation error of 1.61e-7 is also reached with 4 layers. On the other hand, the lowest \(\L2\) error of 0.0031 and the lowest training error of 1.40e-7 are reached with 6 layers. Contrary to the previously discussed hydrodynamic case, the training error and \(\L2\) errors are of lower magnitude for 6 layers, except for the validation error. Looking at <ref>(right), it is observed that the model with 6 layers starts to overfit after the 1500 epochs, yielding a decreasing training error and a stagnating validation error. Hence the model improved in the optimization task which additionally improves the \(\L2\) error. Its generalization ability, measured by the validation error, did not improve and is larger than the validation error reached with 4 layers. This concludes a model with 4 layers constitutes the best-performing.

Qualitatively, the overall training for both rarefaction levels is very stable. Training and validation errors do not diverge excessively and converge early in training. Separation of training and validation error occurs prominently for the hydrodynamic solution.

The width of the two remaining hidden layers is examined in the following. For both the hydrodynamic and the rarefied regime five experiments are conducted, lowering the hidden units of the hidden layers from fifty to ten. Note that the decoder is chosen to be structurally a reflection of the encoder. Therefore only one parameter is changed. Results for the hydrodynamic case \(\hy\) and the rarefied case \(\rare\) are shown in <ref>. Note that the contribution of over-and underfitting is negligible and therefore the training error is omitted. A model with 30 hidden units in encoder and decoder performs best with the hydrodynamic case \(\hy\) and reaches a validation error of \(\num{1.77e-08}\). The corresponding \(\L2\) error is equal to \( \num{1.5e-3}\) with a shrinkage factor of 0.015. Overall, the loss of each experiment with the hydrodynamic case \(\hy\) is quite similar and ranges from \(\num{1.77e-8}\) to \(\num{5.11e-8}\). The \(\L2\) error behaves similarly and is even equal for 50 and 30 layers. A model with 40 hidden units performs best for the rarefied case \(\rare\). The corresponding validation error is \(\num{1.65e-8}\) with \(\L2=\num{1.4e-3}\), which is smaller than for the hydrodynamic case \(\hy\). The shrinkage factor here is 0.125. In all experiments, a model with 10 hidden nodes performs worst. Training and validation errors over 4000 epochs for both experiments can be seen in <ref>.

Results for the variation of the width. Given are the minimum value of the validation error as well as the \(\L2\) error. The minima were reached around the last 50 epochs of the training, the \(\L2\) error is evaluated with the model at the last epoch.

Hidden units     2cValidation error     2c$\L2$     2cShrinkage factor
\(\hy\)   \(\rare\)   \(\hy\)   \(\rare\)   \(\hy\)   \(\rare\)

50     1.91e-8      5.05e-8     0.0015      0.0025     0.06      0.01

40     2.65e-08     1.65e-8     0.0018      0.0014     0.075     0.125

30     1.77e-08     3.40e-8     0.0015      0.0021     0.015     0.0167

20     2.50e-08     5.25e-8     0.0017      0.0027     0.1       0.25

10     5.11e-08     3.97e-7     0.0025      0.0077     0.3       0.5

Next, the mini-batch size is analyzed. Results are displayed in <ref>. Experiments are conducted with mini-batch sizes of 2, 4, 8, 16, 32. The smallest batch size of 2 yields the lowest validation error of \(\num{1.15e-8}\) with corresponding \(\L2 = 0.0012\) at epoch 4956 for the hydrodynamic case \(\hy\). The lowest validation error with \(\num{6.30e-9}\) is achieved for the rarefied case \(\rare\) at epoch 4534 with a batch size of 4. The corresponding \(\L2\) error equals  \(0.001\). Small batch sizes have a regularizing effect on the training and therefore are beneficial to generalization. At the same time, the lower the batch size is, the more unstable is the training as seen in <ref>. The oscillations that begin with batch sizes of 8 and lower, which make the training unstable, can be cured with a lower learning rate as soon as training starts to tremble. Additionally, small batch sizes drastically increase training time, thus a batch size as low as 2 is not used for the experiments. In conclusion, a batch size of 4 is chosen. Furthermore, a reduction of the learning rate from \(\num{1e-4}\) to \(\num{1e-5}\) is applied after the 3000th epoch.

Results for the variation of the batch sizes. Given are the minimum value of validation error as well as the corresponding epoch. Additionally, the \(\L2\) error is given but evaluated with the model at the last epoch.

Batch Size     2cValidation error     2c$\L2$    2cEpoch
\(\hy\)   \(\rare\)   \(\hy\)   \(\rare\)   \(\hy\)   \(\rare\)

32    5.40e-8     2.17e-8     0.0024      0.0017   4998   4992

16    1.95e-8     2.06e-8     0.0015      0.0016   4999   5000

8     2.25e-8     1.03e-8     0.0017      0.0012   4965   4961

4     1.52e-8     6.30e-9     0.0013      0.0010   3956   4534

2     1.15e-8     9.18e-9     0.0012      0.0013   4956   4872

Eight experiments with different activation functions, ReLU, ELU, Tanh, SiLU, and LeakyReLU, are performed. The experiment designs and results are given in <ref> for hidden and code activations. With the hydrodynamic case \(\hy\), combining ELU and ELU for hidden and code activation, respectively, yields the best result for the validation error with 4.44e-9 and a corresponding \(\L2\) error of 0.0008. These values are achieved at the last epoch. For the rarefied case \(\rare\), a combination of ReLU and ReLU for hidden and code activation, respectively, produces a validation error of 7.18e-9 and a corresponding \(\L2\) error of 0.0009. Both are reached close to the last epoch. Note that all models reach their lowest loss at or very close to the last epoch. The reason is the stable training after the 3000th epoch, where the learning rate is lowered to 1e-5 as seen in <ref>. This measure shows in all experiments an immediate success for learning. Both validation and training errors fall at the 3001st epoch and only decrease slightly thereafter. This behavior clearly shows that the updates to the free parameters \(\frepar\) were too big, which prohibitively slowed down or even prevented the learning process. Small updates to \(\frepar\) made all models quickly reach a minimum.

Results for the variation of the activation functions for the hidden-/code layers. Given are the minimum value of validation error as well as the corresponding epoch and the \(\L2\) error. The \(\L2\) error is evaluated with the models saved when the minimum validation error was achieved during training.

Activ. hidden/code     2cValidation error     2c$\L2$    2cEpoch
\(\hy\)   \(\rare\)   \(\hy\)   \(\rare\)   \(\hy\)   \(\rare\)

ReLU/ReLU 	           9.79e-9     7.18e-9     0.0010      0.0009   5000    4998

ELU/ELU                4.44e-9     1.11e-8     0.0008      0.0012   5000    5000

Tanh/Tanh 	           7.83e-9     2.58e-8     0.0011      0.0018   5000    5000

SiLU/SiLU 	           7.69e-9     1.37e-8     0.0011      0.0013   5000    5000

LeakyReLU/LeakyReLU    1.86e-8     9.39e-9     0.0015      0.0010   5000    4997

ELU/Tanh               5.49e-9     1.87e-8     0.0008      0.0014   5000    5000

LeakyReLU/Tanh         1.00e-8     1.42e-8     0.0010      0.0012   4997    4992

ELU/SiLU               8.11e-9     1.93e-8     0.0011      0.0015   5000    5000

The final hyperparameters for both input data are summarized below in <ref>. From the initial models to the final models, the decrease in the validation error gained \(\approx \num{1.5e-7}\) for hydrodynamic case \(\hy\) and \(\approx \num{7.2e-8}\) for rarefied case \(\rare\) which amount to 93% of the initial values for both models.

Five experiments over the different depths with the hydrodynamic case $\hy$ left and the rarefied case $\rare$ right. The number of layers used for every experiment is given. Training and validation loss are shown over 2000 epochs.

Five experiments over different widths with the hydrodynamic case $\hy$ left and the rarefied case $\rare$ right. The number of nodes used for every experiment is given. Training and validation loss are shown over 4000 epochs.

Five experiments over different batch sizes with the hydrodynamic case $\hy$ left and the rarefied case $\rare$ right. The batch size used for every experiment is given. Training and validation loss are shown over 5000 epochs.

Eight experiments with different combinations of activation functions for the hydrodynamic case $\hy$ left and the rarefied case $\rare$ right. Shown are training- and validation error over 5000 epochs.