# Enhancing Neural Network Differential Equation Solvers

Matthew J. H. Wright

(December 2022)

###### 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.

##  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}$.

###  2.1 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)$ (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.

###  2.2 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}$ (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. 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. 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. 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. 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,

$\displaystyle\mathbb{E}[\nabla_{\theta}\mathcal{L}(\theta_{t})|\theta_{t}]$
$\displaystyle=\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]$
$\displaystyle=\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]$
$\displaystyle=\nabla_{\theta}\left[\int_{\Omega}(\nabla^{2}\mathcal{N}(\mathbf{x};\theta_{t})-f(\mathbf{x}))^{2}\nu_{1}(\mathbf{x})\,d\mathbf{x}+\int_{\partial\Omega}(\mathcal{N}(\mathbf{y};\theta_{t})-g(\mathbf{y}))^{2}\nu_{2}(\mathbf{y})\,d\mathbf{y}\right]$
$\displaystyle=\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{a}{2}}+||\nabla\phi||^{\frac{b}{2}}\right)||\nabla\mathcal{N}-\nabla\phi||$
(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

$\displaystyle|\nabla^{2}\mathcal{N}-\nabla^{2}\phi|^{2}$
$\displaystyle\leq\left(||\nabla\mathcal{N}||^{\frac{a}{2}}+||\nabla\phi||^{\frac{b}{2}}\right)^{2}||\nabla\mathcal{N}-\nabla\phi||^{2}$
$\displaystyle\leq\left(||\nabla\mathcal{N}||^{a}+||\nabla\phi||^{b}+2||\nabla\mathcal{N}||^{\frac{a}{2}}||\nabla\phi||^{\frac{b}{2}}\right)||\nabla\mathcal{N}-\nabla\phi||^{2}$
$\displaystyle\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,

$\displaystyle\int_{\Omega}|\nabla^{2}\mathcal{N}(\mathbf{x})-\nabla^{2}\phi(\mathbf{x})|^{2}\,d\nu_{1}(\mathbf{x})$
$\displaystyle\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})$
$\displaystyle\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 Hölder’s inequality [16] for exponents $p,q$ satisfying
$\frac{1}{p}+\frac{1}{q}=1$ and $1\leq p,q\leq\infty$. Furthermore,

$\displaystyle\int_{\Omega}|\nabla^{2}\mathcal{N}(\mathbf{x})-\nabla^{2}\phi(\mathbf{x})|^{2}\,d\nu_{1}(\mathbf{x})$
$\displaystyle\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}}$
$\displaystyle\quad\cdot\left[\int_{\Omega}||\nabla\mathcal{N}(\mathbf{x})-\nabla\phi(\mathbf{x})||^{2q}\,d\nu_{1}(\mathbf{x})\right]^{\frac{1}{q}}$
$\displaystyle\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:

$\displaystyle\mathcal{J}(\mathcal{N})$
$\displaystyle=\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})$
$\displaystyle=\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})$
$\displaystyle\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\limits_{\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).

##  3 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.

###  3.1 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.

###  3.2 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{Bx}),\sin(2\pi\mathbf{Bx})]^{T}$ (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,

$\displaystyle NTK(\gamma(\mathbf{x}),\gamma(\mathbf{x^{\prime}}))$
$\displaystyle=h(\gamma(\mathbf{x})^{T}\gamma(\mathbf{x^{\prime}}))$
$\displaystyle=h\left(\cos(2\pi\mathbf{Bx})\cos(2\pi\mathbf{Bx^{\prime}})+\sin(2\pi\mathbf{Bx})\sin(2\pi\mathbf{Bx^{\prime}})\right)$
$\displaystyle=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.’

###  3.3 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,

$\displaystyle\mathbf{F}[\mathcal{N}]$ $\displaystyle=\nabla^{2}\mathcal{N}-f$
$\displaystyle=\nabla^{2}[\phi-\phi_{\epsilon}]-f$
$\displaystyle=\nabla^{2}\phi-f-\nabla^{2}\phi_{\epsilon}$
$\displaystyle=-\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,

$\displaystyle\mathbf{G}[\mathcal{N}]$
$\displaystyle=\nabla^{2}\mathcal{N}+\sinh\mathcal{N}-f$
$\displaystyle=\nabla^{2}[\phi-\phi_{\epsilon}]+\sinh\mathcal{N}+\sinh\phi-\sinh\phi-f$
$\displaystyle=\nabla^{2}\phi+\sinh\phi-f-\nabla^{2}\phi_{\epsilon}+\sinh\mathcal{N}-\sinh\phi$
$\displaystyle=-\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$ (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}\to\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}\to\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}\to\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}\to\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}\to\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}$.

$\displaystyle\mathbf{F_{0}}[\mathcal{N}_{0}]$
$\displaystyle=\mathbf{A}[\mathcal{N}_{0}]+\mathbf{B}[\mathcal{N}_{0}]+\mathbf{C}$
$\displaystyle=\mathbf{A}[\phi-\phi_{1}]+\mathbf{B}[\mathcal{N}_{0}]+\mathbf{B}[\phi]-\mathbf{B}[\phi]+\mathbf{C}$
$\displaystyle=\mathbf{A}[\phi]+\mathbf{B}[\phi]+\mathbf{C}-\mathbf{A}[\phi_{1}]+\mathbf{B}[\mathcal{N}_{0}]-\mathbf{B}[\phi]$
$\displaystyle=-\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,

$\displaystyle\mathbf{F_{1}}[\mathcal{N}_{1}]$
$\displaystyle=\mathbf{F_{0}}[\mathcal{N}_{0}]+\mathbf{A}[\mathcal{N}_{1}]-\mathbf{B}[\mathcal{N}_{0}]+\mathbf{B}[\mathcal{N}_{0}+\mathcal{N}_{1}]$
$\displaystyle=\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}]$
$\displaystyle=\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}]$
$\displaystyle=-\mathbf{A}[\phi_{2}]+\mathbf{B}[\mathcal{N}_{0}+\mathcal{N}_{1}]-\mathbf{B}[\mathcal{N}_{0}+\phi_{1}]$
$\displaystyle=-\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\textsuperscript{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_{k}}[\phi_{k}]=\mathbf{F_{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$
(6)

######  Remark.

$\mathbf{F_{k}}[\mathcal{N}_{k}]\equiv\mathbf{F_{0}}[\mathcal{N}^{(k)}]$.

######  Corollary 7.

Suppose that $\mathbf{F_{k}}:\mathbb{R}\to\mathbb{R}$ in (6) is continuously
differentiable in some open set containing $\phi_{k}^{*}$, that
$D\mathbf{F_{k}}[\phi_{k}^{*}]$ is invertible, and
$\mathbf{F_{k}}[\phi_{k}^{*}]=0$. Then, there is a neighbourhood of $0$ small
enough such that

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

Furthermore,

$|\mathcal{N}^{(k)}-\phi|=\mathcal{O}\left(|\mathbf{F_{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_{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

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

and therefore,

$|\mathcal{N}^{(k)}-\phi|=\mathcal{O}\left(|\mathbf{F_{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\textsuperscript{th}$ error correction can be defined as

$\mathcal{L}_{k}(\theta^{(k)})=\frac{1}{M}\sum_{i=1}^{M}\left(\mathbf{F_{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}$
(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. 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. 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. 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\textsuperscript{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.

###  4.1 3D Poisson’s Equation

Figure 1: $\phi(x,y,z)=\sin(5x)\sin(5y)\sin(5z)$ at $z=\frac{\pi}{10}$

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}$ (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 2: Numerical solutions $\mathcal{N}^{(0)},\mathcal{N}^{(1)}$ and
$\mathcal{N}^{(2)}$ to (8) Figure 3: Per-epoch loss and relative errors for
numerical solutions to (8)

###  4.2 2D Poisson’s Equation

#### (i)

Figure 4: $\phi(x,y)=\sin(20x)\sin(20y)$

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}$ (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.

Figure 5: Naive attempt at a numerical solution to (9)

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 6: Numerical solutions $\mathcal{N}^{(0)},\mathcal{N}^{(1)}$ and
$\mathcal{N}^{(2)}$ to (9), trained using random Fourier features with
$\Sigma=1$ and $n=256$ Figure 7: Per-epoch loss and relative errors for
numerical solutions to (9)

#### (ii)

Figure 8: $\phi(x,y)=(\pi^{2}-y^{2})\sin(10x)$

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}$ (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 9: Numerical solution $\mathcal{N}^{(0)}$ to (10), trained for $2^{11}$
epochs Figure 10: Numerical solutions $\mathcal{N}_{0}$ and
$\mathcal{N}^{(1)}$ to (10), trained for a total of $2^{11}$ epochs

##  5 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\textsuperscript{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.

## References

  * [1] I.E. Lagaris, A. Likas, and D.I. Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE Transactions on Neural Networks, 9(5):987–1000, 1998.
  * [2] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Physics informed deep learning (part I): Data-driven solutions of nonlinear partial differential equations.  arXiv:1711.10561, 2017.
  * [3] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Physics informed deep learning (part II): Data-driven discovery of nonlinear partial differential equations.  arXiv:1711.10566, 2017.
  * [4] Justin Sirignano and Konstantinos Spiliopoulos.  DGM: A deep learning algorithm for solving partial differential equations.  Journal of Computational Physics, 375:1339–1364, 2018.
  * [5] Anil Ananthaswamy.  Latest neural nets solve world’s hardest equations faster than ever before.  Quanta Magazine, 2021.
  * [6] Kurt Hornik.  Approximation capabilities of multilayer feedforward networks.  Neural Networks, 4(2):251–257, 1991.
  * [7] Andrew R. Barron.  Universal approximation bounds for superpositions of a sigmoidal function.  IEEE Transactions on Information Theory, 39(3):930–945, 1993.
  * [8] Vincent Sitzmann, Julien N.P. Martel, Alexander W. Bergman, David B. Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  In Proceedings of the 34th Conference on Neural Information Processing System, 2020.
  * [9] Matthew Tancik, Pratul P. Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan T. Barron, and Ren Ng.  Fourier features let networks learn high frequency functions in low dimensional domains.  In Proceedings of the 34th Conference on Neural Information Processing Systems, 2020.
  * [10] Akshunna S. Dogra.  Error estimation and correction from within neural network differential equation solvers.  arXiv:2007.04433, 2020.
  * [11] Akshunna S. Dogra et. al. Neural network differential equation solvers allow unsupervised error analysis and correction.  under review, 2023.
  * [12] Marios Mattheakis, David Sondak, Akshunna S. Dogra, and Pavlos Protopapas.  Hamiltonian neural networks for solving equations of motion.  Physical Review E, 105(6), 2022.
  * [13] Akshunna S. Dogra and William T. Redman.  Local error quantification for neural network differential equation solvers.  arXiv:2008.12190, 2021.
  * [14] Akshunna S. Dogra.  Dynamical systems and neural networks.  arXiv:2004.11826, 2020.
  * [15] Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  Neural Networks, 2(5):359–366, 1989.
  * [16] G. H. Hardy, John E. Littlewood, and George Pólya.  Inequalities.  Cambridge University Press, Cambridge, second edition, 1988.
  * [17] Simon S. Du, Jason D. Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai.  Gradient descent finds global minima of deep neural networks.  In Proceedings of the 36th International Conference on Machine Learning, 2019.
  * [18] Rudrajit Das.  Recent Advances in Non-Convex Optimization for Deep Learning.  https://rudrajit15.github.io/posts/2018/09/blog-post-2/.  Accessed: 2 May 2022.
  * [19] Xavier Glorot, Antoine Bordes, and Yoshua Bengio.  Deep sparse rectifier neural networks.  In Proceedings of the 14th International Conference on Artificial Intelligence and Statistics, 2011.
  * [20] Ronen Basri, Meirav Galun, Amnon Geifman, David Jacobs, Yoni Kasten, and Shira Kritchman.  Frequency bias in neural networks for input of non-uniform density.  In Proceedings of the 37th International Conference on Machine Learning, 2020.
  * [21] Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred A. Hamprecht, Yoshua Bengio, and Aaron Courville.  On the spectral bias of neural networks.  In Proceedings of the 36th International Conference on Machine Learning, 2019.
  * [22] Yuanzhi Li and Yingyu Liang.  Learning overparameterized neural networks via stochastic gradient descent on structured data.  In Proceedings of the 32nd Conference on Neural Information Processing Systems, 2018.
  * [23] Ronen Basri, David Jacobs, Yoni Kasten, and Shira Kritchman.  The convergence rate of neural networks for learned functions of different frequencies.  In Proceedings of the 33rd Conference on Neural Information Processing Systems, 2019.
  * [24] Arthur Jacot, Franck Gabriel, and Clement Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  In Proceedings of the 32nd Conference on Neural Information Processing Systems, 2018.
  * [25] Akshunna S. Dogra and William T. Redman.  Optimizing neural networks via koopman operator theory.  Advances in Neural Information Processing Systems 33 (NeurIPS 2020), 2020.
  * [26] Johannes Nicaise.  Lecture notes in real analysis, 2020.
  * [27] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In Proceedings of the 3rd International Conference for Learning Representations, 2015.