Enhancing Neural Network Differential Equation Solvers
We motivate the use of neural networks for the construction of numerical
solutions to differential equations. We prove that there exists a feed-forward
neural network that can arbitrarily minimise an objective function that is zero
at the solution of Poisson's equation, allowing us to guarantee that neural
network solution estimates can get arbitrarily close to the exact solutions. We
also show how these estimates can be appreciably enhanced through various
strategies, in particular through the construction of error correction
networks, for which we propose a general method. We conclude by providing
numerical experiments that attest to the validity of all such strategies for
variants of Poisson's equation.
Matthew J. H. Wright
# Enhancing Neural Network Differential Equation Solvers

###### Abstract

We motivate the use of neural networks for the construction of numerical solutions to differential equations. We prove that there exists a feed-forward neural network that can arbitrarily minimise an objective function that is zero at the solution of Poisson's equation, allowing us to guarantee that neural network solution estimates can get arbitrarily close to the exact solutions. We also show how these estimates can be appreciably enhanced through various strategies, in particular through the construction of error correction networks, for which we propose a general method. We conclude by providing numerical experiments that attest to the validity of all such strategies for variants of Poisson's equation. The source code for this project can be found at [https://github.com/mjhwright/error-correction](https://github.com/mjhwright/error-correction).

## 1 Introduction

Differential equations are among the most ubiquitous problems in contemporary mathematics. In recent years, developments in artificial neural networks have prompted new research into their capacity to be excellent differential equation solvers [1, 2, 3, 4, 5]. They are universal approximators [6]; they can circumvent the curse of dimensionality [7]; and they are continuous. However, practically, their construction and optimisation costs are enough to deter the discerning user.

In this paper, we explain a method by which neural networks can numerically solve differential equations. We further this by providing three strategies that can be targeted to improve the efficacy of the solver. The first two - sinusoidal representation networks [8] and random Fourier features [9] - are well-established in the field of artificial neural networks and machine learning. The third is a novel technique called error correction [10, 11, 12, 13, 14]. We explain how error correction can be implemented recursively, with little modification to the original solver, to give enhanced numerical solutions to differential equations, and we present results that demonstrate this.

This paper is designed to give a flavour of the competence of artificial neural networks in this field, while also highlighting their certain limitations.

## 2 Background

Throughout this paper, we consider differential equations with solution \(\phi:\mathbb{R}^{d}\to\mathbb{R}\). Consequently, our neural network approximation is a function \(\mathcal{N}:\mathbb{R}^{d}\to\mathbb{R}\).

### Universal approximation theorems

The realisation of neural networks' capabilities to learn seemingly any function has brought about numerous universal approximation theorems. These state that, under certain conditions, neural networks are able to approximate any function to arbitrary closeness. We recall one of these theorems by Hornik [6].

First, define the set of all functions represented by a neural network with a single hidden layer of width \(n\) and identity activation on the output layer as
\[\mathscr{A}^{n}(\sigma)=\left\{\mathcal{N}:\mathbb{R}^{d}\to\mathbb{R},\mathcal{N}( \mathbf{x})=\mathbf{W}^{(1)}\left(\sigma\left(\mathbf{W}^{(0)}\mathbf{x}+ \mathbf{b}^{(0)}\right)\right)+b^{(1)}\right\}\]

where \(\mathbf{x}\in\mathbb{R}^{d},\mathbf{W}^{(0)}\in\mathbb{R}^{n\times d},\mathbf{ W}^{(1)}\in\mathbb{R}^{1\times n},\mathbf{b}^{(0)}\in\mathbb{R}^{n},b^{(1)}\in \mathbb{R}\), and \(\sigma:\mathbb{R}\to\mathbb{R}\) is applied element-wise. Then,

\[\mathscr{A}(\sigma)=\bigcup_{n=1}^{\infty}\mathscr{A}^{n}(\sigma) \tag{1}\]

is the set of all such functions with any number of neurons. Define also \(\mathcal{C}^{m}(\mathbb{R}^{d})\) as the space of all functions that, together with their partial derivatives of order \(|\alpha|\leq m\), are continuous on \(\mathbb{R}^{d}\).

**Theorem 1**.: [6] _If \(\sigma\in\mathcal{C}^{m}(\mathbb{R}^{d})\) is nonconstant and bounded, then \(\mathscr{A}(\sigma)\) is uniformly \(m\)-dense on all compact sets of \(\mathcal{C}^{m}(\mathbb{R}^{d})\), i.e. for all \(\phi\in\mathcal{C}^{m}(\mathbb{R}^{d})\), for all compact sets \(\Omega\subset\mathbb{R}^{d}\) and for all \(\epsilon>0\), there exists \(\mathcal{N}\in\mathscr{A}(\sigma)\) such that_

\[\max_{|\alpha|\leq m}\sup_{x\in\Omega}|\partial_{x}^{(\alpha)}\mathcal{N}(x)- \partial_{x}^{(\alpha)}\phi(x)|<\epsilon\]

Theorem 1 illustrates the universal approximation quality for single-layer networks of arbitrary width. Applying results from an earlier paper by Hornik et al. [15], this can be extended to multilayer networks. Crucially, these theorems tell us that neural networks are dense on certain function spaces, but they do not tell us how to train a network to realise this.

### Neural network differential equation solvers

Using neural networks to solve differential equations was introduced in the late 1990s [1], but experienced a modern resurgence through the publication of two papers [2, 3] on physics-informed neural networks. The deep Galerkin method [4] which we describe below is very similar to the method described in [2] only, instead of using experimental data, we train a network on points randomly sampled across the domain of the differential equation.

Consider Poisson's equation with Dirichlet boundary conditions:

\[\begin{cases}\nabla^{2}\phi&=f\text{ in }\Omega\\ \phi&=g\text{ on }\partial\Omega\end{cases} \tag{2}\]

**Lemma 2**.: _Let \(\Omega\subset\mathbb{R}^{d}\) be a smooth, compact domain. Then there exists at most one solution \(\phi\) to (2)._

Proof.: Suppose \(\phi\) and \(\varphi\) both satisfy the conditions of (2) and let \(\omega=\phi-\varphi\). Then \(\omega\) is harmonic in \(\Omega\) and zero on \(\partial\Omega\). Then,

\[\int_{\Omega}\omega(\mathbf{x})\nabla^{2}\omega(\mathbf{x})\,d\mathbf{x}=\int _{\partial\Omega}\omega(\mathbf{x})\delta_{\mathbf{n}}\omega(\mathbf{x})\,d \mathbf{x}-\int_{\Omega}||\nabla\omega(\mathbf{x})||^{2}\,d\mathbf{x}=-\int_{ \Omega}||\nabla\omega(\mathbf{x})||^{2}\,d\mathbf{x}=0\]

and \(\nabla\omega=0\). Thus, \(\omega=0\) and \(\phi=\varphi\). 

We now seek an approximation \(\mathcal{N}\) to \(\phi\). Define the objective function

\[\mathcal{J}(\mathcal{N})=\int_{\Omega}|\nabla^{2}\mathcal{N}(\mathbf{x})-f( \mathbf{x})|^{2}\nu_{1}(\mathbf{x})\,d\mathbf{x}+\int_{\partial\Omega}| \mathcal{N}(\mathbf{x})-g(\mathbf{x})|^{2}\nu_{2}(\mathbf{x})\,d\mathbf{x}\]
for probability distributions \(\nu_{1}\) on \(\Omega\) and \(\nu_{2}\) on \(\partial\Omega\). By uniqueness of \(\phi\), \(\mathcal{J}(\mathcal{N})=0\implies\mathcal{N}=\phi\). However, minimising the objective function directly is impractical. First, we transform the problem into a machine learning framework. Our approximation \(\mathcal{N}=\mathcal{N}(\cdot;\theta)\) becomes a neural network with parameters \(\theta\).

#### Deep Galerkin method

We demonstrate the algorithm for the deep Galerkin method [4] when applied to Poisson's equation (2):

1. Randomly sample points \(\{\mathbf{x}_{i}\}_{i=1}^{M}\) from \(\Omega\) and \(\{\mathbf{y}_{j}\}_{j=1}^{N}\) from \(\partial\Omega\) according to respective probability distributions \(\nu_{1}\) and \(\nu_{2}\), and propagate them through a feed-forward neural network \(\mathcal{N}(\cdot;\theta)\).
2. Calculate the loss: \[\mathcal{L}(\theta)=\frac{1}{M}\sum_{i=1}^{M}\left(\nabla^{2}\mathcal{N}( \mathbf{x}_{i};\theta)-f(\mathbf{x}_{i})\right)^{2}+\frac{1}{N}\sum_{j=1}^{N }\left(\mathcal{N}(\mathbf{y}_{j};\theta)-g(\mathbf{y}_{j})\right)^{2}\]
3. Update parameters \(\theta_{t+1}=\theta_{t}-\eta\nabla_{\theta}\mathcal{L}(\theta_{t})\) with learning rate \(\eta>0\) and \(t\in\mathbb{N}_{0}\).
4. Repeat until \(\nabla_{\theta}\mathcal{L}(\theta_{t})\approx 0\).

This is a minibatch gradient descent implementation, where \(M\) and \(N\) are the size of the minibatches and \(M>N\).

**Lemma 3**.: \(\mathbb{E}[\nabla_{\theta}\mathcal{L}(\theta_{t})|\theta_{t}]=\nabla_{\theta} \mathcal{J}(\mathcal{N}(\cdot;\theta_{t}))\)__

Proof.: Assume \(\mathcal{L}\) sufficiently smooth and bounded to interchange derivatives and integrals. Then,

\[\mathbb{E}[\nabla_{\theta}\mathcal{L}(\theta_{t})|\theta_{t}] =\nabla_{\theta}\left[\frac{1}{M}\sum_{i=1}^{M}\mathbb{E}\left[( \nabla^{2}\mathcal{N}(\mathbf{x}_{i};\theta_{t})-f(\mathbf{x}_{i}))^{2}\right] +\frac{1}{N}\sum_{j=1}^{N}\mathbb{E}\left[(\mathcal{N}(\mathbf{y}_{j};\theta_{ t})-g(\mathbf{y}_{j}))^{2}\right]\right]\] \[=\nabla_{\theta}\left[\frac{1}{M}\sum_{i=1}^{M}\int_{\Omega}( \nabla^{2}\mathcal{N}(\mathbf{x};\theta_{t})-f(\mathbf{x}))^{2}\nu_{1}( \mathbf{x})\,d\mathbf{x}+\frac{1}{N}\sum_{j=1}^{N}\int_{\partial\Omega}( \mathcal{N}(\mathbf{y};\theta_{t})-g(\mathbf{y}))^{2}\nu_{2}(\mathbf{y})\,d \mathbf{y}\right]\] \[=\nabla_{\theta}\mathcal{J}(\mathcal{N}(\cdot;\theta_{t}))\]

Therefore, the \(\nabla_{\theta}\mathcal{L}(\theta_{t})\) are unbiased estimates of \(\nabla_{\theta}\mathcal{J}(\mathcal{N}(\cdot;\theta_{t}))\), and we can assume a step in the descent direction of \(\mathcal{L}\) is also one in \(\mathcal{J}\). Thus, any minimisation of \(\mathcal{L}\) should translate to a local minimisation of \(\mathcal{J}\).

#### Minimisation of \(\mathcal{J}(\mathcal{N})\)

We prove the following theorem, adapted from the original deep Galerkin method paper [4].

**Theorem 4**.: _Let \(\mathscr{A}(\sigma)\) be given by (1), for nonconstant, bounded \(\sigma\), and let \(\Omega\in\mathbb{R}^{d}\) be a compact domain and consider measures \(\nu_{1},\nu_{2}\) whose supports are contained in \(\Omega,\partial\Omega\) respectively. Assume further that \(\nabla^{2}\phi\) is locally Lipschitz with Lipschitz constant that can have at most polynomial growth on \(\nabla\phi\), uniformly with respect to \(x\), i.e._

\[|\nabla^{2}\mathcal{N}-\nabla^{2}\phi|\leq\left(||\nabla\mathcal{N}||^{\frac{ \alpha}{2}}+||\nabla\phi||^{\frac{\alpha}{2}}\right)||\nabla\mathcal{N}-\nabla \phi|| \tag{3}\]

_for some constants \(0\leq a,b<\infty\). Then, for all \(\epsilon>0\), there exists a constant \(\kappa>0\) such that there exists a function \(\mathcal{N}\in\mathscr{A}(\sigma)\) with_

\[\mathcal{J}(\mathcal{N})\leq\kappa\epsilon\]
Proof.: The condition given by (3) implies that

\[|\nabla^{2}\mathcal{N}-\nabla^{2}\phi|^{2} \leq\left(||\nabla\mathcal{N}||^{\frac{p}{2}}+||\nabla\phi||^{ \frac{p}{2}}\right)^{2}||\nabla\mathcal{N}-\nabla\phi||^{2}\] \[\leq\left(||\nabla\mathcal{N}||^{a}+||\nabla\phi||^{b}+2||\nabla \mathcal{N}||^{\frac{p}{2}}||\nabla\phi||^{\frac{p}{2}}\right)||\nabla\mathcal{N }-\nabla\phi||^{2}\] \[\leq 2\left(||\nabla\mathcal{N}||^{a}+||\nabla\phi||^{b}\right)|| \nabla\mathcal{N}-\nabla\phi||^{2}\]

with the last line following from Young's inequality [16]. Then,

\[\int_{\Omega}|\nabla^{2}\mathcal{N}(\mathbf{x})-\nabla^{2}\phi( \mathbf{x})|^{2}\,d\nu_{1}(\mathbf{x}) \leq 2\int_{\Omega}\left(||\nabla\mathcal{N}(\mathbf{x})||^{a}+|| \nabla\phi(\mathbf{x})||^{b}\right)||\nabla\mathcal{N}(\mathbf{x})-\nabla \phi(\mathbf{x})||^{2}\,d\nu_{1}(\mathbf{x})\] \[\leq 2\left[\int_{\Omega}\left(||\nabla\mathcal{N}(\mathbf{x})|| ^{a}+||\nabla\phi(\mathbf{x})||^{b}\right)^{p}\,d\nu_{1}(\mathbf{x})\right]^ {\frac{1}{p}}\left[\int_{\Omega}||\nabla\mathcal{N}(\mathbf{x})-\nabla\phi( \mathbf{x})||^{2q}\,d\nu_{1}(\mathbf{x})\right]^{\frac{1}{q}}\]

if we apply Holder's inequality [16] for exponents \(p,q\) satisfying \(\frac{1}{p}+\frac{1}{q}=1\) and \(1\leq p,q\leq\infty\). Furthermore,

\[\int_{\Omega}|\nabla^{2}\mathcal{N}(\mathbf{x})-\nabla^{2}\phi( \mathbf{x})|^{2}\,d\nu_{1}(\mathbf{x}) \leq K\left[\int_{\Omega}\left(||\nabla\mathcal{N}(\mathbf{x})- \nabla\phi(\mathbf{x})||^{a}+||\nabla\phi(\mathbf{x})||^{\max\{a,b\}}\right)^ {p}\,d\nu_{1}(\mathbf{x})\right]^{\frac{1}{p}}\] \[\quad\cdot\left[\int_{\Omega}||\nabla\mathcal{N}(\mathbf{x})- \nabla\phi(\mathbf{x})||^{2q}\,d\nu_{1}(\mathbf{x})\right]^{\frac{1}{q}}\] \[\leq K(\epsilon^{a}+\sup_{\mathbf{x}\in\Omega}||\nabla\phi( \mathbf{x})||^{\max\{a,b\}})\epsilon^{2}\]

for some constant \(K\). The last line follows from Theorem 1. Applying this result and Theorem 1 again to the objective function \(\mathcal{J}\), we obtain:

\[\mathcal{J}(\mathcal{N}) =\int_{\Omega}|\nabla^{2}\mathcal{N}(\mathbf{x})-f(\mathbf{x})|^{2 }\,d\nu_{1}(\mathbf{x})+\int_{\partial\Omega}|\mathcal{N}(\mathbf{x})-g( \mathbf{x})|^{2}\,d\nu_{2}(\mathbf{x})\] \[=\int_{\Omega}|\nabla^{2}\mathcal{N}(\mathbf{x})-\nabla^{2}\phi( \mathbf{x})|^{2}\,d\nu_{1}(\mathbf{x})+\int_{\partial\Omega}|\mathcal{N}( \mathbf{x})-\phi(\mathbf{x})|^{2}\,d\nu_{2}(\mathbf{x})\] \[\leq K(\epsilon^{a}+\sup_{\mathbf{x}\in\Omega}||\nabla\phi( \mathbf{x})||^{\max\{a,b\}})\epsilon^{2}+\epsilon^{2}\]

Finally, a rescaling of \(\epsilon>0\) yields

\[\mathcal{J}(\mathcal{N})\leq\kappa\epsilon\]

for some constant \(\kappa>0\) which may depend on \(\sup_{\mathbf{x}\in\Omega}||\nabla\phi(\mathbf{x})||\).

Theorem 4 guarantees the existence of a feed-forward neural network \(\mathcal{N}\) that, under relatively relaxed conditions, makes the objective function \(\mathcal{J}(\mathcal{N})\) for Poisson's equation arbitrarily small. However, neural network objective functions are highly non-convex. This means they have numerous minima and, while gradient descent algorithms like the deep Galerkin method are extremely effective at reaching said minima [17], there is no guarantee of achieving the global minimum i.e., in our case, finding the unique solution. Many authors research such ideas in non-convex optimisation [18], but we do not touch on them here, and present only empirical evidence of our solver finding/not finding global minima in the Results section (see **4**).
Methods

We now present three highly accessible methods to enhance the performance of a neural network trained to solve differential equations via the deep Galerkin method.

### Sinusoidal representation networks

Consider a neural network that is trained to approximate a function directly. We need only the first-order derivatives of the activation functions to backpropagate, and thus ReLU seems a natural choice [19]. However, our framework requires a network to learn a function via its derivatives. ReLU networks cannot do this without significant loss of information since they have second derivative zero. They are incapable of accurately modelling a signal's higher-order derivatives.

A recent paper [8] highlighting these limitations proposes something the authors call a sinusoidal representation network or SIREN. This is a neural network that implicitly defines a function, in our case \(\mathcal{N}\), with sinusoidal activations. Thus, while regular feed-forward networks with, say, ReLU activation may be excellent function approximators, a SIREN can further accurately fit derivatives of functions \(\phi\) through its own derivatives. ReLU networks typically cannot, due to their piecewise linear nature. This idea is hidden in Theorem 1 since ReLU is continuous but not differentiable, and so a network \(\mathcal{N}\) with ReLU activation could only achieve

\[\sup_{x}|\partial_{x}^{(\alpha)}\mathcal{N}(x)-\partial_{x}^{(\alpha)}\phi(x )|<\epsilon\]

for \(\alpha=0\). By contrast, \(\sin\in\mathcal{C}^{\infty}\), so the equivalent statement is true for any \(|\alpha|<\infty\).

Evaluating the gradient of a SIREN scales quadratically in the number of layers of the SIREN [8]. So, fitting higher-order derivatives is no easy task. However, for simple differential equations like Poisson's equation, it is computationally feasible, and the authors of [8] provide experimental results that show SIRENs are excellent at modelling first and second-order derivatives of complicated signals, as well as the high-frequency signals themselves.

### Random Fourier features

Recent works [20, 21] have described a spectral bias inherent to neural networks learning functions. They prioritise learning the low-frequency modes of the functions and thus, high frequencies are captured much later in the training procedure.

In many ways, this is a key reason behind the immense success of neural networks. Often, they are over-parameterised, i.e. the number of parameters far exceeds the number of training samples yet, counter-intuitively, they still show remarkable capacity to generalise well [22]. Spectral bias may explain part of this phenomenon because it suggests, if there is a way to fit data effectively with only low frequencies, then a neural network will do just this, without needing to resort to high frequencies that overfit the data.

However, this also means that neural networks struggle to learn high frequency functions. Theoretical results in [23] show that a one-dimensional function of pure frequency \(\omega\), e.g. \(\cos(\omega x)\), is learned in time that scales with \(\omega^{2}\). This is ratified experimentally.

A 2020 paper [9] publishes results on the use of a Fourier feature mapping to effectively overcome this spectral bias, and allow multilayer perceptrons (MLPs) to learn high frequency functions in low-dimensional domains. The authors motivate such work with neural tangent kernel (NTK) theory. NTKs have been shown to model the behaviour of MLPs in the infinite-width limit during training [24]. We do not describe them in detail here, but give a summary of the main idea behind Fourier feature mapping. For two different inputs \(\mathbf{x},\mathbf{x}^{\prime}\) to the MLP, the corresponding NTK can be given by
\[NTK(\mathbf{x},\mathbf{x}^{\prime})=h(\mathbf{x}^{T}\mathbf{x}^{\prime})\]

where \(h\) is some scalar function [9].

The mapping

\[\gamma(\mathbf{x})=[\cos(2\pi\mathbf{B}\mathbf{x}),\sin(2\pi\mathbf{B}\mathbf{x })]^{T} \tag{4}\]

is a Gaussian random Fourier feature mapping for \(\mathbf{x}\in\mathbb{R}^{d}\), where each entry in \(\mathbf{B}\in\mathbb{R}^{n\times d}\) is sampled from a normal distribution with mean zero and variance \(\Sigma^{2}\). Therefore,

\[NTK(\gamma(\mathbf{x}),\gamma(\mathbf{x}^{\prime})) =h(\gamma(\mathbf{x})^{T}\gamma(\mathbf{x}^{\prime}))\] \[=h\left(\cos(2\pi\mathbf{B}\mathbf{x})\cos(2\pi\mathbf{B}\mathbf{ x}^{\prime})+\sin(2\pi\mathbf{B}\mathbf{x})\sin(2\pi\mathbf{B}\mathbf{x}^{ \prime})\right)\] \[=h(\cos(2\pi\mathbf{B}(\mathbf{x}-\mathbf{x}^{\prime})))\]

Crucially, this defines a kernel function with width controlled by the random matrix \(\mathbf{B}\). Kernel functions are used to fit data, and their width directly influences whether they overfit (with high frequencies) or underfit (with low frequencies). So, given that this function characterises the evolution of the MLP during training, we can tune the network towards learning particular frequencies by simply changing \(\Sigma\):

* A small \(\Sigma\) gives a wide kernel that will underfit a high-frequency function.
* A large \(\Sigma\) gives a narrow kernel that will overfit a low-frequency function.

In our framework, \(\Sigma\) is now just another hyperparameter, and we can find the optimal \(\Sigma\) through a simple sweep of values. We choose the value that gives the fastest convergence. The authors of [9] also advise that \(n\), the number of Fourier features, improves performance with size. Of course, there is a computational cost associated with increasing \(n\), so it is best taken 'as small as gives good enough results.'

### Error correction

We introduce the main work of this paper; the novel technique error correction [10, 11, 12, 13, 14] is designed to increase the efficacy of any neural network differential equation solver. This method is general and can be applied to all differential equations, in combination with any such similar strategies, such as Koopman boosting [25] or those presented above. Much of the work here was proposed in [10] and formalised in [11], which the reader should refer to as supplement.

When dealing with neural networks, we bank on the idea that a'small enough' loss implies a 'good enough' accuracy. Now, in many scenarios, this ideology fails because zero loss would represent drastic overfitting. Conveniently, this does not concern us as we want our network to fit the (training) data as accurately as possible. Still, the original problem remains; how can we know how close we are to the true solution \(\phi\)?

It turns out analysis and estimation of the unknown error between \(\phi\) and \(\mathcal{N}\) is possible. Indeed, in [12], the author shows how you can obtain specific bounds on this error, without knowledge of \(\phi\). In this section, we provide a correction method (based on this error) to enhance neural network differential equation solvers, by overcoming performance saturation when the network settles around a local minimum of the loss function. Here, we also make use of differential equation operators which send true solutions to zero. Consider this for Poisson's equation (2):

\[\mathbf{F}[\cdot]=\nabla^{2}[\cdot]-f\]
Define \(\phi_{\epsilon}=\phi-\mathcal{N}\) as the error between the unknown solution \(\phi\) and a fixed approximation \(\mathcal{N}\). Clearly,

\[\mathbf{F}[\mathcal{N}] =\nabla^{2}\mathcal{N}-f\] \[=\nabla^{2}[\phi-\phi_{\epsilon}]-f\] \[=\nabla^{2}\phi-f-\nabla^{2}\phi_{\epsilon}\] \[=-\nabla^{2}\phi_{\epsilon}\]

since \(\mathbf{F}[\phi]=\nabla^{2}\phi-f=0\). Thus, \(\mathbf{F}[\mathcal{N}]+\nabla^{2}\phi_{\epsilon}=0\) and, given that \(\mathbf{F}[\mathcal{N}]\) is completely independent to \(\phi_{\epsilon}\), we have defined a new Poisson's equation. Our general strategy now will be to train a neural network \(\mathcal{N}_{\epsilon}\) to approximate \(\phi_{\epsilon}\) through the conditions of this new differential equation. Then, \(\mathcal{N}+\mathcal{N}_{\epsilon}\approx\mathcal{N}+\phi_{\epsilon}=\phi\).

Before we formalise and evaluate this method, note that it applies also to differential equations with non-linear terms. Consider the Poisson-Boltzmann equation with Dirichlet boundary conditions:

\[\begin{cases}\nabla^{2}\phi+\sinh\phi&=f\text{ in }\Omega\\ \phi&=g\text{ on }\partial\Omega\end{cases}\]

Define the operator

\[\mathbf{G}[\cdot]=\nabla^{2}[\cdot]+\sinh[\cdot]-f\]

and, once again, have \(\phi_{\epsilon}=\phi-\mathcal{N}\). Then,

\[\mathbf{G}[\mathcal{N}] =\nabla^{2}\mathcal{N}+\sinh\mathcal{N}-f\] \[=\nabla^{2}[\phi-\phi_{\epsilon}]+\sinh\mathcal{N}+\sinh\phi- \sinh\phi-f\] \[=\nabla^{2}\phi+\sinh\phi-f-\nabla^{2}\phi_{\epsilon}+\sinh \mathcal{N}-\sinh\phi\] \[=-\nabla^{2}\phi_{\epsilon}+\sinh\mathcal{N}-\sinh(\mathcal{N}+ \phi_{\epsilon})\]

since \(\mathbf{G}[\phi]=\nabla^{2}\phi+\sinh\phi-f=0\). A clever trick of adding and subtracting \(\sinh\phi\) allows the \(\mathbf{G}[\phi]\) term to be removed from the equation. In the last line, we simply seek to keep the equation explicit in \(\mathcal{N}\) and \(\phi_{\epsilon}\).

#### Theoretical results

Now, we formalise this idea of error correction, adapting the approach from [11]. Consider a differential equation over \(\Omega\) in operator form:

\[\mathbf{F_{0}}[\phi]=\mathbf{A}[\phi]+\mathbf{B}[\phi]+\mathbf{C}=0 \tag{5}\]

where \(\mathbf{A}\) represents the terms that depend linearly on \(\phi\), \(\mathbf{B}\) represents those that depend non-linearly on \(\phi\), and \(\mathbf{C}\) is independent of \(\phi\). The solution \(\phi\) may also admit some constraints on the boundary \(\partial\Omega\) but, for now, these are not of interest. Assume also that \(\phi\) is unique.

We first prove a result that follows from the inverse function theorem [26]:

**Theorem 5**.: (Inverse function theorem). _Suppose that \(F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}\) is continuously differentiable in some open set containing \(x^{*}\), and suppose moreover that the Jacobian \(DF(x^{*})\) is invertible. Then there exists open sets \(U,V\subset\mathbb{R}^{n}\) with \(x^{*}\in U\) and \(F(x^{*})\in V\) such that \(F:U\to V\) is a bijection, and \(F^{-1}:V\to U\) is continuously differentiable for all \(y\in V\) with_

\[DF^{-1}(y)=\left[DF(F^{-1}(y))\right]^{-1}\]
**Corollary 6**.: _Suppose that \(\mathbf{F_{0}}:\mathbb{R}\rightarrow\mathbb{R}\) in (5) is continuously differentiable in some open set containing \(\phi^{*}\), that \(D\mathbf{F_{0}}[\phi^{*}]\) is invertible, and \(\mathbf{F_{0}}[\phi^{*}]=0\). Then, there is a neighbourhood of \(0\) small enough such that_

\[\mathbf{F_{0}}[\mathcal{N}]\to 0\implies\mathcal{N}\rightarrow\phi^{*}\]

Proof.: By Theorem 5, choose neighbourhoods \(U,V\subset\mathbb{R}\) with \(\phi^{*}\in U,0\in V\) such that \(\mathbf{F_{0}}:U\to V\) is a bijection and \(\mathbf{F_{0}}^{-1}:V\to U\) is continuous differentiable for all \(y\in V\). For \(\mathcal{N}\in U\), the continuity of \(\mathbf{F_{0}}^{-1}\) implies that

\[\mathbf{F_{0}}[\mathcal{N}]\to 0\implies\mathcal{N}\rightarrow\phi^{*}\]

Thus, assuming we can minimise the loss function for some neural network \(\mathcal{N}\) such that \(\mathbf{F_{0}}[\mathcal{N}]\to 0\) at all points, then \(\mathcal{N}\rightarrow\phi\) at all points. So, let us train such a network \(\mathcal{N}_{0}\) to approximate \(\phi\) via (5). Define also \(\phi_{1}=\phi-\mathcal{N}_{0}\).

\[\mathbf{F_{0}}[\mathcal{N}_{0}] =\mathbf{A}[\mathcal{N}_{0}]+\mathbf{B}[\mathcal{N}_{0}]+\mathbf{C}\] \[=\mathbf{A}[\phi-\phi_{1}]+\mathbf{B}[\mathcal{N}_{0}]+\mathbf{B }[\phi]-\mathbf{B}[\phi]+\mathbf{C}\] \[=\mathbf{A}[\phi]+\mathbf{B}[\phi]+\mathbf{C}-\mathbf{A}[\phi_{1} ]+\mathbf{B}[\mathcal{N}_{0}]-\mathbf{B}[\phi]\] \[=-\mathbf{A}[\phi_{1}]+\mathbf{B}[\mathcal{N}_{0}]-\mathbf{B}[ \mathcal{N}_{0}+\phi_{1}]\]

since \(\mathbf{F_{0}}[\phi]=\mathbf{A}[\phi]+\mathbf{B}[\phi]+\mathbf{C}=0\) by definition. We have defined a new differential equation in operator form:

\[\mathbf{F_{1}}[\phi_{1}]=\mathbf{F_{0}}[\mathcal{N}_{0}]+\mathbf{A}[\phi_{1}]- \mathbf{B}[\mathcal{N}_{0}]+\mathbf{B}[\mathcal{N}_{0}+\phi_{1}]=0\]

\(\phi_{1}\) solves the above equation exactly and, given the uniqueness of \(\phi\), is also unique. Now, train some other neural network \(\mathcal{N}_{1}\) to approximate \(\phi_{1}\), and define \(\phi_{2}=\phi_{1}-\mathcal{N}_{1}\). Once again,

\[\mathbf{F_{1}}[\mathcal{N}_{1}] =\mathbf{F_{0}}[\mathcal{N}_{0}]+\mathbf{A}[\mathcal{N}_{1}]- \mathbf{B}[\mathcal{N}_{0}]+\mathbf{B}[\mathcal{N}_{0}+\mathcal{N}_{1}]\] \[=\mathbf{F_{0}}[\mathcal{N}_{0}]+\mathbf{A}[\phi_{1}-\phi_{2}]- \mathbf{B}[\mathcal{N}_{0}]+\mathbf{B}[\mathcal{N}_{0}+\mathcal{N}_{1}]+ \mathbf{B}[\mathcal{N}_{0}+\phi_{1}]-\mathbf{B}[\mathcal{N}_{0}+\phi_{1}]\] \[=\mathbf{F_{0}}[\mathcal{N}_{0}]+\mathbf{A}[\phi_{1}]-\mathbf{B}[ \mathcal{N}_{0}]+\mathbf{B}[\mathcal{N}_{0}+\phi_{1}]-\mathbf{A}[\phi_{2}]+ \mathbf{B}[\mathcal{N}_{0}+\mathcal{N}_{1}]-\mathbf{B}[\mathcal{N}+\phi_{1}]\] \[=-\mathbf{A}[\phi_{2}]+\mathbf{B}[\mathcal{N}_{0}+\mathcal{N}_{1}] -\mathbf{B}[\mathcal{N}_{0}+\phi_{1}]\] \[=-\mathbf{A}[\phi_{2}]+\mathbf{B}[\mathcal{N}_{0}+\mathcal{N}_{1} ]-\mathbf{B}[\mathcal{N}_{0}+\mathcal{N}_{1}+\phi_{2}]\]

since \(\mathbf{F_{1}}[\phi_{1}]=\mathbf{F_{0}}[\mathcal{N}_{0}]+\mathbf{A}[\phi_{1}]- \mathbf{B}[\mathcal{N}_{0}]+\mathbf{B}[\mathcal{N}_{0}+\phi_{1}]=0\), and \(\phi_{1}=\mathcal{N}_{1}+\phi_{2}\). We define a further differential equation in operator form:

\[\mathbf{F_{2}}[\phi_{2}]=\mathbf{F_{1}}[\mathcal{N}_{1}]+\mathbf{A}[\phi_{2}]- \mathbf{B}[\mathcal{N}_{0}+\mathcal{N}_{1}]+\mathbf{B}[\mathcal{N}_{0}+\mathcal{ N}_{1}+\phi_{2}]=0\]

Now, repeat the process. This algorithm can continue indefinitely, and we summarise the steps below. The idea is that our error-corrected approximation \(\mathcal{N}_{0}+\mathcal{N}_{1}+\mathcal{N}_{2}+...\) will be more accurate than the once-trained approximation \(\mathcal{N}_{0}\). This strategy is not unseen in the field of numerical methods to differential equations, we just apply it here to neural network solvers.

Let us define a recursive differential equation for the \(k^{\text{th}}\) error correction. At this point, we have trained the initial network \(\mathcal{N}_{0}\), and also a further \(k-1\) residual networks \(\mathcal{N}_{1},\mathcal{N}_{2},...,\mathcal{N}_{k-1}\). Our current error-corrected
approximation is \(\mathcal{N}^{(k-1)}=\mathcal{N}_{0}+\mathcal{N}_{1}+\mathcal{N}_{2}+...+\mathcal{ N}_{k-1}\). Define \(\phi_{k}=\phi_{k-1}-\mathcal{N}_{k-1}\). Now, train a new network \(\mathcal{N}_{k}\) to approximate \(\phi_{k}\) through the following differential equation:

\[\mathbf{F}_{\mathbf{k}}[\phi_{k}]=\mathbf{F}_{\mathbf{k-1}}[\mathcal{N}_{k-1}]+ \mathbf{A}[\phi_{k}]-\mathbf{B}[\mathcal{N}^{(k-1)}]+\mathbf{B}[\mathcal{N}^{( k-1)}+\phi_{k}]=0 \tag{6}\]

**Remark**.: \(\mathbf{F}_{\mathbf{k}}[\mathcal{N}_{k}]\equiv\mathbf{F}_{\mathbf{0}}[ \mathcal{N}^{(k)}]\)_._

**Corollary 7**.: _Suppose that \(\mathbf{F}_{\mathbf{k}}:\mathbb{R}\rightarrow\mathbb{R}\) in (6) is continuously differentiable in some open set containing \(\phi_{k}^{*}\), that \(D\mathbf{F}_{\mathbf{k}}[\phi_{k}^{*}]\) is invertible, and \(\mathbf{F}_{\mathbf{k}}[\phi_{k}^{*}]=0\). Then, there is a neighbourhood of \(0\) small enough such that_

\[\mathbf{F}_{\mathbf{k}}[\mathcal{N}_{k}]\to 0\implies\mathcal{N}_{k} \rightarrow\phi_{k}^{*}\]

_Furthermore,_

\[|\mathcal{N}^{(k)}-\phi|=\mathcal{O}\left(|\mathbf{F}_{\mathbf{k}}[\mathcal{N }_{k}]|\right)\]

Proof.: The first result follows analogously from the inverse function theorem as in Corollary 6.

By Theorem 5, \(\mathbf{F_{\mathbf{0}}}^{-1}\) is continuously differentiable on some open set around \(0\). Thus, it is also locally Lipschitz continuous around \(0\), meaning there exists some constant \(\alpha\geq 0\) such that

\[|\mathcal{N}^{(k)}-\phi| =|\mathbf{F_{\mathbf{0}}}^{-1}[\mathbf{F_{\mathbf{0}}}[\mathcal{ N}^{(k)}]-\mathbf{F_{\mathbf{0}}}^{-1}[\mathbf{F_{\mathbf{0}}}[\phi]]|\] \[\leq\alpha|\mathbf{F_{\mathbf{0}}}[\mathcal{N}^{(k)}]-\mathbf{F_ {\mathbf{0}}}[\phi]|\] \[\leq\alpha|\mathbf{F_{\mathbf{0}}}[\mathcal{N}^{(k)}]|\] \[\leq\alpha|\mathbf{F_{\mathbf{k}}}[\mathcal{N}_{k}]|\]

and therefore,

\[|\mathcal{N}^{(k)}-\phi|=\mathcal{O}\left(|\mathbf{F}_{\mathbf{k}}[\mathcal{N }_{k}]|\right)\]

Finally, given Dirichlet boundary conditions \(\phi=g\) on \(\partial\Omega\), any \(\phi_{k}\) is known exactly over \(\partial\Omega\) since \(\phi_{k}=\phi-\mathcal{N}^{(k-1)}\). Thus, the loss function for the \(k^{\text{th}}\) error correction can be defined as

\[\mathcal{L}_{k}(\theta^{(k)})=\frac{1}{M}\sum_{i=1}^{M}\left( \mathbf{F}_{\mathbf{k}}\left[\mathcal{N}_{k}\left(\mathbf{x}_{i};\theta^{(k) }\right)\right]\right)^{2}+\frac{1}{N}\sum_{j=1}^{N}\left(\mathcal{N}_{k} \left(\mathbf{y}_{j};\theta^{(k)}\right)-\phi_{k}(\mathbf{y}_{j})\right)^{2} \tag{7}\]

for some randomly sampled points \(\{\mathbf{x}_{i}\}_{i=1}^{M}\) from \(\Omega\) and \(\{\mathbf{y}_{j}\}_{j=1}^{N}\) from \(\partial\Omega\).

**Algorithm**

The error correction algorithm to order \(K\) proceeds as follows:

1. Train a neural network \(\mathcal{N}_{0}\) to satisfy the conditions of a differential equation given by (5) and constraint conditions. Once the loss has converged, stop training and freeze the parameters of \(\mathcal{N}_{0}\).
2. Initiate and train new neural networks \(\{\mathcal{N}_{k}\}_{k=1}^{K}\) in sequence to satisfy differential equations given by (6), via loss functions (7). Once the loss has converged, stop training, freeze the parameters of \(\mathcal{N}_{k}\), and proceed with \(\mathcal{N}_{k+1}\).
3. The solution to (5) is approximated by \(\mathcal{N}:=\mathcal{N}^{(K)}=\sum\limits_{k=0}^{K}\mathcal{N}_{k}\).

This is given above for Dirichlet boundary conditions, but works generally if you incorporate the constraint conditions into all loss functions.

#### Poisson's equation

For Poisson's equation (2), \(\mathbf{B}\equiv 0\) since the Laplacian is linear, so we can write the \(k^{\text{th}}\) differential equation as

\[\mathbf{F_{k}}[\phi_{k}]=\mathbf{F_{k-1}}[\mathcal{N}_{k-1}]+\nabla^{2}[\phi_ {k}]=0\]

which is a Poisson's equation with our usual \(f=-\mathbf{F_{k-1}}[\mathcal{N}_{k-1}]\). Thus, we can apply Theorem 4 to guarantee that there are neural networks out there that can get very, very close to the \(\phi_{k}\). In the next section, we provide evidence that shows, if we can train just two or three of these networks to reasonably approximate their true solutions, our error-corrected approximation will be a more accurate numerical solution to the original differential equation.

## 4 Results

We present results for a variety of different Poisson's equations (2). Our choice of Poisson's equation is motivated by its immense application in many areas of theoretical physics, including electrostatics and fluid dynamics. It is also the simplest second-order, linear PDE, making for a concise yet insightful demonstration of the power of error correction in neural network differential equation solvers.

To achieve this, we choose the function \(f\) on the RHS to force a particular solution \(\phi\) that we want to capture. For example, \(f(x)=1\) would force the solution \(\phi(x)=\frac{1}{2}x^{2}+c_{1}x+c_{0}\). However, in general, \(f\) can be anything, particularly something which does not admit a closed-form solution to (2), and we do this for ease of visualising \(\phi\).

Knowing the ground truth solution \(\phi\) in closed form also allows us to compute the relative error

\[\frac{\sum\limits_{\mathbf{x}\in S}\left(\phi(\mathbf{x})-\mathcal{N}(\mathbf{ x})\right)^{2}}{\sum\limits_{\mathbf{x}\in S}\phi(\mathbf{x})^{2}}\]

at each epoch (iteration) of the training procedure, so we have an understanding of the success of our solver. It is important to note that, while we know \(\phi\) and the relative error associated with our approximation, the neural network does not, and is solely trained via the loss function.

All neural networks used are SIRENs with 5 hidden layers and 128 hidden units per layer. They are trained on batches of 256, using the stochastic gradient descent variant Adam [27], and learning rates are manually tuned for each case of Poisson's equation. All experiments are run on a 1.8 GHz Dual-Core Intel Core i5 CPU.
### 3D Poisson's Equation

Figure 1 shows the solution to Poisson's equation:

\[\begin{cases}\nabla^{2}\phi&=-75\sin(5x)\sin(5y)\sin(5z)\text{ in }\Omega=[-\pi,\pi] \times[-\pi,\pi]\times[-\pi,\pi]\\ \phi&=0\text{ on }\partial\Omega\end{cases} \tag{8}\]

Figure 2 shows our numerical solutions, with \(\mathcal{N}^{(0)}=\mathcal{N}_{0}\) on the left, \(\mathcal{N}^{(1)}=\mathcal{N}_{0}+\mathcal{N}_{1}\) in the centre, and \(\mathcal{N}^{(2)}=\mathcal{N}_{0}+\mathcal{N}_{1}+\mathcal{N}_{2}\) on the right. We refer to these as Error Correction 0, 1 and 2, respectively.

Visually, all error corrections seem to capture the solution well. Furthermore, each correction decreases the relative error (printed at the bottom of of Figure 2). Error Correction 1 does so significantly, while the improvement in accuracy from Error Correction 2 is marginal.

This is further captured in Figure 3, which is a plot of the loss and relative error per epoch. After finding a local minimum in Error Correction 0, the loss fluctuates erratically until we initialise Error Correction 1. The improvement is truly appreciable, and felt across the trends in relative error too.
Figure 3: Per-epoch loss and relative errors for numerical solutions to (8)
### 2D Poisson's Equation

Figure 4 shows the solution to Poisson's equation:

\[\begin{cases}\nabla^{2}\phi&=-800\sin(5x)\sin(5y)\text{ in }\Omega=[-\pi,\pi] \times[-\pi,\pi]\\ \phi&=0\text{ on }\partial\Omega\end{cases} \tag{9}\]

Due to the highly oscillatory nature of the solution, a neural network will struggle to accurately capture its structure. This is demonstrated in Figure 5, where the approximation cannot account for so many peaks and troughs in the solution.

To obtain a realistic solution, we apply a Gaussian random Fourier feature mapping to the input, before passing it through the network. After a simple sweep of values, we take \(\Sigma=1\) and \(n=256\), as defined in (4). Figures 6 and 7 show similar trends to those in the previous experiment.

Figure 5: Naive attempt at a numerical solution to (9)
Figure 7: Per-epoch loss and relative errors for numerical solutions to (9)
Figure 8 shows the solution to Poisson's equation:

\[\begin{cases}\nabla^{2}\phi&=(100y^{2}-100\pi^{2}-2)\sin(10x)\text{ in }\Omega=[-\pi,\pi]\times[-\pi,\pi]\\ \phi&=0\text{ on }\partial\Omega\end{cases} \tag{10}\]

In Figure 9, we train a neural network \(\mathcal{N}^{(0)}\) to approximate the solution to (10) for \(2^{11}\) epochs, but we save its parameter states after \(2^{10}\) epochs. These define a new network which we call \(\mathcal{N}_{0}\). The fully-trained \(\mathcal{N}^{(0)}\) achieves a reasonable relative error. Roughness is clearly visible in the plot.

In Figure 10, we plot the half-trained \(\mathcal{N}_{0}\) on the left. As expected, it has not yet reached the accuracy of \(\mathcal{N}^{(0)}\). However, we also initiate an error correction \(\mathcal{N}_{1}\), of \(\mathcal{N}_{0}\), that trains for another \(2^{10}\) epochs. Thus, we produce an approximation \(\mathcal{N}^{(1)}=\mathcal{N}_{0}+\mathcal{N}_{1}\) that has also trained for a total of \(2^{11}\) epochs. This is significantly more accurate than \(\mathcal{N}^{(0)}\), and the plot is visibly smoother. Figures 9 and 10 provide a clear exemplification of the immediate fruitfulness of a single error correction.
Figure 10: Numerical solutions \(\mathcal{N}_{0}\) and \(\mathcal{N}^{(1)}\) to (10), trained for a total of \(2^{11}\) epochs

Figure 9: Numerical solution \(\mathcal{N}^{(0)}\) to (10), trained for \(2^{11}\) epochs
Discussion

Our results do not endorse error correction as a tool to marginally reduce error across tens of corrections. Instead, they suggest training a network for half the allotted time, and devoting the other half to a single error correction. This can yield significantly more accurate results.

Error correction is not without cost however. In our implementation, we train correction networks on newly sampled points. This means that to obtain \(\mathbf{F_{k}}[\mathcal{N}_{k}]\), we must first make \(k\) forward passes of the new data through \(\mathcal{N}_{0},\mathcal{N}_{1},...,\mathcal{N}_{k-1}\) and differentiate these to compute \(\mathbf{F_{k-1}}[\mathcal{N}_{k-1}]\). The time complexity of producing a \(k^{\text{th}}\) order approximation \(\mathcal{N}^{(k)}\), assuming the number of epochs \(E\) and batch size \(B\) per correction, and optimisation costs, are kept constant across all corrections, is \(\mathcal{O}\left(EB(k+1)^{2}\right)\). If we instead pass identical batches through each correction network, storing the \(\mathbf{F_{k-1}}[\mathcal{N}_{k-1}]\) in memory, we can have a time complexity of \(\mathcal{O}(EB(k+1))\), however the space complexity would be substantially increased.

## 6 Further Work

Over time, this study of neural network differential equation solvers naturally lent itself to hot-off-the-press topics in machine learning like sinusoidal representation networks [8] and random Fourier features [9], for the simple reason that such concepts are inextricably linked through their applications. Outside of differential equations, neural networks as continuous parameterisations of discrete signals have immense potential in 3D shape representation, but also in image, video and audio representation and reconstruction. These problems may utilise neural networks as function approximators or, as we did, derivative approximators. There is no reason to suggest why the ideas of error correction cannot be employed here, and every reason to further explore the interplay of these techniques when applied to problems in computer vision.
