# Structure-Informed Neural Networks for Boundary Observation Problems

Jakub Horsky, Andrew Wynn

(April 2023)

###### Abstract

We introduce Structure Informed Neural Networks (SINNs), a novel method for
solving boundary observation problems involving PDEs. The SINN methodology is
a data-driven framework for creating approximate solutions to internal
variables on the interior of a domain, given only boundary data. The key idea
is to use neural networks to identify a co-ordinate transformation to a latent
space, upon which a well-posed elliptic system of partial differential
equations is constructed. The use of elliptic systems enables the low-cost
transfer of information from the domain’s boundary to its interior. This
enables approximate solutions to PDE boundary observation problems to be
constructed for generic, and even ill-posed, problems. A further advantage of
the proposed method is its ability to be trained on experimental or numerical
data without any knowledge of the underlying PDE. We demonstrate the ability
of SINNs to accurately solve boundary observation problems by considering two
challenging examples of a non-linear heat equation and boundary observation
for the Navier-Stokes equations.

Key words: Data driven scientific computing, Reduced order modeling, Machine
learning, Partial differential equations, Operator learning

##  1 Introduction

Boundary observation problems aim to discover the value of a physical quantity
inside a domain by using only observations from its boundary. If possible,
this means that potentially complex physical information can be obtained
without the need for invasive internal sensors. Many fundamental problems in
engineering and physics have this form with applications, for example, in
fluid mechanics [4], medical imaging [16], geophysics [15], and thermal
sensing [2].

Typically, the internal physical quantity of interest is linked to the
boundary observations by a partial differential equation (PDE). In many
applications, this can make the problem highly challenging to analyse
analytically and computationally impractical to solve numerically. In fluid
mechanics, for example, the nonlinear Navier-Stokes equations govern the
relation between the internal fluid properties, such as its velocity or
temperature, and boundary data which are convenient to observe experimentally,
such as pressure or shear stress. The well-known complexity of solutions to
such nonlinear PDEs implies that solving boundary observation problems of
practical importance is a significant challenge.

In this paper, we propose a new data-driven methodology, called Structure
Informed Neural Networks (SINNs), for solving boundary observation problems
involving nonlinear PDEs. The idea is to embed an inherently well-posed
structure for boundary observation problems into a data-driven framework with
the aim to enable efficient, low-order, approximate solutions. This is
achieved in a three-stage process, indicated schematically in Figure 1. First,
a neural network encodes both the boundary data and the structure of the
boundary geometry into a simpler latent space of boundary variables.
Information is then passed from the boundary to the interior of the latent
space using an elliptic system [3]. This embeds a general class of well-posed
PDEs into the SINN. Finally, a second neural network is used to decode the
interior latent to physical variables.

The idea of using elliptic systems in a data-driven approach is the main
novelty of this paper. Boundary value problems for elliptic systems were
widely studied in the “golden age” of PDE analysis in the 1950s [11]. Our
motivation for using them now is that they can describe a significant range of
boundary value problems, are numerically tractable to solve, and can be
defined with only small number of parameters. The second major contribution of
this paper is to develop an operator-theoretic framework for embedding
elliptic systems within the classical encoder-decoder structure of neural
network-based reduced order modelling. This underpins the efficient numerical
identification of SINNs, enables a powerful coupling of elliptic systems with
deep neural networks, and opens the door to the data-driven solution of a wide
range of challenging nonlinear boundary observation problems.

The structure-informed neural networks (SINNs) developed here have some
similarities, and take inspiration from, a number of existing data-driven
methods for PDE analysis. For example, Koopman-based modal decomposition
methods [1, 14, 10] possess the same three-stage mapping structure as in
Figure 1; Physics-Inspired Neural Networks (PINNs) [13] use neural networks to
efficiently solve PDEs, including boundary value problems [7]; and Neural
Operators [5] use a kernel-based Neural Networks to construct solution
operators for PDE parameter identification. To enable a full discussion of the
relation and distinction between SINNs and existing methods in §1.3 we must
first define the mathematical structure of the boundary observations problems
we aim to solve and give an overview of the SINN methodology.

Figure 1: Schematic overview of a Structure-Informed Neural Network (SINN)
operator $\mathcal{F}=\delta\circ\mathcal{L}\circ\epsilon^{\partial}$ for
solving boundary observation problems.

###  1.1 Boundary observation problems

Consider a physical domain $\Omega\subset\mathbb{R}^{d}$, with $d=2$ or $3$,
and let $\partial\Omega$ denote its boundary. At each point
$\boldsymbol{x}\in\Omega$, we want to recover the value of an $n$-dimensional
physical variable $\boldsymbol{u}(\boldsymbol{x})\in\mathbb{R}^{n}$. To do
this, we can only make use of boundary data
$\boldsymbol{b}(\boldsymbol{z})\in\mathbb{R}^{n_{\partial}}$ which can be
measured at each point $\boldsymbol{z}\in\partial\Omega$. It will be assumed
that both interior and boundary data are square-integrable functions in the
sense that $\boldsymbol{u}\in X$, where

$X=L^{2}(\Omega,\mathbb{R}^{n})=\left\\{f:\Omega\rightarrow\mathbb{R}^{n}:\int_{\Omega}\|f(\boldsymbol{x})\|^{2}d\boldsymbol{x}<\infty\right\\},$

and that $\boldsymbol{b}\in Y$, where

$Y=L^{2}(\partial\Omega,\mathbb{R}^{n_{\partial}})=\left\\{f:\partial\Omega\rightarrow\mathbb{R}^{n_{\partial}}:\int_{\partial\Omega}\|f(\boldsymbol{z})\|^{2}d\boldsymbol{z}<\infty\right\\}.$

A typical situation in which such data arises is if the internal and boundary
data satisfy a PDE of the form

$\begin{split}\mathcal{L}(\boldsymbol{u},\boldsymbol{\lambda})&=0,\qquad\text{in}\;\Omega\\\
\mathcal{B}(\boldsymbol{u})&=\boldsymbol{b},\qquad\text{on}\;\partial\Omega,\end{split}$
(1)

where $\mathcal{L}$ is a differential operator, $\boldsymbol{\lambda}$ are any
parameters, and $\mathcal{B}$ is an output operator linking the interior to
boundary variables. We do not assume that (1) is a well-posed in the sense
that for any boundary function $\boldsymbol{b}$, there is a unique solution
$\boldsymbol{u}$ satisfying the PDE. Instead, the output operator
$\mathcal{B}$ should be viewed simply as shorthand for the “available
information” which may be observed on the boundary, given that a physical
system is in state $\boldsymbol{u}$ inside the domain.

In this abstract language, the structure-informed neural networks (SINNs) that
will be constructed in this paper are operators

$\begin{split}\mathcal{F}:&Y\rightarrow X\\\
&\boldsymbol{b}\mapsto\boldsymbol{u}\end{split}$ (2)

acting between the function space $Y$ of observable measurements and the space
$X$ of all possible distributions of interior physical values. The fact we
will identify operators is important since it means that a single SINN
$\mathcal{F}$ is able to approximate the internal variables $\boldsymbol{u}$,
given any possible boundary observation $\boldsymbol{b}\in X$. As will be
discussed in §1.3, this operator-based philosophy places SINNs within the
recent class of neural-network-based operator identification methods such as
Neural Operators [5] or DeepONets [9].

A key objective of this paper is to identify SINN operators from data. We
assume that a data ensemble

$\mathcal{U}:=\\{(\boldsymbol{u}_{i},\boldsymbol{b}_{i})\\}_{i=1}^{N_{T}}\subset(X\times
Y)^{N_{T}}$

is available consisting of $N_{T}$ pairs of internal and boundary data, each
arising from a solution to (1). The idea will be to construct a mapping of the
form (2) which optimally fits the data $\mathcal{U}$. Given the infinite-
dimensional nature of the underlying problem and the finite-dimensional nature
of the data $\mathcal{U}$, however, any numerically tractable method must make
an a-priori restriction on the possible forms that the $\mathcal{F}$ can take.

###  1.2 SINN operators

We assume that $\mathcal{F}$ is the composition of three operators

$\mathcal{F}=\delta\circ\mathcal{E}_{\mathcal{A}}\circ\epsilon^{\partial},$
(3)

the structure of which is shown schematically in Figure 1. The first operator
is called a boundary encoder $\epsilon^{\partial}:Y\rightarrow Y_{L}$. This is
a nonlinear operator, defined in terms of a neural network, which maps both
the boundary data and geometry into a boundary latent space,
$Y_{L}:=L^{2}(\partial\Omega,\mathbb{R}^{r})$, where the parameter $r$ governs
the order and complexity of the latent space.

To enable data-driven training we must further restrict the form of the
operator $\epsilon^{\partial}$, and we assume that $\epsilon^{\partial}$ acts
semi-locally in the following sense. Given boundary data $\boldsymbol{b}\in
Y$, the value of $(\epsilon^{\partial}\boldsymbol{b})(\boldsymbol{z})$ at any
$\boldsymbol{z}\in\partial\Omega$ can only depend on the values of
$\boldsymbol{b}$ in a small neighbourhood
$N_{\boldsymbol{z}}\subset\partial\Omega$ of $\boldsymbol{z}$. Practically,
this will be achieved by training a neural network111A formal mathematical
definition of neural networks used in this paper is given in §4.5.
$\mathcal{N}^{\partial}:\left\\{\boldsymbol{b}(\boldsymbol{\xi}):\boldsymbol{\xi}\in
N_{\boldsymbol{z}}\right\\}\mapsto(\epsilon^{\partial}\boldsymbol{b})(\boldsymbol{z})$
As will be described in detail in §2.2, the fact that the input to
$\mathcal{N}_{\partial}$ is defined in terms of a local neighbourhood will
enable the use of a single neural network to be repeatably to build up the
definition $\epsilon^{\partial}:Y\rightarrow Y_{L}$. This enables a wide class
of nonlinear operators to be considered without significantly increasing the
number of optimisation parameters.

The purpose of introducing latent variables is to define a common structure
within which information can be passed from the boundary latent space $Y_{L}$
to an interior latent space $X_{L}=L^{2}(\Omega,\mathbb{R}^{r})$. A SINN
implements this transfer of information by using an elliptic system of PDEs.
An elliptic system is governed by a second-order differential operator

$D_{\mathcal{A}}\boldsymbol{\ell}=\sum_{i,j=1}^{d}A_{ij}\frac{\partial^{2}\boldsymbol{\ell}}{\partial
x_{i}\partial x_{j}},$

where $A_{ij}\in\mathbb{R}^{r\times r}$ are symmetric matrices satisfying the
two conditions: i) that $A_{ij}=A_{ji}$, for any $i,j=1,\dots,r$; and ii) that
the block matrix $\mathcal{A}=(A_{ij})\in\mathbb{R}^{rd\times rd}$ is strictly
positive definite.

Information is passed from the boundary latent space $Y_{L}$ to an interior
latent space $X_{L}$ by solving the following boundary problem:

$\begin{split}D_{\mathcal{A}}\boldsymbol{\ell}&=0,\phantom{(b)^{\partial}}\qquad\text{in}\;\Omega\\\
\boldsymbol{\ell}&=\epsilon^{\partial}(\boldsymbol{b}),\qquad\text{on}\;\partial\Omega.\end{split}$
(4)

The assumption that $\mathcal{A}$ is positive definite is crucial. This
implies that (4) is a strongly elliptic system of PDEs. It then follows, under
appropriate smoothness conditions [3] on the latent boundary data
$\epsilon^{\partial}(\boldsymbol{b})$ and the boundary geometry, that (4) has
a unique solution $\boldsymbol{\ell}\in X_{L}$. We let
$\mathcal{E}_{\mathcal{A}}:Y_{L}\rightarrow X_{L}$ denote the operator which
maps boundary data to internal variables when solving the elliptic boundary
value problem (4). The structure of the SINN mapping (3) is hence specifically
designed to create a latent space in which passage of data from boundary to
the interior is well-posed. This is achieved irrespective of the properties of
the PDE or the observation mapping structure (1) from which the physical data
was sampled.

This third, and final, component of a SINN operator (3) is a decoder

$\begin{split}\delta&:X_{L}\longrightarrow X\\\
&\quad\;\;\;\boldsymbol{\ell}\longmapsto\boldsymbol{u}\end{split}$ (5)

which lifts a distribution of interior latent $\boldsymbol{\ell}\in X_{L}$
back into physical space $\boldsymbol{u}\in X$. Analogous to the boundary
encoder, $\delta$ is assumed to be nonlinear and semi-local. That is, for any
$\boldsymbol{x}\in\Omega$, the value of
$(\delta\boldsymbol{\ell})(\boldsymbol{x})$ must only depend on the values of
$\boldsymbol{\ell}$ in a small neighbourhood $N_{\boldsymbol{x}}\subset\Omega$
of $\boldsymbol{x}$. Again, this can be implemented using a single neural
network $\mathcal{N}:\\{\boldsymbol{\ell}(\boldsymbol{y}):\boldsymbol{y}\in
N_{\boldsymbol{x}}\\}\mapsto\boldsymbol{\ell}(\boldsymbol{x})$ which is
applied repeatably to form the definition of the operator $\delta$, as
described in detail in §2.3.

In summary, a structure-informed neural network (SINN) $\mathcal{F}$ is an
operator of the following form

$\mathcal{F}=\left\\{\footnotesize\begin{array}[]{c}\text{semi-local}\\\
\text{nonlinear NN}\\\ \delta:X_{L}\rightarrow
X\end{array}\right\\}\circ\left\\{\footnotesize\begin{array}[]{c}\text{global}\\\
\text{elliptic system}\\\ \mathcal{E}_{\mathcal{A}}:Y_{L}\rightarrow
X_{L}\end{array}\right\\}\circ\left\\{\footnotesize\begin{array}[]{c}\text{semi-
local}\\\ \text{nonlinear NN}\\\ \epsilon^{\partial}:Y\rightarrow
Y_{L}\end{array}\right\\}$

The semi-local architecture of the encoder and decoder mappings is chosen
specifically so as to restruct the number of degrees of freedom involved in
defining the nonlinear components of the operator. The global transfer of
information from boundary to interior is performed in the latent space via a
well-posed elliptic system. This embeds a natural, yet very general, object
into a SINN which is specifically tailored to the structure of the boundary
observation problems that are our aim to solve. Furthermore, as will be
described in §3, a key advantage of using elliptic systems of PDEs is that
identification of their coefficients can be performed in a computationally-
efficient manner using only local training data. However, once trained, the
resulting elliptic system can then be applied globally to give a SINN solution
to the original boundary observation problem.

In §2 we introduce the concept of a generating function which underpins the
semi-local structure of the encoder and decoder operators, before introducing
these operator formally and deriving their inherited mathematical properties.
The method of training SINNs from data is described in §3 and its numerical
implemention discussed in §4. Implementation of our approach on a pair of
challenging test-cases is given in §5. Before this, we first comment briefly
on the relation between the proposed SINN architecture and other, related,
data-driven approaches to PDE analysis.

###  1.3 Relation of SINNs to existing methods

The use of neural networks to solve PDEs has received much recent interest
with the development of Physics-inspired Neural Networks (PINNs) [13]. In the
context of solving a PDE of the form (1), the idea is to view the solution
$\boldsymbol{u}$ as a mapping
$\mathbb{R}^{d}\ni\boldsymbol{x}\mapsto\boldsymbol{u}(\boldsymbol{x})\in\mathbb{R}^{n}$
and to therefore seek to construct a neural network
$\mathcal{N}_{P}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{n}$ which approximates
the solution. The crucial step is to add so-called physics-inspired
constraints, namely
$\mathcal{L}(\mathcal{N}_{P}(\boldsymbol{x}),\boldsymbol{\lambda})_{|_{\Omega}}=0$
and
$\left[\mathcal{B}(\mathcal{N}_{P}(\boldsymbol{x}))-\boldsymbol{b}\right]_{|_{\partial\Omega}}=0$,
to force the constructed solution to satisfy the underlying PDE.

In contrast to the SINN operators $\mathcal{F}:Y\rightarrow X$ which act
between functions spaces, PINNs are finite-dimensional mappings that directly
attempt to replicate the solution mapping
$\boldsymbol{x}\mapsto\boldsymbol{u}(\boldsymbol{x})$. They require knowledge
of the underlying PDE they seek to solve (i.e., of
$\mathcal{L},\boldsymbol{\lambda}$ and $\mathcal{B}$) and, when applied to
boundary observation problems, must be trained using knowledge of the specific
boundary data $\boldsymbol{b}$. In contrast, SINNs do not require such
information: the identification of operators means that such boundary data is
not required in the SINN methodology.

The three-operator structure of the mapping
$\mathcal{F}=\delta\circ\mathcal{E}_{\mathcal{A}}\circ\epsilon^{\partial}$ is
widely used in a variety of data-driven approaches to low-order modelling. In
Koopman-based modelling, for example, operators with this three-level
structure are used to approximate the time-evolution of chaotic, infinite-
dimensional, dynamical systems. In these approaches, the role of the central
operator $\mathcal{E}_{\mathcal{A}}$ is to model temporal evolution, rather
than the passage of information from a domain’s boundary to its interior as in
this paper. The main distinction between the SINNs and the Koopman methodology
is that, in the latter approach, the latent space is finite dimensional and
the temporal operator a finite-dimensional ODE.

This represents an important distinction with the SINN methodology. To
explain, consider the case of the decoder $\delta$ operator, and assume that
it maps from a finite dimensional latent space, say $\mathbb{R}^{r}$, into the
infinite dimensional space of physical variables
$X=L^{2}(\Omega,\mathbb{R}^{n})$. The mismatch in dimensions between latent
and physical space implies that the decoder must have an inherent method of
translating finite-dimensional latent variables to infinite-dimensional
functions. In Koopman-based approaches, this is typically achieved by
considering a basis of functions $\\{\Phi_{i}\\}_{i=1}^{N}\subset X$ and
letting $\delta$ involve a mapping from the latent space $\mathbb{R}^{r}$ to
the coefficients $\\{\hat{f}_{i}\\}\subset\mathbb{R}^{N}$ of a series
expansion $\sum_{i=1}^{N}\hat{f}_{i}\Phi_{i}\in X$. A major challenge of this
approach is the choice of an appropriate basis $\\{\Phi_{i}\\}$ and attempting
to solve this problem has motivated a range of different Koopman-based methods
[17, 6].

In contrast, in the SINN approach developed here, the use of an elliptic
system $\mathcal{E}_{\mathcal{A}}$ removes the need for assigning or
identifying a set of basis functions, and therefore the imposition of
unnecessary structure on the operator $\mathcal{F}$. Constructing an
appropriate elliptic system only requires identifying the PDE coefficient
matrix $\mathcal{A}$, which potentially offers a significant reduction in
dimension to identifying a set of basis function $\\{\Phi_{i}\\}\subset X$.
This advantage comes at the cost of requiring solution of a PDE, as opposed to
an ODE, as the central component of the model $\mathcal{F}$. However, the SINN
methodology deliberately imposes a well-posed elliptic structure which, in
many cases, enables this PDE to be solved at accurately and at low cost using
existing algorithms. In addition, as will be explained in §3, since our aim is
to identify a PDE, the cost function for SINN training can be chosen to
involve only low-cost, local, solutions to elliptic systems during training.
However, once trained, the identified elliptic systems can then be used to
transfer information across a domain globally.

Finally, the philosophy taken in this paper to identify operators using the
SINN methodology is related to the recent interest in using neural networks to
identify operators between function spaces, such as Neural Operators [5] or
DeepONets [9]. The Neural Operator framework [5] seeks to construct solution
operators $G:\boldsymbol{\lambda}\mapsto\boldsymbol{u}$ which solve PDEs of
the form (1) with Dirichlet boundary conditions $\boldsymbol{b}=0$ using
knowledge of their distributed parameters $\boldsymbol{\lambda}$. In this
approach, $G$ transfers information globally in the domain $\Omega$ using an
iterative sequence of integral operators whose kernels are identified using
neural networks. For practicable computational implementation in model
training, structure needs to be imposed in the integral kernels, such as using
low-rank approximations, Convolutional Neural Networks, Graph Neural Networks
[12], or Fourier Neural Operators [8]. Any such choice of structure is
philosophically similar to need to prescribe a functional basis in the
Koopman-based methodology described previously. Again, the contrast to the
SINN methodology is that by training an elliptic operators, only a relative
small number of coefficients are required to enable global transfer of problem
information, and this is achieved without the need to impose any additional
structure on the operator ansatz. We note, finally, that the DeepONet
methodology [9], which can be viewed as a special case of the Neural Operator
approach, also essentially requires the identification of a functional basis
during training.

##  2 Encoders and Decoders for SINNs

In this section, we give a detailed description of the mathematical structure
of the encoder and decoder operators required to create a SINN. We will
describe three classes of operator: interior encoders, boundary encoders, and
decoders. As indicated in Figure 1, only the boundary encoder and decoder are
required to define a SINN mapping. However, as will be explained in §3,
interior encoders will be required to enable data-driven training.

The semi-local structure of all encoder and decoder operators will be
implemented by defining generating functions (GFs), which act as the building
blocks of the SINN methodology. In each of the follow sections we first
introduce a generating function, use it to define the respective operator,
then comment on the regularity properties inherited by that operator.

###  2.1 Interior Encoders

For the purposes of model training only, we will construct interior encoders
$\epsilon$ which, given any distribution of physical variable
$\boldsymbol{u}\in X$, transforms these into a distribution of latent
variables $\boldsymbol{\ell}=\epsilon\boldsymbol{u}$ on the domain interior.

Interior Encoder GFs: Given a compact set $0\in E\subset\mathbb{R}^{d}$, a
generating function for an interior encoder is any continuous, compact222A
compact mapping is one which maps bounded subsets to relatively compact
subsets., generally nonlinear mapping

$e:L^{2}(E,\mathbb{R}^{n})\longrightarrow\mathbb{R}^{r}.$ (6)

This should be thought of as a mapping

$e:\left\\{\footnotesize\begin{array}[]{c}\text{Local patch of}\\\ \text{
interior
data}\end{array}\right\\}\longmapsto\left\\{\footnotesize\begin{array}[]{c}\text{Latent}\\\
\text{variables}\end{array}\right\\}.$

which will be used to endow the interior encoder $\epsilon$ with the desired
semi-local structure.

Definition of Interior Encoders: For any $\boldsymbol{x}\in\Omega$, define a
local neighbourhood

$E_{\boldsymbol{x}}:=\boldsymbol{x}+E=\\{\boldsymbol{x}+\boldsymbol{y}:\boldsymbol{y}\in
E\\},$

and let
$\Omega_{E}:=\\{\boldsymbol{x}\in\Omega:E_{\boldsymbol{x}}\subset\Omega\\}$ be
the set of points whose neighbourhoods $E_{\boldsymbol{x}}$ are entirely
contained in $\Omega$. These sets are shown in Figure 3.

Next let $\boldsymbol{u}\in X$. For any $\boldsymbol{x}\in\Omega_{E}$, a local
function $\boldsymbol{u}_{\boldsymbol{x}}:E\rightarrow\mathbb{R}^{n}$ can be
defined by

$\boldsymbol{u}_{\boldsymbol{x}}(\boldsymbol{y}):=\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{y}),\qquad\boldsymbol{y}\in
E.$

Given a generating function
$e:L^{2}(E,\mathbb{R}^{n})\rightarrow\mathbb{R}^{r}$, we then define an
interior encoder by

$\left(\epsilon\boldsymbol{u}\right)(\boldsymbol{x}):=e\left(\boldsymbol{u}_{\boldsymbol{x}}\right),\qquad\boldsymbol{x}\in\Omega_{E},$
(7)

This definition should be thought of as mapping the physical data
$\boldsymbol{u}$, viewed as a function in $X=L^{2}(\Omega,\mathbb{R}^{n})$, to
a new function $\epsilon\boldsymbol{u}:\Omega_{E}\rightarrow\mathbb{R}^{r}$.
This allows the latent variables
$\boldsymbol{\ell}(\boldsymbol{x})=(\epsilon\boldsymbol{u})(\boldsymbol{x})$
corresponding to $\boldsymbol{u}$ to be defined on the subdomain $\Omega_{E}$.

It follows trivially from its definition that interior encoders are operators
satisfying $\epsilon:L^{2}(\Omega,\mathbb{R}^{n})\rightarrow
L^{2}(\Omega_{E},\mathbb{R}^{r})$. However, the following result shows that
the latent variable field created using the encoder $\epsilon$ are, in fact,
continuous, uniformly bounded, functions.

######  Lemma 1.

Let $e:L^{2}(E,\mathbb{R}^{n})\rightarrow\mathbb{R}^{r}$ be an interior
encoder generating function and let $\epsilon$ be defined by (7). Then
$\epsilon:L^{2}(\Omega,\mathbb{R}^{n})\rightarrow
C(\Omega_{E},\mathbb{R}^{r})$.

###### Proof.

See Appendix 8.1. ∎

Figure 2: Schematic structure of action of an interior generating function
$e:L^{2}(E)\rightarrow\mathbb{R}^{r}$. At a point
$\boldsymbol{x}\in\Omega_{E}$, the local physical data
$\boldsymbol{u}_{|_{E_{x}}}$ is viewed as a function in $L^{2}(E)$. The
generating function then maps local information to the value of the latent
varibles at $\boldsymbol{x}$ via
$e:\boldsymbol{u}_{|_{E_{\boldsymbol{x}}}}\mapsto\boldsymbol{\ell}(\boldsymbol{x})\in\mathbb{R}^{r}$.

###  2.2 Boundary Encoders

We describe how to construct a boundary encoder $\epsilon^{\partial}$ which,
given any distribution of boundary values $\boldsymbol{b}\in Y$, transforms
these into a distribution of latent variables
$\boldsymbol{\ell}=\epsilon^{\partial}\boldsymbol{b}$ on the boundary
$\partial\Omega$. The construction is analogous to that of the interior
encoder in (7) but with the added complication of including information about
the boundary geometry.

Boundary encoder GFs: Given a fixed, compact, set $0\in
E_{\partial}\subset\mathbb{R}^{d-1}$ containing the origin, a generating
function for the boundary encoder is any continuous, compact, and generally
nonlinear function

$e^{\partial}:L^{2}(E_{\partial},\mathbb{R}^{n_{\partial}})\times
L^{2}(E_{\partial},\mathbb{R}^{d})\rightarrow\mathbb{R}^{r}.$ (8)

This should be understood as a mapping

$e^{\partial}:\left\\{\footnotesize\begin{array}[]{c}\text{Section of}\\\
\text{ boundary
data}\end{array}\right\\}\times\left\\{\footnotesize\begin{array}[]{c}\text{Section
of}\\\ \text{ boundary
geometry}\end{array}\right\\}\longmapsto\left\\{\footnotesize\begin{array}[]{c}\text{Latent}\\\
\text{variables}\end{array}\right\\}$

which will be used repeatably to define a semi-local boundary encoder
operator.

Definition of Boundary Encoders: We assume thoughout that $\partial\Omega$ is
sufficiently regular that a normal vector
$\boldsymbol{n}({\boldsymbol{z}})\in\mathbb{R}^{d}$ and a tangent plane
$T_{\boldsymbol{z}}\subset\mathbb{R}^{d-1}$ exists for every
$\boldsymbol{z}\in\partial\Omega$. Each tangent plane $T_{\boldsymbol{z}}$ is
defined in terms of a local-coordinate system with origin at $\boldsymbol{z}$
and whose basis vectors $(\boldsymbol{e}^{\boldsymbol{z}}_{i})_{i=1}^{d-1}$
are orthogonal to $\boldsymbol{n}({\boldsymbol{z}})$. We assume further that
there exists a ball $B_{R}(\boldsymbol{z})\subset\mathbb{R}^{d}$ of radius $R$
such that the local projection $P_{\boldsymbol{z}}:\partial\Omega\cap
B_{R}(\boldsymbol{z})\rightarrow T_{\boldsymbol{z}}$ from the boundary to
tangent plane $T_{\boldsymbol{z}}$ is one-to-one and, in addition, that there
exists $\tau>0$, independent of $\boldsymbol{z}\in\partial\Omega$, such that

$\\{t_{i}\boldsymbol{e}_{i}^{\boldsymbol{z}}:0\leq t_{i}<\tau\\}\subset
P_{\boldsymbol{z}}(B_{R}(\boldsymbol{z})\cap\partial\Omega)\subset
T_{\boldsymbol{z}},\qquad\boldsymbol{z}\in\partial\Omega$ (9)

and we also assume that

$E_{\partial}\subseteq(0,\tau)^{d-1}.$ (10)

Property (10) implies that we can view $E_{\partial}$ as a subset of the
tangent plane, while (9) then implies that a well-defined, continuous, inverse
$P_{\boldsymbol{z}}^{-1}:\\{t_{i}\boldsymbol{e}_{i}^{\boldsymbol{z}}:\boldsymbol{t}\in
E_{\partial}\\}\rightarrow\partial\Omega$ exists. A schematic illustration of
this construction is shown in Figure 3.

This technical construction allows us, for each
$\boldsymbol{z}\in\partial\Omega$, to define a function
$\boldsymbol{b}_{\boldsymbol{z}}:E_{\partial}\rightarrow\mathbb{R}^{n_{\partial}}$,
which depends on the boundary data local to $\boldsymbol{z}$, by

$\boldsymbol{b}_{\boldsymbol{z}}(\boldsymbol{t}):=\boldsymbol{b}\left(P^{-1}_{\boldsymbol{z}}(t_{i}\boldsymbol{e}_{i}^{\boldsymbol{z}})\right),\qquad\boldsymbol{t}=(t_{i})_{i=1}^{d-1}\in
E_{\partial},$ (11)

Similarly, we can also define a function
$\boldsymbol{n}_{\boldsymbol{z}}:E_{\partial}\rightarrow\mathbb{R}^{d}$ which
describes the boundary geometry local to $\boldsymbol{z}$ by

$\boldsymbol{n}_{\boldsymbol{z}}(\boldsymbol{t}):=\boldsymbol{n}\left(P^{-1}_{\boldsymbol{z}}(t_{i}\boldsymbol{e}_{i}^{\boldsymbol{z}})\right),\qquad\boldsymbol{t}=(t_{i})_{i=1}^{d-1}\in
E_{\partial}.$ (12)

Next, using the boundary generating function $e^{\partial}$, the corresponding
boundary encoder is defined by

$\left(\epsilon^{\partial}\boldsymbol{b}\right)(\boldsymbol{z}):=e^{\partial}(\boldsymbol{b}_{\boldsymbol{z}},\boldsymbol{n}_{\boldsymbol{z}}),\qquad\boldsymbol{z}\in\partial\Omega.$
(13)

Similar to the case of the interior encoder, since $e^{\partial}$ is assumed
to be compact and continuous, an analogous proof to that of Lemma 1 implies
that

$\epsilon^{\partial}:L^{2}(\partial\Omega,\mathbb{R}^{n_{\partial}})\rightarrow
C(\partial\Omega,\mathbb{R}^{r}).$

Hence, the boundary latent variables
$\boldsymbol{\ell}_{|_{\partial\Omega}}=\epsilon^{\partial}\boldsymbol{b}$ are
continuous functions.

Figure 3: Schematic structure of the objects explaining the action of a
boundary generating function $e^{\partial}:L^{2}(E_{\partial})\times
L^{2}(E_{\partial})\rightarrow\mathbb{R}^{r}$. Local boundary data
$\boldsymbol{b}_{\boldsymbol{z}}:\boldsymbol{t}\mapsto\boldsymbol{b}(P_{\boldsymbol{z}}^{-1}\boldsymbol{t})$
and geometry
$\boldsymbol{n}_{\boldsymbol{z}}:\boldsymbol{t}\mapsto\boldsymbol{n}(P_{\boldsymbol{z}}^{-1}\boldsymbol{t})$
functions are considered as elements of $L^{2}(E_{\partial})$. The generating
function $e^{\partial}$ maps local boundary and geometric information to the
value of the boundary latent variable at $\boldsymbol{z}\in\partial\Omega$ via
$e^{\partial}:(\boldsymbol{b}_{\boldsymbol{z}},\boldsymbol{n}_{\boldsymbol{z}})\mapsto\boldsymbol{\ell}(\boldsymbol{z})$.

###  2.3 Decoders

We describe how to construct a decoder mapping $\delta$ such that, given a
distribution of latent variables $\ell\in X_{L}$, one can transform these into
a distribution of physical variables $\boldsymbol{u}=\delta\boldsymbol{\ell}$
on the domain interior.

Decoder GFs: Given a compact, symmetric, set $0\in D\subset\mathbb{R}^{d}$, a
decoder generating function is any continuous, compact, and generally
nonlinear mapping

$d:\mathbb{R}^{r}\rightarrow C(D,\mathbb{R}^{n}).$ (14)

This GF should be thought of as follows: given latent variables
$\boldsymbol{\ell}(\boldsymbol{y})\in\mathbb{R}^{r}$ at a point
$\boldsymbol{y}\in\Omega$, then
$d(\boldsymbol{\ell}(\boldsymbol{y}))(\boldsymbol{x})$ gives a local
prediction of the physical variables
$\boldsymbol{u}(\boldsymbol{x})\in\mathbb{R}^{n}$ for any $\boldsymbol{x}\in
D_{\boldsymbol{y}}=\boldsymbol{y}+D$. The idea is to use this map repeatably
to build up a semi-local decoder operator.

We define decoders in two situations, which we refer to as partition decoders
and averaging decoders.

Definition of Partition Decoders: In this case, it is assumed that there exist
points $\\{\boldsymbol{y}_{i}\\}_{i=1}^{N_{d}}\subset\Omega$ such that the
collections of sets $(\boldsymbol{y}_{i}+D)_{i=1}^{N_{d}}$ forms a disjoint
partition of $\Omega$. Now, let $\boldsymbol{\ell}\in X_{L}$ be a latent
variable distribution and let $\boldsymbol{x}\in\Omega$. Due to the assumed
partition property, there is a unique index $j\in\\{1,\dots,N_{d}\\}$ such
that $\boldsymbol{x}\in\boldsymbol{y}_{j}+D$. Consequently,
$\boldsymbol{x}-\boldsymbol{y}_{j}\in D$ and we define a decoded value
$\boldsymbol{u}_{\boldsymbol{\ell}}(\boldsymbol{x})$ by

$(\delta\boldsymbol{\ell})(\boldsymbol{x}):=d(\boldsymbol{\ell}(\boldsymbol{y}_{j}))(\boldsymbol{x}-\boldsymbol{y}_{j}).$

Consequently, we can view $\delta$ as an operator $\delta:X_{L}\rightarrow X$
and we can create an approximation to the physical variables by letting
$\boldsymbol{u}(\boldsymbol{x})=(\delta\boldsymbol{\ell})(\boldsymbol{x})$.

An advantage of using a partition decoder is that if $D$ is be chosen as a
coarse discretization of $\Omega$, then this can reduce the computational cost
of implementing the decoder. However, there are two potential disadvantages of
this choice. First, requiring the decoder to extrapolate from latent to
physical variables over a large set $D$ may introduce approximation errors to
the solution. Second, while $\delta\boldsymbol{\ell}$ is guaranteed to be
square-integrable (as an element of $X$), there is no guarantee that the
resulting physical solution $\delta\boldsymbol{\ell}$ is smooth, or even
continuous. If such a proprerty is desirable, the it is possible to instead
implement the following notion of an averaging decoder.

Definition of Averaging Decoders: Let $\boldsymbol{\ell}\in X_{L}$ be a latent
variable distribution and fix $\boldsymbol{x}\in\Omega$. Now, for any point
$\boldsymbol{y}$ such that $\boldsymbol{x}\in D_{\boldsymbol{y}}$, it follows
from the definition of decoder GFs that a prediction of the physical variables
at $\boldsymbol{x}$ can be obtained using the function
$d(\boldsymbol{\ell}({\boldsymbol{y}}))$. The idea is to average all such
possible predictions. To simplify the resulting expression, note that since
$D$ is symmetric,

$\boldsymbol{x}\in
D_{\boldsymbol{y}}\Leftrightarrow\boldsymbol{x}-\boldsymbol{y}\in
D\Leftrightarrow\boldsymbol{y}-\boldsymbol{x}\in
D\Leftrightarrow\boldsymbol{y}\in D_{\boldsymbol{x}},$

meaning that $\boldsymbol{x}$ can be predicted from any point
$\boldsymbol{y}\in D_{\boldsymbol{x}}\cap\Omega$, as illustrated schematically
in Figure 4, and that the value of the prediction from the point
$\boldsymbol{y}$ at $\boldsymbol{x}$ is
$d(\boldsymbol{\ell}(\boldsymbol{y}))(\boldsymbol{x}-\boldsymbol{y})$.

Figure 4: Schematic structure of the action of a decoder generating function
$d:\mathbb{R}^{r}\rightarrow C(D)$. For any $\boldsymbol{x}\in\Omega$, the
predicted physical value from a point $\boldsymbol{y}\in
D_{\boldsymbol{x}}\cap\Omega$ is
$d(\boldsymbol{\ell}(\boldsymbol{y}))(\boldsymbol{x}-\boldsymbol{y})$. The
decoder operator $\delta$ averages over all such possible values.

Consequently, given a function $\ell\in C(\Omega,\mathbb{R}^{r})$, we define a
decoder mapping $\delta$ by

$(\delta\boldsymbol{\ell})(\boldsymbol{x}):=\frac{1}{|D_{\boldsymbol{x}}\cap\Omega|}\int_{D_{\boldsymbol{x}}\cap\Omega}d(\boldsymbol{\ell}(\boldsymbol{y}))(\boldsymbol{x}-\boldsymbol{y})\,d\boldsymbol{y},\qquad\boldsymbol{x}\in\Omega.$
(15)

The following lemma shows that continuous generating functions create decoders
which themselves produce continuous functions on the entire domain $\Omega$.

######  Lemma 2.

Let $d:\mathbb{R}^{r}\rightarrow C(D,\mathbb{R}^{n})$ be a decoder generating
function and let $\delta$ be defined by (15). Then
$\delta:C(\Omega,\mathbb{R}^{r})\rightarrow C(\Omega,\mathbb{R}^{n})$.

###### Proof.

See Appendix 8.2. ∎

###  2.4 Elliptic Systems

The final component required to define a SINN operator is an elliptic system.
We simply refer to any symmetric, strictly positive definite, matrix

$\mathcal{A}\in\mathbb{S}_{++}^{(rd)^{2}}$ (16)

as a generating function from which an elliptic operator $D_{\mathcal{A}}$ and
the associated boundary value problem (4) can be defined. This then generates
the solution operator $\mathcal{E}_{\mathcal{A}}:Y_{L}\rightarrow X_{L}$.

##  3 Training generating functions

It is worth summarising the constructions developed in §2. Given localisation
sets $G:=(D,E,E_{\partial})$ and local generating functions
$(d,e,e^{\partial},\mathcal{A})$, one can use (7), (13) and (15) to define a
globalisation mapping

$\mathcal{G}_{G}:(d,e,e^{\partial},\mathcal{A})\mapsto(\delta,\epsilon,\epsilon^{\partial},\mathcal{E}_{\mathcal{A}})$

which outputs an interior encoder $\epsilon$, boundary encoder
$\epsilon^{\partial}$, decoder $\delta$ and elliptic system solution operator
$\mathcal{E}_{\mathcal{A}}$. These components can then be combined to give a
SINN operator
$\mathcal{F}=\delta\circ\mathcal{E}_{\mathcal{A}}\circ\epsilon^{\partial}$ in
(3). The aim now is to use the available data ensemble

$\mathcal{U}=\left(\boldsymbol{u}_{j}(x),\boldsymbol{b}_{j}(\boldsymbol{z})\right)_{j=1}^{N_{T}},\qquad\boldsymbol{x}\in\Omega,\boldsymbol{z}\in\partial\Omega,$
(17)

to obtain identify an optimal generating functions and, consequently, optimal
SINN operators.

###  3.1 The cost function for SINN training

Training will be posed as a minimisation problem, and a schematic for the cost
function to be minimised is given by the four-stage process shown in Figure 5.
In the following, it is assumed that
$(\boldsymbol{u},\boldsymbol{b})\in\mathcal{U}$ is a snapshot selected from
the training data ensemble.

Stage 1: Fix a set $\\{\boldsymbol{p}_{i}\\}_{i=1}^{M}\subset\Omega_{E}$ of
training points. At each training point, it is assumed that a training patch
exists, which is defined as the convex hull333The convex hull of a set of
points is the smallest convex subset containing all such points. of a set of
points $\\{\boldsymbol{q}_{ij}\\}_{j=1}^{N}\in\Omega\cup\partial\Omega$ local
to $\boldsymbol{p}_{i}$, an example of which is shown in Figure 5 (a).
Specifically, for each $i$, we assume that there exist points satisfying

  * $(i)$

$\boldsymbol{q}_{ij}\in\Omega_{E}\cup\partial\Omega$, for each $j=1,\dots,N$;

  * $(ii)$

$\boldsymbol{p}_{i}\in\text{int}(Q_{i})$ where
$Q_{i}=\text{conv}\\{\boldsymbol{q}_{ij}:j=1,\dots,N\\}$;

  * $(iii)$

$\boldsymbol{q}_{ij}\in\partial Q_{i}$, for each $j=1,\dots,N$.

For any generating functions $e,e^{\partial}$, assumption (i) implies that we
can compute the latent variables
$\boldsymbol{\ell}_{i}:=e(\boldsymbol{u}_{\boldsymbol{p}_{i}})$ and

$\boldsymbol{\ell}_{ij}:=\left\\{\begin{array}[]{rcl}e(\boldsymbol{u}_{\boldsymbol{q}_{ij}}),&&\text{if}\;\boldsymbol{q}_{ij}\in\Omega;\\\
e^{\partial}(\boldsymbol{b}_{\boldsymbol{q}_{ij}},\boldsymbol{\eta}_{\boldsymbol{q}_{ij}}),&&\text{if}\;\boldsymbol{q}_{ij}\in\partial\Omega.\end{array}\right.,\qquad
j=1,\dots,N.$

Stage 2: Conditions $(ii)$ and $(iii)$ from Stage 1 imply that the convex hull
$Q_{i}\subset\mathbb{R}^{d}$ is a polytope and that each $\boldsymbol{q}_{ij}$
is an exterior point on its boundary $\partial Q_{i}\subset\mathbb{R}^{d-1}$.
Linear interpolation can then be used to obtain a function
$\boldsymbol{f}_{i}\in C(\partial Q_{i},\mathbb{R}^{r})$ satisfying
$\boldsymbol{f}_{i}(\boldsymbol{q}_{ij})=\boldsymbol{\ell}_{ij}$ for each
$j=1,\dots,N$. This process is indicted in Figure 5 (b).

(a) Training Patch

(b) Latent variables on $\partial Q$

(c) Elliptic extension of latent variables to $Q$.

(d) Decoded physical variables at $\boldsymbol{p}$.

Figure 5: Schematic overview of variables involved in a training run. (a) A
training patch with $Q$, containing central point $\boldsymbol{p}$, formed as
the convex hull of exterior points $\\{\boldsymbol{q}_{i}\\}$. (b) Latent
variables computed to form a function $\boldsymbol{f}:\partial
Q\rightarrow\mathbb{R}^{r}$ on the training patch boundary. (c) Elliptic
extension defines latent variables in $Q$, in particular at
$\hat{\ell}(\boldsymbol{p})$. (d) The decoder allows a prediction of the
physical variables at $\boldsymbol{p}$.

Stage 3: Given a generating matrix $\mathcal{A}\in\mathbb{R}^{dr\times dr}$,
define the associated linear elliptic operator $D_{\mathcal{A}}$, and solve
the boundary value problem

$\begin{split}D_{\mathcal{A}}\hat{\boldsymbol{\ell}}&=0\qquad\text{in}\;Q_{i},\\\
\hat{\boldsymbol{\ell}}_{|_{\partial
Q_{i}}}&=\boldsymbol{f}_{i}\qquad\text{on}\;\partial Q_{i},\end{split}$ (18)

on the training patch $Q_{i}$. From this, a predicted value
$\hat{\boldsymbol{\ell}}_{i}:=\hat{\boldsymbol{\ell}}(\boldsymbol{p}_{i})\in\mathbb{R}^{r}$
can be obtained, as shown in Figure 5 (c). It is then natural to define the
error function

$\Psi_{1}((\boldsymbol{u},\boldsymbol{b}),(e,e^{\partial},\mathcal{A})):=\frac{1}{Mr}\sum_{i=1}^{M}\left\|\boldsymbol{\ell}_{i}-\hat{\boldsymbol{\ell}}_{i}\right\|_{2}^{2}.$

which quantifies the error, averaged over all training points
$\\{\boldsymbol{p}_{i}\\}_{i=1}^{M}$, between the encoded latent variables
computed using the generating function $e$, and their predictions from the
boundary $\partial Q_{i}$ using the elliptic system
$\mathcal{E}_{\mathcal{A}}$.

Stage 4. Given a decoder $d:\mathbb{R}^{r}\rightarrow C(D,\mathbb{R}^{n})$, we
create predictions for the physical variables, in the local sets
$D_{\boldsymbol{p}_{i}}$ of points close to the training points, using

$\hat{\boldsymbol{u}}(\boldsymbol{x}):=d(\hat{\boldsymbol{\ell}}_{i})(\boldsymbol{x}-\boldsymbol{p}_{i}),\qquad\boldsymbol{x}\in
D_{\boldsymbol{p}_{i}}.$

A second error function

$\Psi_{2}((\boldsymbol{u},\boldsymbol{b}),(d,e,e^{\partial},\mathcal{A})):=\frac{1}{M|D|}\sum_{i=1}^{M}\int_{D}|\boldsymbol{u}(\boldsymbol{p}_{i}+\boldsymbol{y})-\hat{\boldsymbol{u}}(\boldsymbol{p}_{i}+\boldsymbol{y})|_{2}^{2}d\boldsymbol{y}.$

then quantifies whether the decoder generating function $d$ is able to
accurately recreate the physical data, averaged across a subset of $\Omega$
local to the chosen training points.

The four-stage process described above allows us, for each data point
$(\boldsymbol{u},\boldsymbol{b})\in\mathcal{U}$ to define two functions
$\Psi_{1},\Psi_{2}$ which quantify the error associated with a given quadruple
of generating functions $\mathcal{X}:=(d,e,e^{\partial},\mathcal{A})$. We
therefore define the ensemble cost function by

$\Psi\left(\mathcal{U},\mathcal{X}\right):=\frac{1}{N_{T}}\sum_{j=1}^{N_{T}}\left[\Psi_{1}((\boldsymbol{u}_{j},\boldsymbol{b}_{j}),\mathcal{X})+\alpha\Psi_{2}((\boldsymbol{u}_{j},\boldsymbol{b}_{j}),\mathcal{X})\right],$
(19)

where $\alpha>0$ is a weighting parameter.

##  4 Numerical Implementation

In §2 it was shown how the generating functions
$(e,e^{\partial},d,\mathcal{A})$ act as building blocks for SINNs. The
generating functions $(e,e^{\partial},d)$ and, however, still infinite
dimensional nonlinear functionals. A feasible approach to SINN training must
therefore parameterise the infinite-dimensional generating functions
$(e,e^{\partial},d)$ by mappings between finite-dimensional spaces. The
assumed semi-local structure of these operators will be used to achieve this
in a simple manner, avoiding the need to impose any restrictive structure on
the resulting approximation.

###  4.1 Interior generating functions

The aim is to create an ensemble of possible mappings
$e:L^{2}(E,\mathbb{R}^{n})\rightarrow\mathbb{R}^{r}$ using only finite
dimensional functionals. The first step is to partition
$E\subset\mathbb{R}^{d}$ into a union of $N_{E}$ subsets,

$E=\bigcup_{i=1}^{N_{E}}E_{i}.$ (20)

Then, given any function $\boldsymbol{u}\in L^{2}(E,\mathbb{R}^{n})$, we can
create a vector $\mathcal{D}\boldsymbol{u}\in\mathbb{R}^{n\times N_{E}}$ by
taking its the average of $\boldsymbol{u}$ over each of the $N_{E}$ subsets.
Consequently, any finite-dimensional mapping $\tilde{e}:\mathbb{R}^{n\times
N_{E}}\rightarrow\mathbb{R}^{r}$ can be used to define an encoder
$e:L^{2}(E,\mathbb{R}^{n})\rightarrow\mathbb{R}^{r}$ by forming the
composition $e=\tilde{e}\circ\mathcal{D}$.

###  4.2 Boundary generating functions

The aim is to create boundary generating function
$e^{\partial}:L^{2}(E_{\partial},\mathbb{R}^{n_{\partial}})\times
L^{2}(E_{\partial},\mathbb{R}^{d})\rightarrow\mathbb{R}^{r}$. Using the same
idea as before, partition $E_{\partial}\subset\mathbb{R}^{d-1}$ into
$N_{E_{\partial}}$ sets

$E_{\partial}=\bigcup_{i=1}^{N_{E_{\partial}}}(E_{\partial})_{i}.$ (21)

Then, for any functions $\boldsymbol{f}\in
L^{2}(E_{\partial},\mathbb{R}^{n_{\partial}})$ and $\boldsymbol{\eta}\in
L^{2}(E_{\partial},\mathbb{R}^{d})$ we can form vectors
$\mathcal{D}^{\partial}\boldsymbol{f}\in\mathbb{R}^{n_{\partial}\times
N_{E_{\partial}}}$ and
$\mathcal{D}^{\partial}\boldsymbol{\eta}\in\mathbb{R}^{d\times
N_{E_{\partial}}}$ by averaging these functions over each subset
$(E_{\partial})_{i}$. Consequently, any finite-dimensional map
$\tilde{e}^{\partial}:\mathbb{R}^{(n_{\partial}+d)\times
N_{E_{\partial}}}\rightarrow\mathbb{R}^{r}$ induces a boundary encoder
$e^{\partial}:L^{2}(E_{\partial},\mathbb{R}^{n_{\partial}})\times
L^{2}(E_{\partial},\mathbb{R}^{d})\rightarrow\mathbb{R}^{r}$ via the
composition $e^{\partial}:=\tilde{e}^{\partial}\circ\mathcal{D}^{\partial}$.

A summary of the constructions developed so-far is given in Table 1. It should
be emphasised that due to the semi-local role of the encoder generating
function $e$, it is not important to impose any particular structure on the
partitions (20) or (21). The only requirement is to form a sufficiently
resolved local approximation to the underlying data.

| Interior Encoder | Boundary Encoder
---|---|---
Semi-local | $\begin{array}[]{rcl}\epsilon:L^{2}(\Omega)&\rightarrow&C(\Omega)\\\ \boldsymbol{u}&\stackrel{{\scriptstyle\eqref{eq:interior_encoder}}}{{\mapsto}}&e(\boldsymbol{u}_{\boldsymbol{x}})\end{array}$ | $\begin{array}[]{rcl}\epsilon^{\partial}:L^{2}(\partial\Omega)\times L^{2}(\partial\Omega)&\rightarrow&C(\partial\Omega)\\\ (\boldsymbol{b},\boldsymbol{n})&\stackrel{{\scriptstyle\eqref{eq:bndry_encoder}}}{{\mapsto}}&e^{\partial}(\boldsymbol{b}_{\boldsymbol{z}},\boldsymbol{n}_{\boldsymbol{z}})\end{array}$
GF | $\begin{array}[]{c}e:L^{2}(E)\rightarrow\mathbb{R}^{r}\\\ e=\tilde{e}\circ\mathcal{D}\end{array}$ | $\begin{array}[]{c}e^{\partial}:L^{2}(E_{\partial})\times L^{2}(E_{\partial})\rightarrow\mathbb{R}^{r}\\\ e^{\partial}=\tilde{e}^{\partial}\circ\mathcal{D}^{\partial}\end{array}$
$\begin{array}[]{r}\text{Finite}\\\ \text{Dim.}\end{array}$ | $\tilde{e}:\mathbb{R}^{n\times N_{E}}\rightarrow\mathbb{R}^{r}$ | $\tilde{e}^{\partial}:\mathbb{R}^{(n_{\partial}+d)\times N_{E_{\partial}}}\rightarrow\mathbb{R}^{r}$
Table 1: Structure of the Interior and Boundary Encoders.

###  4.3 Decoder Generating Functions

We aim to create functions $d:\mathbb{R}^{r}\rightarrow C(D,\mathbb{R}^{n})$.
To do this, choose $N_{D}$ points
$\\{d_{i}\\}_{i=1}^{N_{D}}\subset\mathbb{R}^{d}$ whose convex hull contains
$D\in\mathbb{R}^{d}$, and let $\mathcal{I}:\mathbb{R}^{N_{D}}\rightarrow
C(D,\mathbb{R}^{n})$ be any interpolation operator, such as linear
interpolation, which continuously extends known functional values at the
points $d_{i}$ to the whole of $D$. Then, any mapping
$\tilde{d}:\mathbb{R}^{r}\rightarrow\mathbb{R}^{N_{D}}$ induces a decoder
generating function via $d:=\mathcal{I}\circ\tilde{d}$.

###  4.4 Elliptic system solution on training patches

Effective training of interior and boundary encoders requires a choice of
training patches which provide a good sample of both the domain interior and
the domain boundary. For simplicity, we only describe the case of rectangular
domains $\Omega$. In this case, by considering training patches
$Q=\text{conv}\\{\boldsymbol{q}_{i}\\}\subset\bar{\Omega}$ whose boundary
$\partial Q$ is also rectangular, these may be chosen to either lie entirely
in the domain interior, or to coincide with a portion of the boundary
$\partial\Omega$.

After appropriate interpolation of the latent variables
$\ell(\boldsymbol{q}_{i})$ to the rectangular boundary $\partial Q$, standard
second-order finite difference schemes can be used to solve the required
elliptic system (18) on each training patch to obtain the prediction
$\hat{\boldsymbol{\ell}}(\boldsymbol{p})$ at the training point
$\boldsymbol{p}\in Q$.

###  4.5 Cost function minimisation

Suppose that an ensemble $\mathcal{U}$ of training data of the form (17) is
available. Suppose that the latent space dimension $r$, the partition
dimensions $N_{E},N_{E_{\partial}}$ for the encoders, and the decoder
discretisation dimension $N_{D}$ are all given. The problem of model training
is now reduced to finding mappings $\tilde{e}:\mathbb{R}^{n\times
N_{E}}\rightarrow\mathbb{R}^{r},\tilde{e}^{\partial}:\mathbb{R}^{(n_{\partial}+d)\times
N_{E_{\partial}}}\rightarrow\mathbb{R}^{r},d:\mathbb{R}^{r}\rightarrow\mathbb{R}^{N_{D}}$
and a matrix $\mathcal{A}\in\mathbb{R}^{rd\times rd}$ which minimise that for
which the cost function (19) is minimised.

Each of the generating functions $\tilde{e},\tilde{e}^{\partial}$ and $d$ are
assumed to be fully connected feed-forward neural networks. A neural network
(NN) with $L$ layers is a mapping
$\mathcal{N}:\mathbb{R}^{m_{I}}\rightarrow\mathbb{R}^{m_{O}}$ formed by
repeated composition of a prescribed nonlinear activation function
$\phi:\mathbb{R}\rightarrow\mathbb{R}$ and a series of affine maps
$A_{i}:\boldsymbol{x}\mapsto W_{i}\boldsymbol{x}+\boldsymbol{b}_{i}$, for
$i=1,\dots,L$. Here, $W_{i}\in\mathbb{R}^{r_{i}\times r_{i-1}}$ and
$\boldsymbol{b}_{i}\in\mathbb{R}^{r_{i}}$ are the free parameters of the
neural network, $r_{i}$ are the number of neurons in the $i^{\text{th}}$
layer, and $r_{0}=m_{I},r_{L}=m_{O}$. The output of the $i^{\text{th}}$ layer
of the network is given by
$\mathcal{N}_{i}(\boldsymbol{x}):=W_{i}\phi(\mathcal{N}_{i-1}(\boldsymbol{x}))+\boldsymbol{b}_{i}$,
$i\geq 2$, where $\phi$ acts component-wise. For an input
$\boldsymbol{x}\in\mathbb{R}^{m_{I}}$, and letting
$\mathcal{N}_{1}\boldsymbol{x}=W_{1}\boldsymbol{x}+\boldsymbol{b}_{1}$, the
output of the neural network, after iterating its $L$ layers, is
$\mathcal{N}(\boldsymbol{x})=\mathcal{N}_{L}(\boldsymbol{x})$. In this paper,
we fix

$\phi(x)=\text{ReLU}(x)=\left\\{\begin{array}[]{rcl}x,&&x\geq 0,\\\
0&&x<0,\end{array}\right.$

meaning that the tunable parameters of the considered neural networks are
$\Theta=(W_{i},\boldsymbol{b}_{i})_{i=1}^{L}$. We then write
$\mathcal{N}=\mathcal{N}_{\Theta}$ to emphasise this dependency.

The aim of model training is to find a quadruple
$\mathcal{X}:=(e,e^{\partial},d,\mathcal{A})$ which minimises the modelling
residual $\Psi(\mathcal{U},\mathcal{X})$ defined in (19). To impose the
positive definiteness constraint (18) on $\mathcal{A}$, we use a logarithmic
barrier function and instead solve the optimisation problem

$\begin{split}\min_{\Theta,\mathcal{A}}\quad&\Psi(\mathcal{U},\mathcal{X})-\rho\log{(\det{\mathcal{A}})}\\\
\quad&\Theta=(\Theta_{e},\Theta_{e^{\partial}},\Theta_{d}),\\\
\quad&\mathcal{A}\in\mathbb{R}^{dr\times dr},\\\
\quad&e=\tilde{e}\circ\mathcal{D},\;e^{\partial}:=\tilde{e}^{\partial}\circ\mathcal{D}^{\partial},\;d:=\mathcal{I}\circ\tilde{d},\\\
\;&\tilde{e}=\mathcal{N}_{\Theta_{e}},\;\tilde{e}^{\partial}=\mathcal{N}_{\Theta_{e^{\partial}}},\;d=\mathcal{N}_{\Theta_{d}}\\\
\quad&\rho>0,\alpha>0.\end{split}$ (22)

The input and output dimensions of each neural network in (22) are prescribed
by the encoder and decoder structure described in §4.1, §4.2, and §4.3. The
number of layers and neurons in each neural network is problem dependent and
will be specified in §5 for the particular numerical examples considered. The
$-\log{(\det{\mathcal{A}})}$ term provides a barrier in the sense that
$-\log{(\det{\mathcal{A}})}\rightarrow\infty$ as $\det{A}\rightarrow 0$, which
penalises sign-changes of the eigenvalues of $\mathcal{A}$, and hence promotes
positive definiteness.

The optimisation problem (22) is solved using a stochastic gradient descent
approach. At each iteration, a random set of training patches of the form
described in §4.4 are created. The gradient of cost function (19), dependent
upon the chosen training patches, is then computed using automatic
differentiation an appropriate step in the decision variables is taken. The
Adam gradient-based optimization algorithm implemented in the TensorFlow
package is used to identify a local minimum. Since solutions to the elliptic
system (4) are invariant upon rescaling of $\mathcal{A}$, after each iteration
the elliptic decision variables are updated via
$\mathcal{A}\mapsto\mathcal{A}/(\text{det}(\mathcal{A}))^{(1/(rd)^{2})}$ to
maintain the value of the determinant to be unity. The process of random
training patch selection and gradient-based weight updates is then iterated
until the cost $\Psi$ has converged to a local minimum.

##  5 Numerical Examples

We consider two nonlinear PDEs to demonstrate performance SINNs for solving
boundary observation problems.

###  5.1 A nonlinear heat equation

Let $\Omega=[0,1]\times[0,1]\subset\mathbb{R}^{2}$ and suppose that $u(x,y)$
satisfies the PDE

$\begin{split}\nabla\cdot(e^{u}\nabla u)&=0,\qquad\text{in}\;\Omega,\\\
u&=g,\qquad\text{on}\;\partial\Omega.\end{split}$ (23)

This can be viewed as a the steady solution of a nonlinear diffusion equation
for which the diffusivity, $e^{u(x,y)}$, depends on the local solution
$u(x,y)$. An alternative view is that (23) is equivalent to the nonlinear PDE
$\Delta u=-|\nabla u|^{2}$. The motivation for studying this example is that
the linearised PDE is simply the Laplace equation $\Delta u=0$, meaning that
this example will facilitate a careful comparison of the SINN methodology to
more traditional approaches which directly employ a system’s linearised
dynamics with an encoder/decoder architecture.

Data for training and testing is obtained by first creating $10^{3}$ boundary
functions $g_{i}\in L^{2}(\partial\Omega,\mathbb{R})$, with $N_{T}=900$ of
these used for training and the remaining $100$ boundary functions used for
testing. The boundary functions $g_{i}$ are created as random sums of
sinusoids. To describe this process, for any boundary point
$\boldsymbol{z}=\boldsymbol{z}(x,y)\in\partial\Omega$, let
$\alpha(\boldsymbol{z})$ be the angle, measured anticlockwise, using a co-
ordinate system with origin at the centre of the square domain $\Omega$,
namely

$\alpha(\boldsymbol{z}(x,y))=\text{atan2}\left(x-0.5,y-0.5\right),\qquad\boldsymbol{z}\in\partial\Omega.$

By sampling coefficients $X_{i},Y_{i}\sim N(0,1)$ from standard Normal
distributions, we first let

$\tilde{g}_{i}(\boldsymbol{z})=\sum_{n=1}^{4}\frac{X_{n}\sin(n\alpha(\boldsymbol{z}))+Y_{n}\cos(n\alpha(\boldsymbol{z}))}{n},\qquad\boldsymbol{z}\in\partial\Omega.$

The cosine component creates a random phase shift, while higher-order
sinusoids are moderately attenuated to encourage the lower frequency data.
Each sampled boundary function is then normalized to define the final data
boundary function $g_{i}=\tilde{g}_{i}/(\Delta_{g})$ where $\Delta_{g}$ is
randomly sampled from a triangular distribution with pdf

$f(\Delta_{g})=\begin{cases}\Delta_{g}/8,&\text{if }0\leq\Delta_{g}\leq 4,\\\
0,&\text{otherwise}.\end{cases}$

This involved approach create an ensemble of boundary functions which have
mostly large differences between their largest and smallest value.

####  5.1.1 Numerical solution and SINN implementation

For each boundary data function $g_{i}\in L^{2}(\partial\Omega,\mathbb{R})$,
the PDE (23) is solved on a uniform grid $38\times 38$ grid with a Newton
Linearization Method with finite difference equations to obtain solution data
$u_{i}\in L^{2}(\Omega,\mathbb{R})$. Given the computational domain
discretisation, we view $\Omega$ as the union of $38\times 38$ square elements
$\mathcal{T}=[1/38]\times[1/38]$, with any solution $u(x,y)$ to (23) assumed
to have a single value in each element.

To define interior encoders, we let $E$ be a square, centred at the
$(0,0)\in\mathbb{R}^{2}$, and formed of the union of $N_{E}=(2m_{e}+1)^{2}$
elements $\mathcal{T}$ for some $m_{e}\in\mathbb{N}$. In this way, the value
of any interior latent variable at $(x,y)\in\Omega$ depends only on the
solution values in the $(2m_{e}+1)^{2}$ elements symmetrically surrounding
$(x,y)\in\Omega$.

Boundary latent variables $\ell(\boldsymbol{z})$ are defined at the centre of
each exterior element, using the construction described in Section 4.2, with
$E_{\partial}$ chosen to be a line segment, centred at $0\in\mathbb{R}$, and
formed of $N_{E_{\partial}}=2m_{e}+1$ line segments whose lengths are equal to
the side length of the tile $\mathcal{T}$. Consequently, boundary latent
variables $\ell{(\boldsymbol{z})}$ depend on the boundary values
$g(\boldsymbol{z})$ on the $2m_{e}+1$ tile boundaries symmetrically
surrounding $\boldsymbol{z}\in\partial\Omega$. Note that is has been assumed
for simplicity that $N_{E}=N_{E_{\partial}}^{2}$.

Decoders are created using the construction in §4.3 by letting $D$ be a
square, centred at $(0,0)\in\mathbb{R}^{2}$, form of the union of
$N_{D}=(2m_{d}+1)^{2}$ tiles $\mathcal{T}$. Consequently, decoders seek to use
latent values at a point $(x,y)\in\Omega$ to predict the solution $u(x,y)$ on
$N_{D}$ tiles symmetrically surrounding $(x,y)$. Since the decoder output
shape is square the domain can be fully tiled by the decoder outputs, we use
partition decoders as descibed in §2.3.

As previously described the training was split into batches where each batch
needs to contain samples from both the interior and domain boundary. Since the
domain is square a distinction will also be made between edge training patches
which do not include the corner points:

$(x_{c},y_{c})\in\left\\{(0,0),(0,1),(1,0),(1,1)\right\\}$

and corner training patches which include a corner point. Each batch contained
128 internal samples, 128 edge samples and 32 corner samples which were
equally divided between the 4 sides of the square.

Finally, to implement the symmetric matrix $\mathcal{A}$ as an optimisation
variable, we define matrices $P_{11},P_{12},P_{22}\in\mathbb{R}^{r\times r}$,
let $A_{ij}:=P_{ij}+P_{ij}^{\top}$ and form the block matrix
$\mathcal{A}=(A_{ij})_{ij=1}^{2}\in\mathbb{R}^{4r^{2}}$. The corresponding
elliptic PDE for a latent variable function
$\boldsymbol{\ell}:\Omega\rightarrow\mathbb{R}^{r}$ is then given by

$D_{\mathcal{A}}\boldsymbol{\ell}=A_{11}\boldsymbol{l}_{xx}+2A_{12}\boldsymbol{l}_{xy}+A_{22}\boldsymbol{l}_{yy}=0.$

We consider the performance of SINN models for two types of boundary data. In
the first, for each boundary data function $g_{i}$, we only assume that the
boundary encoder can access pure boundary data via

$\boldsymbol{b}_{i}(\boldsymbol{z})=\left(g_{i}(\boldsymbol{z}),\boldsymbol{n}(\boldsymbol{z})\right),\qquad\boldsymbol{z}\in\partial\Omega,$
(24)

where at the corners of the square domain, the boundary vector is defined
diagonally in an outward pointing manner (e.g.
$(-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$ for the north-west corner).
Following this, we will also consider the case in which for each boundary data
function $g_{i}$, the boundary derivative of its associated solution $u_{i}$
is available to the encoder, by letting

$\boldsymbol{b}_{i}(\boldsymbol{z})=\left(g_{i}(\boldsymbol{z}),\frac{\partial
u_{i}}{\partial\boldsymbol{n}}(\boldsymbol{z}),\boldsymbol{n}(\boldsymbol{z})\right),\qquad\boldsymbol{z}\in\partial\Omega,$
(25)

The code available for generating results is available at:
https://github.com/jh6220/SINNs-for-boundary-observtion-problems.git

####  5.1.2 SINN performance with pure boundary data (24)

Table 2 shows, for different choices of the latent space dimension $r$ and the
encoder and decoder complexities $N_{E},N_{D}$, the mean square error between
the $N_{\text{test}}$ test boundary functions solved by the SINN model
$\mathcal{F}$ and ground truth is

$\mathcal{E}=\frac{1}{N_{\text{test}}}\sum_{i=1}^{N_{\text{test}}}\|u_{i}-\mathcal{F}(\boldsymbol{b}_{i})\|^{2}_{L^{2}(\Omega,\mathbb{R})}$

In the above equation, integrals are interpreted as sums over the $38^{2}$
square elements comprising the domain $\Omega$. The results in Table 2 use
SINNs in which each component (encoder, decoder, boundary encoder) is
parameterised using a neural network with $5$ hidden layers of $60$ nodes.

Table 2: Mean square error $\mathcal{E}$ of SINNs using pure boundary data
$\boldsymbol{b}\in(u_{|\partial\Omega},\boldsymbol{n}_{|\partial\Omega})$) for
the nonlinear heat equation (23).

$r$ | 1 | 2 | 3 | 8 | 2 | 3 | 5 | 8
---|---|---|---|---|---|---|---|---
$N_{E}$ | 1 | 3 | 3 | 3 | 3 | 3 | 3 | 3
$N_{D}$ | 1 | 1 | 1 | 1 | 3 | 3 | 3 | 3
$r/N_{D}^{2}$ | 1 | 2 | 3 | 8 | 2/9 | 1/3 | 5/9 | 8/9
$\mathcal{E}$ | $4.54\text{\times}{10}^{-3}$ | $7.38\text{\times}{10}^{-4}$ | $5.3\text{\times}{10}^{-4}$ | $6.56\text{\times}{10}^{-5}$ | $2.18\text{\times}{10}^{-2}$ | $8.14\text{\times}{10}^{-4}$ | $6.74\text{\times}{10}^{-4}$ | $5.29\text{\times}{10}^{-4}$

It is evident that increasing either the latent dimension $r$ or the encoder
complexity $N_{E}$ reduces the SINN error $\mathcal{E}$. The former allows for
a more complex latent space, while the latter effectively allows the encoder
to access higher order derivatives of the underlying data. For example, when
$N_{E}=3$ an encoder has access to nine local function values and is therefore
has the potential to access approximate second-order derivatives. Conversely,
increasing the decoder dimension $N_{D}$ increases the error $\mathcal{E}$
which occurs due to the choice of partition decoder used for this example. In
this case, a increasing $N_{D}$ corresponds to requiring the decoder to
extrapolate to a larger sets, naturally increasing $\mathcal{E}$. However, if
the number of degrees of freedom, i.e. $r/N_{D}^{2}$, of a trained SINN are
considered it can be seen from the penultimate row of Table 2 that a higher
value of $N_{D}$ can possibly be viewed as computationally advantageous.

To discuss the influence of latent variable dimension $r$, we consider the
trained internal elliptic models $D_{\mathcal{A}}$ in two cases. In the
simplest case $(r=1,N_{E}=1,N_{D}=1)$, the internal elliptic model is

$1.608\,\ell_{xx}+0.001\,\ell_{xy}+1.613\,\ell_{yy}=0,$

which is very close to the linearised PDE $\Delta u=0$. On the other hand,
when more modelling degrees of freedom are available for the case
$(r=3,N_{E}=3,N_{D}=1)$, the trained internal elliptic system

$\left(\begin{smallmatrix}1.46&0.49&-1.14\\\ 0.49&1.57&-0.72\\\
-1.14&-0.72&0.36\end{smallmatrix}\right)\ell_{xx}+\left(\begin{smallmatrix}0.01&-0.03&-0.02\\\
-0.03&0.02&-0.02\\\
-0.02&-0.02&0.06\end{smallmatrix}\right)\ell_{xy}+\left(\begin{smallmatrix}0.92&1.66&-1.30\\\
1.66&0.74&-1.58\\\ -1.30&-1.58&1.69\end{smallmatrix}\right)\ell_{yy}=0$

of the SINN is non-trivial, and the resulting error $\mathcal{E}$ is an order
of magnitude lower than that of the simplest model.

(a) $\epsilon u$

(b)
$\left(\mathcal{E}_{\mathcal{A}}\circ\epsilon^{\partial}\right)\boldsymbol{b}$

(c) $\epsilon
u-\left(\mathcal{E}_{\mathcal{A}}\circ\epsilon^{\partial}\right)\boldsymbol{b}$

(d) $u$

(e) $\mathcal{F}(\boldsymbol{b})$

(f) $\mathcal{F}(\boldsymbol{b})-u$

Figure 6: Indicative example of SINN performance with $(r=1,N_{E}=1,N_{D}=1)$.
(a) shows the encoded latent variable $\epsilon u$; (b) the solved latent
variable
$\left(\mathcal{E}_{\mathcal{A}}\circ\epsilon^{\partial}\right)\boldsymbol{b}$;
(c) latent variable error; (d) original data $u$; (e) SINN solution; (f) SINN
error. The SINN error for this data pair is comparable to the mean value
reported in Table 2

We next discuss the structure of the trained encoder and boundary encoders.
For the simplest case $(r=1,N_{E}=1,N_{D}=1)$, Figure 6 visualises, for a
selected test data pair $(u,\boldsymbol{b})$, the SINN solution
$\mathcal{F}(\boldsymbol{b})$ computed using the boundary data in (e); the
encoded latent variable $\epsilon u$ in (a); and the reconstructed latent
variable $(\mathcal{E}_{\mathcal{A}}\circ\epsilon^{\partial})\boldsymbol{b}$
in (b). The error field of both the SINN and the internal latent variables are
also shown in Figure 6 (c),(f) and these show small, but non-trivial,
discrepancies.

Since the trained PDE is approximately equivalent the linearised PDE, and the
latent variable space has scalar-values $(r=1)$, it is not surprising that the
encoded latent variable $\epsilon u$ and the original data $u$ are
superficially similar. Model accuracy in this case is achieved purely from the
nonlinearity of the boundary encoder and decoder. Indeed, Figure 7 (a–b) shows
two slices of the SINN solution $\mathcal{F}(\boldsymbol{b})$, ground-truth
data $u$, and the solution obtained by extending the boundary data using just
the linearised PDE $\Delta u=0$. It is clear that the nonlinear SINN solution
is substantially more accurate than the linearised model. To understand the
precise way in which the nonlinear structure achieves this increase in
accuracy, Figure 7 (c) shows both the true boundary data $\boldsymbol{b}$ and
the encoded data $\epsilon^{\partial}\boldsymbol{b}$. The encoder
$\epsilon^{\partial}$ behaves assymetricaly in the sense that it attenuates
positive boundary values and amplifies negative ones. The reason for this
behaviour is that the local diffusion coefficient $e^{u(x,y)}$ of the PDE (23)
increases exponentially with the $u(x,y)$. Thus, positive boundary values
imply diffusion on shorter length scales compared to negative boundary values.
The boundary encoder’s behaviour can now be interpreted as
$\epsilon^{\partial}$ reflecting this local nonlinear structure of the
underlying nonlinear diffusion coefficient, and this improves the accuracy of
the SINN operator $\mathcal{F}$. While this behaviour of the boundary encoder
is now interpretable, there is a persistent error if a simple model with $r=1$
is used. As indicated in Table 2, this error can be avoided by increasing the
dimension $r$ of the latent space.

Figure 7: Indicative example of SINN performance for $(r=1,N_{E}=1,N_{D}=1)$.
(a) shows solutions on a slice through the domain at $x=0.5$; (b) the
respective variables on the slice $y=0.5$; (c) Original and encoded boundary
data.

Finally, we seek to understand the geometric properties of the latent
variables for cases in which $r=3$ and $r=5$. While the latent variables do
not have an strict physical meaning, one can use sensitive analysis to extract
their underlying structure. In particular, Figure 8 shows

$\frac{\partial\delta}{\partial\boldsymbol{\ell}}(\bar{\boldsymbol{\ell}})$

where $\bar{\ell}$ is the mean value of the latent variables computed across
the entire data ensemble. The gradient was computed numerically using central
finite difference method with a step size equal to a standard deviation of
each latent dimension computed from the data-set in a similar manner to the
mean.

Figure 8: Latent dimensions interpretation using sensitivity analysis for
models with pure boundary data
$\boldsymbol{b}\in(u_{|\partial\Omega},\boldsymbol{n}_{|\partial\Omega})$).
The figure shows the gradients
$(\nabla\delta)(\bar{\boldsymbol{\ell}})_{i}=\frac{\partial\delta}{\partial\boldsymbol{\ell}_{i}}(\overline{\boldsymbol{\ell}})$
interpolated onto a fine grid

The results of the sensitivity analysis for two with $(r=3,N_{E}=3,N_{D}=3)$
and $(r=5,N_{E}=3,N_{D}=3)$ are shown in Figure 8. It can be observed that the
implies latent variables are spatially coherent. In the case $r=3$, the
structures are approximately orthogonal linear surfaces, while for $r=5$ more
complex, yet still coherent, spatial structures can be observed.

####  5.1.3 SINN performance with extended boundary data (25)

We now consider the case in which extra boundary data, namely the normal
boundary derivatives, are available to the boundary encoder. Table 2 shows
SINN errors $\mathcal{E}$ for a variety of model choices. These follow a
similar trend to the case of pure boundary conditions, although the final row
of Table 5 indicates that extra boundary information gives a consistent
performance improvement, and that this improvement is more pronounced for
higher values of $N_{D}$ and of $r$.

Table 3: SINN errors $\mathcal{E}$ when extra boundary information
($\boldsymbol{b}=(u_{|\partial\Omega},\boldsymbol{n}_{|\partial\Omega},\frac{\partial
u}{\partial\boldsymbol{n}}_{|\partial\Omega})$) is used.

$r$ | 2 | 3 | 8 | 2 | 3 | 5
---|---|---|---|---|---|---
$N_{E}$ | 3 | 3 | 3 | 3 | 3 | 3
$N_{D}$ | 1 | 1 | 1 | 3 | 3 | 3
$\mathcal{E}$ | $6.63\text{\times}{10}^{-4}$ | $2.50\text{\times}{10}^{-4}$ | $1.88\text{\times}{10}^{-5}$ | $7.3\text{\times}{10}^{-3}$ | $3.41\text{\times}{10}^{-4}$ | $1.43\text{\times}{10}^{-4}$
$\mathcal{E}/\mathcal{E}_{\text{pure BCs}}$ | 1.11 | 2.11 | 3.49 | 1.78 | 2.37 | 4.68

Figure 9 shows, for an test function $u$ whose error is indicative of the mean
values presented in Table 5, slices thought the true solution and SINN
solutions. Both slices at $x=0$ and $y=0$ in Figure 9 (a–b) show very good
agreement of the SINN solution $\mathcal{F}(\boldsymbol{b})$ with the true
solution $u$. Consider first the SINN models with $N_{D}=1$, whose decoders
are not required to extrapolate. It can be seen in both error plots of Figure
9 (c–d) that the pointwise error decreases uniformly as $r$ increases. Next,
consider the SINN models with $N_{D}=3$, whose decoders must extrapolate to
two adjacent elements. In this case, while absolute error also decreases with
increasing $r$, the error plots are oscillatory as a result of the
extrapolation error involved with the chosen partition decoder.

Figure 9: An indicative example of SINN performance for the parameters shown
in Table 3. Results are plotted along slices through the domain where $x=0.5$
(left) and $y=0.5$ (right). The upper plots (a-b) show the test data $u$ and
SINN reconstructions $\mathcal{F}(\boldsymbol{b})$; the lower plots (c–d) show
errors $u_{\text{error}}=u-\mathcal{F}(\boldsymbol{b})$.

Finally, we discuss the influence of underlying neural network complexity on
the SINN error. For the case $(r=5,N_{E}=3,N_{D}=3)$, Table 4 shows the
average SINN error $\mathcal{E}$ over the testing ensemble for different
choices of Neural Network dimensions. Three different neural network
structures are considered, in terms of the number of layers and nodes per
layer, with each case applying to the encoder, boundary encoder and decoder.
It is interesting to note that $\mathcal{E}$ decreases with neural network
complexity, suggesting the model over fitting has not occurred for these
parametric values. This highlights a potential benefit of the SINN methodology
in that, due to the use of training patches, significant training information
can be obtained from each training data pair
$(\boldsymbol{u}_{i},\boldsymbol{b}_{i})$. Consequently, the SINN approach
appears robust to overfitting, even when employing only a relatively small
training data ensemble.

Table 4: SINN error $\mathcal{E}$ as a function of the neural network
parameters for the case $(r=5,N_{E}=3,N_{D}=3)$.

$N_{\text{layers}}$ | 4 | 5 | 5
---|---|---|---
$N_{\text{nodes}}$ | 40 | 60 | 200
$N_{\text{parameters}}$ | $5.53\text{\times}{10}^{3}$ | $1.55\text{\times}{10}^{4}$ | $1.64\text{\times}{10}^{5}$
$\mathcal{E}$ | $1.03\text{\times}{10}^{-3}$ | $1.43\text{\times}{10}^{-4}$ | $2.34\text{\times}{10}^{-5}$

###  5.2 Steady laminar fluid flow

Let $\Omega=[0,1]\times[0,1]\subset\mathbb{R}^{2}$ and suppose that
$\boldsymbol{u}:\Omega\rightarrow\mathbb{R}^{2}$ and
$p:\Omega\rightarrow\mathbb{R}$ satisfy the steady, incompressible, Navier-
Stokes equations

$\displaystyle\begin{split}\boldsymbol{u}\cdot\nabla\boldsymbol{u}+\nabla
p&=\nu\Delta\boldsymbol{u},\qquad\text{in}\;\Omega\\\
\nabla\cdot\boldsymbol{u}&=0,\phantom{\Delta u}\qquad\;\text{in}\;\Omega,\\\
\boldsymbol{u}&=\boldsymbol{g},\phantom{\Delta
u}\qquad\;\text{on}\;\partial\Omega.\end{split}$ (26)

Here, $\boldsymbol{u}=(u_{x},u_{y})$ represents the velocity compents of a
fluid contained in the square domain $\Omega$, $p$ is the pressure of the
fluid, and $\nu>0$ is the kinematic viscosity of the fluid. In this example,
$\Omega$ should be thought of a control volume in a larger fluid flow. The
corresponding velocity boundary conditions must, by the divergence theorem and
incompressibility, then satisfy

$\int_{\partial\Omega}\boldsymbol{g}\cdot\boldsymbol{n}\,dS=\int_{\Omega}\nabla\cdot\boldsymbol{u}\,dV=0.$
(27)

The aim of this example is to study whether SINNs can reconstruct the fluid
velocity in the domain $\Omega$, using only boundary velocity data. This
presents a more complex and challenging example than the nonlinear heat
equation in the previous section, in view of the three-dimensional state space
comprising of two velocity components and the pressure. Furthermore, to
increase the challenge presented by this example, we will not seek to exploit
any pressure information, meaning that its influence must be automatically
discovered during SINN training. We note that, although we do not consider
such an example here, the SINN methodology could analogously be applied to a
fluid flow example with no-slip (i.e. Dirichlet) boundary conditions in which
the interior fluid velocity must be recovered using only boundary pressure
data.

A data ensemble is created using a similar method to that in §5.1. Here,
$2\times 10^{3}$ random, sinusoidal, boundary velocity functions
$\boldsymbol{g}_{i}=((u_{x})_{i},(u_{y})_{i})$ are generated, each of which
also satisfies the constraint (27). For each boundary data function, the PDE
(26) was solved using a SIMPLE algorithm on a staggered grid, with the domain
$\Omega$ discredited into $30\times 30$ rectangular cells, with the pressure
$p$ computed at each grid centre, and the velocity components $u_{x},u_{y}$
computed at the mid-point of the cell’s sides. The chosen grid was non-uniform
with greater resolution near the boundaries to aid numerical convergence. The
resulting “ground-truth” solutions were then re-sampled onto a $38\times 38$
uniform grid with the same properties as in §4.2.

The generating functions and training loop was implemented in an equivalent
manner to the example in §5.1. The only difference for this example is that
encoders and boundary encoders take the two velocity components $u_{x},u_{y}$
as inputs. We note again, that pressure information is not available during
training. Here, for brevity, we focus on case in which extra boundary
information
$\boldsymbol{b}=(\boldsymbol{u}_{|\partial\Omega}=\boldsymbol{g},\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{n}}_{|\partial\Omega},\boldsymbol{n})$
is available. Each trained SINN model this section has the same architecture
for generating functions $\epsilon$, $\epsilon^{\partial}$ and uses a
partition decoder $\delta$. The models with $N_{D}=1$ have $5$ hidden layers
of $60$ nodes, and models $N_{D}=3$ have $5$ hidden layers of $200$ nodes to
compensate for the higher dimensional decoder output.

Table 5 shows the testing errors $\mathcal{E}$ of SINN models created with a
selected parameter values. The errors exhibit the same trends as observed for
the nonlinear heat equation in §5.1, with modelling error decreasing with
increasing latent space dimension $r$ and encoder dimension $N_{E}$, and with
errors increasing as the extrapolation dimension $N_{D}$ of the partition
decoder is increased. The final row of Table 5 shows the increase in error if
only standard boundary conditions
$(\boldsymbol{u}_{|_{\partial\Omega}},\boldsymbol{n})$ are available in model
training although, for brevity, we do not discuss these results in detail
here.

Table 5: Mean square testing error $\mathcal{E}$ of SINNs using boundary
information
$\boldsymbol{b}=(\boldsymbol{u}_{|\partial\Omega},\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{n}}_{|\partial\Omega},\boldsymbol{n})$
to solve the PDE (26).

$r$ | 4 | 6 | 6 | 8 | 10
---|---|---|---|---|---
$N_{E}$ | 3 | 3 | 3 | 3 | 3
$N_{D}$ | 1 | 1 | 3 | 3 | 3
$\mathcal{E}$ | $2.32\text{\times}{10}^{-3}$ | $1.69\text{\times}{10}^{-3}$ | $2.61\text{\times}{10}^{-3}$ | $8.32\text{\times}{10}^{-4}$ | $6.75\text{\times}{10}^{-4}$
$\mathcal{E}/\mathcal{E}_{\text{pure BCs}}$ | 2 | 2.15 | 1.79 | 3.61 | 2.65

An indicative visualisation of the internal structure of two of the trained
SINNs is given in Figure 10 for the cases $(r=6,N_{E}=3,N_{D}=1)$ and
$(r=10,N_{E}=3,N_{D}=3)$. The true test function velocity components
$u_{x},u_{y}$ are shown in the top row. Both models give a solution to the
boundary value problem $\mathcal{F}(\boldsymbol{b})$ which is a very good
approximation to the true flow, as is shown in the right-hand column of Figure
10. The main noticeable difference is that the SINN model with more latent
variables $r=10$ is able to more accurately capture the elongated vertical
structure corresponding to values $u_{y}<0$. Both models have coherent, yet
non-trivial, latent variables which are shown in the middle two columns of
Figure 10. The additional degrees of freedom enjoyed by the SINN with $r=10$
allows for more accurate reconstruction of the finer scale flow features than
the simpler model with $r=6$.

Figure 10: An indicative example of SINN performance and latent variables when
solving the boundary value problem (26). The true solution, $\boldsymbol{u}$,
is shown in the top-right corner, and its SINN approximations and error fields
for the two indicated models shown in the right-hand column. The encoded
latent variables $\epsilon\boldsymbol{u}$ and solved latent variables
$(\mathcal{E}_{\mathcal{A}}\circ\epsilon^{\partial})(\boldsymbol{b})$ are
shown in the left-hand and middle columns, respectively.

To look more closely at SINN performance for boundary observation of the
Navier-Stokes PDE (26), Figure 11 shows, for one indicative example, slices
through the SINN solutions $\mathcal{F}(\boldsymbol{b})$ at $x=0.5$ and
$y=0.5$ for all models considered in Table 5. All models exhibit a good
approximation to the main trends of the true data, with approximation error
decreasing with increasing latent space dimension $r$. For this more
challenging example, modelling errors have not yet converged for the parameter
values shown. However, as can be observed from the reconstructed flow fields
shown in Figures (10), all SINN models are able to reconstruct a very close
approximation to the dominant internal vortex structures in the flow domain.
Such performance, if replicated in experimental applications, for example,
would be of great practical use.

Figure 11: An indicative example of SINN performance for the parameters shown
in Table 5. Results are plotted along slides of the domain where $x=0.5$
(left) and $y=0.5$ (right), with the original “ground truth” solution to the
PDE (26) also shown.

Finally, we perform a latent variable sensitivity analysis for the two trained
models with $(r=8,N_{E}=3,N_{D}=3)$ and $(r=10,N_{E}=3,N_{D}=3)$. The decoder
sensitivities
$\frac{\partial\delta}{\partial\boldsymbol{\ell}_{i}}(\overline{\boldsymbol{\ell}})$
each have two components which correspond to the two velocity components in
the $x$ and $y$ directions. These sensitivities are shown in Figure 12. The
model with $r=8$ exhibits approximately planar sensitivities about the mean
latent variable value $\bar{\boldsymbol{\ell}}$, suggesting that the $r=8$
latent variables are being used by the trained model to enable planar
perturbations to the solution in four directions for each of the two velocity
components. Conversely, the model with $r=10$ clearly exhibits nonlinear, yet
spatially coherent, latent-variable sensitivities which appear to enable a
more accurate solution to the underlying boundary observation problem.

Figure 12: Latent variable sensitivities
$(\nabla\delta)(\bar{\boldsymbol{\ell}})$ at the mean ensemble latent variable
value $\bar{\boldsymbol{\ell}}$.

##  6 Discussion

The numerical examples discussed in §5 suggest that SINNs are able to provide
very good approximations to the nonlinear solution operators

$\mathcal{F}:L^{2}(\partial\Omega,\mathbb{R}^{n_{\partial}})\rightarrow
L^{2}(\Omega,\mathbb{R}^{n})$

for nonlinear boundary observation problems of the form (1). The power of the
SINN approach is that, via training only finite-dimensional neural networks,
it provides nonlinear infinite-dimensional operators $\mathcal{F}$ which can
give approximate solutions $\mathcal{F}(\boldsymbol{b})\in
L^{2}(\Omega,\mathbb{R}^{n})$ to a boundary observation problem for any given
boundary data function $\boldsymbol{b}\in
L^{2}(\partial\Omega,\mathbb{R}^{n_{\partial}})$. This represents a step-
change in utility in comparison to data-driven approaches in which model
training, and hence also the trained models, directly depends on a fixed
instance of the boundary data.

From the viewpoint of operator identification, the fact that ensemble errors
in the range of $\mathcal{O}(10^{-3})$ to $\mathcal{O}(10^{-5})$ can be
obtained by SINNs with very few latent variables ($3\leq r\leq 10$) indicates
that the approach has strong potential to be successfully applied to more
complex examples. This is supported by the evidence, discussed in §5.1, that
the semi-local structure of the SINN training algorithm endows the approach
with significant robustness against over fitting. A further advantage of our
data-driven approach is that SINNs can be obtained regardless of whether the
available boundary data renders the underlying PDE boundary observation
problem over- or under-determined. The data-driven operator $\mathcal{F}$
merely attempts to find an optimal approximation to the PDE solution, given
the available training ensemble. We also emphasise that SINN training does not
require knowledge of the underling PDE, meaning that our method can be applied
to experimental data and subsequently used to solve unseen boundary
conditions.

A natural question is to ask whether the approximation error will converge to
zero with increased complexity of the trained SINN operator (e.g. as
$r,N_{E}\rightarrow\infty$, or with the complexity of the underlying neural
networks). Since our aim is identify operators which solve nonlinear boundary
observation problems which may have no closed-form solutions, and in view of
the fact that linear elliptic systems are used as the central non-local
building blocks of SINNs, it is unlikely that such convergence will hold in
general. However, even without such a property, the numerical evidence
presented in this paper suggests that SINNs can provide very good, low-
complexity, approximations to nonlinear boundary observation problems which,
furthermore, capture key physical features of the solution.

Viewing performance from an approximation accuracy philosophy is not out of
line with the motivation for many well-established approaches to the
simulation of complex nonlinear PDEs. For example, in fluid mechanics, if one
numerically solves the Reynolds Averaged Navier Stokes (RANS) equations, there
is no expectation that the solution will agree with a fully resolved direct
numerical simulation (DNS) of the governing Navier-Stokes equations. However,
in many practical cases, a RANS solution may provide sufficiently physical
insights at a substantially reduced computational cost than DNS. A similar
philosophy applies to more accurate, yet still approximate, numerical
approaches such as Large Eddy Simulation (LES). From the perspective of
creating low-cost SINN models, it should be noted that there is technically no
limit to using of significantly larger choice of extrapolation dimension
$N_{D}$ than those used in the numerical examples considered in this paper.
Furthermore, even if a SINN is trained using a finely-resolved spatial grid,
the fact that an elliptic PDE is identified implies that the SINN operator
$\mathcal{F}=\delta\circ\mathcal{E}_{\mathcal{A}}\circ\epsilon^{\partial}$ can
be implemented using an arbitrary resolution, and potentially low-cost,
solution to the central elliptic system
$\mathcal{D}_{\mathcal{A}}\boldsymbol{\ell}=0$.

The boundary encoding and internal decoding can be computationally expensive
if large neural networks are used but unlike the latent elliptic system this
computation is applied to each section of the domain independently and is
trivial to parallelise. Given a discretisation grid of $n_{i}$ internal points
the computational complexity of decoding would scale linearly with
$\mathcal{O}(n_{i})$ and the boundary encoding would scale even more
favourably, since a typical choice of the number of boundary points $n_{b}$ is
lower (e.g., for a 2D domain it may be assumed to scale as
$n_{b}\propto\sqrt{n_{i}}$). On the other hand the latent elliptic system
requires solving a linear system with $n_{i}\times r$ variables. Computational
complexity of this step depends on specific algorithm, with direct methods
such as Cholesky decomposition scaling as $\mathcal{O}((n_{i}r)^{3})$. Since
the SINN method only provides an approximate solution this precision is not
required and so an iterative method could be used which has a smaller per
iteration complexity. This is still significantly higher than the encoding or
decoding step which, at large enough $n_{i}$, would dominate the computational
cost. In summary, SINNs scale very well for large neural networks and, as
shown in the section §5.1, this can significantly increase modeling accuracy.

Finally, we comment on the computational cost of SINN training. A potential
bottleneck is that, for each update to the elliptic system coefficients
$\mathcal{A}$, one must repeatably solve a new elliptic system of PDEs on each
training patch that is used to build up the cost function
$\Psi(\mathcal{U},\mathcal{X})$. The cost of evaluating the cost function can
be controlled by using a fixed number of training patch geometries $Q$, and by
parallelising the elliptic system solutions on each training patch. To give an
example, suppose that each training patch is as shown in Figure 5 and requires
the solution of an elliptic system at $n_{p}$ internal points
$\boldsymbol{p}_{i}$ in the training patch. Solution of this elliptic system
on the training patch domain involves solving a linear system:

$A(\mathcal{A},\boldsymbol{p}_{i},\boldsymbol{q}_{j})\boldsymbol{x}=\boldsymbol{b}(\mathcal{A},\boldsymbol{p}_{i},\boldsymbol{q}_{j},\boldsymbol{\ell}(\boldsymbol{q}_{j}))$

where $A\in\mathbb{R}^{(n_{p}r)\times(n_{p}r)}$ is a symmetric positive
definite matrix which depends linearly on the coefficients of $\mathcal{A}$,
on the boundary points $\boldsymbol{q}_{j}$, and on the interior points
$\boldsymbol{p}_{ij}$. The vector $\boldsymbol{b}\in\mathbb{R}^{n_{p}r}$
depends on $\mathcal{A}$, $\boldsymbol{q}_{j}$, $\boldsymbol{p}_{i}$ and the
encoded boundary values $\boldsymbol{\ell}(\boldsymbol{q}_{i})$. Solving the
above linear system can be performed in two steps: forming the Cholesky
decomposition $A=LL^{\top}$, then using $L$ to solve the linear system via
$\boldsymbol{x}=(LL^{\top})^{-1}\boldsymbol{b}$. If a common training patch
geometry is used, the first step only needs to be computed once per training
iterate, with the matrix $L$ stored in memory. This computation can be
performed in parallel across all training patches required to compute $\Psi$.
In a similar manner, any required evaluations of the encoder and decoder can
also be parallelised. These steps imply that very efficient training of SINNs
is possible.

##  7 Conclusions

We have presented a data-driven method for solving boundary observation
problems which identifies a solution operator which can approximate the PDE
solution for arbitrary boundary data. The constructed models, referred to here
as Structure Informed Neural Networks (SINNs), embed an elliptic system into a
classical encoder/decoder Neural-Network architecture for reduced-order
modelling. The use of elliptic systems, which are well-posed with respect to
the global passage of problem data, enables very efficient model training to
be performed on small patches of the underlying domain. Numerical evidence
suggests that this endows the proposed SINN methodology with significant
robustness to over-fitting.

The methodology presented in this paper can be used to solve boundary
observation problems which are both time-independent and have boundary data
which is known on the entire boundary. Future research will investigate the
possibility of extending the SINN methodology to handle cases in which only
partial boundary data is available for training or testing, the potential for
SINN operators to be embedded in time-dependent algorithms for boundary
observation, and the application of the developed methodology to more complex
domain geometries.

##  8 Appendix

We present the proofs of the regularity results stated in the paper.

###  8.1 Proof of Lemma 1

Regularity of $\epsilon\boldsymbol{u}$: Given
$\boldsymbol{x},\boldsymbol{y}\in\Omega_{E}$, note that

$|(\epsilon\boldsymbol{u})(\boldsymbol{x})-(\epsilon\boldsymbol{u})(\boldsymbol{y})|=|e(\boldsymbol{u}_{\boldsymbol{x}})-e(\boldsymbol{u}_{\boldsymbol{y}})|.$
(28)

Now, if $\boldsymbol{x}\rightarrow\boldsymbol{y}$ in $\Omega_{E}$, then by a
standard approximation argument,
$\|\boldsymbol{u}_{\boldsymbol{x}}-\boldsymbol{u}_{\boldsymbol{y}}\|_{L^{2}(E)}\rightarrow
0$. It then follows from (28) and the assumed continuity of the generating
function $e$ that
$(\epsilon\boldsymbol{u})(\boldsymbol{x})\rightarrow(\epsilon\boldsymbol{u})(\boldsymbol{y})$,
meaning that $(\epsilon\boldsymbol{u}):\Omega_{E}\rightarrow\mathbb{R}^{r}$ is
continuous.

To prove uniform boundedness of $\epsilon\boldsymbol{u}$, note that for for
any $\boldsymbol{u}\in L^{2}(\Omega)$,

$\sup_{x\in\Omega_{E}}\|\boldsymbol{u}_{\boldsymbol{x}}\|_{L^{2}(E)}^{2}=\sup_{x\in\Omega_{E}}\int_{E}|u(\boldsymbol{x}+\boldsymbol{y})|^{2}d\boldsymbol{y}\leq\|\boldsymbol{u}\|^{2}_{L^{2}(\Omega)}$

Since $e:L^{2}(E)\rightarrow\mathbb{R}^{r}$ is compact, it maps bounded
subsets of $L^{2}(E)$ to bounded subsets of $\mathbb{R}^{r}$. Hence,

$\sup_{x\in\Omega_{E}}|(\epsilon\boldsymbol{u})(\boldsymbol{x})|=\sup_{x\in\Omega_{E}}|e(\boldsymbol{u}_{\boldsymbol{x}})|_{2}<\infty.$

Consequently, $\epsilon\boldsymbol{u}\in C(\Omega_{E},\mathbb{R}^{r})$. ∎

###  8.2 Proof of Lemma 2

Let $\boldsymbol{\ell}\in C(\Omega,\mathbb{R}^{r})$ and let $\epsilon>0$. Let
$\boldsymbol{x},\boldsymbol{z}\in\Omega$ and define sets
$D_{\boldsymbol{x}\boldsymbol{y}}:=(D_{\boldsymbol{x}}\cap
D_{\boldsymbol{z}}\cap\Omega)$ and

$D_{\boldsymbol{x}\setminus\boldsymbol{z}}=(D_{\boldsymbol{x}}\cap\Omega)\setminus
D_{\boldsymbol{x}\boldsymbol{y}},\quad
D_{\boldsymbol{z}\setminus\boldsymbol{x}}=(D_{\boldsymbol{z}}\cap\Omega)\setminus
D_{\boldsymbol{x}\boldsymbol{y}}$

and set volumes by

$c_{\boldsymbol{x}}=|D_{\boldsymbol{x}}\cap\Omega|,\quad
c_{\boldsymbol{z}}=|D_{\boldsymbol{z}}\cap\Omega|.$

For convenience, we also let $f(\cdot):=(\delta\boldsymbol{\ell})(\cdot)$ and
$g_{\boldsymbol{y}}(\cdot):=d(\ell(\boldsymbol{y}))(\cdot)$. Then,

$\displaystyle|f(\boldsymbol{x})-f(\boldsymbol{z})|$
$\displaystyle=\left|\frac{1}{c_{\boldsymbol{x}}}\int_{D_{\boldsymbol{x}}}g_{\boldsymbol{y}}(\boldsymbol{x}-\boldsymbol{y})d\boldsymbol{y}-\frac{1}{c_{\boldsymbol{z}}}\int_{D_{\boldsymbol{z}}}g_{\boldsymbol{y}}(\boldsymbol{z}-\boldsymbol{y})d\boldsymbol{y}\right|$
$\displaystyle\leq\underbrace{\frac{1}{c_{\boldsymbol{x}}}\int_{D_{\boldsymbol{x}\setminus\boldsymbol{z}}}|g_{\boldsymbol{y}}(\boldsymbol{x}-\boldsymbol{y})|d\boldsymbol{y}+\frac{1}{c_{\boldsymbol{z}}}\int_{D_{\boldsymbol{z}\setminus\boldsymbol{x}}}|g_{\boldsymbol{y}}(\boldsymbol{z}-\boldsymbol{y})|d\boldsymbol{y}}_{:=I_{1}}$
$\displaystyle\quad+\underbrace{\frac{1}{c_{\boldsymbol{x}}}\int_{D_{\boldsymbol{x}\boldsymbol{z}}}\left|g_{\boldsymbol{y}}(\boldsymbol{x}-\boldsymbol{y})-g_{\boldsymbol{y}}(\boldsymbol{z}-\boldsymbol{y})\right|d\boldsymbol{y}}_{:=I_{2}}$
$\displaystyle\quad+\underbrace{\left|\frac{1}{c_{\boldsymbol{x}}}-\frac{1}{c_{\boldsymbol{z}}}\right|\int_{D_{\boldsymbol{x}{\boldsymbol{z}}}}|g_{\boldsymbol{y}}(\boldsymbol{z}-\boldsymbol{y})|d\boldsymbol{y}}_{:=I_{3}}$

Now, since $\boldsymbol{\ell}\in C(\Omega,\mathbb{R}^{r})$, it follows that
$\boldsymbol{\ell}(\Omega)\subset\mathbb{R}^{r}$ is bounded. Then, using
compactness of the encoder generating function $d$, it follows that
$\\{d(\boldsymbol{\ell})(\boldsymbol{y})\\}_{\boldsymbol{y}\in\Omega}=\\{g_{\boldsymbol{y}}\\}_{\boldsymbol{y}\in\Omega}$
is a bounded subset of $C(\Omega,\mathbb{R}^{n})$. Hence, there exists $K>0$
such that

$\sup_{\boldsymbol{y}\in\Omega}\|g_{\boldsymbol{y}}\|_{C(\Omega,\mathbb{R}^{n})}\leq
K<\infty.$ (29)

Then, since
$D_{\boldsymbol{x}\setminus\boldsymbol{z}},D_{\boldsymbol{z}\setminus\boldsymbol{x}}\rightarrow
0$ and $c_{\boldsymbol{x}}-c_{\boldsymbol{z}}\rightarrow 0$ as
$\boldsymbol{x}\rightarrow\boldsymbol{z}$, it follows that there exists
$\delta_{1}>0$ such that

$I_{1}+I_{3}\leq
K\left(|D_{\boldsymbol{x}\setminus\boldsymbol{z}}|+|D_{\boldsymbol{z}\setminus\boldsymbol{x}}|+|D_{\boldsymbol{x}\boldsymbol{z}}|\left|\frac{1}{c_{\boldsymbol{x}}}-\frac{1}{c_{\boldsymbol{z}}}\right|\right)<\frac{\epsilon}{2}.$

whenever $|\boldsymbol{x}-\boldsymbol{y}|<\delta_{1}$.

Finally, since $d$ is continuous and
$\boldsymbol{\ell}(\bar{\Omega})\subset\mathbb{R}^{d}$ is compact, it follows
that $\mathcal{F}:=\\{g_{\boldsymbol{y}}\\}_{\boldsymbol{y}\in\bar{\Omega}}$
is a compact subset of $C(\Omega,\mathbb{R}^{n})$. Consequently, the set of
functions $\mathcal{F}$ is equicontinuous and, hence, there exists
$\delta_{2}>0$ such that
$|\boldsymbol{g}_{y}(\boldsymbol{x}-\boldsymbol{y})-g_{\boldsymbol{y}}(\boldsymbol{z}-\boldsymbol{y})|<\epsilon/2$,
for any $\boldsymbol{y}\in\Omega$, whenever
$|\boldsymbol{x}-\boldsymbol{z}|<\delta_{2}$. Hence,

$|f(\boldsymbol{x})-f(\boldsymbol{z})|\leq I_{1}+I_{2}+I_{3}\leq\epsilon$

whenever $|\boldsymbol{x}-\boldsymbol{z}|<\min\\{\delta_{1},\delta_{2}\\}$,
meaning that $f=\delta\boldsymbol{\ell}$ is continuous. That
$\delta\boldsymbol{\ell}\in C(\Omega,\mathbb{R}^{n})$ then follows from the
upper bound (29). ∎

## References

  * Bevanda et al. [2021] P. Bevanda, S. Sosnowski, and S. Hirche.  Koopman operator dynamical models: Learning, analysis and control.  _Annual Reviews in Control_ , 52:197–212, 2021.
  * Bryan and Caudill [1996] K. Bryan and L. F. Caudill.  An inverse problem in thermal imaging.  _SIAM Journal on Applied Mathematics_ , 56(3):715–735, 1996.
  * Giaquinta and Martinazzi [2012] M. Giaquinta and L. Martinazzi.  _An Introduction to the Regularity Theory for Elliptic Systems, Harmonic Maps and Minimal Graphs_.  Publications of the Scuola Normale Superiore (PSNS). Springer, 2012.
  * Illingworth et al. [2018] S. J. Illingworth, J. P. Monty, and I. Marusic.  Estimating large-scale structures in wall turbulence using linear models.  _Journal of Fluid Mechanics_ , 842:146–162, 2018.
  * Kovachki et al. [2023] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar.  Neural operator: Learning maps between function spaces with applications to pdes.  _Journal of Machine Learning Research_ , 24(89):1–97, 2023.  URL http://jmlr.org/papers/v24/21-1524.html.
  * Li et al. [2017] Q. Li, F. Dietrich, E. M. Bollt, and 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(10):103111, 2017.
  * Li et al. [2023] X.-A. Li, Z.-Q. J. Xu, and L. Zhang.  Subspace decomposition based dnn algorithm for elliptic type multi-scale pdes.  _Journal of Computational Physics_ , 488:112242, 2023.
  * Li et al. [2021] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar.  Fourier neural operator for parametric partial differential equations.  _arXiv_ , 2010.08895, 2021.
  * Lu et al. [2021] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  _Nat. Mach. Intell._ , 3:218–229, 2021.
  * Lusch et al. [2018] B. Lusch, J. Kutz, and S. Brunton.  Deep learning for universal linear embeddings of nonlinear dynamics.  _Nat. Commun._ , 9:4950, 2018.
  * Morrey Jr. and Nirenberg [1957] C. B. Morrey Jr. and L. Nirenberg.  On the analyticity of the solutions of linear elliptic systems of partial differential equations.  _Communications on Pure and Applied Mathematics_ , 10(2):271–290, 1957.
  * Pilva and Zareei [2022] P. Pilva and A. Zareei.  Learning time-dependent pde solver using message passing graph neural networks.  _arXiv_ , 2204.07651, 2022.
  * Raissi et al. [2019] M. Raissi, P. Perdikaris, and G. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ , 378:686–707, 2019\.
  * Schmid [2022] P. J. Schmid.  Dynamic mode decomposition and its variants.  _Annual Review of Fluid Mechanics_ , 54(1):225–254, 2022.
  * Snieder [1988] R. Snieder.  Large-scale waveform inversions of surface waves for lateral heterogeneity: 1. theory and numerical examples.  _Journal of Geophysical Research: Solid Earth_ , 93(B10):12055–12065, 1988.  doi: https://doi.org/10.1029/JB093iB10p12055.
  * Song et al. [2022] Y. Song, L. Shen, L. Xing, and S. Ermon.  Solving inverse problems in medical imaging with score-based generative models.  In _ICLR-22. 10th International Conference on Learning Representations_ , 2022.  doi: arXiv:2111.08005v2.
  * Wynn et al. [2013] A. Wynn, D. S. Pearson, B. Ganapathisubramani, and P. J. Goulart.  Optimal mode decomposition for unsteady flows.  _Journal of Fluid Mechanics_ , 733:473–503, 2013.