Epistemic Modeling Uncertainty of Rapid Neural Network Ensembles for
  Adaptive Learning
Emulator embedded neural networks, which are a type of physics informed
neural network, leverage multi-fidelity data sources for efficient design
exploration of aerospace engineering systems. Multiple realizations of the
neural network models are trained with different random initializations. The
ensemble of model realizations is used to assess epistemic modeling uncertainty
caused due to lack of training samples. This uncertainty estimation is crucial
information for successful goal-oriented adaptive learning in an aerospace
system design exploration. However, the costs of training the ensemble models
often become prohibitive and pose a computational challenge, especially when
the models are not trained in parallel during adaptive learning. In this work,
a new type of emulator embedded neural network is presented using the rapid
neural network paradigm. Unlike the conventional neural network training that
optimizes the weights and biases of all the network layers by using
gradient-based backpropagation, rapid neural network training adjusts only the
last layer connection weights by applying a linear regression technique. It is
found that the proposed emulator embedded neural network trains
near-instantaneously, typically without loss of prediction accuracy. The
proposed method is demonstrated on multiple analytical examples, as well as an
aerospace flight parameter study of a generic hypersonic vehicle.
Atticus Beachy, Harok Bae, Jose Camberos, Ramana Grandhi
# Epistemic Modeling Uncertainty of Rapid Neural Network Ensembles for Adaptive Learning

###### Abstract

**Emulator embedded neural networks, which are a type of physics informed neural network, leverage multi-fidelity data sources for efficient design exploration of aerospace engineering systems. Multiple realizations of the neural network models are trained with different random initializations. The ensemble of model realizations is used to assess epistemic modeling uncertainty caused due to lack of training samples. This uncertainty estimation is crucial information for successful goal-oriented adaptive learning in an aerospace system design exploration. However, the costs of training the ensemble models often become prohibitive and pose a computational challenge, especially when the models are not trained in parallel during adaptive learning. In this work, a new type of emulator embedded neural network is presented using the rapid neural network paradigm. Unlike the conventional neural network training that optimizes the weights and biases of all the network layers by using gradient-based backpropagation, rapid neural network training adjusts only the last layer connection weights by applying a linear regression technique. It is found that the proposed emulator embedded neural network trains near-instantaneously, typically without loss of prediction accuracy. The proposed method is demonstrated on multiple analytical examples, as well as an aerospace flight parameter study of a generic hypersonic vehicle.**

_Keywords:_

Machine Learning, Neural Network, Multifidelity, Active Learning, Aircraft Design

## 1 Introduction

Design exploration for high-performance aircraft requires computational modeling of multiple unconventional design configurations. Models must capture aerodynamics, structure, and propulsion, as well as interactions between these disciplines. A common design exploration technique is to sample the expensive physics-based models in a design of experiments and then use the sample data to train an inexpensive meta-model. Conventional metamodels include regression models such as ridge regression [1], Lasso [2], Polynomial Chaos Expansion [3], Gaussian process regression (GPR) or kriging [4, 5, 6], and neural networks [7]. However, many simulation evaluations are needed for the design of experiments because of the large number of independent parameters for each design and the complex responses resulting from interactions across multiple disciplines. Because high-fidelity simulations are expensive, the total computational costs can easily become computationally intractable.

Computational cost reduction is often achieved using Multi-Fidelity Methods (MFM) and Active Learning (AL). MFMs work by supplementing High-Fidelity (HF) simulations with less accurate but inexpensive
Low-Fidelity (LF) simulations. AL involves intelligent generation of training data in an iterative process: rebuilding the metamodel after each HF sample is added, and then using the metamodel to determine the best HF sample to add next.

When performing Multi-Fidelity (MF) modeling [8, 9], the usual strategy is to generate many affordable Low-Fidelity (LF) samples to capture the design space and correct them using a small number of expensive High-Fidelity (HF) samples. MF modeling is a more cost-effective way of training an accurate surrogate model than using a single fidelity, which will suffer from sparseness with only HF samples and inaccuracy with only LF samples. One method for generating LF data is Model Order Reduction (MOR). Nachar et al. [10] used this technique to generate LF data for a multi-fidelity kriging model. Reduced order models [11, 12, 13] are constructed by approximating the low-dimensional manifold on which the solutions lie. The manifold can be approximated linearly using Proper Orthogonal Decomposition, or nonlinearly using various methods such as an autoencoder neural network [14]. This approximation enables vastly decreasing the model degrees of freedom, which in turn reduces the computational costs of a simulation. Typically, reduced order models are used to give final predictions instead of being used as LF models in a multi-fidelity surrogate. However, this limits the aggressiveness with which the model can be reduced without introducing unacceptable errors into the final predictions. Even after very aggressive reduction, reduced order models can still provide valuable information about trends when used as LF models. LF models can also be constructed by coarsening the mesh, simplifying the physics, or utilizing historical data from a similar problem.

When performing AL, location-specific epistemic uncertainty information is critical for determining where to add additional samples. Kriging is popular in large part because it returns this uncertainty. Cokriging [15, 16] is a popular type of MFM which extends kriging to multiple fidelities. It typically performs well but has difficulties if the LF function is not well correlated with the HF function. It also cannot incorporate more than a single LF function unless they fall into a strict hierarchy known beforehand. As an alternative to cokriging, a localized-Galerkin kriging approach has been developed which can combine multiple nonhierarchical LF functions and enable adaptive learning [17, 18]. However, kriging and kriging-based methods have fundamental numerical limitations. Fitting a kriging model requires optimization of the hyperparameters \(\theta\), where a different \(\theta\) parameter exists for each dimension. Additionally, each evaluation of the loss function requires inverting the covariance matrix, an operation with computational complexity on the order of \(O(d^{3})\), where \(d\) is the number of dimensions. This makes kriging-based methods poorly suited for modeling functions above a couple of dozen dimensions, especially if the number of training samples is high.

In this work, we use Neural Networks (NNs) to avoid the limitations of kriging. A common NN architecture, feed forward fully-connected, is shown in Fig. 1. A neural network is a structure composed of layers of neurons with weighted connections to the neurons of other layers. Each neuron takes as input the weighted sum of neurons in the previous layer, plus the neuron's bias term. This input is then transformed by the neuron's activation function, and output in turn to the neurons of the next layer.

Figure 1: Illustration of a fully connected multi-layer NN. A single hidden layer with sufficiently many neurons will enable a NN to model any continuous function.
Training a NN involves adjusting the many weights and biases to output accurate predictions on training data. Typically, weight optimization is performed using a gradient descent algorithm such as Adam, with the gradients calculated via backpropagation.

The more weights and biases a NN has, the higher the degree of flexibility. This allows NNs to function as universal approximators. Specifically, the Universal Approximation Theorem (UAT) states that a feed-forward neural network with a single hidden layer can approximate any continuous function to any degree of accuracy over a bounded region if sufficiently many neurons with nonpolynomial activation functions are included in the hidden layer [19, 20, 21]. However, the UAT does not put upper bounds on the error of a neural network with a limited number of neurons. Therefore, in practice, the number of layers and neurons in each layer is typically selected based on experience or trial and error.

Unfortunately, using conventional NNs for aerospace engineering design studies faces several practical drawbacks. In aerospace applications, physical experiments or high-fidelity simulations are typically expensive, meaning HF data will be sparse. Therefore, design studies often require extrapolating beyond the bounds of the HF data, which interpolators such as NNs are poorly suited to.

One approach for mitigating these practical problems is to use Physics Informed Neural Networks (PINNs) [22, 23]. These combine physics information with a NN model. Typically, the physics information takes the form of differential equations, such as governing equations and boundary conditions. PINNs have various structures, such as NNs modifying parameters of the differential equations to improve accuracy, or NNs adding corrections and detail on top of a simplified model that does not capture the full physics. The physics information constrains the model, alleviating overfitting and increasing the accuracy of extrapolations beyond the HF data points.

PINNs can also use physics models other than differential equations, such as first-principle models, data-driven models, and expert knowledge models. PINNs can also include multiple fidelities of data, for instance, by training a NN on LF training data and using the output as an input to the PINN trained on HF data [24, 25, 26].

Recently, the authors proposed the Emulator Embedded Neural Network (E2NN) [27, 28] a generic framework for combining any mix of physics-based models, experimental data, and any other information sources. Instead of using an LF model as an input to the main NN, E2NN embeds emulator models into its hidden layer neurons intrusively. A neuron with an LF model embedded into it is called an emulator. The E2NN is trained to transform and mingle the information from multiple emulators within the NN by finding the optimal connection weights and biases. This has an effect like that of a standard PINN architecture, reducing overfitting and enabling accurate extrapolations beyond sparse HF data points. E2NN performed well in the tests of the reference papers, handling non-stationary function behaviors, capturing a high dimensional Rosenbrock function with low-quality LF models, and successfully modeling stress in a Generic Hypersonic Vehicle (GHV) engineering example.

To enable AL, it is necessary to capture epistemic prediction uncertainty for a data acquisition metric. Using a Bayesian NN [29, 30] is one option for obtaining epistemic learning uncertainty. In a Bayesian NN, the weights are assigned probability distributions rather than fixed values. Therefore, the prediction output is not a point estimate but a Gaussian distribution. Bayesian NNs are trained through backpropagation, but they require significantly higher computational costs than conventional NNs to optimize the probability density distribution of weights and biases.

Ensemble methods [31, 32] combine multiple models to improve accuracy and provide uncertainty information. Accuracy improves more when the errors of the models are less correlated. In the extreme case where \(M\) models have equal and entirely uncorrelated errors, averaging the models decreases error by a factor of \(1/M\)[33]. Error correlation decreases when models are less similar, which can be achieved by training on different subsets of data or using models with different underlying assumptions [34].

Uteva et al. [35] tested an active learning scheme of sampling where two different GPR models disagreed the most. Because the GPR models were trained on different sets of data, they had different \(\theta\)-hyperparameters and thus different behavior. This method outperformed sampling the location of maximum predicted variance of a single GPR model. Lin et al. [36] combined two NNs with different architectures into an ensemble, and added adaptive samples at the locations of maximum disagreement. Christiano et al. [37] combined three neural networks in an ensemble and added additional training samples where the variance of the predictions was maximized.

While existing ensemble methods can output a location of maximum variance for adaptive learning, they
do not output a predictive probability distribution like GPR does. Such a predictive distribution is highly useful, offering additional insight into the model and enabling the use of probabilistic acquisition functions such as Expected Improvement. In this work, we present a formal statistical treatment to extract such a probability distribution. Specifically, we estimate the epistemic learning uncertainty by combining multiple E2NN models and calculating a Bayesian predictive distribution.

Training a single E2NN model typically requires several minutes to reach convergence, which means updating the entire ensemble can be time-consuming. Therefore, a method combining Neural Networks and linear regression is explored to alleviate the cost of neural network training. Early exploration of the method was performed by Schmidt et al. [38]. The method combines two basic steps: first, creating a neural network with random weights and biases, and second, setting the last layer of weights using a linear regression technique such as ordinary least squares or ridge regression. These steps are far cheaper than backpropagation-based neural network training, enabling accelerated retraining of the E2NN ensemble during active learning.

The next section, Section 2, contains a brief review of E2NN. Section 3 covers the proposed approach of adaptive and rapid learning, which is followed by Section 4 showing fundamental mathematical demonstrations and a practical aerospace engineering example involving predicting the aerodynamic performance of the generic hypersonic vehicle.

## 2 Review of Emulator Embedded Neural Networks

The E2NN has low fidelity models, called emulators, embedded in specific neurons of the neural network's architecture. The emulators contribute to regularization and preconditioning for increased accuracy. The emulators can take the form of any low-cost information source, such as reduced/decomposed physics-based models, legacy equations, data from a related problem, models with coarsened mesh or missing physics, etc.

An emulator embedded in the last hidden layer can only be scaled before being added to the response, while an emulator embedded in the second-to-last hidden layer can be transformed through any functional mapping (by the Universal Approximation Theorem [19, 20, 21]) before being added to the response. The simplest solution when selecting the architecture is to embed the emulator in all hidden layers and allow the E2NN training to select which instances are useful. The flexibility of E2NN in selecting between an arbitrary number of LF models, each embedded in multiple hidden layers, enables wide applicability to problems with LF models of high, low, or unknown accuracy. The architecture of E2NN is illustrated in Fig. 2.

Figure 2: Architecture of an Emulator Embedded Neural Network.
## 3 Proposed Approach: Adaptive Learning with Non-Deterministic Emulator Embedded Neural Network

The main technical contribution of this work is a novel method for combining predictions of an ensemble of models to approximate epistemic modeling uncertainty. This uncertainty information lowers training data costs by enabling Adaptive Learning (AL). For optimization problems, the Expected Improvement (EI) metric identifies adaptive sampling locations of maximum information gain. The EI assessment requires a probability density function, which is estimated by training an ensemble of E2NN models and estimating the underlying distribution of predictions. Greater disagreement among individual E2NN realizations is typically due to a lack of training samples and implies greater epistemic uncertainty.

A second contribution involved speeding up E2NN training hundreds of times using Rapid Neural Network (RaNN) training, allowing fast model updating during AL iterations. Ensemble uncertainty estimation is discussed in Subsection 3.1, and AL using EI is explained in Subsection 3.2 Finally, the RaNN methodology is summarized in Subsection 3.3, and practical considerations for avoiding numerical errors are discussed in Subsection 3.4.

### Proposed Approach for Assessing Epistemic Modeling Uncertainty of E2NN

The epistemic modeling uncertainty is estimated using an ensemble of E2NN models. The E2NN model predictions agree at the training points, but the predictions between points depend on the random initializations of the weights and biases. In other words, the E2NN model predictions are samples drawn from a stochastic prediction function in the design space. A higher magnitude of disagreement implies greater epistemic modeling uncertainty. This suggests a useful assumption: model the epistemic uncertainty as the aleatoric uncertainty of the E2NN model predictions. Specifically, assume the two pdfs are approximately equal at each prediction point \(x\):

\[p(y_{true}(x))\approx p(y_{E2NN}(x)) \tag{1}\]

However, finding the exact aleatoric distribution of \(y_{E2NN}(x)\) requires training an infinite number of E2NN models. Instead, we use a finite number of model predictions at a point \(x\) as data \(D(x)\) to construct a posterior predictive distribution (ppd) for the E2NN predictions. This ppd is then used as an estimate of the epistemic pdf of the true function.

\[p(y_{true}(x))\approx p(y_{E2NN}(x)|D(x)) \tag{2}\]

Finding the ppd requires introducing a second assumption: The aleatoric E2NN predictions are normally distributed at each point \(x\)

\[y_{E2NN}(x)\sim N(\mu,\sigma) \tag{3}\]

The process of combining multiple E2NN model predictions to estimate the epistemic uncertainty is illustrated in Fig. 3.

For the general case of finding a posterior predictive distribution for a normal distribution from which we have \(n\) iid data points \(D\) but no additional information, we include a proof in Appendix A. First, we need a prior and a likelihood. The prior for the mean and variance is a normal-inverse-chi-squared distribution (\(NI\chi^{2}\))

\[p(\mu,\sigma^{2})=NI\chi^{2}(\mu_{0},\kappa_{0},\nu_{0},\sigma_{0}^{2})=N(\mu |\mu_{0},\sigma^{2}/\kappa_{0})\cdot\chi^{2}(\sigma^{2}|\nu_{0},\sigma_{0}^{2}) \tag{4}\]

Here \(\mu_{0}\) is the prior mean and \(\kappa_{0}\) is the strength of the prior mean, while \(\sigma_{0}^{2}\) is the prior variance and \(\nu_{0}\) is the strength of the prior variance. The likelihood of the \(n\) iid data points is the product of the likelihoods of the individual data points. By Bayes' rule, the joint posterior distribution for \(\mu\) and \(\sigma^{2}\) is proportional to the prior times the likelihood. Selecting the constants \(\kappa_{0}\), \(\sigma_{0}\) and \(\nu_{0}\) for an uninformative prior and integrating to marginalize out \(\mu\) and \(\sigma\) in the posterior yields the predictive posterior distribution. Given \(n\) iid samples, this is given by a t-distribution
\[p(y|D)=t_{n-1}\left(\bar{y},\frac{1+n}{n}s^{2}\right) \tag{5}\]

where \(\bar{y}\) is the sample mean and \(s^{2}\) is the sample variance.

\[s^{2}=\frac{1}{n-1}\sum_{i}(y_{i}-\bar{y})^{2} \tag{6}\]

The final pdf of \(y\) given the data \(D\) is a t-distribution instead of a normal distribution because it combines both Bayesian epistemic and aleatoric uncertainty. The epistemic component of the uncertainty can be reduced by increasing the number of samples. As the number of samples or ensemble models \(n\) approaches \(\infty\), the pdf of Eq. 5 will approach a normal distribution with the correct mean and standard deviation. For \(n=2\) samples (\(\nu=1\)) this yields a Cauchy, or Lorentzian, distribution, which has tails so heavy that the variance \(\sigma^{2}\) and mean \(\mu\) are undefined [39]. For \(n=3\) samples (\(\nu=2\)) the mean \(\mu\) is defined but the variance \(\sigma^{2}\) is infinite. Therefore, more than 3 ensemble models should always be used to estimate the mean when an uninformative prior is used.

Substituting Eq. 5 into Eq. 1 yields the final estimate of the epistemic modeling uncertainty distribution,

\[p(y_{true}(x))=t_{n-1}\left(\bar{y}_{E2NN}(x),\frac{1+n}{n}\cdot s_{E2NN}(x) \right)^{2} \tag{7}\]

where \(\bar{y}_{E2NN}(x)\) is the sample mean prediction of the \(n\) E2NN models in the ensemble,

\[\bar{y}_{E2NN}(x)=\frac{1}{n}\sum_{i}y_{E2NN_{i}}(x) \tag{8}\]

and \(s_{E2NN}(x)\) is the sample variance of E2NN model predictions in the ensemble.

\[s_{E2NN}(x)^{2}=\frac{1}{n-1}\sum_{i}(y_{E2NN_{i}}(x)-\bar{y}_{E2NN}(x))^{2} \tag{9}\]

Use of an uninformative prior results in conservative and robust estimation of the epistemic modeling uncertainty.

Figure 3: Illustration of using multiple E2NN model realizations to estimate the underlying aleatoric probability distribution. This estimate is used as an approximation of the epistemic modeling uncertainty.
### Adaptive Learning (AL) Using Expected Improvement

Adaptive learning allows data collection to be focused in areas of interest and areas with high uncertainty, reducing the number of samples required for design exploration and alleviating computational costs. Typically, an acquisition function is used to measure the value of information gained from adding data at a new location. The overall process requires a few simple steps:

1. Generate HF responses from an initial design of experiments.
2. Use sample data to build an ensemble of E2NN models.
3. Maximize the acquisition function using an optimization technique.
4. If the maximum acquisition value is above tolerance, add a training sample at the location and go to step 2. Otherwise, stop because the optimization has converged.

The criteria for where to add additional data depends on the design exploration goals. A common goal is global optimization. Global minimization seeks to find the minimum of a function in \(d\)-dimensional space.

\[x_{opt}=\operatorname*{argmin}_{x\in\mathbb{R}^{d}}y(x) \tag{10}\]

For this objective, Expected Improvement (EI) is applicable [4, 18]. Informally, EI at a design point is the amount by which the design will, on average, improve over the current optimum. Formally, this is calculated by integrating the product of the degree and likelihood of every level of improvement. Likelihoods are given by the predictive probability distribution. Any new sample worse than the current optimum yields an improvement of 0. The general expression for EI is given by Eq. 11 and illustrated in Fig. 4.

\[EI(x)=\int_{-\infty}^{\infty}pdf(y)\cdot max(y_{current\ opt}-y,0)\cdot dy \tag{11}\]

The EI value for a Gaussian predicted probability distribution is expressed as the closed-form expression,

\[EI(x)=\left(f_{min}-\hat{y}(x)\cdot\Phi\left(\frac{f_{min}-\hat{y}(x)}{\sigma_ {z}(x)}\right)+\sigma_{z}(x)\cdot\phi\left(\frac{f_{min}-\hat{y}(x)}{\sigma_ {z}(x)}\right)\right) \tag{12}\]

where \(\phi(\cdot)\) and \(\Phi(\cdot)\) are the pdf and cdf of the standard normal distribution, respectively; \(\hat{y}(x)\) and \(\sigma_{z}(x)\) are the mean and standard deviation of the predictive probability distribution, respectively; and \(f_{min}\) is the current optimum. The current optimum can be defined as either the best sample point found so far, or the best mean prediction of the current surrogate. We use the former definition. The two definitions approach each other as the adaptive learning converges on the optimum.

Figure 4: Illustration of Expected Improvement calculation for adaptive learning.
Unlike a kriging model, an E2NN ensemble returns a Student's t-distribution instead of a Gaussian distribution. The resulting formulation of EI in this case is [40]

\[E[I(x)]=(f_{min}-\mu)\Phi_{t}(z)+\frac{\nu}{\nu-1}\left(1+\frac{z^{2}}{\nu} \right)\sigma\phi_{t}(z) \tag{13}\]

where the t-score of the best point found so far is

\[z=\frac{f_{min}-\mu}{\sigma} \tag{14}\]

and \(\phi_{t}(\cdot)\) and \(\Phi_{t}(\cdot)\) are the pdf and cdf of the standard Student's t-distribution, respectively. Also, \(\mu\), \(\sigma\), and \(\nu\) are the mean, scale factor, and degrees of freedom of the predictive t-distribution.

### Training of E2NNs As Rapid Neural Networks

Typically, NNs are trained using auto-differentiation and backpropagation. The main problem in NN training is saddle points, which look like local minima because all gradients are zero and the objective function is increasing along almost all direction vectors. However, in a tiny fraction of possible directions the objective function decreases. If a NN has millions of weights, determining which direction to follow to escape from the saddle point is difficult. Optimizers such as Adam use momentum and other techniques to reduce the chance of getting stuck at a saddle point. However, a single E2NN model still takes significant computational time to reach convergence, requiring minutes to hours for a moderately sized neural network. Depending on the number of samples and the dimensionality of the problem, training an ensemble of E2NNs for adaptive learning will introduce significant computational costs.

To increase speed, the E2NN models can be trained as Rapid Neural Networks (RaNNs). Essentially, this involves initializing a neural network with random weights and biases in all layers and then training only the last layer connections. The last layer's weights and bias are trained by formulating and solving a linear regression problem such as ridge regression, skipping the iterative training process entirely. This accelerates training multiple orders of magnitude. These models are sometimes referred to as extreme learning machines, a term coined by Dr. Guang-Bin Huang, although the idea existed previously.

An early example of RaNN by Schmidt et al. in 1992 utilized a feed-forward network with one hidden layer [38]. The weights connecting the input to the first layer were set randomly, and the weights connecting the hidden layer to the output were computed by minimizing the sum of squared errors. The optimal weights are analytically computable using standard matrix operations, resulting in very fast training. Saunders et al. [41] used ridge regression instead of least squares regression when computing the weights connecting the hidden layer to the output.

Huang et al. later developed equivalent methods [42, 43], and demonstrated that like conventional NNs, RaNNs are universal approximators [44], with the ability to capture any continuous function over a bounded region to arbitrary accuracy if a sufficient number of neurons exist in the hidden layer.

However, despite being universal approximators, RaNNs require many more hidden layer neurons to accurately approximate a complex function than NNs trained with backpropagation and gradient descent. This is because the hidden layer neurons are essentially functions that are linearly combined to yield the NN prediction. Backpropagation intelligently selects these functions, while RaNN relies on randomly chosen functions and therefore requires more of them to construct a good fit. If fewer neurons and higher robustness are desired for the final model, an ensemble of E2NN models can be trained using backpropagation after the adaptive learning has converged.

In this work, we propose applying the random initialization and linear regression techniques to E2NN models instead of standard NNs. Multiple realizations of E2NN with RaNN training are combined into an ensemble, enabling uncertainty estimation and adaptive learning. The proposed framework of adaptive learning is demonstrated in multiple example problems in Section 4.

### Practical Considerations for Avoiding Large Numerical Errors

Setting the last layer weights using linear regression sometimes causes numerical issues when capturing highly nonlinear functions, even when enough neurons are included. If ridge regression or LASSO is used,
the E2NN model will not interpolate the training data points. If no regularization is used, the weights become very large to force interpolation. Large positive values are added to large negative values, resulting in round-off error of the E2NN prediction. The resulting fit is not smooth, but jitters up and down randomly.

Numerical stability can be improved by using a Fourier activation function such as \(\sin(x)\). This is reminiscent of a Fourier series, which can capture even non-continuous function to arbitrary accuracy. In fact, a NN with Fourier activation functions and a single hidden layer with \(n\) neurons can capture the first \(n\) terms of a Fourier series. Each neuron computes a different term of the Fourier Series, with the first layer of weights controlling frequency, the bias terms controlling offset, and the second layer weights controlling amplitudes. However, when using rapid training only the last layer weights are optimized.

As points are added to a highly nonlinear function, interpolation becomes more difficult and numerical instability is introduced despite the Fourier activation. This is counteracted by increasing the frequency of the Fourier activation function, which enables more rapid changes of the fit. For smooth or benign functions, the Swish activation function tends to outperform Fourier. Therefore, we include some E2NN models using Swish and some using Fourier within the ensemble. Any model with any weights of magnitude above the tolerance of 100 are considered unstable and dropped from the ensemble. Additionally, any model with NRMSE on the training data above the tolerance of 0.001 is dropped from the ensemble, where NRMSE is defined as

\[NRMSE=\sqrt{\frac{\sum_{i=1}^{N}(\hat{y}_{i}-Y_{i})^{2}}{\sum_{i=1}^{N}(\bar{ Y}-Y_{i})^{2}}} \tag{15}\]

for predictions \(\hat{y}_{i}\) and \(N\) training samples \(Y_{i}\). If over half of the Fourier models are dropped from the ensemble, all Fourier models have their frequencies increased and are retrained. These changes eliminate noisy and unstable models from the ensemble, as well as modifying models to remain stable as new points are added to the training data.

## 4 Numerical Experiments

In aerospace applications, the costs of HF sample generation, i.e., computational fluid dynamics simulation, aeroelasticity analysis, coupled aerothermal analysis, etc., are typically far higher than the costs of generating LF samples and training NN models. Therefore, in the following examples, we compare the cost of various prediction models in terms of the number of HF samples required, rather than the computer wall-clock or GPU time. We assume that enough LF samples are collected to train an accurate metamodel, which is used to cheaply compute the emulator activations whenever the neural network makes predictions.

In the following examples, we use a fully connected feed-forward neural network architecture. All LF functions are embedded as emulators in all hidden layers. The input variables are scaled to \([-1,1]\). The weights are initialized using the Glorot normal distribution. Biases are initialized differently for Swish and Fourier activation functions. Fourier biases are initialized uniformly between \([0,2\pi]\) to set the functions to a random phase or offset, while Swish biases are initialized uniformly on the region \([-4,4]\).

In total, each ensemble contains 16 E2NN models with a variety of architectures and activation functions. Identical models make the same assumptions about underlying function behavior when deciding how to interpolate between points. If these assumptions are incorrect, the error will affect the whole ensemble. Therefore, including dissimilar models results in more robust predictions.

Four activation functions and two different architectures yield 8 unique NN models. Each unique model is included twice for a total of 16 E2NN models in the ensemble. The four activation functions are \(\text{swish}(x)\), \(\sin(\text{scale }x)\), \(\sin(1.1\text{ }scale\ x)\), and \(\sin(1.2\text{ }scale\ x)\). The two architectures include a small network and a large network. The small network has a single hidden layer with \(2n\) neurons, where \(n\) is the number of training samples. This means the number of neurons is dynamically increased as new training samples are added. The large network has two hidden layers, where the first hidden layer has 200 neurons, and the second hidden layer has 5000 neurons. Having most neurons in the second hidden layer enables more of the NN weights to be adjusted by linear regression. The large and small NNs use different scale terms for the Fourier activation functions. The large NN scale term is increased whenever more than half the large Fourier NNs have bad fits, and the small NN scale term is increased whenever over half the small Fourier NNs have bad fits. Because numerical instability is already corrected for, we do not use ridge regression. Instead,
we perform unregularized linear regression using the Moore-Penrose inverse with a numerically stabilized \(\Sigma\) term.

### One-dimensional analytic example with a linearly deviated LF model

An optimization problem with the following form is considered.

\[x_{opt}=\operatorname*{argmin}_{x\in[0,1]}y_{HF}(x) \tag{16}\]

Here \(x\) is a design variable on the interval \([0,1]\). The high-fidelity function \(y_{HF}(x)\) and its low-fidelity counterpart \(y_{LF}(x)\) are given in Eqs. 17 and 18. These functions have been used previously in the literature when discussing MF modeling methods [17, 45].

\[y_{HF}(x)=(6x-2)^{2}\sin(12x-4) \tag{17}\]

\[y_{LF}(x)=0.5y_{HF}(x)+10(x-0.5)-5 \tag{18}\]

As shown in Fig. 4(a), the initial fit uses three HF samples at \(x=[0,0.5,1]\) and an LF function which is linearly deviated from the HF function. The 16 E2NN models used in the ensemble are shown in Fig. 4(b). Three ensemble models are outliers, significantly overestimating the function value near the optimum. Two of these models nearly overlap, looking like a single model. All three of these inferior models are small NNs with Fourier activation functions.

The mean and 95% probability range of the predictive t-distribution are both shown in Fig. 5(a). From this t-distribution, the Expected Improvement is calculated in Fig. 5(b). A new sample is added at the location of maximum EI.

The true optimum is \(x_{opt}=0.7572\), \(y_{opt}=-6.0207\). As shown in Fig. 6(a), the first adaptive sample at \(x=0.7545\) lands very near the optimum, and the retrained ensemble's mean prediction is highly accurate.

Figure 5: Initial problem and fitting of the E2NN model.
After the first adaptive sample, the maximum EI is still above tolerance. Therefore, an additional sample is added as shown in Fig. 6(b). The second adaptive sample at \(x=0.7571\) is only \(10^{-4}\) from the true optimum. After the second adaptive sample, the maximum EI value falls below tolerance, and the adaptive sampling converges.

### Two-dimensional analytical example

The proposed ensemble method is compared with the popular kriging method for minimization of the following two-dimensional function.

\[y_{HF}(x_{1},x_{2})=\sin(21(x_{1}-0.9)^{4})\cos(2(x_{1}-0.9))+(x_{1}-0.7)/2+2*x_ {2}^{2}\sin(x_{1}x_{2}) \tag{19}\]

This nonstationary test function was introduced in [46, 47][49, 50]. The kriging method uses only HF training samples during optimization, but E2NN makes use of the following LF function.

\[y_{LF}(x_{1},x_{2})=\frac{y_{HF}(x_{1},x_{2})-2+x_{1}+x_{2}}{1+0.25x_{1}+0.5x_{2}} \tag{20}\]

Figure 6: Initial model and expected improvement.

Figure 7: Iterative model as adaptive samples are added.
The independent variables \(x_{1}\) and \(x_{2}\) are constrained to the intervals \(x_{1}\in[0.05,1.05]\), \(x_{2}\in[0,1]\). The LF function exhibits nonlinear deviation from the HF function as shown in Fig. 8.

Eight training samples are selected to build the initial model using Latin hypercube sampling. The initial E2NN ensemble prediction is shown in Fig. 8(a). The resulting EI is shown in Fig. 8(b), and the ensemble prediction after a new sample is added is shown in Fig. 8(c). The initial fit is excellent, with only a small difference between the mean prediction and HF function. After the first adaptive sample is added, the optimum is accurately pinpointed.

To compare the performance of adaptive learning with single fidelity kriging, EGO with a kriging model is run with the same initial set of 8 points. The best kriging sample is shown in Fig. 9(a). After adding a sample near the location of the optimum, the kriging model still does not capture the underlying trend of the HF model, as shown in Fig. 9(c). The E2NN ensemble maintains higher accuracy over the design space by leveraging the LF emulator.

For the E2NN ensemble, the algorithm adds one more sample further up the valley that the optimum lies in, and then terminates because the expected improvement converges below tolerance. The final fit is shown in Fig. 10(a). The kriging model adds 29 samples before it converges, and still doesn't find the exact optimum, as shown in Fig. 10(b).

Figure 8: Comparison of HF and LF functions for a nonstationary test function.

Figure 9: Adaptive sampling of E2NN ensemble.
### Three-dimensional CFD example using a Hypersonic Vehicle Wing

This example explores modeling the Lift to Drag Ratio (\(CL/CD\)) of a wing of the Generic Hypersonic Vehicle (GHV) given various flight conditions. The GHV was developed at the Wright-Patterson Air Force Base to allow researchers outside the Air Force Research Laboratory to perform hypersonic modeling studies. The parametric geometry of the wing used was developed by researchers at the Air Force Institute of Technology [48] and is shown in Fig. 12.

For this example, we studied the maximum lift-to-drag ratio of the GHV wing design with respect to three operational condition variables: Mach number (normalized as \(x_{1}\)), Angle of Attack (normalized as \(x_{2}\)) and Altitude (normalized as \(x_{3}\)). The Mach number ranges from \([1.2,4.0]\), while the angle of attack ranges from \([-5^{\circ},8^{\circ}]\) and the altitude ranges from \([0,50\text{ km}]\). The speed of sound decreases with altitude, so the same Mach number denotes a lower speed at higher altitude. Atmospheric properties were calculated using the scikit-aero python library based on the U.S. 1976 Standard Atmosphere.

FUN3D [49] performed the CFD calculations. We used Reynolds Averaged Navier Stokes (RANS) with a mesh of \(272,007\) tetrahedral cells for the HF model, and Euler with a mesh of \(29,643\) tetrahedral cells for the LF model. To enable rapid calling of the LF model during each NN evaluation, \(300\) evaluations of the LF model were performed and used to train a GPR model. This GPR model was then used as the LF function. The HF and LF meshes are compared in Fig. 13.
the Angle of Attack (\(x_{2}\)) is the most influential variable, followed by Mach number (\(x_{1}\)), with Altitude (\(x_{3}\)) contributing little effect. The models show similar trends, except for Mach number which is linear according to the LF model and quadratic according to the HF model. The captured viscous effects and finer mesh enable the HF model to capture more complexity resulting from the underlying physics.

The problem is formulated as minimizing the negative of the lift-to-drag ratio rather than maximizing the lift-to-drag ratio, following optimization convention. Both the ensemble and GPR models are initialized with 10 HF samples selected using Latin Hypercube Sampling. Five different optimization runs are completed for each method, using the same 5 random LHS initializations to ensure a fair comparison.

The optimization convergence as samples are added is shown in Fig. 15. Both models start with the same initial optimum values because they share the same initial design of experiments. However, the E2NN ensemble model improves much more quickly than GPR, and converges to a better optimum. E2NN requires only 11 HF samples to reach a better optimum on average than GPR reached after 51 HF samples. In some cases, the GPR model converges prematurely before finding the optimum solution. E2NN much more consistently finds the optimum of \(-13\), which corresponds to a lift-to-drag ratio of 13.

Figure 12: Illustration of GHV wing used in CFD analysis.

Figure 13: Comparison of HF and LF meshes used in CFD analysis.
## 5 Conclusions

In this research, we present a novel framework of adaptive machine learning for engineering design exploration using an ensemble of Emulator Embedded Neural Networks (E2NN) trained as rapid neural networks. This approach allows for the assessment of modeling uncertainty and expedites learning based on the design study's objectives. The proposed framework is successfully demonstrated with multiple analytical examples and an aerospace problem utilizing CFD analysis for a hypersonic vehicle. The performance of the proposed E2NN framework is highlighted when compared to a single-fidelity kriging optimization. E2NN exhibited superior robustness and reduced computational costs in converging to the true optimum.

The central contribution in this work is a novel technique for approximating epistemic modeling uncertainty using an ensemble of E2NN models. The ensemble predictions are analyzed using Bayesian statistics to produce a t-distribution estimate of the true function response. The uncertainty estimation methodology allows for active learning, which reduces the training data required through efficient goal-oriented design exploration. This is a sequential process of intelligently adding more data based on the ensemble, and then rebuilding the ensemble. The training of the ensemble is drastically accelerated by applying the Rapid Neural Network (RaNN) methodology, enabling individual E2NNs to be trained near-instantaneously. This

Figure 14: Comparison of HF and LF CFD models.

Figure 15: Comparison of E2NN and GPR convergence towards optimum \(-CL/CD\). Average values across all 5 runs are shown as thick lines, while the individual run values are shown as thin lines.
speedup is enabled by setting some neural network weights to random values and setting others using fast techniques such as linear regression and ridge regression.

The essential components of the proposed framework are 1) the inclusion of emulators, 2) the ensemble uncertainty estimation, 3) active learning, and 4) the RaNN methodology. These components work together to make the overall methodology feasible and effective. For instance, the active learning would not be possible without the uncertainty estimation. Emulators make the ensemble more robust by stabilizing individual E2NN fits, resulting in fewer defective or outlier fits and better uncertainty estimation. The inclusion of the emulators and the use of adaptive sampling reduces the cost of training data generation more than either method could individually. Finally, the RaNN methodology accelerates training of individual E2NN models by hundreds to thousands of times, preventing the many re-trainings required by the ensemble and adaptive sampling methodologies from introducing significant computational costs. All these techniques combine synergistically to create a robust and efficient framework for engineering design exploration.

**Declaration of Competing Interests**

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Research Data**

No data was used for the research described in this article. The computer codes used for the numerical examples are available at [https://github.com/AtticusBeachy/multi-fidelity-nn-ensemble-examples](https://github.com/AtticusBeachy/multi-fidelity-nn-ensemble-examples).

**Acknowledgement**

This research was sponsored by the Air Force Research Laboratory (AFRL) and Strategic Council for Higher Education (SOCHE) under agreement FA8650-19-2-9300. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied by the Strategic Council for Higher Education and AFRL or the U.S. Government.