# DeepGreen: Deep Learning of Green’s Functions
for Nonlinear Boundary Value Problems

Craig R. Gin1111Co-first authors. , Daniel E. Shea2∗, Steven L. Brunton3, and
J. Nathan Kutz4
1Department of Population Health and Pathobiology, North Carolina State
University,
Raleigh, NC 27695, United States
2Department of Materials Science and Engineering, University of Washington,
Seattle, WA 98195, United States
3Department of Mechanical Engineering, University of Washington,
Seattle, WA 98195, United States
4Department of Applied Mathematics, University of Washington,
Seattle, WA 98195, United States

∗Co-first authors and corresponding authors. Email<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>
###### Abstract

Boundary value problems (BVPs) play a central role in the mathematical
analysis of constrained physical systems subjected to external forces.
Consequently, BVPs frequently emerge in nearly every engineering discipline
and span problem domains including fluid mechanics, electromagnetics, quantum
mechanics, and elasticity. The fundamental solution, or Green’s function, is a
leading method for solving linear BVPs that enables facile computation of new
solutions to systems under any external forcing. However, fundamental Green’s
function solutions for nonlinear BVPs are not feasible since linear
superposition no longer holds. In this work, we propose a flexible deep
learning approach to solve nonlinear BVPs using a dual-autoencoder
architecture. The autoencoders discover an invertible coordinate transform
that linearizes the nonlinear BVP and identifies both a linear operator $L$
and Green’s function $G$ which can be used to solve new nonlinear BVPs. We
find that the method succeeds on a variety of nonlinear systems including
nonlinear Helmholtz and Sturm–Liouville problems, nonlinear elasticity, and a
2D nonlinear Poisson equation. The method merges the strengths of the
universal approximation capabilities of deep learning with the physics
knowledge of Green’s functions to yield a flexible tool for identifying
fundamental solutions to a variety of nonlinear systems.

Keywords: Deep learning, Green’s function, Nonlinearity, Koopman operator
theory

##  1 Introduction

Boundary value problems (BVPs) are ubiquitous in the sciences (?). From
elasticity to quantum electronics, BVPs have been fundamental in the
development and engineering design of numerous transformative technologies of
the 20th century. Historically, the formulation of many canonical problems in
physics and engineering result in linear BVPs: from Fourier formulating the
heat equation in 1822 (?) to more modern applications such as designing chip
architectures in the semi-conductor industry (?, ?). Much of our theoretical
understanding of BVPs comes from the construction of the fundamental solution
of the BVP, commonly known as the Green’s function (?). The Green’s function
solution relies on a common property of many BVPs: linearity. Specifically,
general solutions rely on linear superposition to hold, thus limiting their
usefulness in many modern applications where BVPs are often heterogeneous and
nonlinear. By leveraging modern deep learning, we are able to learn
linearizing transformations of BVPs that render nonlinear BVPs linear so that
we can construct the Green’s function solution. Our deep learning of Green’s
functions, DeepGreen, provides a transformative architecture for modern
solutions of nonlinear BVPs.

DeepGreen is inspired by recent works which use deep neural networks (DNNs) to
discover advantageous coordinate transformations for dynamical systems (?, ?,
?, ?, ?, ?, ?, ?, ?, ?). The universal approximation properties of DNNs (?, ?)
are ideal for learning coordinate transformations that linearize nonlinear
BVPs, ODEs and PDEs. Specifically, such linearizing transforms fall broadly
under the umbrella of Koopman operator theory (?), which has a modern
interpretation in terms of dynamical systems theory (?, ?, ?, ?). There are
only limited cases in which Koopman operators can only be constructed
explicitly (?). However Dynamic Mode Decomposition (DMD) (?) provides a
numerical algorithm for approximating the Koopman operator (?), with many
recent extensions that improve on the DMD approximation (?). More recently,
neural networks have been used to construct Koopman embeddings (?, ?, ?, ?, ?,
?, ?, ?). This is an alternative to enriching the observables of DMD (?, ?, ?,
?, ?, ?, ?). Thus, neural networks have emerged as a highly effective
mathematical tool for approximating complex data (?, ?) with a linear model.
DNNs have been used in this context to discover time-stepping algorithms for
complex systems (?, ?, ?, ?, ?). Moreover, DNNs have been used to approximate
constitutive models of BVPs (?).

DeepGreen leverages the success of DNNs for dynamical systems to discover
coordinate transformations that linearize nonlinear BVPs so that the Green’s
function solution can be recovered. This allows for the discovery of the
fundamental solutions for nonlinear BVPs, opening many opportunities for the
engineering and physical sciences. DeepGreen exploits physics-informed
learning by using autoenconders (AEs) to take data from the original high-
dimensional input space to the new coordinates at the intrinsic rank of the
underlying physics (?, ?, ?). The architecture also leverages the success of
Deep Residual Networks (DRN) (?) which enables our approach to efficiently
handle near-identity coordinate transformations (?).

The Green’s function constructs the solution to a BVP for any given forcing by
superposition. Specifically, consider the classical linear BVP (?)

$\begin{array}[]{ll}L[v({\mathbf{x}})]=f({\mathbf{x}})\end{array}$ (1)

where $L$ is a linear differential operator, $f$ is a forcing,
${\mathbf{x}}\in{\Omega}$ is the spatial coordinate, and ${\Omega}$ is an open
set. The boundary conditions $Bv({\mathbf{x}})=0$ are imposed on
$\partial{\Omega}$ with a linear operator $B$. The fundamental solution is
constructed by considering the adjoint equation

$\begin{array}[]{ll}L^{\dagger}[G({\mathbf{x},\boldsymbol{\xi}})]=\delta({\mathbf{x}}-{\boldsymbol{\xi}})\end{array}$
(2)

where $L^{\dagger}$ is the adjoint operator (along with its associated
boundary conditions) and $\delta({\mathbf{x}}-{\bf\xi})$ is the Dirac delta
function. Taking the inner product of (1) with respect to the Green’s function
gives the fundamental solution

$v({\bf x})=(f({\boldsymbol{\xi}}),G(\boldsymbol{\xi},{\bf
x}))=\int\limits_{{\Omega}}G(\mathbf{\xi},{\bf
x})f(\boldsymbol{\xi})d\boldsymbol{\xi},$ (3)

which is valid for any forcing $f({\mathbf{x}})$. Thus once the Green’s
function is computed, the solution for arbitrary forcing functions can be
easily extracted from integration. This integration represents a superposition
of a continuum of delta function forcings that are used to represent
$f({\mathbf{x}})$.

In many modern applications, nonlinearity plays a fundamental role so that the
BVP is of the form

$\begin{array}[]{ll}N[u({\mathbf{x}})]=F({\mathbf{x}})\end{array}$ (4)

where $N[\cdot]$ is a nonlinear differential operator. For this case, the
principle of linear superposition no longer holds and the notion of a
fundamental solution is lost. However, modern deep learning algorithms allow
us the flexibility of learning a coordinate transformation (and their
inverses) of the form

$\displaystyle v$ $\displaystyle=\boldsymbol{\psi}(u),$ (5a) $\displaystyle f$
$\displaystyle=\boldsymbol{\phi}(F),$ (5b)

such that $v$ and $f$ satisfy the linear BVP (1) for which we generated the
fundamental solution (3). This gives a nonlinear fundamental solution through
use of this deep learning transformation.

DeepGreen is a supervised learning algorithm which is ultimately a high-
dimensional interpolation problem (?) for learning the coordinate
transformations $\boldsymbol{\psi}(u)$ and $\boldsymbol{\phi}(F)$. DeepGreen
is enabled by a physics-informed deep autoencoder coordinate transformation
which establishes superposition for nonlinear BVPs, thus enabling a Koopman
BVP framework. The learned Green’s function enables accurate construction of
solutions with new forcing functions in the same way as a linear BVP. We
demonstrate the DeepGreen method on a variety of nonlinear boundary value
problems, including a nonlinear 2D Poisson problem, showing that such an
architecture can be used in many modern and diverse applications in aerospace,
electromagnetics, elasticity, materials, and chemical reactors.

Figure 1: DeepGreen solves nonlinear BVPs by identifying the Green’s Function
of the nonlinear problem using a deep learning approach with a dual
autoencoder architecture. An nonhomogenous linear BVP can be solved using the
Green’s function approach, but a nonlinear BVP cannot. DeepGreen transforms a
nonlinear BVP to a linear BVP, solves the linearized BVP, and then inverse
transforms the linear solution to solve the nonlinear BVP.

##  2 Deep Autoencoders for Linearizing BVPs

Deep AEs have been used to linearize dynamical systems, which are initial
value problems. We extend this idea to BVPs. To be precise, we consider BVPs
of the form

$\displaystyle N[u(\mathbf{x})]$ $\displaystyle=F(\mathbf{x}),$
$\displaystyle\mathbf{x}$ $\displaystyle\in\Omega,$ (6a) $\displaystyle
B[u(\mathbf{x})]$ $\displaystyle=0,$ $\displaystyle\mathbf{x}$
$\displaystyle\in\partial\Omega,$ (6b)

where $\Omega$ is a simply connected open set in $\mathbb{R}^{n}$ with
boundary $\partial\Omega$, $N$ is a nonlinear differential operator,
$F(\mathbf{x})$ is the nonhomogeneous forcing function, $B$ is a boundary
condition, and $u(\mathbf{x})$ is the solution to the BVP. We wish to find a
pair of coordinate transformations of the form (5) such that $v$ and $f$
satisfy a linear BVP

$\displaystyle L[v(\widehat{\mathbf{x}})]$
$\displaystyle=f(\widehat{\mathbf{x}}),$ $\displaystyle\widehat{\mathbf{x}}$
$\displaystyle\in\widehat{\Omega},$ (7a)
$\displaystyle\widehat{B}[v(\widehat{\mathbf{x}})]$ $\displaystyle=0,$
$\displaystyle\widehat{\mathbf{x}}$
$\displaystyle\in\partial\widehat{\Omega},$ (7b)

where $L$ is a linear differential operator, $\widehat{\mathbf{x}}$ is the
spatial coordinate in the transformed domain $\widehat{\Omega}$ with boundary
$\partial\widehat{\Omega}$. Because $L$ is linear, there is a Green’s function
$G(\widehat{\mathbf{x}},\boldsymbol{\xi})$ such that the solution $v$ to the
BVP (7) can be obtained through convolution of the Green’s function and
transformed forcing function

$v(\widehat{\mathbf{x}})=\int\limits_{\widehat{\Omega}}G(\mathbf{\xi},\widehat{\mathbf{x}})f(\boldsymbol{\xi})d\boldsymbol{\xi}.$
(8)

The coordinate transformation along with the Green’s function of the
linearized BVP provide the analog of a Green’s function for the nonlinear BVP
(6). In particular, for a forcing function $F(\mathbf{x})$, the transformed
forcing function is $f=\boldsymbol{\phi}(F)$. The solution to the linearized
BVP can be obtained using the Green’s function $v=\int
G(\boldsymbol{\xi},\widehat{\mathbf{x}})f(\boldsymbol{\xi})d\boldsymbol{\xi}$.
Then the solution to the nonlinear BVP (6) is obtained by inverting the
coordinate transformation $u=\boldsymbol{\psi}^{-1}(v)$ to obtain the solution
to the nonlinear BVP, $u(\mathbf{x})$.

The question that remains is how to discover the appropriate coordinate
transformations $\boldsymbol{\psi}$ and $\boldsymbol{\phi}$. We leverage the
universal approximation properties of neural networks in order to learn these
transformations. In order to use neural networks, we first need to discretize
the BVP. Let $\mathbf{u}$ be a spatial discretization of $u(\mathbf{x})$ and
$\mathbf{F}$ be a discretization of $F(\mathbf{x})$. Then the discretized
version of the BVP (6) is

$\displaystyle\mathbf{N}[\mathbf{u}]$ $\displaystyle=\mathbf{F},$ (9a)
$\displaystyle\mathbf{B}[\mathbf{u}]$ $\displaystyle=\mathbf{0}.$ (9b)

Neural networks $\boldsymbol{\psi}_{u}$ and $\boldsymbol{\phi}_{F}$ are used
to transform $\mathbf{u}$ and $\mathbf{F}$ to the latent space vectors
$\mathbf{v}$ and $\mathbf{f}$

$\displaystyle\mathbf{v}$ $\displaystyle=\boldsymbol{\psi}_{u}(\mathbf{u}),$
(10a) $\displaystyle\mathbf{f}$
$\displaystyle=\boldsymbol{\phi}_{F}(\mathbf{F}),$ (10b)

where $\mathbf{v}$ and $\mathbf{f}$ satisfy the linear equation

$\mathbf{L}\mathbf{v}=\mathbf{f},$ (11)

for some matrix $\mathbf{L}$, which is also learned. In order to learn
invertible transforms $\boldsymbol{\psi}_{u}$ and $\boldsymbol{\phi}_{F}$, we
construct the problem as a pair of autoencoder networks.

In this construction, the transforms $\boldsymbol{\psi}_{u}$ and
$\boldsymbol{\phi}_{F}$ are the encoders and the transform inverses are the
decoders. The network architecture and loss functions are shown in Figure 2.

Figure 2: DeepGreen architecture. Two autoencoders learn invertible coordinate
transformations that linearize a nonlinear boundary value problem. The latent
space is constrained to exhibit properties of a linear system, including
linear superposition, which enables discovery of a Green’s function for
nonlinear boundary value problems.

The neural network is trained using numerous and diverse solutions to the
nonlinear BVP (9), which can be obtained with many different forcings
$\mathbf{F}_{k}$. Consider a dataset comprised of pairs of discretized
solutions and forcing functions
$\\{\mathbf{u}_{k},\mathbf{F}_{k}\\}_{k=1}^{N}$. The loss function for
training the network is the sum of six losses, each of which enforces a
desired condition. The loss functions can be split into three categories:

  1. 1.

Autoencoder losses: We wish to learn invertible coordinate transformations
given by equations (10a) and (10b). In order to do so, we use two
autoencoders. The autoencoder for $\mathbf{u}$ consists of an encoder
$\boldsymbol{\psi}_{u}$ which performs the transformation (10a) and a decoder
$\boldsymbol{\psi}_{u}^{-1}$ which inverts the transformation. In order to
enforce that the encoder and decoder are inverses, we use the autoencoder loss

$\mathcal{L}_{1}=\frac{1}{N}\sum_{k=1}^{N}\frac{\left\lVert\mathbf{u}_{k}-\boldsymbol{\psi}_{u}^{-1}\circ\boldsymbol{\psi}_{u}(\mathbf{u}_{k})\right\rVert_{2}^{2}}{\left\lVert\mathbf{u}_{k}\right\rVert_{2}^{2}}.$
(12)

Similarly, there is an autoencoder for $\mathbf{F}$ where the encoder
$\boldsymbol{\phi}_{F}$ performs the transformation (10b). This transformation
also has an inverse enforced by the associated autoencoder loss function

$\mathcal{L}_{2}=\frac{1}{N}\sum_{k=1}^{N}\frac{\left\lVert\mathbf{F}_{k}-\boldsymbol{\phi}_{F}^{-1}\circ\boldsymbol{\phi}_{F}(\mathbf{F}_{k})\right\rVert_{2}^{2}}{\left\lVert\mathbf{F}_{k}\right\rVert_{2}^{2}}.$
(13)

  2. 2.

Linearity losses: In the transformed coordinate system, we wish for the BVP to
be linear so that the operator can be represented by a matrix $\mathbf{L}$.
The matrix $\mathbf{L}$ and the encoded vectors $\mathbf{v}$ and $\mathbf{f}$
should satisfy equation (11). This is enforced with the linear operator loss

$\mathcal{L}_{3}=\frac{1}{N}\sum_{k=1}^{N}\frac{\left\lVert\mathbf{f}_{k}-\mathbf{L}\mathbf{v}_{k}\right\rVert_{2}^{2}}{\left\lVert\mathbf{f}_{k}\right\rVert_{2}^{2}}.$
(14)

The major advantage of working with a linear operator is that linear
superposition holds. We use a linear superposition loss in order to further
enforce the linearity of the operator in the latent space

$\mathcal{L}_{4}=\frac{1}{N^{2}}\sum_{j=1}^{N}\sum_{i=1}^{N}\frac{\left\lVert(\mathbf{f}_{i}+\mathbf{f}_{j})-\mathbf{L}(\mathbf{v}_{i}+\mathbf{v}_{j})\right\rVert_{2}^{2}}{\left\lVert\mathbf{f}_{i}+\mathbf{f}_{j}\right\rVert_{2}^{2}}.$
(15)

  3. 3.

Cross-mapping losses: The losses described above are theoretically sufficient
to find coordinate transformations for $\mathbf{u}$ and $\mathbf{F}$ as well
as a linear operator $\mathbf{L}$. However, in practice the two autoencoders
were not capable of generating the Green’s function solution. To rectify this,
we add two “cross-mapping” loss functions that incorporate parts of both
autoencoders. The first cross-mapping loss enforces the following mapping from
$\mathbf{u}$ to $\mathbf{F}$. First, one of the solutions from the dataset
$\mathbf{u}_{k}$ is encoded with $\boldsymbol{\psi}_{u}$. This is an
approximation for $\mathbf{v}_{k}$. This is then multiplied by the matrix
$\mathbf{L}$, giving an approximation of $\mathbf{f}_{k}$. Then the result is
decoded with $\boldsymbol{\phi}_{F}^{-1}$. This gives an approximation of
$\mathbf{F}_{k}$. The $\mathbf{u}$ to $\mathbf{F}$ cross-mapping loss is given
by the formula

$\mathcal{L}_{5}=\frac{1}{N}\sum_{k=1}^{N}\frac{\left\lVert\mathbf{F}_{k}-\boldsymbol{\phi}_{F}^{-1}\circ\mathbf{L}\circ\boldsymbol{\psi}_{u}(\mathbf{u}_{k})\right\rVert_{2}^{2}}{\left\lVert\mathbf{F}_{k}\right\rVert_{2}^{2}}.$
(16)

We can similarly define a cross-mapping from $\mathbf{F}$ to $\mathbf{u}$. For
a forcing function $\mathbf{F}_{k}$ from the dataset, it is encoded with
$\boldsymbol{\phi}_{F}$, multiplied by the Green’s function
($\mathbf{G}=\mathbf{L}^{-1}$), and then decoded with
$\boldsymbol{\psi}_{u}^{-1}$ to give an approximation of $\mathbf{u}_{k}$. The
$\mathbf{F}$ to $\mathbf{u}$ cross-mapping loss is

$\mathcal{L}_{6}=\frac{1}{N}\sum_{k=1}^{N}\frac{\left\lVert\mathbf{u}_{k}-\boldsymbol{\psi}_{u}^{-1}\circ\mathbf{L}^{-1}\circ\boldsymbol{\phi}_{F}(\mathbf{F}_{k})\right\rVert_{2}^{2}}{\left\lVert\mathbf{u}_{k}\right\rVert_{2}^{2}}.$
(17)

Note that this final loss function gives the best indication of the
performance of the network to solve the nonlinear BVP (9) using the Green’s
function. The strategy for solving (9) for a given discrete forcing function
$\mathbf{F}$ is to encode the forcing function to obtain
$\mathbf{f}=\boldsymbol{\phi}_{F}(\mathbf{F})$, apply the Green’s function as
in equation (8) to obtain $\mathbf{v}$, and then decode this function to get
the solution $\mathbf{u}=\boldsymbol{\psi}_{u}^{-1}(\mathbf{v})$. The discrete
version of the convolution with the Green’s function given in equation (8) is
multiplication by the matrix $\mathbf{L}^{-1}$.

For the encoders $\boldsymbol{\phi}$ and $\boldsymbol{\psi}$ and decoders
$\boldsymbol{\phi}^{-1}$ and $\boldsymbol{\psi}^{-1}$, we use a residual
neural network (ResNet) architecture (?). The ResNet architecture has been
successful in learning coordinate transformations for physical systems (?) and
is motivated by near-identity transformations in physics. The linear operator
$\mathbf{L}$ is constrained to be a real symmetric matrix and therefore is
self-adjoint. Additionally, $\mathbf{L}$ is initialized as the identity
matrix. Therefore, $\mathbf{L}$ is strictly diagonally dominant for at least
the early parts of training which guarantees $\mathbf{L}$ is invertible and
well-conditioned. For more information on the network architecture and
training procedure, see Appendix B.

##  3 Results

The DeepGreen architecture, which is highlighted in Fig. 2 and whose detailed
loss functions are discussed in the last section, is demonstrated on a number
of canonical nonlinear BVPs. The first three BVPs are one-dimensional systems
and the final one is a two-dimensional system. The nonlinearities in these
problems do not allow for a fundamental solution, thus recourse is typically
made to numerical computations to achieve a solution. DeepGreen, however, can
produce a fundamental solution which can then be used for any new forcing of
the BVP.

###  3.1 Cubic Helmholtz

The architecture and methodology is best illustrated using a basic example
problem. The example problem uses a nonhomogeneous second-order nonlinear
Sturm–Liouville model with constant coefficients and a cubic nonlinearity,
thus making it a cubic Helmholtz equation. The differential equation is given
by

$\displaystyle u^{\prime\prime}+\alpha u+\epsilon u^{3}=F(x),$ (18a)
$\displaystyle u(0)=u(2\pi)=0,$ (18b)

where $u=u(x)$ is the solution when the system is forced with $F(x)$ with
$x\in(0,2\pi)$, $\alpha=-1$ and $\epsilon=-0.3$. The notation
$u^{\prime\prime}$ denotes $\frac{d^{2}}{dx^{2}}u(x)$. The dataset contains
discretized solutions and forcings,
$\\{\mathbf{u}_{k},\mathbf{F}_{k}\\}_{k=1}^{N}$. The forcing functions used
for training are cosine and Gaussian functions; details of data generation and
the forcing functions are provided in Appendix A.

Figure 3: Learning curve. This is a typical learning curve for the DeepGreen
architecture. The vertical dashed line indicates where the training procedure
transitions from autoencoders-only (only $\mathcal{L}_{1}$ and
$\mathcal{L}_{2}$) to a full-network training procedure (all losses). Figure
4: Latent space representations $\mathbf{v}_{k}$ and $\mathbf{f}_{k}$. The
autoencoder transformation $\boldsymbol{\psi}_{u}$ encodes $\mathbf{u}_{k}$ to
the latent space, producing the vector $\mathbf{v}_{k}$ (orange). The forcing
vector $\mathbf{F}_{k}$ is transformed by $\boldsymbol{\psi}_{F}$ to the
encoded vector $\mathbf{f}_{v}$ (blue). Figure 5: Visualized operator and
Green’s function. Discovered Green’s function $\mathbf{G}=\mathbf{L}^{-1}$ and
corresponding linear operator $\mathbf{L}$.

The data is divided into three groups: training, validation, and test. The
training and validation sets are used for training the model. The test set is
used to evaluate the results. The training set contains $N_{train}=8906$
vector pairs $\mathbf{u}_{k}$ and $\mathbf{F}_{k}$. The validation set
contains $N_{validation}=2227$ and test set contains $N_{test}=1238$.

####  3.1.1 Training the Model

The autoencoders used in this example are constructed with fully connected
layers. In both autoencoders, a ResNet-like identity skip connection connects
the input layer to the layer before dimension reduction in the encoder, and
the first full-dimension layer in the decoder with the final output layer (see
Figure 14).

The model is trained in a two-step procedure. First, the autoencoders are
trained, without connection in the latent space, to condition the networks as
autoencoders. In this first phase, only the autoencoder loss functions listed
in Figure 2 are active ($\mathcal{L}_{1}$ and $\mathcal{L}_{2}$). After a set
number of epochs, the latent spaces are connected by an invertible matrix
operator, $\mathbf{L}$, and the remaining 4 loss functions in Figure 2 become
active ($\mathcal{L}_{3}$–$\mathcal{L}_{6}$). In the final phase of training,
the autoencoder learns to encode a latent space representation of the system
where properties associated with linear systems hold true, such as linear
superposition.

Figure 3 shows a typical training loss curve. The vertical dashed line
indicates the transition between the two training phases. The models in this
work are trained for 75 epochs in the first autoencoder-only phase and 2750
epochs in the final phase. The first-phase epoch count was tuned empirically
based on final model performance. The final phase epoch count was selected for
practical reasons; the training curve tended to flatten around 2750 epochs in
all of our tested systems. The autoencoder latent spaces are critically
important. The latent space is the transformed vector space where linear
properties (e.g. superposition) are enforced which enables the solution of
nonlinear problems. In the one-dimensional problems, the latent spaces vectors
$\mathbf{v}$ and $\mathbf{f}$ are in $\mathbb{R}^{20}$.

The latent spaces did not have any obvious physical interpretation, and
qualitatively appeared similar to the representations shown in Figure 4. We
trained 100 models to check the consistency in the learned model and latent
space representations, but discovered the latent spaces varied considerably
(see Appendix C). This implies the existence of an infinity of solutions to
the coordinate transform problem, which indicates further constraints could be
placed on the model.

Figure 6: Model predictions on test data. The top row shows the true solution
$\mathbf{u}_{k}(x)$ and the solution predicted by the network given the
forcing $\mathbf{F}_{k}(x)$ using the Green’s function $\mathbf{G}$. The
bottom row shows the true forcing function $\mathbf{F}_{k}(x)$ compared to the
forcing computed by applying the operator $\mathbf{L}$ to the solution
$\mathbf{u}_{k}$. Three columns show the best, mean, and worst case samples as
evaluated by the sum of normalized $\ell 2$ reconstruction errors.

Despite lacking obvious physical interpretations, the latent space enables
discovery of an invertible operator $\mathbf{L}$ which described the linear
system $\mathbf{L}[\mathbf{v}_{k}]=\mathbf{f}_{k}$. The operator matrix
$\mathbf{L}$ can be inverted to yield the Green’s function matrix
$\mathbf{G}$, which allows computation of solutions to the linearized system
$\mathbf{v}_{k}=\mathbf{G}[\mathbf{f}_{k}]$. An example of the operator
$\mathbf{L}$ and its inverse $\mathbf{G}$ are shown in Figure 5. The operator
and Green’s function shown in Figure 5 display an important prominent feature
seen in all of the results: a diagonally-dominant structure. We initialize the
operator as an identity matrix, but the initialization had little impact on
the diagonally-dominant form of the learned operator and Green’s function
matrices (see Appendix C). The diagonally-dominant operators indicate that the
deep learning network tends to discover a coordinate transform yielding a
nearly-orthonormal basis, which mirrors the common approach of diagonalization
in spectral theory for Hermitian operators. Furthermore, diagonally-dominant
matrices guarantee favorable properties for this application such as being
well-conditioned and non-singular.

We emphasize that training parameters and model construction choices used in
this work were not extensively optimized. We expect the model performance can
be improved in a myriad of ways including extending training times, optimizing
model architecture, modifying the size of the latent spaces, restricting the
form of the operator, and applying additional constraints to the model.
However, these topics are not the main scope of the present work; our focus is
to illustrate the use of autoencoders as a coordinate transform for finding
solutions to nonlinear BVPs.

####  3.1.2 Evaluating the Model

The goal for this model is to find a Green’s function $\mathbf{G}$ for
computing solutions $\mathbf{u}_{k}$ to a nonlinear BVP governed by (6) for a
given forcing function $\mathbf{F}_{k}$. Similarly, we can estimate the
forcing term, $\mathbf{F}_{k}$, given the solution $\mathbf{u}_{k}$. The model
is consequently evaluated by its ability to use the learned Green’s function
and operator for predicting solutions and forcings, respectively, for new
problems from a withheld test data set.

Recall the original model is trained on data where the forcing function is a
cosine or Gaussian function. As shown in Figure 6, the model performs well on
withheld test data where the forcing functions are cosine or Gaussian
functions, producing a cumulative loss around $10^{-4}$. The solutions
$\mathbf{u}_{k}$ and forcing $\mathbf{F}_{k}$ are depicted for the best, mean,
and worst samples scored by cumulative loss.

It’s important to note the test data used in Figure 6 is similar to the
training and validation data. Because ML models typically work extremely well
in interpolation problems, it is reasonable to expect the model to perform
well on this test data set.

As an interesting test to demonstrate the ability of the model to extrapolate,
we prepared a separate set of test data
$\\{\mathbf{u}_{k},\mathbf{F}_{k}\\}_{k=1}^{N}$ containing solutions where
$\mathbf{F}_{k}$ are cubic polynomial forcing functions. This type of data was
not present in training, and provides some insight into the generality of the
learned linear operator and Green’s function matrices. Figure 7 shows examples
of how the model performs on these cubic polynomial-type forcing functions.
Similar to Figure 6, the best, mean, and worst samples are shown as graded by
overall loss.

Figure 7: Model predictions on cubic Helmholtz forced system. The top row
shows the true solution $\mathbf{u}_{k}(x)$ and the solution predicted by the
network given the forcing $\mathbf{F}_{k}(x)$ using the Green’s function
$\mathbf{G}$. The bottom row shows the true forcing function
$\mathbf{F}_{k}(x)$ compared to the forcing computed by applying the operator
$\mathbf{L}$ to the solution $\mathbf{u}_{k}$. Three columns show the best,
mean, and worst case samples as evaluated by the sum of normalized $\ell 2$
reconstruction errors. Figure 8: Model performance summary. Distribution of
loss values are shown for every sample in the test data set. Model loss
functions are minimized during training, making them a natural metric to use
for summarizing performance.

Figures 6 and 7 provide some qualitative insight into the model’s performance
on specific instances selected from the pool of evaluated data. A quantitative
perspective of the model’s performance is presented in Figure 8. This box plot
shows statistics (median value, $Q_{1}$, $Q_{3}$, and range) for four of the
loss functions evaluated on the similar (cosine and Gaussian) test data. Note
the superposition loss function is _not_ scored in this plot because the
superposition loss function can only be evaluated within a single batch, and
the loss depends on batch size and composition.

In conclusion, the DeepGreen architecture enables discovery of invertible,
linearizing transformations that facilitate identification of a linear
operator and Green’s function to solve nonlinear BVPs. It is tested on data
similar and dissimilar to the training data, and evaluated on the loss
functions that guide the training procedure. The discovered operator and
Green’s function take on a surprisingly diagonally-dominant structure, which
hints at the model’s preference to learn an optimal basis. The model appears
to extrapolate beyond the test data, suggesting that the learned operator is
somewhat general to the system.

###  3.2 Nonlinear Sturm–Liouville and Biharmonic Operators

In addition to the example system described above, the approach was applied to
two other one-dimensional systems. We used the same training procedure and
forcing functions that were described in Section 3.1. The first is a system
governed by the nonlinear Sturm–Liouville equation

$\displaystyle[-p(x)u^{\prime}]^{\prime}+q(x)(u+\epsilon u^{3})=F(x),$
$\displaystyle u(0)=u(2\pi)=0,$

where $\epsilon=0.4$ controls the extent of nonlinearity, and $p(x)$ and
$q(x)$ are spatially-varying coefficients

$\displaystyle p(x)=0.5\sin(x)-3,$ $\displaystyle q(x)=0.6\sin(x)-2,$

with $x\in[0,2\pi]$. The final one-dimensional system is a biharmonic operator
with an added cubic nonlinearity

$\displaystyle[-pu^{\prime\prime}]^{\prime\prime}+q(u+\epsilon u^{3})=F(x),$
$\displaystyle u(0)=u(2\pi)=u^{\prime}(0)=u^{\prime}(2\pi)=0,$

where $p=-4$ and $q=2$ are the coefficients and $\epsilon=0.4$ controls the
nonlinearity. As in the prior example, the forcing functions in the training
data are cosine and Gaussian functions, which are described further in
Appendix A.

Results for all the one-dimensional models, including the cubic Helmholtz
example from Section 3.1, are presented in Table 1. Model performance is
quantitatively summarized by box plots and the Green’s function matrix is
shown for each model.

Nonlinear cubic Helmholtz (constant coefficients) $u^{\prime\prime}+\alpha u+\epsilon u^{3}=F$ $u(0)=u(L)=0$ |  Nonlinear Sturm–Liouville (varying $p(x)$, $q(x)$) $[-pu^{\prime}]^{\prime}\\!+\\!qu\\!+\\!\alpha qu^{3}\\!=\\!F$ $u(0)=u(L)=0$ |  Nonlinear Biharmonic operator (constant coefficients) $-pu^{\prime\prime\prime\prime}\\!+\\!qu\\!+\\!\alpha qu^{3}\\!=\\!F$ $u(0)=u(L)=0$
---|---|---
|  |
Table 1: Summary of results for three one-dimensional models. The models are
provided with the Green’s function learned by DeepGreen. A summary box plot
shows the relative losses $\mathcal{L}_{1}$, $\mathcal{L}_{2}$,
$\mathcal{L}_{3}$, $\mathcal{L}_{5}$, and $\mathcal{L}_{6}$ for all three
model systems.

Importantly, the learned Green’s function matrices consistently exhibit
diagonally-dominant structure. The losses for the nonlinear cubic Helmholtz
equation and the nonlinear Sturm–Liouville equation are similar which
indicates that spatially-varying coefficients do not make the problem
significantly more difficult for the DeepGreen architecture. In contrast, the
loss for the nonlinear biharmonic equation are about an order of magnitude
higher than the other two systems. This result implies the fourth-order
problem is more difficult than the second-order problems. Also of note is that
the linear operator loss $\mathcal{L}_{3}$ is consistently the highest loss
across all models. Therefore, it is easier for DeepGreen to find invertible
transformations for the solutions and forcing functions than it is to find a
linear operator that connects the two latent spaces.

###  3.3 Nonlinear Poisson Equation

Figure 9: Model predictions for the (a) best and (b) worst examples from test
data with Gaussian and cosine forcings. In both (a) and (b), the top row shows
the true solution $\mathbf{u}(\mathbf{x})$, the predicted solution using the
Green’s function, and the difference between the true and predicted solution.
The bottom row shows the true forcing function $\mathbf{F}(\mathbf{x})$, the
predicted forcing function, and the difference between the true and predicted
forces. In order to account for the difference in scale between
$\mathbf{u}(\mathbf{x})$ and $\mathbf{F}(\mathbf{x})$, the differences are
scaled by the infinity norm of the true solution or forcing function
($\text{Difference}=(\text{True}-\text{Predicted})/||\text{True}||_{\infty}$).

We also tested our method on a two-dimensional system. The two-dimensional
model is a nonlinear version of the Poisson equation with Dirichlet boundary
conditions

$\displaystyle-\nabla\cdot\left[(1+u^{2})\nabla u\right]=F(\mathbf{x}),$
$\displaystyle\mathbf{x}$ $\displaystyle\in\Omega,$ (19a) $\displaystyle u=0,$
$\displaystyle\mathbf{x}$ $\displaystyle\in\partial\Omega,$ (19b)

where $\Omega:=(0,2\pi)\times(0,2\pi)$. Similar to the one-dimensional models,
the forcing functions used to train the model are cosine and Gaussian
functions, the details of which are provided in Appendix A. The sizes of the
data sets are also similar to the one-dimensional data sets. The training data
contains $N_{train}=9806$ vector pairs $\mathbf{u}_{k}$ and $\mathbf{F}_{k}$,
the validation data contains $N_{validation}=2452$, and the test data contains
$N_{test}=1363$.

The network architecture of the encoders and decoders for the two-dimensional
example differs from the one-dimensional examples. Instead of fully connected
layers, convolutional layers were used in the encoders and decoders. However,
we still use a ResNet architecture. Additionally, the latent space vectors are
in $\mathbb{R}^{200}$. Full details on the network architecture can be found
in Appendix B. Note that the method proposed for discovering Green’s functions
allows for any network architecture to be used for the encoders and decoders.
For the one-dimensional example, similar results were obtained using fully
connected and convolutional layers. However, the convolutional architecture
was better in the two-layer case and also allowed for a more manageable number
of parameters for the wider network that resulted from discretizing the two-
dimensional space.

Figure 10: Model predictions for the (a) best and (b) worst examples from test
data with cubic polynomial forcings. In both (a) and (b), the top row shows
the true solution $\mathbf{u}(\mathbf{x})$, the predicted solution using the
Green’s function, and the difference between the true and predicted solution.
The bottom row shows the true forcing function $\mathbf{F}(\mathbf{x})$, the
predicted forcing function, and the difference between the true and predicted
forces. In order to account for the difference in scale between
$\mathbf{u}(\mathbf{x})$ and $\mathbf{F}(\mathbf{x})$, the differences are
scaled by the infinity norm of the true solution or forcing function
($\text{Difference}=(\text{True}-\text{Predicted})/||\text{True}||_{\infty}$).

The operator and Green’s function for the two-dimensional model are similar to
those displayed in shown in Figure 5. The diagonal dominance is even more
prevalent in this case than the one-dimensional example. The model was
evaluated on test data containing cosine and Gaussian forcing functions.
Figure 9a shows the true solution $\mathbf{u}(x)$ and forcing function
$\mathbf{F}(x)$ as well as the network predictions for the example from the
test data for which the model performed the best (i.e. the smallest value of
the loss). The difference between the true and predicted functions is shown in
the right column of Figure 9a and is scaled by the infinity norm of the true
solution or forcing functions. Figure 9b shows similar results but for the
worst example from the test data. In both cases, the model gives a
qualitatively correct solution for both $\mathbf{u}(x)$ and $\mathbf{F}(x)$.
Unsurprisingly, the network struggles most on highly localized forcing
functions and has the highest error in the region where the forcing occurs.

The model was also evaluated on test data that has cubic polynomial forcing
functions, a type of forcing function not found in the training data. The best
and worst examples are shown in Figure 10. Although the model does not perform
as well for test data which is not similar to the training data, the
qualitative features of the predicted solutions are still consistent with the
true solutions. Figure 11 shows a box plot of the model’s performance on the
similar (cosine and Gaussian forcing test data). The results are similar to
the one-dimensional results, and, in fact, better than the biharmonic operator
model.

Figure 11: Two-dimensional Poisson model performance summary. Distribution of
loss values are shown for every sample in the test data set. Model loss
functions are minimized during training, making them a natural metric to use
for summarizing performance.

##  4 Conclusion

We have leveraged the expressive capabilities of deep learning to discover
linearizing coordinates for nonlinear BVPs, thus allowing for the construction
of the fundamental solution or nonlinear Green’s function. Much like the
Koopman operator for time-dependent problems, the linearizing transformation
provides a framework whereby the fundamental solution of the linear operator
can be constructed and used for any arbitrary forcing. This provides a broadly
applicable mathematical architecture for constructing solutions for nonlinear
BVPs, which typically rely on numerical methods to achieve solutions. Our
DeepGreen architecture can achieve solutions for arbitrary forcings by simply
computing the convolution of the forcing with the Green’s function in the
linearized coordinates.

Given the critical role that BVPs play in the mathematical analysis of
constrained physical systems subjected to external forces, the DeepGreen
architecture can be broadly applied in nearly every engineering discipline
since BVPs are prevalent in diverse problem domains including fluid mechanics,
electromagnetics, quantum mechanics, and elasticity. Importantly, DeepGreen
provides a bridge between a classic and widely used solution technique to
nonlinear BVP problems which generically do not have principled techniques for
achieving solutions aside from brute-force computation. DeepGreen establishes
this bridge by providing a transformation which allows linear superposition to
hold. DeepGreen is a flexible, data-driven, deep learning approach to solving
nonlinear boundary value problems (BVPs) using a dual-autoencoder
architecture. The autoencoders discover an invertible coordinate transform
that linearizes the nonlinear BVP and identifies both a linear operator $L$
and Green’s function $G$ which can be used to solve new nonlinear BVPs. We
demonstrated that the method succeeds on a variety of nonlinear systems
including nonlinear Helmholtz and Sturm–Liouville problems, nonlinear
elasticity, and a 2D nonlinear Poisson equation. The method merges the
strengths of the universal approximation capabilities of deep learning with
the physics knowledge of Green’s functions to yield a flexible tool for
identifying fundamental solutions to a variety of nonlinear systems.

## Acknowledgments

SLB is grateful for funding support from the Army Research Office (ARO
W911NF-17-1-0306). JNK acknowledges support from the Air Force Office of
Scientific Research (FA9550-19-1-0011).

## References

  * 1. I. Stakgold, Boundary Value Problems of Mathematical Physics: 2-Volume Set, vol. 29 (Siam, 2000).
  * 2. J.-B. J. Fourier, Théorie analytique de la chaleur (Chez Firmin Didot, 1822).
  * 3. J. D. Jackson, Classical electrodynamics (John Wiley & Sons, 2007).
  * 4. A. Yariv, Quantum Electronics (John Wiley & Sons, 1989).
  * 5. I. Stakgold, M. J. Holst, Green’s functions and boundary value problems, vol. 99 (John Wiley & Sons, 2011).
  * 6. B. Lusch, J. N. Kutz, S. L. Brunton, Deep learning for universal linear embeddings of nonlinear dynamics.  Nature communications 9, 4950 (2018).
  * 7. K. Champion, B. Lusch, J. N. Kutz, S. L. Brunton, Data-driven discovery of coordinates and governing equations.  Proceedings of the National Academy of Sciences 116, 22445–22451 (2019).
  * 8. C. Wehmeyer, F. Noé, Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics.  The Journal of Chemical Physics 148, 241703 (2017).
  * 9. A. Mardt, L. Pasquali, H. Wu, F. Noé, VAMPnets: Deep learning of molecular kinetics.  Nature Communications 9 (2018).
  * 10. N. Takeishi, Y. Kawahara, T. Yairi, Advances in Neural Information Processing Systems (2017), pp. 1130–1140.
  * 11. E. Yeung, S. Kundu, N. Hodas, 2019 American Control Conference (ACC) (2019), pp. 4832–4839.
  * 12. S. E. Otto, C. W. Rowley, Linearly-recurrent autoencoder networks for learning dynamics.  SIAM Journal on Applied Dynamical Systems 18, 558–593 (2019).
  * 13. Q. Li, F. Dietrich, E. M. Bollt, I. G. Kevrekidis, Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the koopman operator.  Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 103111 (2017).
  * 14. C. J. Dsilva, R. Talmon, R. R. Coifman, I. G. Kevrekidis, Parsimonious representation of nonlinear dynamical systems through manifold learning: A chemotaxis case study.  Applied and Computational Harmonic Analysis 44, 759–773 (2018).
  * 15. C. Gin, B. Lusch, J. N. Kutz, S. L. Brunton, Deep learning models for global coordinate transformations that linearize PDEs.  To appear in the European Journal of Applied Mathematics (2020). Preprint Available: arXiv:1911.02710.
  * 16. G. Cybenko, Approximation by superpositions of a sigmoidal function.  Mathematics of Control, Signals, and Systems (MCSS) 2, 303–314 (1989).
  * 17. K. Hornik, M. Stinchcombe, H. White, Universal approximation of an unknown mapping and its derivatives using multilayer feedforward networks.  Neural networks 3, 551–560 (1990).
  * 18. B. O. Koopman, Hamiltonian systems and transformation in Hilbert space.  Proceedings of the National Academy of Sciences 17, 315-318 (1931).
  * 19. I. Mezić, A. Banaszuk, Comparison of systems with complex behavior.  Physica D: Nonlinear Phenomena 197, 101 - 133 (2004).
  * 20. I. Mezić, Spectral properties of dynamical systems, model reduction and decompositions.  Nonlinear Dynamics 41, 309–325 (2005).
  * 21. M. Budišić, I. Mezić, Geometry of the ergodic quotient reveals coherent structures in flows.  Physica D: Nonlinear Phenomena 241, 1255–1269 (2012).
  * 22. I. Mezić, Analysis of fluid flows via spectral properties of the Koopman operator.  Annual Review of Fluid Mechanics 45, 357–378 (2013).
  * 23. S. L. Brunton, B. W. Brunton, J. L. Proctor, J. N. Kutz, Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control.  PLOS ONE 11, 1–19 (2016).
  * 24. P. J. Schmid, Dynamic mode decomposition of numerical and experimental data.  Journal of Fluid Mechanics 656, 5–28 (2010).
  * 25. C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, D. Henningson, Spectral analysis of nonlinear flows.  J. Fluid Mech. 645, 115–127 (2009).
  * 26. J. N. Kutz, S. L. Brunton, B. W. Brunton, J. L. Proctor, Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems (SIAM, 2016).
  * 27. F. Noé, F. Nüske, A variational approach to modeling slow processes in stochastic dynamical systems.  Multiscale Modeling & Simulation 11, 635–655 (2013).
  * 28. F. Nüske, B. G. Keller, G. Pérez-Hernández, A. S. Mey, F. Noé, Variational approach to molecular kinetics.  Journal of chemical theory and computation 10, 1739–1752 (2014).
  * 29. M. O. Williams, I. G. Kevrekidis, C. W. Rowley, A data-driven approximation of the Koopman operator: Extending dynamic mode decomposition.  Journal of Nonlinear Science 25, 1307–1346 (2015).
  * 30. M. O. Williams, C. W. Rowley, I. G. Kevrekidis, A kernel-based method for data-driven Koopman spectral analysis.  Journal of Computational Dynamics 2, 247–265 (2015).
  * 31. S. Klus, F. Nüske, P. Koltai, H. Wu, I. Kevrekidis, C. Schütte, F. Noé, Data-driven model reduction and transfer operator approximation.  Journal of Nonlinear Science 28, 985–1010 (2018).
  * 32. J. N. Kutz, J. L. Proctor, S. L. Brunton, Applied Koopman theory for partial differential equations and data-driven modeling of spatio-temporal systems.  Complexity 2018, 1–16 (2018).
  * 33. J. Page, R. R. Kerswell, Koopman analysis of burgers equation.  Physical Review Fluids 3, 071901 (2018).
  * 34. C. Bishop, Pattern Recognition and Machine Learning (Springer, 2006).
  * 35. I. Goodfellow, Y. Bengio, A. Courville, Deep Learning (MIT Press, 2016).
  * 36. R. Rico-Martinez, I. Kevrekidis, K. Krischer, Nonlinear system identification using neural networks: dynamics and instabilities.  Neural networks for chemical engineers pp. 409–442 (1995).
  * 37. R. Gonzalez-Garcia, R. Rico-Martinez, I. Kevrekidis, Identification of distributed parameter systems: A neural net based approach.  Computers & chemical engineering 22, S965–S968 (1998).
  * 38. S. H. Rudy, J. N. Kutz, S. L. Brunton, Deep learning of dynamics and signal-noise decomposition with time-stepping constraints.  Journal of Computational Physics 396, 483–506 (2019).
  * 39. H. Lange, S. L. Brunton, N. Kutz, From fourier to koopman: Spectral methods for long-term time series prediction.  arXiv preprint arXiv:2004.00574 (2020).
  * 40. Y. Liu, J. N. Kutz, S. L. Brunton, Hierarchical deep learning of multiscale differential equation time-steppers.  arXiv preprint arXiv:2008.09768 (2020).
  * 41. D. Z. Huang, K. Xu, C. Farhat, E. Darve, Learning constitutive relations from indirect observations using deep neural networks.  Journal of Computational Physics 416, 109491 (2020).
  * 42. S. Pan, K. Duraisamy, Physics-informed probabilistic learning of linear embeddings of nonlinear dynamics with guaranteed stability.  SIAM Journal on Applied Dynamical Systems 19, 480–509 (2020).
  * 43. K. He, X. Zhang, S. Ren, J. Sun, Proceedings of the IEEE conference on computer vision and pattern recognition (2016), pp. 770–778.
  * 44. S. Mallat, Understanding deep convolutional networks.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374, 20150203 (2016).
  * 45. A. Logg, G. N. Wells, Dolfin: Automated finite element computing.  ACM Transactions on Mathematical Software 37 (2010).
  * 46. M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, G. N. Wells, The fenics project version 1.5.  Archive of Numerical Software 3 (2015).
  * 47. A. Logg, K.-A. Mardal, G. N. Wells, et al., Automated Solution of Differential Equations by the Finite Element Method (Springer, 2012).

##  Appendix A Data Generation

###  A.1 1D Problems

The data for all of the one-dimensional systems are created using the same
method and forcing functions. Each solution is computed on an evenly-spaced
128-point grid using MATLAB’s bvp5c solver with a relative error tolerance of
$10^{-8}$ and an absolute error tolerance of $10^{-10}$. The forcing functions
$\mathbf{F}_{k}(x)$ are designed to yield a variety of solutions
$\mathbf{u}_{k}$ such that $||\mathbf{u}_{k}||_{2}\simeq 1$.

The training data consists of two types of systems: Gaussian-forced and
cosine-forced systems. The Gaussian-forced systems have forcing functions of
the form

$F_{k}(x)=a\exp\left(\frac{-(x-b)^{2}}{2c^{2}}\right),$

where $a\in\\{-25,-20,-15,-10,-5,5,10,15,20,25\\}$,
$b\in\\{0,2\pi/19,4\pi/19,\dots,2\pi\\}$, and
$c\in\\{0.1,0.3,0.5,\dots,4.9\\}$. The cosine forcing functions are of the
form

$F_{k}(x)=\alpha\cos(\beta x),$

where $\alpha\in\\{1,1.1,1.2,\dots,10\\}$ and
$\beta\in\\{1,1.05,1.10,\dots,5\\}$. This gives a total of 5000 Gaussian-
forced solutions and 7371 cosine-forced solutions. For the cubic Helmholtz
equation and the nonlinear Sturm–Liouville equation with spatially-varying
coefficients, all of the 12371 solutions were within the error tolerance.
However, there were 97 solutions of the nonlinear biharmonic equation that did
not meet the error tolerance and were therefore discarded. Of the remaining
data, 10% are randomly chosen and withheld as test data, 80% are used as
training data, and 20% are used as validation data.

In order to test the ability of the network to generalize, we also have
another test data set that consists of solutions with cubic forcing functions
of the form

$F_{i}(x)=\gamma(x-\pi)^{3},$

where $\gamma\in\\{0.01,0.03,0.05,\dots,0.29\\}$, and cubic forcing functions
of the form

$F_{i}(x)=\gamma(x-\pi)^{3}+\zeta(x-\pi)^{2}+\boldsymbol{\psi},$

where $\gamma\in\\{0.01,0.03,0.05,\dots,0.29\\}$,
$\zeta\in\\{0.01,0.03,0.05,\dots,0.49\\}$, and
$\boldsymbol{\psi}\in\\{-5,-4,-3,\dots,5\\}$. There are a total of 4140
solutions with cubic forcing functions.

###  A.2 2D Problem

The two-dimensional data satisfies the nonlinear Poisson equation (19). The
solutions are computed with a finite element method using the DOLFIN library
(?) of the FEniCS Project (?, ?). The forcing functions are similar to the
one-dimensional data in that there are Gaussian and cosine forcing functions
along with a separate data set of cubic polynomial forcing functions used to
test the ability of the network to generalize. The Gaussian forcing functions
are of the form

$F_{k}(x,y)=a\exp\left(\frac{-(x-b_{x})^{2}-(y-b_{y})^{2}}{2c^{2}}\right),$

where $a\in\\{-25,-20,-15,-10,-5,5,10,15,20,25\\}$,
$b_{x},b_{y}\in\\{\pi/3,2\pi/3,\pi,4\pi/3,5\pi/3\\}$, and
$c\in\\{0.1,0.3,0.5,\dots,4.9\\}$. The cosine forcing functions are of the
form

$F_{k}(x,y)=\alpha\cos(\beta_{x}x)\cos(\beta_{y}y),$

where $\alpha\in\\{1,1.1,1.2,\dots,10\\}$ and
$\beta_{x},\beta_{y}\in\\{1,1.5,2,\dots,5\\}$. The cubic forcing functions are
of the form

$F_{i}(x,y)=\gamma_{x}(x-\pi)^{3}+\gamma_{y}(y-\pi)^{3},$

where $\gamma_{x},\gamma_{y}\in\\{0.01+0.28k/3|k=0,1,2,3\\}$, and cubic
forcing functions of the form

$\begin{split}F_{i}(x,y)=&\gamma_{x}(x-\pi)^{3}+\gamma_{y}(y-\pi)^{3}\\\
&+\zeta_{x}(x-\pi)^{2}+\zeta_{y}(y-\pi)^{2}+\boldsymbol{\psi},\end{split}$

where $\gamma_{x},\gamma_{y}\in\\{0.01+0.28k/3|k=0,1,2,3\\}$,
$\zeta_{x},\zeta_{y}\in\\{0.01,0.07,0.13,0.19,0.25\\}$, and
$\boldsymbol{\psi}\in\\{-5,-4,-3,\dots,5\\}$. There are 6250 solutions with
Gaussian forcing functions, 7371 solutions with cosine forcing functions, and
4416 solutions with cubic forcing functions.

Figure 12: Encoder network architecture for the two-dimensional data. All
convolutional layers use $4\times 4$ kernels with stride size 1, zero-padding,
and ReLU activation functions. All pooling layers are average pooling layers
with pool size 2 and stride size 2. Figure 13: Decoder network architecture
for the two-dimensional data. All transposed convolutional layers use $4\times
4$ kernels with stride size 2, zero-padding, and ReLU activation functions
except for the last layer which has stride size 1. Figure 14: Layer-by-layer
autoencoder architecture for 1D problems.

##  Appendix B Neural Network Implementation Details

The model training procedure is kept constant for all of the examples in this
work. The networks are optimized with an Adam optimizer ($\beta_{1}=0.9$,
$\beta_{2}=0.999$). Every numerical experiment starts by training a set of 20
models for a ‘small’ number of epochs. Each of the 20 models has a randomly
selected learning rate for the Adam optimizer, uniformly selected between
$10^{-2}$ and $10^{-5}$. The initial training period consists of two phases:
autoencoder-only (75 epochs) and full model (250 epochs). The autoencoder-only
phase only enforces the autoencoder losses $\mathcal{L}_{1}$ and
$\mathcal{L}_{2}$ during backpropagation. A checkpoint algorithm is used to
keep track of the model with the lowest overall loss throughout the training
procedure. At the end of the initial period, the best model is selected and
the others are eliminated. The best model is trained for an additional 2500
epochs.

There are two network architectures in this work. The architectures depicted
in Figures 12 and 13 are applied to the two-dimensional nonlinear Poisson BVP.
The architecture depicted in Figure 14, is applied to one-dimensional
problems.

The two architectures have a few training variables in common. Both models use
variance scaling initialization, $\ell_{2}$ regularization
($\lambda=10^{-6}$), and ReLu activation functions for fully connected (1D
architecture) and convolutional (2D architecture) layers. Notably, the two
layers immediately before and after the latent space do not have activation
functions. A normalized mean squared error loss function is used for all of
the loss functions, as described in Section 2. The models are trained in
batches of 64 samples.

The 2D architecture utilizes convolutional layers and pooling layers, as shown
in Figures 12 and 13. All convolutional layers use a kernel size of $4\times
4$. There are differences between the convolutional layers in the encoder and
the convolutional layers in the decoder. The encoder convolutional layers use
a stride size of $1\times 1$ and an increasing number of filters
($8,16,32,64$). The deconvolutional layers use a stride size of $2\times 2$
with decreasing filter size ($32,16,8$). Pooling layers are similar for both
the encoder and decoder with a stride size of $2\times 2$ and a pool size of
$2\times 2$.

##  Appendix C Additional Results

The repeatability of the results and models learned by the DeepGreen
architecture are interesting to study from the perspective of operator
convergence and latent space representations. In both cases, we aim to
investigate the convergence of the model parameters to determine if the
learned latent spaces and operators are unique or non-unique.

###  C.1 Operator Initialization

We repeat the training procedure for DeepGreen with three different
initialization approaches for the operator $L$. Again, we train with data from
the example nonlinear cubic Helmholtz model. This experiment focuses on
comparing the initial values of the operator $L$ with the final values of the
operator at the end of training to determine if the DeepGreen approach tends
to converge to a specific operator construction. The results in Figure 15 show
the initial and final operator for identity-initialized, randomly initialized,
and Toeplitz-initialized operator matrices. Impressively, the result shows
that the network tends to learn operators with diagonal dominance for all of
the tested initialization strategies. This approach, which DeepGreen appears
to prefer, draws strong parallels to the coordinate diagonalization approach
commonly used in physics.

Figure 15: Initial vs learned operators for an operator matrix $L$ for
different initial conditions. The top row shows identity matrix
initialization, the middle row shows random initialization (He normal), and
the bottom row shows a Toeplitz gradient initialization.

###  C.2 Latent Space Analysis

We repeat the training procedure for the example system, the nonlinear cubic
Helmholtz model, a total of one hundred times. A single sample was selected
from the training data and the latent space representation, $\mathbf{v}_{i}$
and $\mathbf{f}_{i}$, of the input vectors $\mathbf{u}_{i}$ and
$\mathbf{F}_{i}$ are computed. Statistics for the latent space representations
are presented in Figure 16. It is evident that the latent space vectors are
not identical between runs, and that the values in the vector do not follow
any particular statistical distribution. This information implies that the
learned weights in the model, and the learned latent space representations,
vary for each training instance and do not appear to converge to a single
representation.

Figure 16: Statistics of latent space values of a single sample over 100
experimental runs.

###  C.3 Residual network architecture

All of the autoencoders used in this work use a residual network (ResNet)
architecture. In order to demonstrate the advantage of the ResNet
architecture, we trained six models using the DeepGreen architecture for each
of the four systems. Three of the models use the ResNet skip connections,
while three do not use the ResNet architecture.

For the two simplest systems, the nonlinear cubic Helmholtz equation and the
nonlinear Sturm–Liouville equation, the difference between the models with and
without the ResNet skip connections was negligible. For the nonlinear cubic
Helmholtz equation, the mean validation loss for the non-ResNet models was
$2.7\times 10^{-3}$ and the median validation loss was $2.4\times 10^{-3}$.
Using the ResNet architecture resulted in a mean validation loss of $3.5\times
10^{-3}$ and a median validation loss of $8.8\times 10^{-4}$. The ResNet
architecture resulted in a lower median validation loss but a higher mean due
to one of the three models performing much more poorly than the other two. The
results for the nonlinear Sturm–Liouville system are analogous. With a non-
ResNet architecture, the mean validation loss was $4.5\times 10^{-}3$ and the
median validation loss was $4.0\times 10^{-3}$. With a ResNet architecture,
the mean validation loss was $5.7\times 10^{-}3$ and the median validation
loss was $3.1\times 10^{-3}$. Therefore, the use of the ResNet architecture
produced similar results to a non-ResNet architecture for these two simple
systems.

For the two systems that had larger losses — the nonlinear biharmonic equation
in 1D and the 2D nonlinear Poisson equation — the ResNet architecture was
clearly superior to a non-ResNet architecture. For the nonlinear biharmonic
equation, the ResNet architecture yields a mean validation loss of $2.5\times
10^{-2}$ and median validation loss of $2.8\times 10^{-2}$ for the three
models, compared with $3.8\times 10^{-2}$ and $4.0\times 10^{-2}$,
respectively, for the non-ResNet architecture. Therefore, the ResNet
architecture performed better in terms of both the mean and median loss. The
ResNet architecture is absolutely vital for the nonlinear Poisson system.
Without the ResNet architecture, the model essentially did not converge. Both
the mean and median validation losses were $1.9\times 10^{0}$. In contrast,
the ResNet architecture had a mean validation loss of $1.8\times 10^{-2}$ and
a median of $1.9\times 10^{-3}$.