Deep learning observables in computational fluid dynamics

TL;DR

This paper introduces a deep learning approach to efficiently approximate the input parameters to observable map in computational fluid dynamics, significantly reducing computational costs for uncertainty quantification and related tasks.

Contribution

It develops a neural network-based method that accurately predicts the input-to-observable map with few training samples and combines it with Monte Carlo methods for faster uncertainty propagation.

Findings

Neural networks achieve low prediction errors with minimal training samples.

The combined deep learning and Monte Carlo methods outperform traditional approaches by one to two orders of magnitude.

Numerical experiments validate the efficiency and accuracy of the proposed approach.

Abstract

Many large scale problems in computational fluid dynamics such as uncertainty quantification, Bayesian inversion, data assimilation and PDE constrained optimization are considered very challenging computationally as they require a large number of expensive (forward) numerical solutions of the corresponding PDEs. We propose a machine learning algorithm, based on deep artificial neural networks, that predicts the underlying \emph{input parameters to observable} map from a few training samples (computed realizations of this map). By a judicious combination of theoretical arguments and empirical observations, we find suitable network architectures and training hyperparameters that result in robust and efficient neural network approximations of the parameters to observable map. Numerical experiments are presented to demonstrate low prediction errors for the trained network networks, even when the network has been trained with a few samples, at a computational cost which is several orders of magnitude lower than the underlying PDE solver. Moreover, we combine the proposed deep learning algorithm with Monte Carlo (MC) and Quasi-Monte Carlo (QMC) methods to efficiently compute uncertainty propagation for nonlinear PDEs. Under the assumption that the underlying neural networks generalize well, we prove that the deep learning MC and QMC algorithms are guaranteed to be faster than the baseline (quasi-) Monte Carlo methods. Numerical experiments demonstrating one to two orders of magnitude speed up over baseline QMC and MC algorithms, for the intricate problem of computing probability distributions of the observable, are also presented.

Tables

Table 1. Table 1: Relative percentage Prediction errors ( 2.21 ) with p = 2 𝑝 2 p=2 , for the trained neural networks in the sum of sines experiment ( 2.23 ).

Map	Error (Mean) in $\%$	Error (Std.) in $\%$
${\sum }_{i=1}^6\sin \left(4\pi x_i\right)$	$133.95$	$417.83$
${\sum }_{i=1}^6\frac{\sin \left(4\pi x_i\right)}{i^3}$	$43.66$	$41.26$

Table 2. Table 2: Hyperparameters in the training algorithm

1.	Number of hidden layers $\left(K-1\right)$
2.	Number of units in $k$-th layer $\left(d_k\right)$
3.	Exponent $p$ in the loss function (2.11).
4.	Exponent $q$ in the regularization term in (2.12)
5.	Value of regularization parameter $\lambda$ in (2.12)
6.	Choice of optimization algorithm (optimizer) – either standard SGD or ADAM.
7.	Initial guess $\overline{\theta }$ in the SGD method (2.13).

Table 3. Table 3: Reference neural network architecture for the flow past airfoil problem.

Layer	Width (Number of Neurons)	Number of parameters
Hidden Layer 1	12	84
Hidden Layer 2	12	156
Hidden Layer 3	10	130
Hidden Layer 4	12	132
Hidden Layer 5	10	130
Hidden Layer 6	12	156
Hidden Layer 7	10	130
Hidden Layer 8	10	110
Hidden Layer 9	12	132
Output Layer	1	13
		1149

Table 4. Table 4: The hyperparameter configurations that correspond to the best performing network (one with least mean prediction error) for the two observables of the flow past airfoils.

Obs	Opt	Loss	$L^1$-reg	$L^2$-reg	Selection	BP. Err mean (std)
Lift	ADAM	MSE	$0.0$	$7.8\times {10}^{-6}$	wass-train	$0.786$ ($0.010$)
Drag	ADAM	MAE	$0.0$	$7.8\times {10}^{-6}$	mean-train	$1.847$ ($0.022$)

Table 5. Table 5: Computational times (cost) of (single) sample generation, network training and (single) network evaluation, for the flow past airfoils problem. Each sample was generated using a parallelized finite-volume solver on 16 Intel(R) Xeon(R) Gold 5118 @ 2.30GHz processor cores. The entire run took 1500 1500 1500 secs on the cluster and translates into a total wall clock time of 1500 × 16 = 24000 1500 16 24000 1500\times 16=24000 secs. The training and evaluation of the neural networks were all performed on a Intel(R) Core(TM) i7-8700K CPU @ 3.70GHz machine. The training and evaluation times are approximations of the average runtimes over all hyperparameter configurations.

	Time (in secs)
Sample generation	$24000$
Training (Lift)	$700$
Evaluation (${ℒ}_1$)	$9\times {10}^{-6}$
Training (Drag)	$840$
Evaluation (${ℒ}_2$)	${10}^{-5}$

Table 6. Table 6: Sensitivity of the best performing networks (listed in table 4 ) to retrainings i.e starting values for ADAM for flow past airfoil. All errors are relative mean L 2 superscript 𝐿 2 L^{2} prediction error ( 2.21 ) in percentage and we list the minimum, maximum, mean and standard deviation of the error over 5 5 5 retrainings.

Obs	Err (min)	Err (max)	Err (mean)	Err (std)
Lift	$0.786$	$1.336$	$1.004$	$0.187$
Drag	$1.847$	$8.016$	$3.841$	$2.332$

Table 7. Table 7: Performance of trained neural networks with respect to the variation of network size for the flow past airfoils problem. The rows represent variation with respect to Width (number of neurons per hidden layer) and the columns represent variation with respect to Depth (number of hidden layers). For each entry of the Width-Depth matrix, we tabulate the vector of relative percentage mean prediction error for the best performing networks. The components of each vector represent error in [ Lift , Drag ] Lift Drag \left[{\rm Lift},{\rm Drag}\right] and ¬ \neg is used to indicate that the training procedure failed for this particular configuration.

Width/Depth	$6$	$12$	$24$
$4$	$[1.21,4.37]$	$[0.83,1.74]$	$[1.09,2.48]$
$8$	$[0.84,3.66]$	$[0.86,2.02]$	$[0.81,1.69]$
$16$	$[0.80,\neg ]$	$[\neg ,74.72]$	$[\neg ,5.11]$

Table 8. Table 8: Real speedups ( 5.6 ), for the lift and the drag in the flow past airfoils problem, comparing DLQMC with baseline QMC algorithm.

Observable	Speedup
Lift	$6.64$
Drag	$8.56$

Table 9. Table 9: Real speedups, for the lift and the drag in the flow past airfoils problem, comparing DLQMC with the baseline MC algorithm.

Observable	Speedup over MC
Lift	$246.02$
Drag	$179.54$

Table 10. Table 10: The hyperparameter configurations that correspond the best performing network (one with least mean prediction error) for the two observables of the flow past airfoils, trained on random (Monte Carlo) training set.

Obs	Opt	Loss	$L^1$-reg	$L^2$-reg	Selection	BP. Err mean
Lift	SGD	MAE	$0.0$	$7.8\times {10}^{-7}$	mean-train	$8.487$
Drag	ADAM	MAE	$7.8\times {10}^{-7}$	$0.0$	train	$20.25$

Table 11. Table 11: The hyperparameter configurations that correspond to the best performing network (one with least mean prediction error) for the three observables of the stochastic shock tube problem.

Obs	Opt	Loss	$L^1$-reg	$L^2$-reg	Selection	BP. Err mean
${ℒ}_1$	ADAM	MSE	$7.8\times {10}^{-6}$	$0$	wass-train	$2.212$
${ℒ}_2$	ADAM	MSE	$7.8\times {10}^{-6}$	$0$	wass-train	$6.575$
${ℒ}_3$	ADAM	MSE	$7.8\times {10}^{-5}$	$0$	mean-train	$5.467$

Equations

$${\mathbf{U}}_t+{\mathrm{div}}_x(\mathbf{F}(\mathbf{U}))=\nu {\mathrm{div}}_x(\mathbf{D}(\mathbf{U}){\nabla }_x\mathbf{U}),(x,t)\in D\subset {\mathbb{R}}^{d_s}\times {\mathbb{R}}_{+},$$

$$L(\mathbf{U})=\underset{D\times {\mathbb{R}}_{+}}{\int }\psi (x,t)g(\mathbf{U}(x,t))dxdt.$$

$${\partial }_t\mathbf{U}(t,x,y)+{\mathrm{div}}_x(\mathbf{F}(y,\mathbf{U}))$$

$$\mathbf{U}(0,x,y)$$

$$L_b\mathbf{U}(t,x,y)$$

$$L_g(y,\mathbf{U}):=\underset{0}{\overset{T}{\int }}\underset{D_y}{\int }\psi (x,t)g(\mathbf{U}(t,x,y))dxdt,\mathrm{for}\mu a.ey\in Y.$$

$$\mathcal{L}:y\in Y\to \mathcal{L}(y)=L_g(y,\mathbf{U}),$$

$$\mathcal{L}:y\in Y\to \mathcal{L}(y)=L_g(y,{\mathbf{U}}^{\Delta }),$$

$$\Vert \mathcal{L}(y)-{\mathcal{L}}^{\Delta }(y){\Vert }_{L_{\mu }^p(Y)}<\epsilon ,$$

$$\Vert f{\Vert }_{L_{\mu }^p(Y)}:={\left({\int }_Y|f(y){|}^{p}d\mu (y)\right)}^{\frac{1}{p}},$$

$${\mathcal{L}}_{\theta }(z)=C_K\odot \sigma \odot C_{K-1}\dots \dots \dots \odot \sigma \odot C_2\odot \sigma \odot C_1(z).$$

$$\sigma (z)=\max (z,0).$$

$$C_k{z}_k=W_k{z}_k+b_k,\mathrm{for}{W}_k\in {\mathbb{R}}^{d_{k+1}\times d_k},z_k\in {\mathbb{R}}_{d_k},b_k\in {\mathbb{R}}^{d_{k+1}}.$$

$$\theta =\{W_k,b_k\},{\theta }_W=\{W_k\}\forall 1⩽k⩽K,$$

$$M=\sum _{k=1}^{K-1}(d_k+1)d_{k+1}.$$

$$J(\theta ):=\sum _{y\in \mathcal{S}}|{\mathcal{L}}^{\Delta }(y)-{\mathcal{L}}_{\theta }(y){|}^p,$$

$${\theta }^{\ast }=\arg \underset{\theta \in \Theta }{\min }\left(J(\theta )+\lambda \mathcal{R}(\theta )\right).$$

$${\theta }_{r+1}={\theta }_r-{\eta }_r{\nabla }_{\theta }\left(J({\theta }_k)+\lambda \mathcal{R}({\theta }_k)\right),$$

$${\nabla }_{\theta }J({\theta }_k)\approx {\nabla }_{\theta }\left(\sum _{y\in {\hat{\mathcal{S}}}_q}|{\mathcal{L}}^{\Delta }(y)-{\mathcal{L}}_{\theta }(y){|}^p\right),$$

$$\Vert f-NN(f){\Vert }_{L^p}⩽\epsilon ,1⩽p⩽\infty .$$

$$\Vert \mathcal{L}{\Vert }_{W^{k,p}(Y)}⩽C\Vert \mathbf{U}{\Vert }_{W^{k,p}(D\times (0,T)\times Y)},$$

$$\Vert \mathcal{L}{\Vert }_{BV(Y)}⩽C.$$

$$\Vert \mathcal{L}{\Vert }_{W^{1,\infty }(Y)}⩽C.$$

$$\Vert {\mathcal{L}}^{\Delta }{\Vert }_{W^{1,\infty }(Y)}⩽C.$$

$$|\mathcal{L}(y)-\mathcal{L}(y^{\ast })|$$

$$⩽\Vert \psi {\Vert }_{\infty }\Vert g{\Vert }_{\mathrm{Lip}}\underset{0}{\overset{T}{\int }}\Vert u_0(.,y)-u_0(.,y^{\ast }){\Vert }_{1}dt,\mathrm{by}{L}^1-\mathrm{contractivity}$$

$$⩽\Vert \psi {\Vert }_{\infty }\Vert g{\Vert }_{\mathrm{Lip}}T\Vert u_0{\Vert }_{W^{1,\infty }(Y,L^1(D))}⩽C$$

$${\mathcal{E}}_G({\theta }^{\ast }):={\left(\underset{Y}{\int }|{\mathcal{L}}^{\Delta }(y)-{\mathcal{L}}_{{\theta }^{\ast }}(y){|}^{p}d\mu (y)\right)}^{\frac{1}{p}},1⩽p⩽\infty .$$

$${\mathcal{E}}_P({\theta }^{\ast }):={\left(\frac{1}{\#(\mathcal{T})}\sum _{y\in \mathcal{T}}|{\mathcal{L}}^{\Delta }(y)-{\mathcal{L}}_{{\theta }^{\ast }}(y){|}^p\right)}^{\frac{1}{p}}$$

$${\mathcal{E}}_G({\theta }^{\ast })⩽\sqrt{\frac{U}{N}},$$

$${\mathcal{L}}^1(y)=\sum _{i=1}^6\sin (4\pi y_i),{\mathcal{L}}^2(y)=\sum _{i=1}^6\frac{1}{i^3}\sin (4\pi y_i).$$

$${\mathcal{L}}^{\Delta }(z)={\mathcal{L}}_{\hat{\theta }}(z),\forall z\in \mathcal{S}.$$

Full text

Moreover, we combine the proposed deep learning algorithm with Monte Carlo (MC) and Quasi-Monte Carlo (QMC) methods to efficiently compute uncertainty propagation for nonlinear PDEs. Under the assumption that the underlying neural networks generalize well, we prove that the deep learning MC and QMC algorithms are guaranteed to be faster than the baseline (quasi-) Monte Carlo methods. Numerical experiments demonstrating one to two orders of magnitude speed up over baseline QMC and MC algorithms, for the intricate problem of computing probability distributions of the observable, are also presented.

1 Introduction

Many interesting fluid flows are modeled by so-called convection-diffusion equations [29] i.e, nonlinear partial differential equations (PDEs) of the generic form,

[TABLE]

with ${\bf U}\in\mathbb{R}^{m}$ denoting the vector of unknowns, ${\bf F}=(F_{i})_{1\leqslant i\leqslant d_{s}}$ the flux vector, ${\bf D}=(D_{ij})_{1\leqslant i,j\leqslant d_{s}}$ the diffusion matrix and $\nu$ a small scale parameter representing kinematic viscosity.

Prototypical examples for (1.1) include the compressible Euler equations of gas dynamics, shallow water equations of oceanography and the magnetohydrodynamics (MHD) equations of plasma physics. These PDEs are hyperbolic systems of conservation laws i.e, special cases of (1.1) with $\nu=0$ and with the flux Jacobian $\partial_{{\bf U}}({\bf F}\cdot n)$ having real eigenvalues for all normal vectors $n$. Another important example for (1.1) is provided by the incompressible Navier-Stokes equations, where $0<\nu<<1$ and the flux function is non-local on account of the divergence-free constraint on the velocity field.

It is well known that solutions to convection-diffusion equations (1.1) can be very complicated. These solutions might include singularities such as shock waves and contact discontinuities in the case of hyperbolic systems of conservation laws. For small values of $\nu$, these solutions can be generically unstable, even turbulent, and contain structures with a large range of spatio-temporal scales [29].

Numerical schemes play a key role in the study of fluid flows and a large variety of robust and efficient numerical methods have been designed to approximate them. These include (conservative) finite difference [31], finite volume [21, 24], discontinuous Galerkin (DG) finite element [24] and spectral (viscosity) methods [50]. These methods have been extremely successful in practice and are widely used in science and engineering today.

The exponential increase in computational power in the last decades provides us with the opportunity to solve very challenging large scale problems in computational fluid dynamics, such as uncertainty quantification (UQ) [4, 19], (Bayesian) inverse problems [48] and real-time optimal control, design and PDE constrained (shape) optimization [6, 51]. In such problems, one is not always interested in computing the whole solution field ${\bf U}$ of (1.1). Rather and in analogy with experimental measurements, one is interested in computing the so-called observables (functionals or quantities of interest) for the solution ${\bf U}$ of (1.1). These can be expressed in the generic form,

[TABLE]

Here $\psi:D\times\mathbb{R}_{+}\rightarrow\mathbb{R}$ and $g:\mathbb{R}^{m}\rightarrow\mathbb{R}$ are suitable test functions. Prototypical examples of such observables (functionals) are provided by body forces, such as the lift and drag in aerodynamic simulations, and by the runup height (at certain probe points) in simulations of tsunamis.

Moreover in practice, one is interested not just in a single value, but rather in the statistics of such observables. Typical statistical quantities of interest are the mean, variance, higher-moments and probability density functions (pdfs) of the observable (1.2). Such statistical quantities quantify uncertainty in the solution, propagated from possible uncertainty in the underlying inputs i.e, fluxes, diffusion coefficients, initial and boundary data. They might also stem from the presence of a statistical spread in design parameters such as those describing the geometry of the computational domain. It is customary to represent the resulting solution field as ${\bf U}={\bf U}(y)$, with $y\in Y\subset\mathbb{R}^{d}$, a possibly very high dimensional parameter space. Thus, the goal of many CFD simulations is to compute (statistics of) the so-called parameters to observable map $y\mapsto L({\bf U}(y))$.

Calculating a single realization of this parameters to observable map might require a very expensive forward solve of (1.1) with a CFD solver and a quadrature to calculate (1.2). However, in UQ, Bayesian inversion or optimal design and control, one needs to evaluate a large number of instances of this map, necessitating a very high computational cost, even on state of the art HPC systems. As a concrete example, we consider a rather simple yet prototypical situation. We are interested in computing statistics of body forces such as lift and drag for a model two-dimensional RAE2822 airfoil [25] (see figure 2). A single forward solve of the underlying compressible Euler equations on this airfoil geometry on a mesh of high resolution (see figure 2), with a state of the art high-resolution finite volume solver [43], takes approximately $7$ wall-clock hours on a HPC cluster (see table 5). This cost scales up in three space dimensions proportionately. However, a typical UQ problem such as Bayesian inversion with a state of the art MCMC algorithm, for instance the Random walk Metropolis-Hastings algorithm [48], might need upto $10^{5}-10^{6}$ such realizations, rendering even a two-dimensional Bayesian inversion infeasible ! This example illustrates the fact that the high-computational cost of the parameters to observable map makes problems such as forward UQ, Bayesian inversion, and optimal control and design very challenging. Hence, we need numerical methods which allow ultrafast (several order of magnitude faster than state of the art) computations of the underlying parameters to observable map.

Machine learning, in the form of artificial neural networks (ANNs), has become extremely popular in computer science in recent years. This term is applied to methods that aim to approximate functions with layers of units (neurons), connected by (affine) linear operations between units and nonlinear activations within units, [20] and references therein. Deep learning, i.e an artificial neural network with a large number of intermediate (hidden) layers has proven extremely successful at diverse tasks, for instance in image processing, computer vision, text and speech recognition, game intelligence and more recently in protein folding [15], see [30] and references therein for more applications of deep learning. A key element in deep learning is the training of tunable parameters in the underlying neural network by (approximately) minimizing suitable loss functions. The resulting (non-convex) optimization problem, on a very high dimensional underlying space, is customarily solved with variants of the stochastic gradient descent method [45].

Deep learning is being increasingly used in the context of numerical solution of partial differential equations. Given that neural networks are very powerful universal function approximators [10, 26, 3, 36, 54], it is natural to consider the space of neural networks as an ansatz space for approximating solutions of PDEs. First proposed in [28] on an underlying collocation approach, it has been successfully used recently in different contexts in [40, 41, 35, 12, 13, 23] and references therein. Given spatial and temporal locations as inputs, the (deep) neural networks, proposed by these authors, approximate the solution of the underlying PDE, by outputting function values. This approach appears to work quite well for problems with high regularity (smoothness) of the underlying solutions (see [46, 23]) and/or if the solution of the underlying PDE possesses a representation formula in terms of integrals [13, 23]. Given their low regularity, it is unclear if solutions of (1.1), realized as functions of space and time, can be efficiently learned by deep neural networks [49].

Several papers applying deep learning techniques in the context of CFD advocate embedding deep learning modules within existing CFD codes to increase their efficiency [49, 42, 33]. In this paper, we adopt a different approach and propose to use deep neural networks to learn and predict the input parameters to observable map. Our algorithm will be based on fully connected networks (multi-layer preceptrons) which output values of the observable (1.2) for different input parameters $y\in Y\subset\mathbb{R}^{d}$ with $d>>1$. The network will be trained on data, generated from a few, say ${\mathcal{O}}(100)$, samples i.e, realizations of (1.2) with expensive CFD solvers.

However, the task of designing deep neural networks that will approximate the parameters to observable map with reasonable accuracy is far from straightforward on account of the following issues,

• •

Approximation results for deep neural networks [54, 5, 39], stipulate a network size of ${\mathcal{O}}\left(\varepsilon^{\frac{-d}{s}}\right)$ for attaining an error of ${\mathcal{O}}(\varepsilon)$, with $s$ denoting the Sobolev regularity of the underlying map. However, and as explained in section 2.3, the parameters to observable map for most CFD problems is atmost Lipschitz continuous. Thus, we might require unrealistically large networks for approximating the underlying maps in this context.

• •

Another relevant metric for estimating errors with deep neural network is the so-called generalization error of the trained network, (see section 2.3), that measures the errors of the network in predicting unseen data. The best available bounds for this error scale (inversely) with the square-root of the number of training samples. However, in this context, our training samples are generated from expensive CFD simulations. Hence, we might end up with large prediction errors in this data poor regime.

Hence, it is quite challenging to find deep neural networks that can accurately approximate maps of low Sobolev regularity in a data poor regime. Moreover, there are several hyperparameters that need to be specified in this framework, for instance size of the networks, choice of which variant of the stochastic gradient algorithm that one uses, choice of loss functions and regularizations thereof etc. A priori, the prediction errors can be sensitive to these choices. We propose an ensemble training procedure to search the hyperparameter space systematically and identify network architectures that are efficient in our context. By a combination of theoretical considerations and rigorous empirical experimentation, we provide a recipe for finding appropriate deep neural networks to compute the parameters to observable map i.e, those network architectures which ensure a low generalization error, even for relatively few training samples.

A second aim of this paper is to employ these trained neural networks, in conjunction with (Quasi-)Monte Carlo methods to compute statistical quantities of interest, particularly the underlying probability distributions of the observable and to demonstrate that these deep learning (quasi)-Monte Carlo methods significantly outperform the baseline algorithms.

The rest of the paper is organized as follows: the deep learning algorithm is presented in section 2. The deep learning (Quasi-)Monte Carlo algorithm for forward UQ is presented in section 3 and some details for the implementation of both sets of algorithms are provided in section 4. In section 5, we present numerical experiments illustrating the proposed algorithms and the results of the paper are summarized and discussed in section 6.

2 Deep learning algorithm

2.1 The problem

We consider the following very general form of a parameterized convection-diffusion equation,

[TABLE]

Here, $Y$ is the parameter space and without loss of generality, we assume it to be $Y=[0,1]^{d}$, for some $d\in{\mathbb{N}}$.

The spatial domain is labeled as $y\rightarrow D(y)\subset\mathbb{R}^{d_{s}}$ and ${\bf U}:[0,T]\times D\times Y\rightarrow\mathbb{R}^{m}$ is the vector of unknowns. The flux vector is denoted as ${\bf F}=(F_{i})_{1\leqslant i\leqslant d_{s}}:\mathbb{R}^{m}\times Y\rightarrow\mathbb{R}^{m}$ and ${\bf D}=(D_{ij})_{1\leqslant i,j\leqslant d_{s}}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ is the diffusion matrix.

The operator $L_{b}$ is a boundary operator that imposes boundary conditions for the PDE, for instance the no-slip boundary condition for incompressible Navier-Stokes equations or characteristic boundary conditions for hyperbolic systems of conservation laws. Additional conditions such as hyperbolicity for the flux function ${\bf F}$ and positive-definiteness of the diffusion matrix ${\bf D}$ might also be imposed, depending on the specific problem that is being considered.

We remark that the parameterized PDE (2.1) will arise in the context of (both forward and inverse) UQ, when the underlying convection-diffusion equation (1.1) contains uncertainties in the domain, the flux and diffusion coefficients and in the initial and boundary data. This input uncertainty propagates into the solution. Following [4] and references therein, it is customary to model such random inputs and the resulting solution uncertainties by random fields. Consequently, one can parameterize the probability space on which the random field is defined in terms of a parameter space $Y\subset\mathbb{R}^{d}$, for instance by expressing random fields in terms of (truncated) Karhunen-Loeve expansions. By normalizing the resulting random variables, one may assume $Y=[0,1]^{d}$, with possibly a large value of the parameter dimension $d$. Moreover, there exists a measure, $\mu\in{\rm Prob}(Y)$, with respect to which the data from the underlying parameter space is drawn.

We point out that the above framework is also relevant in problems of optimal control, design and PDE constrained optimization. In these problems, the parameter space $Y$ represents the set of design or control parameters.

For the parameterized PDE (2.1), we aim to compute observables of the following general form,

[TABLE]

Here, $\psi\in L^{1}_{{\rm loc}}(D_{y}\times(0,T))$ is a test function and $g\in C^{s}(\mathbb{R}^{m})$, for $s\geqslant 1$. Most interesting observables encountered in experiments, such as the lift and the drag, can be cast in this general form.

For fixed functions $\psi,g$, we also define the parameters to observable map:

[TABLE]

with $L_{g}$ being defined by (2.2).

We also assume that there exist suitable numerical schemes for approximating the convection-diffusion equation (2.1) for every parameter vector $y\in Y$. These schemes could be of the finite difference, finite volume, DG or spectral type, depending on the problem and on the baseline CFD code. Hence for any mesh parameter (grid size, time step) $\Delta$, we are assuming that for any parameter vector $y\in Y$, a high-resolution approximate solution ${\bf U}^{\Delta}(y)\approx{\bf U}(y)$ is available. Hence, there exists an approximation to the input to observable map $\mathcal{L}$ of the form,

[TABLE]

with the integrals in (2.2) being approximated to high accuracy by quadratures. Therefore, the original input parameters to observable map $\mathcal{L}$ is approximated by $\mathcal{L}^{\Delta}$ to very high accuracy i.e, for every value of a tolerance $\varepsilon>0$, there exists a $\Delta<<1$, such that

[TABLE]

for some $1\leqslant p\leqslant\infty$ and weighted norm,

[TABLE]

for $1\leqslant p<\infty$. The $L^{\infty}_{\mu}$ norm is analogously defined.

2.2 Deep learning the parameters to observable map

As stated earlier, it can be very expensive to compute the map $\mathcal{L}^{\Delta}(y)$ for each single realization, $y\in Y$, as a high-resolution CFD solver, possibly entailing a very large number of degrees of freedom, needs to be used. We propose instead, to learn this map by deep neural networks. This process entails the following steps,

2.2.1 Training set.

As is customary in supervised learning [20] and references therein, we need to generate or obtain data to train the network. To this end, we select a set of parameters $\mathcal{S}=\{y_{i}\}_{1\leqslant i\leqslant N}$, with each $y_{i}\in Y$. The points in $\mathcal{S}$ can be chosen as,

• (i.)

randomly from $Y$, independently and identically distributed with the underlying probability distribution $\mu$.

• (ii.)

from a suitable set of quadrature points in $Y$, for instance the so-called low discrepancy sequences that arise in Quasi-Monte Carlo (QMC) quadrature algorithms [8]. Examples of such sequences include Sobol or Halton QMC quadrature points [8].

We emphasize that the generation of QMC points is very cheap, particularly for Sobol or Halton sequences. Moreover, these points are better spread out over the parameter space than a random selection of points and might provide more detailed information about it [8]. Hence, a priori QMC points might be a better choice for sampling the data. One can also replace QMC points with other hierarchical algorithms such as nodes of sparse grids (Smolyak quadrature points) [7] to form the set $\mathcal{S}$.

Once the training parameter set $\mathcal{S}$ is chosen, we perform a set of high-resolution CFD simulations to obtain $\mathcal{L}^{\Delta}(y)$, for all $y\in\mathcal{S}$. As each high-resolution CFD simulation could be very expensive, we will require that $N=\#(\mathcal{S})$ will not be very large. It will typically be of at most ${\mathcal{O}}(100)$.

2.2.2 Neural network.

Given an input vector, a feedforward neural network (also termed as a multi-layer perceptron), shown in figure 1, consists of layer of units (neurons) which compose of either affine-linear maps between units (in successive layers) or scalar non-linear activation functions within units, culminating in an output [20]. In our framework, for any input vector $z\in Y$, we represent an artificial neural network as,

[TABLE]

Here, $\odot$ refers to the composition of functions and $\sigma$ is a scalar (non-linear) activation function. A large variety of activation functions have been considered in the machine learning literature [20]. A very popular choice is the ReLU function,

[TABLE]

When, $z\in\mathbb{R}^{p}$ for some $p>1$, then the output of the ReLU function in (2.7) is evaluated componentwise.

Moreover, for any $1\leqslant k\leqslant K$, we define

[TABLE]

For consistency of notation, we set $d_{1}=d$ and $d_{K}=1$.

Thus in the terminology of machine learning (see also figure 1), our neural network (2.6) consists of an input layer, an output layer and $(K-1)$ hidden layers for some $1<K\in\mathbb{N}$. The $k$-th hidden layer (with $d_{k}$ neurons) is given an input vector $z_{k}\in\mathbb{R}^{d_{k}}$ and transforms it first by an affine linear map $C_{k}$ (2.8) and then by a ReLU (or another) nonlinear (component wise) activation $\sigma$ (2.7). Although the neural network consists of composition of very elementary functions, its complexity and ability to learn very general functions arises from the interactions between large number of hidden layers [20].

A straightforward addition shows that our network contains $\left(d+1+\sum\limits_{k=2}^{K-1}d_{k}\right)$ neurons.

We denote,

[TABLE]

to be the concatenated set of (tunable) weights for our network. It is straightforward to check that $\theta\in\Theta\subset\mathbb{R}^{M}$ with

[TABLE]

Thus, depending on the dimensions of the input parameter vector and the number (depth) and size (width) of the hidden layers, our proposed neural network can contain a large number of weights. Moreover, the neural network explicitly depends on the choice of the weight vector $\theta\in\Theta$, justifying the notation in (2.6).

Although a variety of network architectures, such as convolutional neural networks or recurrent neural networks, have been proposed in the machine learning literature, [20] and references therein, we will restrict ourselves to fully connected architectures i.e, we do not a priori assume any sparsity structure for our set $\Theta$.

2.2.3 Loss functions and optimization.

For any $y\in\mathcal{S}$, we have already evaluated $\mathcal{L}^{\Delta}(y)$ from the high-resolution CFD simulation. One can readily compute the output of the neural network $\mathcal{L}_{\theta}(y)$ for any weight vector $\theta\in\Theta$. We define the so-called loss function or mismatch function, as

[TABLE]

for some $1\leqslant p<\infty$.

The goal of the training process in machine learning is to find the weight vector $\theta\in\Theta$, for which the loss function (2.11) is minimized. The resulting optimization (minimization) problem might lead to searching a minimum of a non-convex loss function. So, it is not uncommon in machine learning [20] to regularize the minimization problem i.e we seek to find,

[TABLE]

Here, $\mathcal{R}:\Theta\mapsto\mathbb{R}$ is a regularization (penalization) term. A popular choice is to set $\mathbb{R}(\theta)=\|\theta_{W}\|^{q}_{q}$ for either $q=1$ (to induce sparsity) or $q=2$. The parameter $0<\lambda<<1$ balances the regularization term with actual loss $J$ (2.11).

The above minimization problem amounts to finding a minimum of a possibly non-convex function over a subset of $\mathbb{R}^{M}$ for very large $M$. We can approximate the solutions to this minimization problem iteratively, either by a full batch gradient descent algorithm or by a mini-batch stochastic gradient descent (SGD) algorithm. A variety of SGD algorithms have been proposed in the literature and are heavily used in machine learning, see [45] for a survey. A generic step in a (stochastic) gradient method is of the form:

[TABLE]

with $\eta_{r}$ being the learning rate. The stochasticity arises in approximating the gradient in (2.13) by,

[TABLE]

and analogously for the gradient of the regularization term in (2.13). Here $\hat{\mathcal{S}}_{q}\subset\mathcal{S}$ refers to the $q$-th batch, with the batches being shuffled randomly. Moreover, the SGD methods are initialized with a starting value $\theta_{0}=\bar{\theta}\in\Theta$. A widely used variant of the SGD method is the so-called ADAM algorithm [27].

For notational simplicity, we denote the (approximate, local) minimum weight vector in (2.12) as $\theta^{\ast}$ and the underlying deep neural network $\mathcal{L}_{\theta^{\ast}}$ will be our neural network surrogate for the parameters to observable map $\mathcal{L}$ (2.4). The algorithm for computing this neural network is summarized below,

Algorithm 2.1**.**

Deep learning of parameters to observable map*.*

• Inputs:

Parameterized PDE (2.1), Observable (2.2), high-resolution numerical method for solving (2.1) and calculating (2.2).

• Goal:

Find neural network $\mathcal{L}_{\theta^{\ast}}$ for approximating the parameters to observable map $\mathcal{L}$ (2.4).

• Step $1$:

Choose the training set $\mathcal{S}$ and evaluate $\mathcal{L}^{\Delta}(y)$ for all $y\in\mathcal{S}$ by high-resolution CFD simulations.

• Step $2$:

For an initial value of the weight vector $\overline{\theta}\in\Theta$, evaluate the neural network $\mathcal{L}_{\overline{\theta}}$ (2.6), the loss function (2.12) and its gradients to initialize the (stochastic) gradient descent algorithm.

• Step $3$:

Run a stochastic gradient descent algorithm of form (2.13) till an approximate local minimum $\theta^{\ast}$ of (2.12) is reached. The map $\mathcal{L}^{\ast}=\mathcal{L}_{\theta^{\ast}}$ is the desired neural network approximating the parameters to observable map $\mathcal{L}$.

2.3 Theory

2.3.1 Approximation with deep neural networks.

Universal approximation theorems [3, 26, 10] prove that neural networks can approximate any continuous (in fact any measurable) function. However, these universality results are not constructive. More recent papers such as [54, 39] provide the following constructive approximation result:

Theorem 2.2**.**

Let $f\in W^{s,p}\left([0,1]^{d},\mathbb{R}\right)$ for some $1\leqslant p\leqslant\infty$, such that $\|f\|_{W^{s,p}}\leqslant 1$, then for every $\varepsilon>0$, there exists a neural network $NN(f)$ of the form (2.6) with the ReLU activation function, with ${\mathcal{O}}\left(1+\log\left(\frac{1}{\varepsilon}\right)\right)$ layers and of size (number of weights) ${\mathcal{O}}\left(\varepsilon^{-\frac{d}{s}}\left(1+\log\left(\frac{1}{\varepsilon}\right)\right)\right)$ such that

[TABLE]

Given the above approximation result, it is essential to investigate regularity of the underlying parameters to observeble map $\mathcal{L}$ in order to estimate the size of neural networks, needed to approximate them.

2.3.2 Regularity of the parameters to observable map (2.3)

For simplicity, we assume that the underlying domain $D_{y}=D\subset B_{R}$ for a. e. $y\in Y$ and $B_{R}$ is a ball of radius $R$. Moreover, we assume that $\psi,g^{k}\in L^{\infty}$, with $g^{k}$ denoting the $k$-th derivative of the map $g$ in (2.2), for some $k\geqslant 1$. Then, it is straightforward to obtain the following upper bound,

[TABLE]

with the constant $C$ depending only on the initial data, $\Psi$ and $g$. Thus, regularity in the parameters to observable map is predicated on the space-time and parametric regularity of the underlying solution field ${\bf U}$.

Unfortunately, it is well-known that one cannot expect such space-time or parametric regularity for solutions of generic convection-diffusion equations (1.1). In particular, even for the simple case of scalar one-dimensional conservation laws, at best we can show that the solution ${\bf U}\in L^{\infty}((0,T);BV(D\times Y))$. Therefore, one can show readily that

[TABLE]

It can be rightly argued that deriving regularity estimates on $\mathcal{L}$ in terms of the underlying solution field ${\bf U}$ can lead to rather pessimistic results and might overestimate possible cancellations. This is indeed the case, at least for scalar conservation laws with random initial data i.e, (2.1) with $m=1,\nu=0$ and ${\bf F}(y,{\bf U})={\bf F}({\bf U})$. In this case, we have the following theorem,

Theorem 2.3**.**

Consider the scalar conservation law, i.e parameterized convection-diffusion PDE (2.1) with $m=1$ and $\nu=0$, in domain $D=\mathbb{R}^{d_{s}}$ and time period $[0,T]$ and assume that the initial data $u_{0}\in W^{1,\infty}\left(Y;L^{1}(D)\right)$ and support of $u_{0}$ is compact. Assume, furthermore that $\psi\in L^{\infty}(D\times(0,T))$ and $g\in W^{1,\infty}(\mathbb{R})$. Then, the parameters to observable map $\mathcal{L}$, defined in (2.3) satisfies,

[TABLE]

for some constant $C$, depending on the initial data $u_{0}$, $\psi$ and $g$. Moreover, if the numerical solution $u^{\Delta}$ is generated by a monotone numerical method, then the approximate parameters to observable map $\mathcal{L}^{\Delta}$ (2.4) satisfies

[TABLE]

The proof of this theorem is a consequence of the $L^{1}$ stability (contractivity) of the solution operator for scalar conservation laws [11]. In fact, for any $y,y^{\ast}\in Y$, a straightforward calculation using the definition (2.3) and the Lipschitz regularity of $g$ yields,

[TABLE]

The above proof also makes it clear that bounds such as (2.18) and (2.19) will hold for systems of conservation laws and the incompressible Navier-Stokes equations as long as there is some stability of the field with respect to the input parameter vector. On the other hand, we cannot expect any such bounds on the higher parametric derivatives of $\mathcal{L}$, due to the lack of differentiability of the underlying solution field.

Given a bound such as (2.18) or (2.19), we can illustrate the difficulty of approximating the map $\mathcal{L}$ by considering a prototypical problem, namely that of learning the lift and the drag of the RAE2822 airfoil (see section 5.1). For this problem, the underlying parameter space is six-dimensional. Assuming that $s=1$ in theorem 2.2 and requiring that the approximation error is at most one percent relative error i.e $\varepsilon=10^{-2}$ in (2.15), yields a neural network of size ${\mathcal{O}}(10^{12})$ tunable parameters and at least $6$ layers. Such a large network is clearly unreasonable as it will be very difficult to train and expensive to evaluate.

2.3.3 Trained networks and generalization error

It can be argued that approximation theoretic results can severely overestimate errors with trained networks. Rather, the relevant measure is the so-called *generalization error *i.e, once a (local) minimum $\theta^{\ast}$ of (2.12) has been computed, we are interested in the following error,

[TABLE]

In practice, one estimates the generalization error on a so-called test set i.e, $\mathcal{T}\subset Y$, with $\#(\mathcal{T})>>N=\#(\mathcal{S})$. For instance, $\mathcal{T}$ could consist of i.i.d random points in $Y$, drawn from the underlying distribution $\mu$. In this case, the generalization error (2.20) can be estimated by the prediction error,

[TABLE]

It turns out that one can estimate the generalization error (2.20) by using tools from machine learning theory [47] as,

[TABLE]

with $N=\#(\mathcal{S})$ being the number of training samples. The numerator $U$ is typically estimated in terms of the Rademacher complexity or the Vapnik-Chervonenkis (VC) dimension of the network, see [47] and references therein for definitions and estimates.

Unfortunately, such bounds on $U$ are very pessimistic in practice and over estimate the generalization error by several (tens of) orders of magnitude [2, 55, 38]. Recent papers such as [2] provide far sharper bounds by estimating the compression i.e number of effective parameters to total parameters in a network.

Nevertheless, even if we make a stringent requirement of $U\sim{\mathcal{O}}(1)$, we might end up with unacceptably high prediction errors of $30-100\%$ for the ${\mathcal{O}}(10^{2}-10^{3})$ training samples that we can afford to generate in practice.

Summarizing, theoretical considerations outlined above indicate the challenges of finding deep neural networks to accurately learn maps of low regularity in a data poor regime. We illustrate these difficulties with a simple numerical example.

2.3.4 An illustrative numerical experiment.

In this experiment we consider the parameter space $Y=[0,1]^{6}$ the following two maps,

[TABLE]

Both maps are infinitely differentiable but with high amplitudes of the derivatives.

We define a neural network of the form (2.6) and with the architecture specified in table 3. In order to generate the training set, we select the first $128$ Sobol points in $[0,1]^{6}$ and sample the functions (2.23) at these points. The training is performed with ADAM algorithm and hyperparameters defined in table 4 (first line). The performance of the neural network is ascertained by choosing the first $8192$ Sobol points in $[0,1]^{6}$ as the test set $\mathcal{T}$. The results of evaluating the trained networks on the test set is shown in table 1 where we present the relative prediction error (2.21) as a percentage and also present the standard deviation (on the test set) of the prediction error. Clearly, the errors are unacceptably high, particularly for the unscaled sum of sines $\mathcal{L}^{1}$ in (2.23). Although still high, the errors reduce by an almost order of magnitude for the scaled sum of sines $\mathcal{L}^{2}$. This is to be expected as the derivatives for this map are smaller. This example illustrates the difficulty to approximating (even very regular) functions, in moderate to high dimensions, by neural networks when the training set is relatively small.

2.3.5 Pragmatic choice of network size

We estimate network size based on the following very deep but easy to prove fact about training neural networks [55],

Lemma 2.4**.**

For the map $\mathcal{L}^{\Delta}$ and for the training set $\mathcal{S}$, with $\#(\mathcal{S})=N$, there exists a weight vector $\hat{\theta}\in\Theta$, and a resulting neural network of the form $\mathcal{L}_{\theta}$, with $\Theta\subset\mathbb{R}^{M}$ and $M={\mathcal{O}}(d+N)$, such that the following holds,

[TABLE]

The above lemma implies that there exists a neural network with weight vector $\hat{\theta}$ of size ${\mathcal{O}}(d+N)$ such that the training error defined by,

[TABLE]

Thus, it is reasonable to expect that a network of size ${\mathcal{O}}(d+N)$ can be trained by a gradient descent method to achieve a very low training error. It is customary to monitor the generalization capacity of a trained neural network by computing a so-called validation set $\mathcal{V}\subset Y$ with $\mathcal{V}\cap\mathcal{S}=\emptyset$, and evaluating the so-called validation loss,

[TABLE]

A low validation loss is observed to correlate with low generalization errors. In order to generate the validation set, one can set aside a small proportion, say $10-20\%$ of the training set as the validation set.

2.4 Hyperparameters and Ensemble training.

Given the theoretical challenges of determining neural networks that can provide low prediction errors in our data poor regime, a suitable choice of the hyperparameters in the training process becomes imperative. To this end, we device a simple yet effective ensemble training algorithm. To this end, we consider the set of hyperparameters, listed in table 2. A few comments regarding the listed hyperparameters are in order. We need to choose the number of layers and width of each layer such that the overall network size, given by (2.10) is ${\mathcal{O}}(N+d)$. The exponents $p,q$ in the loss function and regularization terms (2.12) usually take the values of $1$ or $2$ and the regularization parameter $\lambda$ is required to be small. The starting value $\bar{\theta}$ is chosen randomly.

For each choice of the hyperparameters, we realize a single sample in the hyperparameter space. Then, the machine learning algorithm 2.1 is run with this sample and the resulting loss function minimized. Further details of this procedure are provided in section 4. We remark that this procedure has many similarities to the active learning procedure proposed recently in [56].

3 Uncertainty quantification in CFD with deep learning.

We will employ the trained deep neural networks that approximate the parameters to observable map $\mathcal{L}^{\Delta}$, to efficiently quantify uncertainty in the underlying map. We are interested in computing the entire measure (probability distribution) of the observable. To this end, we assume that the underlying probability distribution on the input parameters to the problem (2.1) is given by $\mu\in{\rm Prob}(Y)$. In the context of forward UQ, we are interested in how this initial measure is changed by the parameters to observable map $\mathcal{L}$ (or its high-resolution numerical surrogate $\mathcal{L}^{\Delta}$). Hence, we consider the push forward measure $\hat{\mu}^{\Delta}\in{\rm Prob}(\mathbb{R})$ given by

[TABLE]

for any $\mu$-measurable function $h:\mathbb{R}\mapsto\mathbb{R}$. It should be emphasized that any statistical moment of interest can be computed by integrating an appropriate test function with respect to this measure $\hat{\mu}^{\Delta}$. For instance, letting $h(w)=w$ in (3.1) yields the mean of the observable. Similarly, the variance can be computed from the mean and letting $h(w)=w^{2}$ in (3.1). The task of efficiently computing this probability distribution is significantly harder than just estimating the mean and variance of the underlying map.

The simplest algorithm for computing this measure $\hat{\mu}^{\Delta}$ is to use the Monte Carlo algorithm ([8] and references therein). It consists of selecting $J$ samples i.e points $\{y_{j}\}$ with each $y_{j}\in Y$ and $1\leqslant j\leqslant J$, that are independent and identically distributed (iid). Then, the Monte Carlo approximation of the push-forward measure is given by,

[TABLE]

It can be shown that $\hat{\mu}_{mc}\approx\hat{\mu}^{\Delta}$ with an error estimate that scales (inversely) as the square root of the number of samples $J$. Hence, the Monte Carlo algorithm can be very expensive in practice.

3.1 Quasi-Monte Carlo (QMC) methods

3.1.1 Baseline QMC algorithm.

Quasi-Monte Carlo (QMC) methods are deterministic quadrature rules for computing integrals, [8] and references therein. The key idea underlying QMC is to choose a (deterministic) set of points on the domain of integration, designed to achieve some measure of equidistribution. Thus, the QMC algorithm is expected to approximate integrals with higher accuracy than MC methods, whose nodes are randomly chosen [8].

For simplicity, we set the underlying measure $\mu$ to be the scaled Lebesgue measure on $Y$. Our aim is to approximate the push forward measure $\hat{\mu}^{\Delta}$ with the QMC algorithm. To this end, we choose a set of points $\mathcal{J}_{q}=\{y_{j}\}\subset Y$ with $1\leqslant j\leqslant J_{q}=\#(\mathcal{J}_{q})$. Then the measure $\hat{\mu}^{\Delta}$ can be approximated by,

[TABLE]

for any measurable function $h:Y\to\mathbb{R}$.

We want to estimate the difference between the measures $\hat{\mu}^{\Delta}$ and $\hat{\mu}_{qmc}$ and need some metric on the space of probability measures, to do so. One such metric is the so-called Wasserstein metric [53]. To define this metric, we need the following,

Definition 3.1**.**

Given two measures, $\nu,\sigma\in{\rm Prob}(\mathbb{R})$, a transport plan $\pi\in{\rm Prob}(\mathbb{R}^{2})$ is a probability measure on the product space such that the following holds for all measurable functions $F,G$,

[TABLE]

The set of transport plans is denoted by $\Pi(\nu,\sigma)$

Then the $1$-Wasserstein distance [53] between $\hat{\mu}^{\Delta}$ and $\hat{\mu}_{qmc}$ is defined as,

[TABLE]

It can be proved using the Koksma-Hlawka inequality and some elementary optimal transport techniques that the error in the Wasserstein metric behaves as,

[TABLE]

Here, $V^{\Delta}$ is some measure of variation of $\mathcal{L}^{\Delta}$. It can be bounded above by

[TABLE]

the Hardy-Krause variation of the function $\mathcal{L}^{\Delta}$. This upper bound requires some regularity for the underlying map in terms of mixed partial derivatives. However, it is well-known that the bound (3.7) is a large overestimate of the integration error [8].

The $\mathcal{D}(\mathcal{J}_{q})$ in (3.6) is the so-called discrepancy of the point sequence $\mathcal{J}_{q}$, defined formally in [8] and references therein. It measures how equally is the point sequence $\mathcal{J}_{q}$ spread out in $Y$. The whole objective in the development of quasi-Monte Carlo (QMC) methods is to design low discrepancy sequences. In particular, there exists many popular QMC sequences such as Sobol, Halton or Niederreiter [8] that satisfy the following bound on discrepancy,

[TABLE]

Thus, the integration error with QMC quadrature based on these rules behaves log-linearly with respect to the number of quadrature points. For simplicity, we assume that the dimension is moderate enough to replace the logarithms in (3.6) by a suitable power leading to

[TABLE]

for some $1/2<\alpha\leqslant 1$. Thus, for moderate dimensions, we see that the QMC algorithm outperforms the baseline MC algorithm.

Requiring that the Wasserstein distance is of size ${\mathcal{O}}(\varepsilon)$ for some tolerance $\varepsilon>0$, entails choosing $J_{q}\sim\frac{V^{\Delta}}{\varepsilon^{\frac{1}{\alpha}}}$ and leads to a cost of

[TABLE]

with $\mathcal{C}$ being the computational cost of a single CFD forward solve. This can be very expensive, particularly for smaller values of $\alpha$, as the cost of each forward solve is very high.

3.1.2 Deep learning quasi-Monte Carlo (DLQMC)

In order to accelerate the baseline QMC method, we describe an algorithm that combines the QMC method with a deep neural network approximation of the underlying parameters to observable map.

Algorithm 3.2**.**

Deep learning Quasi-Monte Carlo (DLQMC).**

• Inputs:

Parameterized PDE (2.1), Observables (2.2), high-resolution numerical method for solving (2.1) and calculating (2.2).

• Goal:

Approximate the push-forward measure $\hat{\mu}^{\Delta}$.

• Step 1:

For $J_{L}\in\mathbb{N}$, select $\mathcal{J}_{L}=\{y_{j}\}\subset Y$, with $1\leqslant j\leqslant J_{L}$, as the first $J_{L}$ points from a Quasi-Monte Carlo low-discrepancy sequence such as Sobol, Halton etc.

• Step 2:

For some $N<<J_{L}$, denote $\mathcal{S}=\{y_{j}\}\subset\mathcal{J}_{L}$, with $1\leqslant j\leqslant N$ i.e, first $N$ points in $\mathcal{J}_{L}$, and evaluate the map $h(\mathcal{L}^{\Delta}(y_{j}))$, for all $y_{j}\in\mathcal{J}_{d}$ with the high-resolution CFD simulation. As a matter of fact, any consecutive set of $N$ QMC points can serve as the training set.

• Step 3:

With $\mathcal{S}$ as the training set, compute an optimal neural network $\mathcal{L}^{\ast}$ for the parameters to observable map $\mathcal{L}^{\Delta}$, by using algorithm 2.1.

• Step 4:

Approximate the measure $\hat{\mu}^{\Delta}$ by

[TABLE]

The complexity of the DLQMC algorithm is analyzed in the following theorem.

Theorem 3.3**.**

For any given tolerance $\varepsilon>0$ and under the assumption that the generalization error (2.20) is of size ${\mathcal{O}}(\varepsilon)$ and the baseline QMC estimator (3.3) follows the error estimate (3.9), the speedup $\Sigma_{dlqmc}$, defined as the ratio of the cost of baseline QMC algorithm and the cost of the DLQMC algorithm 3.2 for approximating the push-forward measure $\hat{\mu}^{\Delta}$ to an error of size ${\mathcal{O}}(\varepsilon)$ in the $1$-Wasserstein metric, with the DLQMC algorithm 3.2, satisfies

[TABLE]

Here, $\mathcal{C},\mathcal{C}_{\ast}$ are the computational cost of computing $\mathcal{L}^{\Delta}(y)$ (high resolution CFD simulation) and $\mathcal{L}^{\ast}(y)$ (deep neural network) for any $y\in Y$ and $V^{\ast}=V(\mathcal{L}^{\ast})$ is an estimate on the variation that arises in a QMC estimate such as (3.9).

Proof.

We define the measure $\hat{\mu}^{\ast}=\mathcal{L}^{\ast}\#\mu\in{\rm Prob}(\mathbb{R})$ and estimate,

[TABLE]

We claim that assuming that the generalization error $\mathcal{E}_{G}\sim\varepsilon$ implies that $W_{1}\left(\hat{\mu}^{\Delta},\hat{\mu}^{\ast}\right)\sim\varepsilon$. To see this, we define a transport plan $\pi^{\ast}\in{\rm Prob}(\mathbb{R}^{2})$ by $\pi^{\ast}=\hat{\mu}^{\Delta}\otimes\hat{\mu}^{\ast}$. Then,

[TABLE]

This provides an upper bound on the Wasserstein distance (3.5) and proves the claim. Similarly, we use an estimate of the type (3.9) to obtain $W_{1}\left(\hat{\mu}^{\ast},\hat{\mu}^{\ast}_{qmc}\right)\sim\left(\frac{V^{\ast}}{J_{L}}\right)^{\alpha}$. Requiring this error to be ${\mathcal{O}}(\varepsilon)$ yields $J_{L}\sim\frac{V^{\ast}}{\varepsilon^{\frac{1}{\alpha}}}$ and leads to the following estimate on the total cost of the DLQMC algorithm,

[TABLE]

with $\mathcal{C},\mathcal{C}_{\ast}$ being the computational cost of a single CFD forward solve and evaluation of the deep neural network $\mathcal{L}^{\ast}$, respectively.

Dividing (3.13) with (3.10) and using $J_{q}\sim\frac{V^{\Delta}}{\varepsilon^{\frac{1}{\alpha}}}$ leads to the speedup estimate (3.12).

∎

Remark 3.4**.**

In particular, we know that the cost $\mathcal{C}={\mathcal{O}}\left(\Delta^{-\frac{1}{d+1}}\right)$ is very high for a small mesh size $\Delta$. On the other hand the cost $\mathcal{C}_{\ast}={\mathcal{O}}(M)$ with $M$ being the number of neurons (2.10) in the network. This cost is expected to be much lower than $\mathcal{C}$. Moreover, it is reasonable to expect that $V^{\ast}\approx V^{\Delta}$. As long as we can guarantee that $N<<J_{q}$, it follows from (3.12) that the speedup $\Sigma_{dlqmc}$, which measures the gain of using the deep learning based DLQMC algorithm 3.2 over the baseline Quasi-Monte Carlo method, can be quite substantial.

4 Implementation of ensemble training

In the following numerical experiments, we use fully connected neural networks with ReLU activation functions. Essential hyperparameters are listed in table 2. For each experiment, we use a reference network architecture i.e number of hidden layers and width per layer, such that the total network size is ${\mathcal{O}}(d+N)$. This ensures consistency with the prescriptions of lemma 2.4.

We use two choices of the exponent $p$ in the loss function (2.11) i.e either $p=1$ or $p=2$, denoted as mean absolute error (MAE) or mean square error (MSE), respectively. Similarly the exponent $q$ in the regularization term (2.12) is either $q=1$ or $q=2$. In order to keep the ensemble size tractable, we chose four different values of the regularization parameter $\lambda$ in (2.12) i.e, $\lambda=7.8\times 10^{-5},7.8\times 10^{-6},7.8\times 10^{-7},0$, corresponding to different orders of magnitude for the regularization parameter. For the optimizer of the loss function, we use either a standard SGD algorithm or ADAM. Both algorithms are used in the full batch mode. This is reasonable as the number of training samples are rather low in our case.

A key hyperparameter is the starting value $\bar{\theta}$ in the SGD or ADAM algorithm (2.13). We remark that the loss function (2.12) is non-convex and we can only ensure that the gradient descent algorithms converge to a local minimum. Therefore, the choice of the initial guess is quite essential as different initial starting values might lead to convergence to local minima with very different loss profiles and generalization properties. Hence, it is customary in machine learning to retrain i.e, start with not one, but several starting values and run the SGD or ADAM algorithm in parallel. Then, one has to select one of the resulting trained networks. The ideal choice would be to select the network with the best generalization error. However, this quantity is not accessible with the data at hand. Hence, we rely on the following surrogates for predicting the best generalization error,

Best trained network. (train): We choose the network that leads to the lowest training error (2.25), once training has been terminated.

• 2.

Best validated network. (val): We choose the network that leads to the lowest validation loss (2.26), once training has been terminated. One disadvantage of this approach is to sacrifice some training data for the validation set.

• 3.

Best trained network wrt mean. (mean-train): We choose the network that has minimized the error in the mean of the training set i.e,

[TABLE]

• 4.

Best trained network wrt Wasserstein. (wass-train): We choose the network that has minimized the error in the Wasserstein metric with respect to the training set $\mathcal{S}$:

[TABLE]

We note that Wasserstein metric for measures on $\mathbb{R}$ can be very efficiently computed with the Hungarian algorithm [37].

Another hyperparameter (property) is whether we choose randomly distributed points (Monte Carlo) or more uniformly distributed points, such as quasi-Monte Carlo (QMC) quadrature points, to select the training set $\mathcal{S}$. We use both sets of points in the numerical experiments below in order to ascertain which choice is to be preferred in practice. Other important hyperparameters in machine learning are the learning rate in (2.13) and the number of training epochs. Unless otherwise stated, we keep them fixed to the values $\eta_{r}==0.01,\forall r$, for both the SGD and ADAM algorithms and to $500000$ epochs, respectively.

4.1 Code

All the following numerical experiments use finite volume schemes as the underlying CFD solver. The machine learning and ensemble training runs are performed by a collection of Python scripts and Jupyter notebooks, utilizing Keras [57] and Tensorflow [58] for machine learning and deep neural networks. The ensemble runs over hyperparameters are parallelized as standalone processes, where each process trains the network for one choice of hyperparameters. Jupyter notebooks for all the numerical experiments can be downloaded from https://github.com/kjetil-lye/learning_airfoils, under the MIT license. A lot of care is taken to ensure reproducability of the numerical results. All plots are labelled with the git commit SHA code, that produced the plot, in the upper left corner of the plot.

5 Numerical results

In the following numerical experiments, we consider the compressible Euler equations of gas dynamics as the underlying PDE (1.1). In two dimensions, these equations are

[TABLE]

where $\rho,\mathbf{u}=(u,v)$ and $p$ denote the fluid density, velocity and pressure, respectively. The quantity $E$ represents the total energy per unit volume

[TABLE]

where $\gamma=c_{p}/c_{v}$ is the ratio of specific heats, chosen as $\gamma=1.4$ for our simulations. Additional important variables associated with the flow include the speed of sound $a=\sqrt{\gamma p/\rho}$, the Mach number $M=|\mathbf{u}|/a$ and the fluid temperature $T$ evaluated using the ideal gas law $p=\rho RT$, where $R$ is the ideal gas constant.

5.1 Flow past an airfoil

5.1.1 Problem description.

We consider flow past an RAE 2822 airfoil, whose shape is shown in Figure 2 (Left). The flow past this airfoil is a benchmark problem for UQ in aerodynamics [25]. The upper and lower surface of the airfoil are denoted at $(x,S^{U}(x))$ and $(x,S^{L}(x))$, with $x\in[0,1]$. We consider the following parametrized (free-stream) initial conditions and perturbed airfoil shapes

[TABLE]

where $\alpha$ is the angle of attack, the gas constant $R=1$, $\varepsilon_{1}=0.1$, $\varepsilon_{2}=0.2$, $y\in Y=[0,1]^{6}$ is a parameter and $G_{k}(y)=2y_{k}-1$ for $k=1,...,6$.

The total drag and lift experienced on the airfoil surface $\hat{S}(y)=\hat{S}^{L}(y)\bigcup\hat{S}^{U}(y)$ are chosen as the observables for this experiment. These are expressed as

[TABLE]

where $\mathcal{K}(y)=\rho^{\infty}(y)|\mathbf{u}^{\infty}(y)|^{2}/2$ is the free-stream kinetic energy, while $\psi_{l}(y)=[-\sin(\alpha(y)),\ \cos(\alpha(y))]$ and $\psi_{d}(y)=[\cos(\alpha(y)),\ \sin(\alpha(y))]$.

The domain is discretized using triangles to form the primary mesh, while the control volumes are the voronoi dual cells formed by joining the circumcenters of the triangles to the face mid-points. A zoomed view of the primary and dual meshes are shown in Figure 2. The solutions are approximated using a second-order finite volume scheme, with the mid-point rule used as the quadrature to approximate the cell-boundary integrals. The numerical flux is chosen as the kinetic energy preserving entropy stable scheme [43], with the solutions reconstructed at the cell-interfaces using a linear reconstruction with the minmod limiter . For details on flux expression and the reconstruction procedure, we refer interested readers to [43]. The scheme is integrated using the SSP-RK3 time integrator. The contour plots for the Mach number with two different realizations of the parameter $y\in[0,1]^{6}$ are shown in Figure 3. It should be emphasized that the stochastic perturbations in (5.2) are rather strong, representing on an average a $10-20$ percent standard deviation (over mean). Hence, there is a significant statistical spread in the resulting solutions. As shown in figure 3, one of the configurations lead to a transsonic shock around the airfoil while another leads to a shock-free flow around the airfoil.

5.1.2 Generation of training data.

In order to generate the training, validation and test sets, we denote $\mathcal{Q}^{sob}$ as the first $1001$ Sobol points on the domain $[0,1]^{6}$. For each point $y_{j}\in\mathcal{Q}^{sob}$, the maps or (rather their numerical surrogate) $\mathcal{L}_{1,2}^{\Delta}(y_{j})$ is then generated from high-resolution numerical approximations ${\bf U}^{\Delta}(y)$ and the corresponding lift and drag are calculated by a numerical quadrature . This forms the test set $\mathcal{T}$. For this numerical experiment, we take the first $128$ Sobol points to generate the training set $\mathcal{S}$ and the next $128$ points to generate the validation set $\mathcal{V}$.

5.1.3 Results of the Ensemble training procedure.

To set up the ensemble training procedure of the previous section, we select a fully connected network, with number of layers and size of each layer listed in table 3. Our network has $1149$ tunable parameters and this choice is clearly consistent with the prescriptions of Lemma 2.4. Varying the hyperparameters as described in section 4 led to a total of $114$ configurations and each was retrained $5$ times with different random initial choices of the weights and the biases. Thus, we trained a total of $570$ networks in parallel.

For each configuration, we select the retraining that minimizes the set selection criteria for the configuration. Once the network with best retraining is selected, it is evaluated on the whole test set $\mathcal{T}=\mathcal{Q}^{sob}$ and the percentage relative $L^{2}$ prediction error (2.21) computed. Histograms, describing the probability distribution over all hyperparameter samples for the lift and the drag are shown in figure 4. From this ensemble, we choose the hyperparameter configuration that led to the smallest prediction error and list them in table 4. As seen from this table, there is a difference in the best performing network corresponding to lift and to drag. However, both networks use ADAM as the optimizer and regularize with the $L^{2}$-norm of the weights in the loss function (2.12).

We plot the training (and validation) loss, for these best performing networks with respect to the number of epochs of the ADAM algorithm in figure 5 (Top) and observe a three orders (for lift) and two orders (for drag) reduction in the underlying loss function during training.

From figure 5 (Bottom) and table 4, we note that the prediction errors for the best performing networks are very small i.e, less than $1\%$ relative error for the lift and less than $2\%$ relative error for the drag, even for the relatively small number ($128$) of training samples. This is orders of magnitude less than the sum of sines regression problem (table 1). These low prediction errors are particularly impressive when viewed in term of the theory presented in section 2.3. A crude upper bound on the generalization error is given by $\sqrt{\frac{M}{N}}$, with $M,N$ being the number of parameters in the network and number of training samples respectively. Setting these numbers for this problem yields a relative error of approximately $300\%$. However, our prediction errors are $100-200$ times (two orders of magnitude) less, demonstrating significant compression in the trained networks [2]. These results are further contextualized, when considered with reference to the computational costs, shown in table 5. In this table, we present the costs of generating a single sample, the average training time (over a subset of hyperparameter configurations) and the time to evaluate a trained network. Each sample is generated on the EULER high-performance cluster of ETH-Zurich. In this case, the training time is a very small fraction of the time to generate a single sample, let alone the whole training set whereas the computational cost of evaluating the neural networks is almost negligible as it is nine orders of magnitude less than evaluating a single training sample.

A priori, it is unclear if the underlying problem has some special features that makes the compression, and hence high accuracy possible. In particular, learning an affine map is relatively straightforward for a neural network. To test this possibility, we approximate the underlying parameters to observable map, for both lift and drag, with classical linear least squares algorithm i.e, $\mathcal{L}^{\text{lsq}}_{j}:\mathbb{R}^{6}\to\mathbb{R}$ be affine map defined by

[TABLE]

The resulting convex optimization problem is solved by the standard gradient descent. The errors with the linear least squares algorithm are shown in figure 4. From this figure we can see that the prediction error with the linear least squares is approximately $4\%$ for lift and $36\%$ for drag. Hence, the best performing neural networks are approximately $5$ times more efficient for lift and $20$ times more efficient for drag, than the linear least squares algorithm. The gain over linear least squares implicitly measures how nonlinear the underlying map is. Clearly in this example, the drag is more nonlinear than the lift and thus more difficult to approximate.

5.1.4 Network sensitivity to hyperparameters.

We summarize the results of the sensitivity of network performance to hyperparameters below,

• •

Overall sensitivity. The overall sensitivity to hyperparameters is depicted as histograms, representing the distribution, over samples (hyperparameter configurations), shown in figure 4. As seen from the figure, the prediction errors with different hyperparameters are spread over a significant range for each observable. There are a few outliers which perform very poorly, particularly for the drag. On the other hand, a large number of hyperparameter configurations concentrate around the best performing network. Moreover all (most) of the hyperparameter choices led to better prediction errors for the lift (drag), than the linear least squares algorithm.

• •

Choice of optimizer. The difference between ADAM and the standard SGD algorithm is shown in figure 6. As in the Sod shock tube problem, ADAM is clearly superior to the standard SGD algorithm. Moreover, there are a few outliers with the ADAM algorithm that perform very poorly for the drag. We call these outliers as the bad set of ADAM and they correspond to the configurations with MSE loss function and $L^{1}$ regularization, or $L^{2}$ regularization, with a regularization parameter $\lambda>10^{-5}$, or MAE loss function without any regularization. The difference between the bad set of ADAM and its complement i.e, the good set of ADAM is depicted in figure 7.

• •

Choice of loss function. The difference between $L^{1}$ mean absolute error (MAE) and $L^{2}$ root mean square error (MSE) i.e, exponent $p=1$ or $p=2$ in the loss function (2.11) is plotted in the form of a histogram over the error distributions in figure 8. The difference in performance, with respect to choice of loss function is minor. However, there is an outlier, based on MSE, for the drag.

• •

Choice of type of regularization. The difference between choosing an $L^{1}$ or $L^{2}$ regularization term in (2.12) is shown in figure 9. Again, the differences are minor but the $L^{2}$ regularization seems to be slightly better, particularly for the drag.

• •

Value of regularization parameter. The variation of network performance with respect to the value of the regularization parameter $\lambda$ is shown in figure 11. We observe from this figure that a little amount of regularization works better than no regularization. However, the variation in prediction errors wrt respect to the non-zero values of $\lambda$ is minor.

• •

Choice of selection criteria. The distribution of prediction errors with respect to the choice of selection criteria for retraining, namely train,val,mean-train,wass-train is shown in figure 10. There are very minor differences between mean-train and wass-train. On the other hand, both these criteria are clearly superior to the other two selection criteria.

• •

Sensitivity to retrainings. The variation of performance with respect to retrainings i.e, different random initial starting values for ADAM, is provided in table 6. In this table, we list the minimum, maximum, mean and standard deviation of the relative prediction error, over $5$ retrainings, for the best performing networks for each observable. As seen from this table, the sensitivity to retraining is rather low for the lift. On the other hand, it is more significant for the drag.

• •

Variation of Network size. We consider the dependence of performance with respect to variation of network size by varying the width and depth for the best performing networks. The resulting prediction errors are shown in table 7 where a $3\times 3$ matrix for the errors (rows represent width and columns depth) is tabulated. All the network sizes are consistent with the requirements of lemma 2.4. As observed from the table, the sensitivity to network size is greater for the drag than for the lift. It appears that increasing the width for constant depth results in lower errors for the drag. On the other hand, training clearly failed for some of the larger networks as the number of training samples appears to be not large enough to train these larger networks.

5.1.5 Network performance with respect to number of training samples.

All the above results were obtained with $N=128$ training samples. We have repeated the entire ensemble training with different samples sizes i.e $N=32,64,256$ also and present a summary of the results in figure 12, where we plot the mean (percentage relative) $L^{2}$ prediction error with respect to the number of samples. For each sample size, three different errors are plotted, namely the minimum (maximum) prediction errors and the mean prediction error over a subset of the hyperparameter space, that omits the bad set of ADAM as the optimizer. Moreover, we only consider errors with selected (optimal) retrainings in each case. As seen from the figure, the prediction error decreases with sample size $N$ as predicted by the theory. The decay seems of the of form $\left(\frac{U_{i}}{N}\right)^{\alpha_{i}}$ with $\alpha_{1}=0.81$ and $U_{1}=0.77$ for the lift and $\alpha_{2}=0.92$ and $U_{2}=3.42$ for the drag. In both cases, the decay of error with respect to samples, is at a higher rate than that predicted by the theory in (2.22). Moreover, the constants are much lower than the number of parameters in the network. Both facts indicate a very high degree of compressibility (at least in this possible pre-asymptotic regime) for the trained networks and explains the unexpectedly low prediction errors.

5.1.6 UQ with deep learning

Next, we approximate the underlying probability distributions (measures) (3.1), with respect to the lift and the drag in the flow past airfoils problem with the DLQMC algorithm 3.2. A reference measure is computed from the whole test set of $1001$ samples and the corresponding histograms are shown in figure 13 to visualize the probability distributions. In the same figure, we also plot corresponding histograms, computed with the best performing networks for each observable (listed in table 4) and observe very good qualitative agreement between the outputs of the DLQMC algorithm and the reference measure.

In order to quantify the gain in computational efficiency with the DLQMC algorithm over the baseline QMC algorithm in approximating probability distributions, we will compute the speedup. To this end, we observe from figure 14 that the baseline QMC algorithm converges to the reference probability measure with respect to the Wasserstein distance (3.9), at a rate of $\alpha=0.81$, for both the lift and drag. Next, from table 5, we see that once the training samples have been generated, the cost of training the network and performing evaluations with the trained network is essentially free. Consequently, if the number of training samples is $N$, we compute a raw speed-up (for each observable) given by,

[TABLE]

Here for $i={\rm Lift},~{}{\rm Drag}$, $\hat{\mu}_{i}^{\Delta}$ is the reference measure (probability distribution) computed from the test set $\mathcal{T}$, $\hat{\mu}_{i,qmc,N}$ is the measure (3.3) with the baseline QMC algorithm and $N$ QMC points, and $\hat{\mu}^{\ast}_{i,qmc}$ is the measure computed with the DLQMC algorithm (3.11). The real speed up is then given by,

[TABLE]

in order to compensate for the rate of convergence of the baseline QMC algorithm. In other words, our real speed up compares the costs of computing errors in the Wasserstein metric, that are of the same magnitude, with the DLQMC algorithm vis a vis the baseline QMC algorithm. Henceforth, we denote $\sigma_{dlqmc}^{real}$ as the speedup and list the speedups (for each observable) with the DLQMC algorithm for $N=128$ training samples in table 8. We observe from this table that we obtain speedups ranging between half an order to an order of magnitude for the lift and the drag. The speedups for drag are slightly better than that for lift.

In figure 15, we plot the speedup of the DLQMC algorithm (over the baseline QMC algorithm) as a function of the number of training samples in the deep learning algorithm 2.1. As seen from the figure, the best speed ups are obtained in case of $64$ training samples for the lift and $128$ training samples for the drag. This is on account of a complex interaction between the prediction error which decreases (rather fast) with the number of samples and the fact that the errors with the baseline QMC algorithm also decay with an increase in the number of samples.

5.1.7 Comparison with Monte Carlo algorithms

We have computed a Monte Carlo approximation of the probability distributions of the lift and the drag by randomly selecting $N=128$ points from the parameter domain $Y=[0,1]^{6}$ and computing the probability measure $\hat{\mu}_{mc}=\frac{1}{N}\sum_{j=1}^{N}\delta_{\mathcal{L}^{\Delta}(y_{j})}$. We compute the Wasserstein error $W_{1}(\hat{\mu}_{mc},\hat{\mu}^{\Delta})$, with respect to the reference measure $\hat{\mu}^{\Delta}$ and divide it with the error obtained with the DLQMC algorithm to obtain a raw speedup of the DLQMC algorithm over MC. As MC errors converge as a square root of the number of samples, the real speedup of DLQMC over MC is calculated by squaring the raw speedup. We present this real speedup over MC in table 9 and observe from that the speedups of the DLQMC algorithm over the baseline MC algorithm are very high and amount to at least two orders of magnitude. This is not surprising as the baseline QMC algorithm significantly outperforms the MC algorithm for this problem. Given that the DLQMC algorithm was an order of magnitude faster than the baseline QMC algorithm, the cumulative effect leads to a two orders of magnitude gain over the baseline MC algorithm.

5.1.8 Choice of training set.

For the sake of comparision with our choice of Sobol points to constitute the training set $\mathcal{S}$ in this problem, we chose $N=128$ random points in $Y$ as the training set. With this training set, we repeated the ensemble training procedure and found the best performing networks. The relative percentage mean $L^{2}$-prediction error (2.21) (with respect to a test set of $320$ random points) was computed and is presented in table 10. As seen from this table, the prediction errors with respect to the best performing networks are considerably (an order of magnitude) higher than the corresponding errors obtained for the QMC training points (compare with table 4). Hence, at least for this problem, Sobol points are a much better choice for the training set than random points. An intuitive reason for this could lie in the fact that the Sobol points are better distributed in the parameter domain than random points. So on an average, a network trained on them will generalize better to unseen data, than it would if trained on random points.

5.2 A stochastic shock tube problem

As a second numerical example, we consider a shock tube problem with the one-dimensional version of the compressible Euler equations (5.1).

5.2.1 Problem description

The underlying computational domain is $[-5,5]$ and initial conditions are prescibed in terms of left state $(\rho_{L},\ u_{L},\ p_{L})$ and a right state $(\rho_{R},\ u_{R},\ p_{R})$ at the initial discontinuity $x_{0}$. These states are defined as,

[TABLE]

$y\in Y=[0,1]^{6}$ is a parameter and $G_{k}(y)=2y_{k}-1$ for $k=1,...,6$. We use the Lebesgue measure on $Y$ as the underlying measure.

The observables for this problem are chosen as the average integral of the density over fixed intervals

[TABLE]

where $I_{1}=[-1.5,-0.5],I_{2}=[0.8,1.8],I_{3}=[2,3]$ are the underlying regions of interest.

A priori, this problem appears to be simpler as it is only one-dimensional when compared to the two-dimensional airfoil. However, the difficulty stems from the large variance of the initial data. As seen from (5.7), the ratio of the standard deviation to the mean, for the underlying variables, is very high and can be as high as $100\%$. This is almost an order of magnitude larger than the corresponding initial stochasticity of the flow past airfoil problem and we expect this large initial variance to be propagated into the solution. This is indeed borne out by the results shown in figure 16, where we plot the mean (and standard deviation) of the density, both initially and at time $T_{f}=1.5$. Moreover, the initial stochastic inputs span a wide range of possible scenarios. Particular choices of the parameters i.e, $(y_{1},y_{2},y_{3},y_{5},y_{6})=(7/9,1/24,0,1/32,1/13,0)$ and $(y_{1},y_{2},y_{3},y_{5},y_{6})=(0.1611,2/3,0.698,0.8213,0.8015,0)$, correspond to the classical Sod shock tube and Lax shock tube problems [16]. Thus, we can think of this stochastic shock tube problem as more or less generic (scaled) Riemann problem for the Euler equations.

In a recent paper [32], the authors have analyzed the role of the underlying variance in the generalization errors for deep neural networks trained to regress functions. In particular, high variance can lead to large generalization errors. In this sense, this shock tube problem might be harder to predict with neural networks than the flow past airfoil.

5.2.2 Generation of training data.

In order to generate the training and test sets, we denote $\mathcal{Q}^{sob}$ as the first $16256$ Sobol points on the domain $[0,1]^{6}$. The training set $\mathcal{S}$ consists of the first $128$ of these Sobol points whereas as the test set is $\mathcal{T}=\mathcal{Q}^{sob}$. For each point $y_{j}\in\mathcal{Q}^{sob}$, the maps or (rather their numerical surrogate) $\mathcal{L}_{1,2,3}^{\Delta}(y_{j})$ are generated from high-resolution numerical approximations ${\bf U}^{\Delta}(y)$ obtained using a second-order finite volume scheme. The domain is discretized into a uniform mesh of $K=2048$ disjoint intervals, with the HLL numerical flux and with the left and right cell-interface values obtained using a second-order WENO reconstruction. Time marching is performed using the second-order Runge-Kutta method with a CFL number of 0.475. Once the finite volume solution is computed, the corresponding observables (5.8) at each $y$ are approximated as

[TABLE]

where $x_{i}$ is the barycenter of the cell $\Omega_{i}$.

5.2.3 Results.

We follow the ensemble training procedure outlined in section 4. The network architecture, specified in table 3 is used together with the $114$ hyperparameter configurations of the previous section, with $5$ retrainings for each hyperparameter configuration. The best performing networks for each observable are identified as the ones with lowest prediction error (2.21) and these networks are presented in table 11. From the table, we observe very slight differences for the best performing networks for the three observables. However, these networks are indeed different from the best performing networks for the flow past airfoil (compare with table 4). Moreover, the predictions errors are larger than the flow past airfoil case, ranging between $2\%$ for the observable $\mathcal{L}_{1}$ to approximately $6\%$ for the other two observables. This larger error is expected as the variance of the underlying maps is significantly higher. Nevertheless, a satisfactory prediction accuracy is obtained for only $128$ training samples.

Results of the ensemble training procedure are depicted in the form of histograms for the prediction errors in figure 17. From this figure, we observe that although there are some outliers, most of the hyperparameter configurations resulted in errors comparable to the best performing networks (see table 11). Furthermore, a large majority of the hyperparameters resulted in substantially (a factor of $3-4$) smaller prediction errors, when compared to the linear least squares regression (5.4).

A sensitivity study (not presented here) indicates very similar robustness results as for the flow past the airfoil. In figure 18, we plot the prediction error with respect to the number of samples and observe that the mean error (over hyperparameter configurations) decays with an exponent of $0.5$ for the observables $\mathcal{L}_{2,3}$. On the other hand, the decay is a bit slower for $\mathcal{L}_{1}$. However, the low value of the constants still indicates significant compression allowing us to approximate this rather intricate problem with accurate neural networks.

Finally, we compute the underlying probability distributions for each observable with the DLQMC algorithm 3.2. In figure 19, we plot the maximum and mean real speedup (5.6), over hyperparameter configurations, for each observable, with respect to the number of training samples. From this figure, we observe maximum speedups of about $16$ for the observables $\mathcal{L}_{1,2}$ and $5$ for the observable $\mathcal{L}_{3}$, indicating that the DLQMC algorithm is significantly superior to the baseline QMC algorithm, even for this problem.

6 Discussion

Problems in CFD such as UQ, Bayesian Inversion, optimal design etc are of the many query type i.e, solving them requires evaluating a large number of realizations of the underlying input parameters to observable map (2.3). Each realization involves a call to an expensive CFD solver. The cumulative total cost of the many-query simulation can be prohibitively expensive.

In this paper, we have tried to harness the power of machine learning by proposing using deep fully connected neural networks, of the form (2.6), to learn and predict the underlying parameters to observable map. However, a priori, this constitutes a formidable challenge. In section 2.3, we argue by a combination of theoretical considerations and numerical experiments that the crux of this challenge it to find neural networks to learn maps of low regularity in a data poor regime. We are in this regime as the many parameters to observable maps (2.3) in fluid dynamics, can be at best, Lipschitz continuous. Moreover, given the cost, we can only afford to compute a few training samples.

We overcome these difficulties with the following novel ideas, first, we chose to focus on learning observables rather than the full solution of the parametrized PDE (2.1). Observables are more regular than fields (see section 2.3) and might have less underlying variance and are thus, easier to learn, see [32] for a more recent theoretical justification. Second, we propose using low-discrepany sequences to generate the training set $\mathcal{S}$. The equi-distribution property of these points might ensure lower generalization errors than those resulting from training data generated at random points, see the forthcoming paper [34] for further details. Finally, we propose an ensemble training procedure (see section 4) to systematically scan the hyperparameter space in order to winnow down the best performing networks as well as to test sensitivity of the trained networks to hyperparameter choices.

We presented two prototypical numerical experiments for the compressible Euler equations in order to demonstrate the efficacy of our approach. We observed from the numerical experiments that the trained networks were indeed able to provide low prediction errors, despite being trained on very few samples. Moreover, a large number of hyperparameter configurations were close in performance to the best performing networks, indicating a significant amount of robustness to most hyperparameters. More crucially, the obtained generalization errors followed the theoretical predictions of section 2.3 and indicating a high degree of compression, explaining why the networks generalized well. It is essential to state that these low prediction errors are obtained even when the cost of evaluating the neural network is several orders of magnitude lower than the full CFD solve.

As a concrete application of our deep learning algorithm 2.1, we combined it with a Quasi-Monte Carlo (QMC) method to propose a deep learning QMC algorithm 3.2, for performing forward UQ, by approximating the underlying probability distribution (3.1) of the observable. In theorem 3.3, we prove that the DLQMC algorithm is guaranteed to out-perform the baseline QMC algorithm as long as the underlying neural network generalizes well. This is indeed borne out in the numerical experiments, where we observed about an order of magnitude speedup of the DLQMC algorithm over the corresponding QMC algorithm and two orders of magnitude speedup over a Monte Carlo algorithm.

Thus, we have demonstrated the viability of using neural networks to learn observables and to solve challenging problems of the many-query type. Our approach can be compared to a standard model order reduction algorithm [52]. Here, generating the training data and training the networks is the offline step whereas evaluating the network constitutes the online step. Our results compare favorably with attempts to use Model order reduction for hyperbolic problems [1, 9], particularly when it comes to the low costs of training and evaluation. On the other hand, typical MOR algorithms provide a surrogate for the solution field, where we only predict the observable of interest in the approach presented here.

Our results can be extended in many directions; we can adapt this setting to learn the full solution field, instead of the observables. However, this might be harder in practice as the field is even less regular than the observable. Our approach is entirely data driven, in the sense that we do not use any specific information about the underlying parameterized PDE (2.1). Thus, the approach can be readily extended to other forms of (2.1), for instance the Navier-Stokes equations or other elliptic and parabolic PDEs. Finally, we plan to apply the trained neural networks in the context of Bayesian inverse problems and shape optimization under uncertainty, in forthcoming papers.

Acknowledgements

The research of SM is partially support by ERC Consolidator grant (CoG) NN 770880 COMANFLO. A large proportion of computations for this paper were performed on the ETH compute cluster EULER.