Physics-Consistency Condition for Infinite Neural Networks and Experimental Characterization ^†

Abstract

It has previously been shown that prior physics knowledge can be incorporated into the structure of an artificial neural network via neural activation functions based on (i) the correspondence under the infinite-width limit between neural networks and Gaussian processes if the central limit theorem holds and (ii) the construction of physics-consistent Gaussian process kernels, i.e., specialized covariance functions that ensure that the Gaussian process fulfills a priori some linear (differential) equation. Such regression models can be useful in many-query problems, e.g., inverse problems, uncertainty quantification or optimization, when a single forward solution or likelihood evaluation is costly. Based on a small set of training data, the learned model or “surrogate” can then be used as a fast approximator. The bottleneck is then for the surrogate to also learn efficiently and effectively from small data sets while at the same time ensuring physically consistent predictions. Based on this, we will further explore the properties of so-constructed neural networks. In particular, we will characterize (i) generalization behavior and (ii) the approximation quality or Gaussianity as a function of network width and discuss (iii) extensions from shallow to deep NNs.

1. Introduction

The impressive achievements of neural networks (NNs) in computer science have sparked a growing interest also in physics and scientific computing. Here, we explore their utilization as surrogate models in computationally intensive tasks like simulation-based optimization or uncertainty quantification [1,2]. For this, the surrogates’ predictions should be consistent with physics, which typically cannot be achieved with data alone as such models often generalize poorly, i.e., they learn the physics implicitly and well only in regimes where (enough) data are available.

Here, we pursue a previously proposed approach [3] as follows: (i) linear differential equations can be inherently integrated into a Gaussian process (GP) model by carefully selecting the covariance function or kernel [4,5,6,7,8]. The GP defined by the so-constructed kernel may be referred to as a “physics-consistent GP” [6]. (ii) Neal [9] established that Bayesian NNs converge toward a GP in an infinite-width limit using the central limit theorem. By combining these two insights, we infer that at least some laws of physics can also be a priori embedded into a suitable infinite NN model by choosing neural architecture and activation functions in alignment with the corresponding GP’s kernel, i.e., the corresponding infinite NN has the same physics-consistency property as the physics-consistent GP. Note that this approach fundamentally differs from “physics-informed” learning machines [10], where the residual of the equation (of the physical law) is included into the loss function independent of NN architecture. Departing from the foundations introduced in [3], we here consider generalization outside of the training data regime, and Gaussianity, i.e., convergence of the infinite-width limit approximation.

2. Background

This section provides a concise overview of the relevant literature on Gaussian processes (GPs), neural networks (NNs), and the infinite-width correspondence between the two. Additionally, related works on physically inspired GPs and NNs are discussed.

2.1. Physics-Consistent Gaussian Processes

Gaussian process regression, also known as “Kriging”, has been extensively studied for several decades [11]. GPs have become one of the most popular classes of machine learning [12] and Bayesian models [13], and found wide application in the field of “Uncertainty Quantification” [14,15]. The defining element of a Gaussian process is its covariance function, also referred to as the kernel. Consequently, significant efforts have been dedicated to the design and selection of appropriate kernel functions [16].

As a matter of fact, a linear (differential) operator applied to a GP results in another GP with a modified kernel [5,6,7,17]. This allows one to derive a “physics-consistency condition” for the kernel, allowing the construction of kernels to define a GP that satisfies physical laws expressed in terms of said differential equations a priori, i.e., before training. Examples include divergence-free fields [18] and the Helmholtz equation [6]. More general linear constraints and boundary conditions can be considered in an analogous manner [17,19,20,21]. Note that it is also possible to construct physics-based or “weakly” physics-informed GPs through training [22] given some (non necessarily physical) base kernel, or through Green’s functions [23,24].

While this approach was generally believed to be limited to linear equations, it has recently been extended even toward non-linear partial differential equations (PDEs) [25].

2.2. Infinite Neural Networks

Neural networks (NNs) are arguably one of the most popular classes of machine learning models. Similar to GPs and kernels, the design and selection of activation functions play a crucial role in NNs. Numerous activation functions have been proposed [26,27], but few principled approaches are known on how to choose an appropriate neural activation.

The infinite-width correspondence relates to the fact that Bayesian NNs converge toward a GP under the infinite-width limit [9]. This infinite-width correspondence has been used to compute kernels for GPs corresponding to specific neural activations and architectures [28,29,30,31,32], which may be regarded as a form of the “kernel trick”. However, the inverse approach, i.e., computing activations from kernels, has received little attention. Here, we aim to take this route due to the beneficial circumstance that we can make GPs physics-consistent. Note that the term “physics-informed neural network” relates to training based on a physics-loss term composed of the differential equation residual, and was popularized by Raissi et al. [33], although the idea was around earlier [34,35].

Quite importantly, it has recently been shown that a very large class of NN architectures converge toward a GP [36], including deep NNs [29], convolutional NNs [37], recurrent NNs [38] and transformers [39]. Furthermore, the infinite-width limit provides valuable insight into the training dynamics of NNs as well [40]. Note that random Fourier features [41] also share the underlying mechanism of the central limit theorem. As for interactions arising from the breaking of CLT through either a finite number of neurons or non-independence of neurons, i.e., the emergence of a non-Gaussian field, see [42].

3. Infinite Neural Networks Converging to Physics-Consistent Gaussian Processes

In this section, we give a brief review of the notions introduced previously in [3]. We introduce two sets of functions (function spaces) $V,W:X\to Y\subseteq {\mathbb{R}}^Q$. From the first set V, we denote two particular elements (i.e., functions from this space) as $g\left(x\right)\in V,x\mapsto g\left(x\right)$ and $f\left(x\right)\in V,x\mapsto f\left(x\right)$, $x\in X\subseteq {\mathbb{R}}^D$, where $D,Q$ denote the input/output dimension. We will later see that GPs are distributions for functions; hence, we need sets of functions accordingly. We aim to learn these functions from data related to some physical law expressed in terms of a linear operator $\widehat{\mathcal{O}}:V\to W$, $\widehat{\mathcal{O}}\circ f\left(x\right)\mapsto \tilde{f}\left(x\right)$. In words, an operator maps a function from the set (function space) V to another function from the set (function space) W. Linearity requires $\widehat{\mathcal{O}}\circ (\alpha ·f+\beta ·g)=\alpha ·\widehat{\mathcal{O}}\circ f+\beta ·\widehat{\mathcal{O}}\circ g$ for some constants $\alpha ,\beta \in \mathbb{R}$. Let the physical law be expressed in the form $\widehat{\mathcal{O}}f\left(x\right)=0$, where, from now on, we omit the ∘-sign for simpler notation, e.g., if $\widehat{\mathcal{O}}$ were the D-dimensional Laplace-operator $\widehat{\mathcal{O}}\triangleq {\nabla }^2:={\sum }_{k=1}^D\frac{{\partial }^2}{\partial x_k^2}$, $x_k$ denoting the $k-th$ element of D-dimensional vector $x={(x_1,\dots ,x_D)}^T$, then V would be at least the set of (continuous) functions twice-differentiable and $\widehat{\mathcal{O}}:C^2\left(X\right)\to C^0\left(X\right)\supset C^2\left(X\right)$, i.e., $W\subset V$.

In the following, g will denote a GP and f will denote an NN.

3.1. Physics-Consistent Gaussian Processes

We will introduce GPs along the lines of [7]. A GP is a distribution on the set of functions ${\mathbb{R}}^D\to {\mathbb{R}}^Q$ such that the function values $g\left(x^{\left(1\right)}\right),\cdots ,g\left(x^{\left(n\right)}\right)$ at $x^{\left(i\right)}\in {\mathbb{R}}^D$ have a joint Gaussian distribution. Note the consistent marginalization property that follows from this definition. The distribution is specified by a mean function $\mu :{\mathbb{R}}^D\to {\mathbb{R}}^Q:x\mapsto \langle g\left(x\right)\rangle$ and a positive semi-definite covariance function $k:{\mathbb{R}}^D\oplus {\mathbb{R}}^D\mapsto {\mathbb{R}}_{⪰0}^{Q\times Q}:(x,x^{\prime })\mapsto 〈\left(g\left(x\right)-\mu \left(x\right)\right){\left(g\left(x^{\prime }\right)-\mu \left(x^{\prime }\right)\right)}^T〉$, where $\langle ·\rangle$ denotes expectations. k is also known as kernel, and the common notation is $g\left(x\right)\sim \mathcal{GP}\left(\mu \left(x\right),k(x,x^{\prime })\right)$. For simplicity, we assume g to be a zero-mean GP,

$$\begin{array}{c}g\left(x\right)\sim \mathcal{GP}\left(0,k(x,x^{\prime })\right)\hfill \end{array}$$

(1)

We further require that g be consistent with a physical law expressed in terms of a linear differential operator ${\widehat{\mathcal{O}}}_x$, i.e., differential with regard to x, i.e., that it admits the following equation:

$$\begin{array}{c}{\widehat{\mathcal{O}}}_{x}g\left(x\right)=0,\hfill \end{array}$$

(2)

In other words, we use a GP as an ansatz for solving Equation (2). As mentioned before, this is a differential operator defining a differential equation. Then, due to the linearity of $\widehat{\mathcal{O}}$,

$$\begin{array}{c}{\widehat{\mathcal{O}}}_{x}g\left(x\right)\sim \mathcal{GP}\left(0,{\widehat{\mathcal{O}}}_x{\widehat{\mathcal{O}}}_{x^{\prime }}^{†}k(x,x^{\prime })\right)\hfill \end{array}$$

(3)

where † denotes the adjoint, acting on the second argument of the kernel function $x^{\prime }$. We have introduced the subscript x to clarify which argument of the kernel the differential operator is acting on. At this point, we restrict ourselves to linear operators, but recent extensions to non-linear differential PDEs have been fairly successfully [25]. The above relationship can be used to formulate a number of equivalent physics-consistency conditions for the GP (or its kernel, respectively), e.g., from [5]:

$$\begin{array}{cc}{\widehat{\mathcal{O}}}_x{\widehat{\mathcal{O}}}_{x^{\prime }}k(x,x^{\prime }){|}_{x=x^{\prime }}\hfill & \stackrel{!}{=}0.\hfill \end{array}$$

(4)

The solution for k in Equation (4) then defines a GP that satisfies Equation (3) a priori, i.e., before training. In other words, these physics-consistency conditions can be used to find and design physics-consistent kernels. For this, Mercer’s theorem can be especially useful, in that it states that a kernel function (i.e., an element from a reproducing kernel Hilbert space) can always be expressed in terms of a set of basis functions $\left\{{\varphi }_i\right\}$ and vice versa [43]. Formally, this means that

$$\begin{array}{cc}k(x,x^{\prime })\hfill & =\underset{P\to \infty }{lim}\sum _{i,j}^P{\varphi }_i\left(x\right)M_{ij}{\varphi }_j\left(x^{\prime }\right)\hfill \end{array}$$

(5)

with a suitable prior covariance matrix M for the function weights. Those basis functions could, e.g., be constructed from Green’s functions as proposed by [6]. Exploiting linearity and requiring uniform convergence, the substitution of Equation (5) into Equation (4) would lead the operators to act on the basis functions in the sum directly. Note that [7] introduced a symbolic, algebraic algorithm program for constructing such kernels.

3.2. Infinite Neural Networks

Let us now consider an NN f with a single hidden layer with N neurons with non-linear but bounded activation functions h,

$$\begin{array}{c}f\left(x\right)=\sum _{k=1}^N{\mathrm{v}}_k{h}_k\left({\mathrm{w}}_k·x+b_k\right),\hfill \end{array}$$

(6)

where the index $k=1,\dots ,N$ sums over the N neurons in the single hidden layer. ${\mathrm{w}}_k$ are the input-to-latent weight matrices, b are bias terms and ${\mathrm{v}}_k$ are the latent-to-output weights. Under the infinite-width limit, $N\to \infty$, Equation (6) converges to a Gaussian process, i.e., a Gaussian distribution $\forall x$, using the central limit theorem if the activation functions $h_k$ are bounded and weights and biases are i.i.d. [9]. Figure 1 shows an illustration of the correspondence between GPs and infinite NNs.

The same finding holds for a considerably larger class of models [36]. Then, the NN converges to a GP with kernel $k(x,x^{\prime })=〈f\left(x\right)f\left(x^{\prime }\right)〉$, i.e., the covariance of NN predictions must satisfy the covariance of the GP, which is a sensible requirement. This gives rise to the following condition for equivalence between a certain GP and a certain NN:

$$\begin{array}{cc}k(x,x^{\prime })\hfill & \stackrel{!}{=}〈f\left(x\right)f\left(x^{\prime }\right)〉\hfill \\ \hfill & \propto \underset{N\to \infty }{lim}\int \sum _{k=1}^N\sum _{k^{\prime }=1}^N{\mathrm{v}}_k{\mathrm{v}}_{k^{\prime }}{h}_k({\mathrm{w}}_{k}x+b_k)h_{k^{\prime }}({\mathrm{w}}_{k^{\prime }}{x}^{\prime }+b_{k^{\prime }})p\left(\theta \right)d{V}_{\theta }\hfill \end{array}$$

(7)

where we used the linearity of the expectation and $p\left(\theta \right)$ is the probability density function of the weights and biases $\theta =\{\mathrm{v},\mathrm{w},b\}$. Note that we exchanged limit and integration, which is (typically) permissible if convergence is uniform and the limit function is integrable [44]. This integral can be solved numerically and, in some cases, analytically [28]. Note that Equation (7) is essentially determined by the covariance of the activations.

3.3. Physics-Consistency Condition for Neural Networks

We could now take the route of finding a physics-consistent kernel k on the left-hand side of Equation (7) and then try to find “physics-consistent” activations, i.e., activations such that the infinite NN becomes a physics-consistent GP. This approach may be regarded as an “inverse kernel trick”, especially if the Mercer representation is used.

An alternative route would be to find physics-consistent activations directly. Substituting the NN-GP consistency condition Equation (7) into the physics-consistency condition for the GP Equation (4) (or variants thereof) allows us to formulate a physics-consistency condition for the NN directly. Reciting uniform convergence and integrable limits as above used for Equation (7) [44], we switch the order of limit, expectation and differentiation:

$$\begin{array}{c}\left(\underset{N\to \infty }{lim}\int p\left(\theta \right)d{V}_{\theta }\sum _{k,k^{\prime }=1}^N{\mathrm{v}}_k{\mathrm{v}}_{k^{\prime }}\left[\widehat{\mathcal{O}}{h}_k\right]({\mathrm{w}}_{k}x+b_k)·\left[\widehat{\mathcal{O}}{h}_{k^{\prime }}\right]({\mathrm{w}}_{k^{\prime }}{x}^{\prime }+b_{k^{\prime }})\right){|}_{x=x^{\prime }}\stackrel{!}{=}0.\hfill \end{array}$$

(8)

where the chain rule has to be applied, giving rise to differentiated activations $x\mapsto \left[\widehat{\mathcal{O}}{h}_k\right]\left(x\right)$ applied to their original argument $({\mathrm{w}}_{k}x+b_k)$, where, if $h_k\left(x\right)\in V$, then $\left[\widehat{\mathcal{O}}{h}_k\right]\in W$. This gives rise to the problem of finding pairs of probability density functions $p\left(\theta \right)$ and activations h such that the physics-consistency condition for the NN Equation (8) is satisfied. This implies that the choice of activation and prior cannot be simply treated separately. The choice of prior $p\left(\theta \right)$ and activation h is then reminiscent of the choice of the prior covariance matrix M and basis functions $\varphi$ in the Mercer representation of the GP kernel in Equation (5).

Through equating Equation (8) l.h.s. with the kernel’s Mercer representation, Equation (5), and recalling that Green’s functions can admit kernels (i.e., applying the kernel trick to Green’s functions), it is easy to see through the choices $M_{ij}={\sigma }_{{\mathrm{v}}_i}{\delta }_{ij}$, Gaussian priors and the limit $p\left({\mathrm{w}}_k\right)={\displaystyle \underset{{\sigma }_{{\mathrm{w}}_k}\to 0}{lim}}\mathcal{N}(0,{\sigma }_{{\mathrm{w}}_k}^2\Sigma )$ (where $\Sigma$ is diagonal) that $h\to \varphi$. This means that the activations become the Green’s functions [3].

3.4. Training and Regularization with Physics-Information

Given pairs of input–output data $D={\{x^{\left(d\right)},y^{\left(d\right)}\}}_{d=1}^{N_d}$, where $y^{\left(d\right)}=f\left(x^{\left(d\right)}\right)+\epsilon$ are noisy observations with noise $\epsilon$, we ought to choose an optimization criterion or loss function for the learning, e.g., here,

$$\begin{array}{cc}\hfill & {\theta }^{*}=arg\underset{\theta }{min}\left[\sum _{d=1}^{N_d}{\left(y^{\left(d\right)}-f\left(x^{\left(d\right)}\right)\right)}^2+\lambda \sum _{p=1}^{N_p}∥\left({\widehat{\mathcal{O}}}_{x^{}}f\left(x^{}\right)\right){|}_{x=x^{\left(p\right)}}∥\right],\hfill \end{array}$$

(9)

which essentially consists of the likelihood and a prior defined via the residual of Equation (3). The last term denotes the differential to be applied to f and then evaluated at pivot points $x=x^{\left(p\right)}$ that may be distinct from the data points $x^{\left(d\right)}$. In more detail, for the second term (the “prior”), a network f is substituted as an ansatz function into Equation (3) and evaluated at pivot points (or sometimes called collocation points) $x^{\left(p\right)}$. This evaluation generally yields non-zero values for a generic NN (as opposed to the zero right-hand side in Equation (3), i.e., the equation at this point is not actually satisfied). The sum of these (absolute) values may be called the residual. In the physics-informed training approach, this residual is to be minimized by adjusting the function parameters (in an NN, the weights and biases). This minimization can be performed conjointly with minimizing the data misfit (i.e., the first term in Equation (9)), where the trade-off between the physics-residual and the data misfit term is controlled by a tunable regularization parameter $\lambda$. The choice of $\lambda$ and the set ${\left\{x^{\left(p\right)}\right\}}_{p=1}^{N_p}$ is generally non-trivial [45].

4. Experiments

4.1. Generalization

We consider a Helmholtz equation

$$({\nabla }^2-c^2)f\left(x\right)=0$$

(10)

where $c=6.0$ is a constant parameter related to the wavenumber. In the one-dimensional case, the solution to this equation is a sinusoidal function. The physics-consistent GP kernel for this physical law in 1D has the form $k(x,x^{\prime })=cos\left(\alpha (x-x^{\prime })\right)$ [46]. As derived in [3], one particular activation that defines an infinite NN with the same physics-consistent covariance is a sinusoidal activation function, $h(·):=sin(·)$.

We used the function $y=Asin(\omega x+\varphi )$ to generate the training data. A noise of normal distribution $\mathcal{N}(0,0.{04}^2)$ was added to the training sets. Next, we tried to learn an NN function $f\left(x\right)$ from this data, and compared the performance of three networks: all three networks are single-layer with the structure introduced before with $N=100$ neurons in the hidden layer. We denote as “vanilla” a simple NN with an $\mathrm{ReLU}(·)$ activation, and training with $L2$-loss. We denote as “physics-informed” the same network where the training is additionally regularized through a physics-loss term as in Equation (9). For the physics-informed neural network, the loss function specifically is $\mathcal{L}={\mathcal{L}}_d+{\mathcal{L}}_p={|y-\widehat{y}|}^2+\lambda |{\nabla }^{2}y+c^{2}y|$. Note that in the physics-consistent NN, there are no training difficulties due to higher-order derivatives and the training is as straightforward and fast as with the vanilla NN.

Note the distinction between training error, validation error and testing error. The training error is the prediction error on the training sets. Additional validation sets are used for early stopping to avoid overfitting, where training is stopped if there is no improvement in validation errors for a number of epochs. The validation sets are distributed within the same range/domain as the training sets, and test sets are distributed in a different range. The test sets are typically used to investigate the generalization behavior of an NN.

One reason for the wide-spread popularity of physics-informed NNs (PINNs) is that they typically generalize better than NNs trained on the data alone. PINNs achieve this by introducing a weak constraint that penalizes un-physical behavior. In our physics-consistent approach, in contrast, we enforce a hard constraint instead. In particular, this avoids the inconsistency of the multiple training objectives or the Pareto front as discussed in [47].

We compare the generalization behavior of the three NNs in Figure 2.

4.2. Convergence to Gaussian as a Function of Network Width

The Gaussianity of an NN is only exact under the infinite-width limit, whereas, on a computer, NNs are always finite. This begs the question of how good the approximation is as a function of NN width, or how fast the NN converges to a GP. Here, we use a single-hidden-layer neural network with prior distributions all i.i.d. normal with zero-mean. The variances are ${\sigma }_{\mathrm{w}}=5$, ${\sigma }_b=5$, ${\sigma }_{\mathrm{v}}=1$, ${\sigma }_a=0.1$.

For pedagogical illustration, we show scatter plots in Figure 3 of predictions collected from simple Monte Carlo samples from the prior on the weights and biases, $p\left(\theta \right)$ similar to [9] with a single-layer NN and $\tanh (·)$-activations. $NN(-0.2)$ and $NN\left(0.4\right)$ denote the outputs of the network corresponding to the inputs $x=-0.2$ and $x=0.4$, respectively.

The NN’s limit kernel is defined by Equation (7). For some NNs, this kernel can be computed analytically, while more often it has to be computed numerically, e.g., through Monte Carlo integration. For a single-layer NN with ReLU-activation, we compare the analytical kernel to the estimated kernel in Figure 4, which shows that such estimates can converge quickly and easily with only few neurons and samples.

This, however, is only an estimate of the second moments, and does not tell “how Gaussian” the distribution itself is. One suitable measure for this would be the Kullback–Leibler divergence, measuring the (dis-)similarity between two probability distributions. By defining a Gaussian distribution $Q=\mathcal{N}(0,K_Q)$ defined by the computed or estimated covariance matrix K, e.g., from Figure 4, the Kullback–Leibler divergence allows us to measure the dissimilarity between the distribution of NN predictions $P\left(\left(f\right(x)\right)$, i.e., the histograms of NN realizations collected from prior samples fed through the NN as in the scatter plots in Figure 3, and the (Gaussian) limit distribution Q. Note that alternative measures would be possible, such as the Wasserstein distance.

The Kullback–Leibler divergence is defined as:

$${\mathcal{D}}_{KL}(P\parallel Q)=\sum _{x\in \mathcal{X}}P\left(f\right)log\left(\frac{P\left(f\right)}{Q\left(f\right)}\right)$$

(11)

which can be estimated, e.g., through Monte Carlo. We plot the KL-divergence between the finite NN and its infinite-limit Gaussian as a function of the finite network’s width in Figure 5. The result confirms the intuition from Figure 3, and the finite NN behaves similar to a Gaussian already with just ${10}^2$ neurons.

5. Conclusions

We have derived a physics-consistency condition for NNs by exploiting the physics-consistency condition for GPs along with infinite NNs converging to GPs using the central limit theorem. We quantify experimentally that finite NNs with only order ∼10${}^2$ neurons already behave like GPs, i.e., that the approximation is sensible. We further show that engineering NNs to be physics-consistent can have advantages over both purely data-driven approaches and physics-based approaches with a similar intent, e.g., PINNs, in terms of convergence rate and generalization.