Surrogate models serve as efficient approximations for expensive high-fidelity simulations, with the goal of providing significant computational savings while maintaining solution accuracy. Data-driven surrogate modeling approaches like Gaussian process regression (GPR)(9; 52; 46), polynomial chaos expansion (PCE) (16; 57; 39; 58), response surface methods(19; 18), manifold-based approaches (17; 26), and deep neural networks (DNNs) (13; 41)
are gaining popularity as a means of speeding up the engineering design/analysis process. For theoretical validity, training and test data sets for surrogate model construction must have the same support. However, in many practical uses, test data are drawn from outside the compact support of the training data – resulting in the use of the surrogate model for prediction on cases for which it has not been trained (i.e. extrapolation rather than interpolation). This may occur, for example, when test data are drawn from a different probability distribution than the training data (so-called out-of-distribution, OOD, data) or when input are degraded by added noise. The ability to construct surrogate models that generalize beyond the immediate support of their training data and into noisy data regimes is particularly challenging and is expressly connected to the parameterization of the surrogate. Resolving this challenge is of broad interest to the physics-based modeling community.
GPR and PCE-based surrogates are particularly favorable because they involve relatively few hyperparameters and thus are easy to implement. However, they fail to handle physical domains with complex geometry, unseen noisy input data, and suffer from the curse-of-dimensionality. Manifold-based approaches, on the other hand, can handle high-dimensional inputs(30). A recent study on manifold-based PCE (m-PCE) (27) compares several manifold learning methods and provides insights into their effectiveness for the construction of accurate and generalizable surrogates for physics-based models. The key benefits of these manifold-based methods, and m-PCE in particular, are their predictive accuracy, robustness to noisy input data, limited number of tunable hyperparameters, and computational efficiency during training. Alternatively, DNNs have proven to be effective surrogate models across a broad spectrum of applications (45; 50; 8; 25; 35; 20). The pivotal work of (28) kicked off an intense period of research to train larger networks with more hidden units in the search of higher test performance. Despite the fact that the high complexity of such over-parameterized models allows for faultless data interpolation, they also frequently achieve low generalization error (29).
Conventional thinking in machine learning holds that using models with increasing capacity will lead to overfitting to the training data. As a result, the capacity of the model is generally controlled either by limiting the size of the model (number of parameters) or by introducing additional explicit regularization terms to avoid overfitting. However, for DNNs, expanding the model size only improves the generalization error, even when the networks are trained without any regularization term (38; 15; 43; 12). Over-parameterization refers to the case where the number of trainable parameters of a model/network is larger than the number of training observations. While over-parameterized DNNs are the state-of-the-art in machine learning, their robustness and ability to generalize even to noisy test data is what makes this architecture stand out. Optimization is usually easier in over-parameterized models (43). The idea is that over-parameterization changes the loss function, resulting in a large number of locally optimum solutions, making it easier for local search algorithms to locate an near-optimal solution. However, it has been shown that the computational cost of training an over-parameterized model grows at least as a fourth-order polynomial with respect to performance, Computation = (Performance), where performance refers to an error measure, for example, root mean squared error, and this may even be an underestimation (51). On the other side of the spectrum, under-parameterized models are characterized by limited expressivity, but are often preferred due to their simplicity, interpretability, computational efficiency, and generalizability for smooth functions. In general, when small training sets are available, adding more features yields improvements in the training process (2).
Nevertheless, model complexity and efficiency cannot be measured solely by the total number of trainable parameters since the expressivity of these parameters should be accounted for as well. For example, in neural network models increasing the depth of the network by adding new layers might lead to higher expressivity than increasing the size of existing layers even if the total number of parameters is the same. In addition, the ratio between the number of data features and number of data should also be considered. Therefore it is imperative to understand how all the aforementioned affect the robustness of the surrogate model (6).
In this work, we present a systematic study comparing under- and over-parameterized manifold-based PCE methods with over-parameterized DNNs. Through this study, we have tried to answer basic questions about how parameterization of these models relates to their ability to generalize to OOD and noisy data. To do so, we specifically investigated the influence of various factors such as model complexity, label noise and dataset size in the construction of surrogate models. From the corpus of the available models, we study under- and over-parameterization performance experiments using m-PCE with kernel-PCA and over-parameterization experiments on two state-of-the-art neural operators, the Deep Operator Network (DeepONet) (35) and the Fourier Neural Operator (FNO) (31). The comparative study includes a recent extension of the DeepONet, namely POD-DeepONet (36), while for FNO we have employed the version from https://github.com/zongyi-li/fourier_neural_operator. Additionally, we introduce some enhancements for both manifold-based surrogates and DeepONet to allow fast training and convergence and better generalization for problems with highly nonlinear and nonsmooth solutions. The new developments proposed in the current work are the following:
We perform manifold-learning for both the model inputs and model outputs for construction of both over- and under-parameterized m-PCE surrogates.
We introduce fully-trainable weight parameters for the DeepONet, referred to as DeepONet with self-adaptivity, to handle highly nonlinear and non-smooth data.
We have considered the Brusselator diffusion-reaction dynamical system to investigate the relative performance of the chosen surrogate models. All codes will become available on Github, https://github.com/katiana22/surrogate-overparameterization
, upon publication of our paper. Construction of PCE surrogates is performed using the open-source UQpy (Uncertainty Quantification with Python)(40) package. The paper is organized as follows. In Section 2, we describe the DeepONet and FNO architectures, present the manifold-based PCE models and introduce the proposed enhancements. In Section 3, we present the Brusselator system and the data generation process. In Section 4, we compare the performance of the studied models for the Brusselator diffusion reaction system. Finally, we summarize our observations and provide concluding remarks in Section 5.
2 Machine Learned Approximation Methods
Consider an analytical or a computational model,, and corresponding output quantities of interest (QoIs), where represent spatio-temporal coordinates. That is, model performs the mapping where the dimensionality of both the inputs and the outputs, is high, e.g. . Here, the high-dimensional inputs may represent random fields and/or processes, such as spatially or temporally varying coefficients, and the corresponding QoIs represent physical quantities, which can also vary in both space and time. Our objective is to approximate the mapping, from a training dataset of input-output pairs where and , and achieve the lowest possible predictive error on a test dataset. In this section, we outline two classes of approaches to learn this approximation: operator learning and manifold learning based surrogates. We then systematically test these approximations for generalization error on out-of-distribution (OOD) and noisy input data for a set of nonlinear PDEs – the Brusselator reaction-diffusion equations – in subsequent sections.
2.1 Operator learning
Neural operators learn nonlinear mappings between infinite dimensional function spaces on bounded domains, providing a unique simulation framework for real-time prediction of multi-dimensional complex dynamics. Once trained, such models are discretization invariant, which means they share the same network parameters across different parameterizations of the underlying functional data. In recent times, a lot of research is being carried out in the domain of operator regression (20; 32). To study the benefits of over-parameterization of DNNs, we have implemented two operator networks that have shown promising results so far, the DeepONet (35) and the Fourier neural operator (FNO) (31). Although the original DeepONet architecture proposed in (35) has shown remarkable success, several extensions have been proposed in (36) to modify its implementation and produce efficient and robust architectures. In this section, we first provide an overview of the DeepONet framework put forth in (35). This is followed by a proposed extension of the DeepONet, DeepONet with self-adaptivity, to enhance the prediction capability of the network for non-smooth problems. We then briefly review the POD-DeepONet and FNO, which we use in this work.
The DeepONet architecture consists of two DNNs: one encodes the input function at fixed sensor points (branch net), while another encodes the information related to the spatio-temporal coordinates of the output function (trunk net). The goal of the DeepONet is to learn the solution operator, that approximates , and can be evaluated at continuous spatio-temporal coordinates, (input to the trunk net). The output of the DeepONet for a specified input vector, , is a scalar-valued function of expressed as , where includes the trainable parameters (weights, , and biases, ) of the networks.
In general, the input to the branch net is flexible, i.e., it can take the shape of the physical domain, the initial or boundary conditions, constant or variable coefficients, source terms, etc. But it must be discretized to represent the input function, in a finite-dimensional space using a finite number of points, termed as sensors. We specifically evaluate at fixed sensor locations to obtain the pointwise evaluations, , which are used as input to the branch net. The trunk net takes as input the spatial and the temporal coordinates, e.g. , at which the solution operator is evaluated to compute the loss function. The solution operator for an input realization, , can be expressed as:
where are outputs of the branch net and are outputs of the trunk net. Conventionally, the trainable parameters of the DeepONet, represented by in Eq. (1), are obtained by minimizing a loss function, which is expressed as:
where and denote the residual loss and the initial condition loss, respectively.
The DeepONet model provides a flexible paradigm that does not limit the branch and trunk networks to any particular architecture. For an equispaced discretization of the input function (used here), a Convolutional Neural Network (CNN) could be used for the branch net architecture. For a sparse representation of the input function, one could also use a Feedforward Neural Network (FNN). A standard practice is to use a FNN for the trunk network to take advantage of the low dimensions of the evaluation points,. A schematic representation of the standard DeepONet used to approximate the Brusselator dynamics described below is shown in Figure 1, where the branch net is considered as a CNN, and the trunk net has a FNN architecture.
DeepONet aims to infer a continuous latent function that arises as the solution to a system of nonlinear PDEs. Neural networks employ back-propagation algorithm to tune the network parameters, while trying to reflect the best fit solution to the training data. Approximating steep gradients is often challenging for neural networks. Efforts in the literature aim to impose such constraints in a soft manner by appropriately penalizing the loss function of CNN approximations. However, manually manoeuvring for the optimal penalizing parameter is a tedious task and often leads to unstable and erroneous predictions. We propose a simple solution to this problem, DeepONet with self-adaptivity, where the penalty parameters are trainable, so the neural network learns by itself which regions of the solution are difficult and is forced to focus on them.
A. DeepONet with Self-adaptivity
The basic idea behind self-adaptivity is to make the penalty parameters increase where the corresponding loss is higher, which is accomplished by training the network to simultaneously minimize the losses and maximize the value of the penalty parameters. Discontinuities or non-smooth features in the solution lead to a non-differentiable loss function at these locations. Hence, even though the overall error is reduced during the training, there is significant error near the discontinuity. The conventional approach to take care of such specific local minimization issues is to introduce constant penalty terms to force the network to satisfy the conditions. Accordingly, the loss function would be of the form:
where and are manually chosen constant parameters, and and are location-specific losses corresponding to the location of the discontinuity. The optimal values for and differ widely for different PDEs. Considering that in a DeepONet we are trying to find the solution operator for multiple PDEs by training the network just once, obtaining an optimal penalty parameter is a challenge.
In this section, we introduce a simple solution to tune these parameters: make these parameters trainable. These hyper-parameters are updated by back-propagation together with the network weights. The new loss function is therefore expressed as:
where are self-adaptive hyperparameters. Typically, in a neural network, we minimize the loss function with respect to the network parameters, . However, in this approach we additionally maximize the loss function with respect to the trainable hyperparameters using a gradient descent/ascent procedure. The proposed approach shares the same concept of introducing trainable hyperparameters as in (37).
In this work, we use an existing extension of the standard DeepONet proposed in (36) referred to as the POD-DeepONet. The standard DeepONet employs the trunk net to learn the basis of the output function from the data. However, in POD-DeepONet, the basis functions are computed by performing proper orthogonal decomposition (POD) on the training data (after the mean has been excluded), and using this basis in place of the trunk net. A deep neural network is employed in the branch net to learn the POD basis coefficients such that the output can be written as:
where is the mean function of all computed from the training dataset, and are the precomputed POD modes of .
2.1.2 Fourier Neural Operator
The Fourier Neural Operator (FNO) is based on parameterizing the integral kernel in the Fourier space. The concept was initially proposed in (31). FNO in its continuous form can be conceptualized as a DeepONet with a particular network architecture of the branch net and the trunk net expressed by discrete trigonometric basis functions. Neural networks are typically trained to approximate functions that are described in the Euclidean space, the conventional graph with spatial coordinate axes. However, in FNO, the network parameters, inputs to the network and the output from the network are specified in the Fourier space, a distinct type of graph used to represent frequencies.
FNO employs evaluations restricted on an equispaced mesh to discretize both the input and output functions, where the mesh and the domain must be the same. Simplistically, the input random field, , the evaluation locations, , and the output function, , are defined on the same domain with the same equispaced discretization. The input to the network is mapped on a linear layer, converting the input to a higher dimension, followed by passing it through four layers of integral operators and later projecting it back to the original dimension. Within each integral operator operation, there is a additional linear and nonlinear transform where the Fourier layers act as low-pass spectral filters. In practice, the retention of more Fourier modes (maximum is half the resolution+1) has helped in achieving better accuracy. Details on the implementation of FNO can be found in (31).
2.2 Manifold learning based surrogates
Manifold learning based surrogates are used to project high-dimensional input-output data from a mathematical model onto a low-dimensional manifold where interpolation is performed for rapid solution prediction. The manifold-based polynomial chaos expansion or manifold PCE (m-PCE), specifically aims to construct a PCE surrogate model to map the reduced input data to the corresponding QoIs (also potentially in a reduced space), thus mitigating the curse of dimensionality. A detailed survey on the use of various dimension reduction (DR) techniques for the construction of m-PCE surrogates was conducted by Kontolati et al. (27).
In the following sections, we introduce the standard approach for constructing PCE surrogates and then present the k-PCA PCE, a specific implementation of m-PCE
that leverages kernel-Principal Component Analysis (k-PCA) for dimension reduction. Despite the robust performance ofm-PCE as demonstrated in (27), in certain cases where QoIs are very high dimensional, the accuracy of m-PCE surrogates may diminish. In the following, we propose an extension of the m-PCE framework in which a manifold projection based DR is performed on both high-dimensional inputs and high-dimensional solutions.
2.2.1 Standard PCE
We assume the model , where is a
-variate random variable defined on the probability space
and characterized by the joint probability density function (PDF), where is the image space, the sample space, the set of events, and the probability measure. By further assuming that satisfies the Doob-Dynkin lemma (5), its output is a random variable dependent on . In the following, we consider for simplicity a scalar output although the extension to multivariate outputs is straightforward by applying the PCE approximation element-wise. The PCE is a spectral approximation of the form
where are scalar coefficients and are multivariate polynomials that are orthonormal with respect to the joint PDF , such that
where denotes the Kronecker delta. Depending on the PDF , the orthonormal polynomials can be chosen according to the Wiener-Askey scheme (57) or constructed numerically(55; 49). Since is assumed to consist of independent random variables , the joint PDF is given as
where is the marginal PDF of random variable . Likewise, the multivariate orthogonal polynomials are constructed as
where are univariate polynomials of degree and are orthonormal with respect to the univariate PDF , such that
The multi-index is equivalent to the multivariate polynomial degree and uniquely associated to the single index employed in Eq. (6), which can now be written in the equivalent form
where is a multi-index set with cardinality . The choice of the multi-index set plays a central role in the construction of the PCE, as it defines which polynomials and corresponding coefficients form the PCE. The most common choice, as well as the one employed in this work, is that of a total-degree multi-index set, such that includes all multi-indices that satisfy , . In that case, the size of the PCE basis is , such that it scales polynomially with the input dimension and the maximum degree .
Several approaches have been proposed for computing the coefficients , including pseudo-spectral projection (11; 10; 56), interpolation (7; 33), and, most commonly, regression (4; 34; 22; 14; 21; 23; 53). We employ a regression approach in which the PCE coefficients are determined by solving the penalized least squares problem (47)
where is a penalty factor, a penalty function acting on the vector of PCE coefficients , and is an experimental design (ED) containing realizations of with corresponding model outputs . Common choices for the penalty function are the and norms, in which cases problem (12
2.2.2 k-PCA PCE
Manifold PCE employs manifold learning techniques for DR to identify lower-dimensional embeddings for the construction of efficient and accurate PCE surrogates. DR allows us to neglect redundant features and noise, therefore avoiding overfitting and tedious training. The method is based on a two-step approach: 1) the lower-dimensional manifold is discovered from the input and/or output data and 2) a PCE surrogate is constructed to approximate the map between reduced inputs and outputs based on a limited set of training data. In general, any DR method can be employed, including linear methods like principal component analysis (PCA) and nonlinear methods such as independent component analysis (ICA), diffusion maps (DMAPs), non-negative matrix factorization (NMF), autoencoders etc. In this work, we employ k-PCA for DR.
In the original development of m-PCE(27), DR is performed only on the high-dimensional input. In this work, we extend the framework to include DR on the output leveraging the method developed in (26). We note that, in cases where the output data undergo a reduction, the DR method must possess an inverse transformation to map predicted outputs in the low-dimensional space to the physically interpretable ambient space. This proposed framework is illustrated in Figure 2 and described in detail in the following where first we describe the original m-PCE with DR only on the input and then describe the proposed methods with DR on both input and output. For completeness, we then present the k-PCA method employed for DR in this work.
Original m-PCE: DR on input only
In the original m-PCE, the method begins with an input training data set where . We assume that the intrinsic dimension of the data set . Therefore, the input data lie on or near a manifold of dimensionality that is embedded in the -dimensional space. We identify a low-dimensional representation of the original dataset , by constructing a mapping , where and represents the -dimensional vector of parameters associated with the reduction method. The DR method is learned by identifying the parameters that minimize an error measure, , as
For example, when the inverse transformation is available and the objective is to reduce the information loss after the data reduction, the error measure can be represented by the mean-square reconstruction error (30).
Using the reduced representation of the input, , we then construct a PCE surrogate model as described in Eq. (11). Given that is not generally scalar-valued, PCEs are constructed component-wise as . For prediction purposes at a new point , the point is first projected onto the low-dimensional manifold by and the PCE is used on this projected point to predict the solution. For additional details, see (27).
Proposed m-PCE: DR on both input and output
Commonly, in physics-based models, QoIs are evaluated at many points in space and time resulting in very high-dimensional () output which standard surrogate modeling techniques cannot handle without loss in predictive accuracy. Assuming that the intrinsic dimension of this output data , we propose to use DR on both the input and output data as a preprocessing step prior to the construction of the PCE surrogate. With respect to the input data, the method proceeds exactly as in the previous section. However, we introduce a second mapping , where and that serves to reduce the dimension of the output data . Here, the DR technique employed must possess an inverse map , which will act as a decoder for projecting low-dimensional outputs back to the ambient space.
Once the two embeddings have been identified, a PCE surrogate is constructed to map the reduced inputs to the outputs as . For prediction purposes, the reduced input at a new point is used to predict the reduced solution , which is then project back to the ambient space by . To evaluate the accuracy of the proposed approach, data generated by the PCE surrogate are transformed back to the original space , and compared with the ground truth. This proposed m-PCE framework is presented in Figure 2.
DR with k-PCA:
In this work, we use k-PCA to identify the lower-dimensional embeddings of both the inputs and the outputs. In this section, we present the method generically and we use to denote the original dimension of the data, which may refer to either or .
where with . Standard PCA is then performed in this higher-dimensional space. The above transformation is obtained implicitly through the use of a kernel function which replaces the scalar product using the relation , thus is never calculated explicitly. This process is known as ‘kernel trick’ and allows the transformation of datasets to very high-dimensional spaces, therefore allowing the encoding of highly nonlinear manifolds without explicit knowledge of suitable feature functions (3). Since data in the featured space are not guaranteed to have zero-mean, the kernel matrix is centralized as
where represents an matrix with values and is the kernel matrix.
The k-PCA formulation does not compute the principal components (PCs) directly, but rather it projects the data onto the PCs. The projection of a mapped data point onto the -th PC in is computed as
where again is the dot product obtained from the entries of . The coefficients
, are determined by solving the eigenvalue problem
are the eigenvalues and eigenvectors ofrespectively and is the number of data points. Application of the method is dependent on the selection of a positive semi-definite kernel function, of which there are many options including linear, polynomial, and Gaussian kernels.
3 Brusselator reaction-diffusion system
As a prototypical physico-chemical system, we consider the Brusselator diffusion-reaction system introduced by Ilya Prigogine in the 1970s (44), which describes an autocatalytic chemical reaction in which a reactant substance interacts with another substance to increase its production rate (1). The system is selected because it is mathematically described by a set of PDEs that exhibit nonlinear features that are particularly difficult for approximate models to replicate. The Brusselator model is characterized by the following reactions
where are positive parameters representing the reaction rate constant. In Eq. (18), a reactant, is converted to a final product , in four steps with the help of four additional species, and . We consider that and are in vast excess and thus can be modeled at constant concentration. The 2D rate equations becomes:
with the given initial conditions
where are the spatial coordinates, represent the diffusion coefficients, are constant concentrations, and represent the concentrations of reactant species ..
3.1 Data generation
We aim to learn the mapping from initial concentration to the evolved concentration , where (as shown in Figure 3(a)). The initial concentration is modeled as a Gaussian random field
where and are the mean and covariance functions respectively. For simplicity, we set , while the covariance matrix is given by the squared exponential kernel as
where , are the correlation length scales along the and spatial directions, respectively. To generate realizations of the input stochastic field, we employ the truncated Karhunen-Loéve expansion.
We consider two cases corresponding to different initial concentrations of reactant as illustrated in Figure 3. In Case I, the concentrations and both approach a stable equilibrium concentration. In Case II, the system reaches a limit cycle which causes periodic oscillations in the concentrations and . In both cases, we develop three sets of data: a training/test data set with specified random field parameters and two OOD data sets with different random field parameters. These parameter sets are presented in Table 1. Datasets OOD and OOD correspond to lower and higher length scales respectively, compared to the reference training data, to test the extrapolation accuracy of the trained models.
In all cases, the initial field is kept constant at . The diffusion coefficients are set equal to and . The simulation takes place in a square domain , discretized with grid points and solved with finite differences (FD) in the time interval for . To form the training data set, we collect solution snapshots at (Case I) and (Case II) equally spaced points in time. We denote the input random field with , where , and corresponding model outputs , where , , and is the number of training data.
3.2 Surrogate Approximations and their Parameterizations
For the data sets described in the previous sections, we develop surrogate models using the DeepONet, FNO, and m-PCE methods. In this section, we briefly discuss some important aspects of these surrogates.
For the DeepONet, we consider different architectures for the two cases. In Case I, we observe a smooth and decaying oscillation, which is encoded into the network simply by modifying the time input of the trunk network, as a trigonometric function, to obtain appropriate basis functions. However, for Case II, we observe that even though the function exhibits periodic oscillations, there are sudden spikes towards the end of each cycle. In this case, we use the proposed modified DeepONet with self-adaptivity (SA-DeepONet). The DeepONet implemented for this case uses a FNN architecture in both the branch and the trunk net. The spikes in the solution field observed toward the end of each cycle lead to a non-differentiable loss function, as seen in Figure 3, which necessitates self-adaptivity. In FNO, we have used a equispaced grid discretization for the spatial and the temporal domain as input to the integral-operator based framework.
For m-PCE surrogate modeling, we investigate both under-parameterized (U m-PCE) and over-parameterized (O m-PCE) versions of the model. In the U m-PCE, the total number of trainable PCE parameters is less than the number of available training data points (i.e. ) where
Meanwhile, for the O m-PCE we have . Note that depends on the PCE order as well as the dimensions of the input/output latent spaces. Thus, increasing any of these quantities will increase model complexity. We further note that we use standard PCE learning methods with full index sets and low-order () polynomials. We do not consider sparse PCE implementations that use, for example Least Angle Regression (4), to reduce the basis set and promote sparsity in the expansion. This is consistent with the DNN architecture we use where, likewise, no attempt is made to explore alternative or more compact architectures.
4 Numerical results
Here, we explore the performance of manifold PCE, DeepONet, and FNO on the Brusselator system. To evaluate the performance of the model, we compute the
relative error of predictions, and we report the mean and standard deviation of this metric based on five independent training trials. Below, we present the results for the two cases of the model discussed in Section3.1.
4.1 Case I
Initially, we train all models with data and test their accuracy on a test data sets of data. In Table 2, we show a brief description of each model (second column, while the DeepONet architectures are shown in Table 8). The total number of trainable parameters for all models is shown in the third column and the relative error on the test dataset in the fourth column. In addition, we show the relative error for the two out-of-distribution data sets (OOD and OOD from Table 1). Finally, we evaluate the model performance on a noisy dataset which consists of the original test dataset with added uncorrelated Gaussian noise (last column).
|Method||Description||of params||Relative error|
|Test data||OOD||OOD||noise data|
|DeepONet||(see Table 8)|
For the under-parameterized m-PCE surrogate, a relatively low dimensional embedding of the input, and output data, has been chosen to reduce the number of trainable parameters, which results in a considerably higher relative error compared to its over-parameterized counterpart. k-PCA111The python package scikit-learn (42) was used for the implementation of k-PCA.
has been used as the DR method, with a radial basis function (‘rbf’) and a polynomial (‘poly’) kernel for the dimensionality reduction of the input and output data, respectively. In both under- and over-parameterized models, the embedding dimension has been optimized manually using a grid search. Among all models, the DeepONet exhibits the best performance with the smallest relative error. However, we note that all models perform comparably and reasonably well on the test and OOD data, pointing to their generalizability for this case. We further observe that all models perform better on OOD, which has a larger lengthscale than the training data set, than they do for OOD. This follows intuitively because a shorter lengthscale for the input random field increases its intrinsic dimension and makes the mapping more complex. When added noise is introduced to the training data, the m-PCE and DeepONet results are consistent with their noise-free errors, while the FNO error increases substantially. We show the performance of both DeepONet and m-PCE for higher levels of noise ( and ) in Table 3. Moreover, to examine the stability of the non-linear mappings we learnt, we added noise to the inputs of the testing data and we show these results in Table 4. The results in both tables suggest that there is no amplification of inference error even for high levels of input noise.
To investigate further the influence of over/under parameterization in the considered model, we systematically study the m-PCE and DeepONet for increasing number of parameters. In Figure 4, we show the relative error for interpolation (test data) and extrapolation (OOD data) for the over-parameterized m-PCE and DeepONet as a function of the number of trainable model parameters. We observe that for m-PCE (left plot) the error is minimized for trainable parameters and subsequently increases as the number of parameters increases. Thus, in m-PCE identifying the optimal number of parameters is crucial. For the DeepONet (right plot), we observe that the testing error continuously decreases as the number of model parameters increases, which implies that simply increasing network size improves performance, albeit with diminishing return.
A closer look at the DeepONet architecture affirms the observation from Figure 4, with some interesting caveats. In Table 3, the model parameters in the DeepONet are modified by altering the width of the FNN in the trunk net and also the depth of the fully connected layers beyond the convolutional layers in the branch net. Overall, we see moving down the table that as network complexity (number of parameters) decreases, the error increases. However, for very large networks having convolutional layers with sizes , the error increases over the less complex networks with convolutional layers. However, using a regularizer can improve the accuracy of the network. When a regularization with is used, the accuracy rather slightly deteriorated. However, when , the accuracy of the network improves. This shows that with enormous over-parameterization, DNNs might require some appropriate regularization to avoid slight over-fitting.
|Branch net||Trunk net||of params||in L||Relative error ()|
|Model description||Relative error ()|
|Case I||U. m-PCE||,|
|Case II||U. m-PCE||,|
|DeepONet||(see Table 8)|
|POD-DeepONet||(see Table 8)|
|SA-DeepONet||(see Table 8)|
|Training data||Relative error|
|Under. m-PCE||Over. m-PCE||DeepONet||FNO|
Another way to look at over/under-parameterization of the models is by increasing the training data set size. As we increase the training data set size, the models become closer to being under-parameterized. In Table 5, the relative error is shown for five different training data set sizes. For better visualization, we also plot the relative error in Figure 5(a). We observe that, for all models, the error decreases as the training set size increases. However, the over-parameterized models (m-PCE, DeepONet, and FNO) have in general a comparable performance, while the under-parameterized m-PCE sees only minimal improvement. We further find that in cases where data are sparse (), m-PCE slightly outperforms the other models. However, for , DeepONet exhibits the best performance, with the lowest relative error for data, with the FNO showing similar accuracy. We further recognize that the over-parameterized m-PCE, in fact, becomes under-parameterized after grows beyond 22,320. The DeepONet and FNO, on the other hand, are still vastly over-parameterized even in the largest data set cases. These results show that the over-parameterized models have inbuilt flexibility to continue learning as data set size increases, while the under-parameterized models reach their training limit. However, the under-parameterized models are superior for small training datasets.
For illustration purposes, the relative errors for the U m-PCE, O m-PCE, and DeepONet are shown for 5 time snapshots from a given initial condition in Figure 6. Despite the high variability in the input fields and the complexity of the model solution, m-PCE and DeepONet are able to accurately predict the model response for all time steps with small local errors that are observed in areas where there is a rapid change in the model response.
Finally, in terms of computational cost, training the m-PCE is significantly faster. For data the cost of training the m-PCE (including both the manifold learning and surrogate construction) was on a GHz quad-core Intel i7 CPU. Meanwhile, the DeepONet took and the FNO for training, both trained on an NVIDIA A6000 GPU. For the DeepONet and FNO, the batch size is 200. We therefore see that the m-PCE training is and times faster than the DeepONet and FNO, respectively. Nonetheless, after the training procedure, all models are capable of making predictions in a fraction of a second.
4.2 Case II
Similar to Section 4.1, we test the model performance for training data for Case II. In this case, the system is enters a limit cycle and thus after sufficient time we observe periodic oscillations. The sharp drop in the model response, as shown in Section 3.1, is particularly challenging for the surrogates to learn. To address this challenge, we additionally employ POD-DeepONet, and the proposed SA-DeepONet. The predictive accuracy of all models are shown in Table 6. We observe that for Case II, the difference between m-PCE and the DeepONets is more significant than in Case I. SA-DeepONet, in particular, achieves a very small relative error and shows very good performance for both in- and out-of-distribution data, and also when Gaussian noise is added to the input random fields. Consistent with the results of the previous section, the over-parameterized m-PCE achieves significantly higher predictive accuracy than its under-parameterized counterpart, and the error does not significantly increase when tested on OOD and noisy data. On the other hand, the prediction accuracy of FNO is comparable for out-of-distribution datasets, however, for noisy inputs ( noise), the network fails to approximate reasonably, with a very high relative error of . The performance of all surrogates models for higher levels of noise ( and ) is shown in Table 4. Compared to DeepONet models, mPCE shows more significant deterioration in performance, however, within acceptable limits.
|Test data||OOD data||OOD data||noise data|
|DeepONet||(see Table 8)|
|POD-DeepONet||(see Table 8)|
|SA-DeepONet||(see Table 8)|
|# of data||Relative error|
|U. m-PCE||O. m-PCE||DeepONet||POD-DeepONet||SA-DeepONet||FNO|
The relative errors for Case II follow the same observations as made in Case I. The surrogate model approximates the field more accurately for OOD over OOD. The intrinsically non-smooth nature of the solution in this cases poses a fundamental problem in developing a surrogate model. Most surrogate models, such as m-PCE, typically demand smoothness in the system response, and therefore we observe higher relative errors in both the under and over-parameterized m-PCE. However, using the DeepONet, the locations of sharp features can be manually selected and penalized along with the regular loss function, which improves optimization of the model. The SA-DeepONet, meanwhile, makes this task easier by automatically selecting and updating the penalty parameters, to give the most accurate solution.
As displayed in Table 7 and also shown in Figure 5(b), we once again see that the dataset size affects significantly the results for all models. Similar to Case I, the error continuously decreases with increasing training set size for all models. However, in this case due to the non-smooth solution, the DeepONet provides superior performance even for small data sets and the SA-DeepONet shows the best overall performance. We further see that the DeepONet models converge very rapidly and do not show significant improvement for much larger data sets.
Once again for illustration, we show representative relative error plots for m-PCE, standard DeepONet, POD-DeepONet and SA-DeepONet for a given random realization of the input stochastic field in Figure 7. Both the under- and over-parameterized m-PCE provide a good approximation of the evolving field, however, when the sharp transition caused by temporal oscillations occurs (), the error grows substantially. The DeepONet models show a greater overall consistency between the prediction and ground truth, while they also handle more effectively the non-smoothness in the model response, leading to a smaller overall test error.
Although m-PCE performed satisfactorily in both dataset cases, we note that there are many possibilities for improvement of the method which were not explored here. Only the k-PCA method was employed in this study for lower-dimensional embedding of both input and output quantities, but the flexibility of m-PCE allows any suitable manifold learning technique to be employed. The choice of which method to use should be based on the complexity and nonlinearity of the model. Exploring other manifold learning methods may, in fact, improve the error for this problem significantly. In addition, improvements in the PCE model construction, e.g., via adaptive refinement strategies such as Least Angle Regression (4) or through adaptive refinement with multi-element PCE (55) could improve its treatment of non-smooth models. These warrant further investigation and are increasingly motivated by the fact that the m-PCE is significantly less computationally expensive, with training time on the order of seconds on a CPU versus minutes or hours for the DeepONet.
In this study, we have investigated and compared the performance of under-and over-parameterized manifold-based surrogates with over-parameterized deep neural operators, specifically DeepONet and FNO, for approximating complex, nonlinear, and high-dimensional models. To compare the surrogate models, we have considered the Brusselator reaction-diffusion system, which is modeled by a 2D time-dependent PDE, where the initial field of one of the species is modeled as a Gaussian random field. The surrogate models at hand are investigated in terms of relative error for a relatively smooth and a highly nonsmooth dataset, and tested on three regression tasks: interpolation, extrapolation, and robustness to noise. Our observations can be summarized as follows:
Over-parameterization of models lead to the construction of more expressive surrogates, achieving higher (often significantly) generalization accuracy than their under-parameterized counterparts (See Tables 2 and 6 fourth column). DeepONet is less sensitive to the total number of trainable parameters, as the extrapolation performance is not significantly affected when this number is varied. On the other hand, m-PCE is more sensitive to the total number of trainable parameters, whether this number is controlled by the input or output embedding dimensionality or the maximum polynomial degree. In cases where the Brusselator model exhibits a relatively smooth response (Case I), we found that the over-parameterized surrogate models in general show comparable performance, while achieving a very small generalization error for both in- and out-of-distribution inputs as well as noisy inputs (See Table 2).
When highly non-smooth dynamics are considered, DeepONet and its extensions (POD-DeepONet and weight self-adaptivity DeepONet) exhibited very high generalization accuracy. While all three DeepONet variants (standard, POD, and SA) performed very well in Case II in both interpolation and extrapolation, the SA-DeepONet outperformed the rest and is found to be very robust to noise. (See Table 6). On the other hand, m-PCE performed well overall but failed to approximate accurately the system response near the sharp transition. This challenge could be addressed by employing enhanced PCE approaches such as multi-element PCE (54).
In general, all surrogate models (under- and over-parameterized) generalize the solution for OOD data having length scales higher than the training data more accurately than for OOD data having lower length scales, which results in higher variability in the input fields and thus model response (See tables 2 and 6 fifth and sixth columns).
The robustness of the surrogate models is evaluated by testing the predictions with noisy inputs. Additionally, we observe the consistency of the surrogate models to make appropriate approximations under the circumstance of noisy input data. We interpolated smoothly between the cases of no noise and inducing , and noise to observe a very small deterioration of the generalization error as the noise level is increased. We found that both DeepONet and m-PCE surrogates can capture the remaining signal in the data while also fitting the noisy input data (see Tables 2, 3, 5,4), while FNO failed to perform well in such cases.
The effect of dataset size has significant impact on the accuracy of the over-parameterized surrogate models. For smooth response (Case I), the under-parameterized model is not significantly improved by the increase in dataset size, while for all the over-parameterized models, the accuracy improves as the training dataset size increases. DeepONet exhibits the highest level of accuracy for large datasets, however, when small datasets are considered, over-parameterized m-PCE outperforms DeepONet (see Figure 5). Considering non-smooth dynamics (Case II), under-parameterized models fail to capture the response accurately while the over-parameterized models exhibit similar behaviour as discussed for Case I (see Table 7).
DeepONet and its proposed extensions exhibit the highest generalization accuracy overall for varying degrees of complexity in model response and provide an expressive surrogate modeling method with the ability to learn the nonlinear operator, thus allowing for predicting model response for new unseen input functions, boundary conditions, domain geometries etc, at high training cost.
The m-PCE provides a powerful and flexible surrogate modeling approach, which combines high predictive accuracy and robustness to noise with an extremely low training computational cost compared to all other methods. This makes it particularly attractive for tasks involving active learning methods that iteratively retrain the surrogate. Finally, significant research potential exists to improve m-PCE methods to achieve better generalization for cases of non-smooth, high-dimensional models.
For KK and MDS, this material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research under Award Number DE-SC0020428. SG and GEK would like to acknowledge support by the DOE project PhILMs (Award Number DE-SC0019453) and the OSD/AFOSR MURI grant FA9550-20-1-0358.
-  (2019) Numerical modeling of three dimensional brusselator reaction diffusion system. AIP Advances 9 (1), pp. 015205. Cited by: §3.
Reconciling modern machine-learning practice and the classical bias–variance trade-off. Proceedings of the National Academy of Sciences 116 (32), pp. 15849–15854. Cited by: §1.
-  (2006) Pattern recognition and machine learning. springer. Cited by: §2.2.2.
-  (2011) Adaptive sparse polynomial chaos expansion based on least angle regression. Journal of Computational Physics 230 (6), pp. 2345–2367. Cited by: §2.2.1, §3.2, §4.2.
-  (2005) Functional analysis for probability and stochastic processes: an introduction. Cambridge University Press. Cited by: §2.2.1.
-  (2021) A universal law of robustness via isoperimetry. Advances in Neural Information Processing Systems 34. Cited by: §1.
-  (2013) Efficient basis change for sparse-grid interpolating polynomials with application to t-cell sensitivity analysis. Computational Biology Journal 2013. Cited by: §2.2.1.
Generative deep neural networks for inverse materials design using backpropagation and active learning. Advanced Science 7 (5), pp. 1902607. Cited by: §1.
Uncertainty propagation using infinite mixture of gaussian processes and variational bayesian inference. Journal of Computational Physics 284, pp. 291–333. Cited by: §1.
-  (2013) Adaptive smolyak pseudospectral approximations. SIAM Journal on Scientific Computing 35 (6), pp. A2643–A2670. Cited by: §2.2.1.
-  (2012) Sparse pseudospectral approximation method. Computer Methods in Applied Mechanics and Engineering 229, pp. 1–12. Cited by: §2.2.1.
-  (2021) A farewell to the bias-variance tradeoff? an overview of the theory of overparameterized machine learning. arXiv preprint arXiv:2109.02355. Cited by: §1.
-  (2021) DeepONet prediction of linear instability waves in high-speed boundary layers. arXiv preprint arXiv:2105.08697. Cited by: §1.
-  (2018) Sparse polynomial chaos expansions via compressed sensing and d-optimal design. Computer Methods in Applied Mechanics and Engineering 336, pp. 640–666. Cited by: §2.2.1.
-  (2018) On the power of over-parametrization in neural networks with quadratic activation. In International Conference on Machine Learning, pp. 1329–1338. Cited by: §1.
-  (1990) Polynomial chaos in stochastic finite elements. Journal of Applied Mechanics. Cited by: §1.
-  (2020) Data-driven surrogates for high dimensional models using gaussian process regression on the grassmann manifold. Computer Methods in Applied Mechanics and Engineering 370, pp. 113269. Cited by: §1.
-  (2013) Adaptive response surface method in structural response approximation under uncertainty. In International Conference on Structural Engineering and Mechanics, pp. 194–202. Cited by: §1.
-  (2016) Reliability analysis of structures by iterative improved response surface method. Structural Safety 60, pp. 56–66. Cited by: §1.
-  (2022) A physics-informed variational deeponet for predicting crack path in quasi-brittle materials. Computer Methods in Applied Mechanics and Engineering 391, pp. 114587. Cited by: §1, §2.1.
-  (2018) Least squares polynomial chaos expansion: a review of sampling strategies. Computer Methods in Applied Mechanics and Engineering 332, pp. 382–407. Cited by: §2.2.1.
-  (2018) Basis adaptive sample efficient polynomial chaos (base-pc). Journal of Computational Physics 371, pp. 20–49. Cited by: §2.2.1.
-  (2020) An adaptive polynomial chaos expansion for high-dimensional reliability analysis. Structural and Multidisciplinary Optimization 62 (4), pp. 2051–2067. Cited by: §2.2.1.
Kernel pca for novelty detection. Pattern Recognition 40 (3), pp. 863–874. Cited by: §2.2.2.
-  (2021) Deep neural networks for the evaluation and design of photonic devices. Nature Reviews Materials 6 (8), pp. 679–700. Cited by: §1.
-  (2021) Manifold learning-based polynomial chaos expansions for high-dimensional surrogate models. arXiv preprint arXiv:2107.09814. Cited by: §1, §2.2.2.
A survey of unsupervised learning methods for high-dimensional uncertainty quantification in black-box-type problems. arXiv preprint arXiv:2202.04648. Cited by: §1, §2.2.2, §2.2.2, §2.2, §2.2.
-  (2012) Imagenet classification with deep convolutional neural networks. Advances in neural information processing systems 25. Cited by: §1.
-  (2021) . arXiv preprint arXiv:2102.09618. Cited by: §1.
-  (2020) Extending classical surrogate modeling to high dimensions through supervised dimensionality reduction: a data-driven approach. International Journal for Uncertainty Quantification 10 (1). Cited by: §1, §2.2.2.
Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08895. Cited by: §1, §2.1.2, §2.1.2, §2.1.
-  (2021) Operator learning for predicting multiscale bubble growth dynamics. The Journal of Chemical Physics 154 (10), pp. 104118. Cited by: §2.1.
-  (2019) Adaptive sparse polynomial chaos expansions via leja interpolation. arXiv preprint arXiv:1911.08312. Cited by: §2.2.1.
-  (2020) Robust adaptive least squares polynomial chaos expansions in high-frequency applications. International Journal of Numerical Modelling: Electronic Networks, Devices and Fields 33 (6), pp. e2725. Cited by: §2.2.1.
-  (2021) Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature Machine Intelligence 3 (3), pp. 218–229. Cited by: §1, §1, §2.1.
-  (2021) A comprehensive and fair comparison of two neural operators (with practical extensions) based on fair data. arXiv preprint arXiv:2111.05512. Cited by: §1, §2.1.1, §2.1.
-  (2020) Self-adaptive physics-informed neural networks using a soft attention mechanism. arXiv preprint arXiv:2009.04544. Cited by: §2.1.1.
-  (2018) The role of over-parametrization in generalization of neural networks. In International Conference on Learning Representations, Cited by: §1.
-  (2012) Data-driven uncertainty quantification using the arbitrary polynomial chaos expansion. Reliability Engineering & System Safety 106, pp. 179–190. Cited by: §1.
-  (2020) UQpy: a general purpose python package and development environment for uncertainty quantification. Journal of Computational Science 47, pp. 101204. Cited by: §1.
-  (2021) Bayesian neural networks for uncertainty quantification in data-driven materials modeling. Computer Methods in Applied Mechanics and Engineering 386, pp. 114079. Cited by: §1.
-  (2011) Scikit-learn: machine learning in Python. Journal of Machine Learning Research 12, pp. 2825–2830. Cited by: footnote 1.
-  (2020) Theoretical issues in deep networks. Proceedings of the National Academy of Sciences 117 (48), pp. 30039–30045. Cited by: §1.
-  (1978) Time, structure, and fluctuations. Science 201 (4358), pp. 777–785. Cited by: §3.
-  (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378, pp. 686–707. Cited by: §1.
-  (2018) Numerical gaussian processes for time-dependent and nonlinear partial differential equations. SIAM Journal on Scientific Computing 40 (1), pp. A172–A198. Cited by: §1.
Notes on regularized least squares.
Technical Report MIT-CSAILTR-2007-025, Computer Science and Artificial Intelligence Laboratory, MIT. Cited by: §2.2.1.
-  (1997) Kernel principal component analysis. In International conference on artificial neural networks, pp. 583–588. Cited by: §2.2.2.
-  (2004) Physical systems with random uncertainties: chaos representations with arbitrary probability measure. SIAM Journal on Scientific Computing 26 (2), pp. 395–410. Cited by: §2.2.1.
-  (2020) Physics-informed deep neural networks for learning parameters and constitutive relationships in subsurface flow problems. Water Resources Research 56 (5), pp. e2019WR026731. Cited by: §1.
-  (2020) The computational limits of deep learning. arXiv preprint arXiv:2007.05558. Cited by: §1.
-  (2016) Gaussian processes with built-in dimensionality reduction: applications to high-dimensional uncertainty propagation. Journal of Computational Physics 321, pp. 191–223. Cited by: §1.
-  (2019) Compressive sensing adaptation for polynomial chaos expansions. Journal of Computational Physics 380, pp. 29–47. Cited by: §2.2.1.
-  (2005) An adaptive multi-element generalized polynomial chaos method for stochastic differential equations. Journal of Computational Physics 209 (2), pp. 617–642. Cited by: item 2.
-  (2006) Multi-element generalized polynomial chaos for arbitrary probability measures. SIAM Journal on Scientific Computing 28 (3), pp. 901–928. Cited by: §2.2.1, §4.2.
-  (2016) Sparse pseudo spectral projection methods with directional adaptation for uncertainty quantification. Journal of Scientific Computing 68 (2), pp. 596–623. Cited by: §2.2.1.
-  (2002) The wiener–askey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing 24 (2), pp. 619–644. Cited by: §1, §2.2.1.
-  (2015) Adaptive multi-element polynomial chaos with discrete measure: algorithms and application to spdes. Applied Numerical Mathematics 90, pp. 91–110. Cited by: §1.
Appendix A Model Architectures
|DeepONet (Case I)||CNN||ReLU|
|DeepONet (Case II)||CNN||ReLU|
|POD-DeepONet (Case II)||CNN||modes||ReLU|
|SA-DeepONet (Case II)||ReLU|