Structure-Informed Neural Networks for Boundary Observation Problems
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.
Jakub Horsky, Andrew Wynn
# Structure-Informed Neural Networks for Boundary Observation Problems

###### 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 non-linear 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 SS1.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.

### 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 \(\mathbf{x}\in\Omega\), we want to recover the value of an \(n\)-dimensional physical variable \(\mathbf{u}(\mathbf{x})\in\mathbb{R}^{n}\). To do this, we can only make use of boundary data \(\mathbf{b}(\mathbf{z})\in\mathbb{R}^{n_{\partial}}\) which can be measured at each point \(\mathbf{z}\in\partial\Omega\). It will be assumed that both interior and boundary data are square-integrable functions
in the sense that \(\mathbf{u}\in X\), where

\[X=L^{2}(\Omega,\mathbb{R}^{n})=\left\{f:\Omega\to\mathbb{R}^{n}:\int_{\Omega}\|f( \mathbf{x})\|^{2}d\mathbf{x}<\infty\right\},\]

and that \(\mathbf{b}\in Y\), where

\[Y=L^{2}(\partial\Omega,\mathbb{R}^{n_{\partial}})=\left\{f:\partial\Omega\to \mathbb{R}^{n_{\partial}}:\int_{\partial\Omega}\|f(\mathbf{z})\|^{2}d\mathbf{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}(\mathbf{u},\mathbf{\lambda})&=0, \qquad\text{in }\Omega\\ \mathcal{B}(\mathbf{u})&=\mathbf{b},\qquad\text{on }\partial\Omega, \end{split} \tag{1}\]

where \(\mathcal{L}\) is a differential operator, \(\mathbf{\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 \(\mathbf{b}\), there is a unique solution \(\mathbf{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 \(\mathbf{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\to X\\ \mathbf{b}&\mapsto\mathbf{u}\end{split} \tag{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 \(\mathbf{u}\), given any possible boundary observation \(\mathbf{b}\in X\). As will be discussed in SS1.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}:=\{(\mathbf{u}_{i},\mathbf{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.

### SINN operators

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

\[\mathcal{F}=\delta\circ\mathcal{E}_{\mathcal{A}}\circ\epsilon^{\partial}, \tag{3}\]

the structure of which is shown schematically in Figure 1. The first operator is called a _boundary encoder_\(\epsilon^{\partial}:Y\to 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 \(\mathbf{b}\in Y\), the value of \((\epsilon^{\partial}\mathbf{b})(\mathbf{z})\) at any \(\mathbf{z}\in\partial\Omega\) can only depend on the values of \(\mathbf{b}\) in a small neighbourhood \(N_{\mathbf{z}}\subset\partial\Omega\) of \(\mathbf{z}\). Practically, this will be achieved by training a neural network1\(\mathcal{N}^{\partial}:\{\mathbf{b}(\mathbf{\xi}):\mathbf{\xi}\in N_{\mathbf{z}}\}\mapsto( \epsilon^{\partial}\mathbf{b})(\mathbf{z})\) As will be described in detail in SS2.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\to Y_{L}\). This enables a wide class of nonlinear operators to be considered without significantly increasing the number of optimisation parameters.

Footnote 1: A formal mathematical definition of neural networks used in this paper is given in §4.5.

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}}\mathbf{\ell}=\sum_{i,j=1}^{d}A_{ij}\frac{\partial^{2}\mathbf{\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,\ldots,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,\qquad\quad\text{in }\Omega\\ \boldsymbol{\ell}&=\epsilon^{\partial}(\boldsymbol{b }),\qquad\text{on }\partial\Omega.\end{split} \tag{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}\to 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\\ \boldsymbol{\ell}&\longmapsto\boldsymbol{u}\end{split} \tag{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 SS2.3.

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

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

The semi-local architecture of the encoder and decoder mappings is chosen specifically so as to restrict 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 SS3, 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 SS2 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 SS3 and its numerical implemention discussed in SS4. Implementation of our approach on a pair of challenging test-cases is given in SS5. Before this, we first comment briefly on the relation between the proposed SINN architecture and other, related, data-driven approaches to PDE analysis.

### 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 \(\mathbf{u}\) as a mapping \(\mathbb{R}^{d}\ni\mathbf{x}\mapsto\mathbf{u}(\mathbf{x})\in\mathbb{R}^{n}\) and to therefore seek to construct a neural network \(\mathcal{N}_{P}:\mathbb{R}^{d}\to\mathbb{R}^{n}\) which approximates the solution. The crucial step is to add so-called physics-inspired constraints, namely \(\mathcal{L}(\mathcal{N}_{P}(\mathbf{x}),\mathbf{\lambda})_{|_{\Omega}}=0\) and \(\left[\mathcal{B}(\mathcal{N}_{P}(\mathbf{x}))-\mathbf{b}\right]_{|_{\partial\Omega}}=0\), to force the constructed solution to satisfy the underlying PDE.

In contrast to the SINN operators \(\mathcal{F}:Y\to X\) which act between functions spaces, PINNs are finite-dimensional mappings that directly attempt to replicate the solution mapping \(\mathbf{x}\mapsto\mathbf{u}(\mathbf{x})\). They require knowledge of the underlying PDE they seek to solve (i.e., of \(\mathcal{L},\mathbf{\lambda}\) and \(\mathcal{B}\)) and, when applied to boundary observation problems, must be trained using knowledge of the specific boundary data \(\mathbf{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 ba
sis 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 SS3, 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 SS3, 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.

### Interior Encoders

For the purposes of model training only, we will construct interior encoders \(\epsilon\) which, given any distribution of physical variable \(\mathbf{u}\in X\), transforms these into a distribution of latent variables \(\mathbf{\ell}=\epsilon\mathbf{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, compact2, generally nonlinear mapping

Footnote 2: A compact mapping is one which maps bounded subsets to relatively compact subsets.

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

This should be thought of as a mapping

\[e:\left\{\begin{array}{c}\text{Local patch of}\\ \text{ interior data}\end{array}\right\}\longmapsto\left\{\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 \(\mathbf{x}\in\Omega\), define a local neighbourhood

\[E_{\mathbf{x}}:=\mathbf{x}+E=\{\mathbf{x}+\mathbf{y}:\mathbf{y}\in E\},\]

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

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

\[\mathbf{u}_{\mathbf{x}}(\mathbf{y}):=\mathbf{u}(\mathbf{x}+\mathbf{y}),\qquad\mathbf{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\mathbf{u}\right)(\mathbf{x}):=e\left(\mathbf{u}_{\mathbf{x}}\right),\qquad\mathbf{ x}\in\Omega_{E}, \tag{7}\]

This definition should be thought of as mapping the physical data \(\mathbf{u}\), viewed as a _function_ in \(X=L^{2}(\Omega,\mathbb{R}^{n})\), to a new function \(\epsilon\mathbf{u}:\Omega_{E}\rightarrow\mathbb{R}^{r}\). This allows the latent variables \(\mathbf{\ell}(\mathbf{x})=(\epsilon\mathbf{u})(\mathbf{x})\) corresponding to \(\mathbf{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})\to 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})\to\mathbb{R}^{r}\) be an interior encoder generating function and let \(\epsilon\) be defined by (7). Then \(\epsilon:L^{2}(\Omega,\mathbb{R}^{n})\to C(\Omega_{E},\mathbb{R}^{r})\)._

Proof.: See Appendix 8.1. 

### 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})\to\mathbb{R}^{r}. \tag{8}\]

This should be understood as a mapping

\[e^{\partial}:\left\{\begin{array}{c}\text{Section of}\\ \text{boundary data}\end{array}\right\}\times\left\{\begin{array}{c}\text{ Section of}\\ \text{boundary geometry}\end{array}\right\}\longmapsto\left\{\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 \(\mathbf{n}(\mathbf{z})\in\mathbb{R}^{d}\) and a tangent plane \(T_{\mathbf{z}}\subset\mathbb{R}^{d-1}\) exists for every \(\mathbf{z}\in\partial\Omega\). Each tangent plane \(T_{\mathbf{z}}\) is defined in terms of a local-coordinate system with origin at \(\mathbf{z}\) and whose basis vectors \((\mathbf{e}_{i}^{\mathbf{z}})_{i=1}^{d-1}\) are orthogonal to \(\mathbf{n}(\mathbf{z})\). We assume further that there exists a ball \(B_{R}(\mathbf{z})\subset\mathbb{R}^{d}\) of radius \(R\) such that the local projection \(P_{\mathbf{z}}:\partial\Omega\cap B_{R}(\mathbf{z})\to T_{\mathbf{z}}\) from the boundary to tangent plane \(T_{\mathbf{z}}\) is one-to-one and, in addition, that there exists \(\tau>0\), independent of \(\mathbf{z}\in\partial\Omega\), such that

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

and we also assume that

\[E_{\partial}\subseteq(0,\tau)^{d-1}. \tag{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_{\mathbf{z}}^{-1}:\{t_{i}\mathbf{e}_{i}^{\mathbf{z}}:\mathbf{t}\in E_{\partial}\}\to\partial\Omega\) exists. A schematic illustration of this construction is shown in Figure 3.

This technical construction allows us, for each \(\mathbf{z}\in\partial\Omega\), to define a function \(\mathbf{b_{z}}:E_{\partial}\to\mathbb{R}^{n_{\partial}}\), which depends on the boundary data local to \(\mathbf{z}\), by

\[\mathbf{b_{z}}(\mathbf{t}):=\mathbf{b}\left(P_{\mathbf{z}}^{-1}(t_{i}\mathbf{e}_{i}^{\mathbf{z}}) \right),\qquad\mathbf{t}=(t_{i})_{i=1}^{d-1}\in E_{\partial}, \tag{11}\]

Similarly, we can also define a function \(\mathbf{n_{z}}:E_{\partial}\to\mathbb{R}^{d}\) which describes the boundary geometry local to \(\mathbf{z}\) by

\[\mathbf{n_{z}}(\mathbf{t}):=\mathbf{n}\left(P_{\mathbf{z}}^{-1}(t_{i}\mathbf{e}_{i}^{\mathbf{z}}) \right),\qquad\mathbf{t}=(t_{i})_{i=1}^{d-1}\in E_{\partial}. \tag{12}\]

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

\[\left(\epsilon^{\partial}\mathbf{b}\right)(\mathbf{z}):=e^{\partial}(\mathbf{b_{z}},\mathbf{ n_{z}}),\qquad\mathbf{z}\in\partial\Omega. \tag{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}})\to C( \partial\Omega,\mathbb{R}^{r}).\]

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

### 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 \(\mathbf{u}=\delta\mathbf{\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}\to C(D,\mathbb{R}^{n}). \tag{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,\ldots,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}\to 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 property 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})\).

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. \tag{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}\to C(D,\mathbb{R}^{n})\) be a decoder generating function and let \(\delta\) be defined by (15). Then \(\delta:C(\Omega,\mathbb{R}^{r})\to C(\Omega,\mathbb{R}^{n})\)._

Proof.: See Appendix 8.2.
### 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}} \tag{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}\to X_{L}\).

## 3 Training generating functions

It is worth summarising the constructions developed in SS2. 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(\mathbf{u}_{j}(x),\mathbf{b}_{j}(\mathbf{z})\right)_{j=1}^{N_{T}},\qquad \mathbf{x}\in\Omega,\mathbf{z}\in\partial\Omega, \tag{17}\]

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

### 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 \((\mathbf{u},\mathbf{b})\in\mathcal{U}\) is a snapshot selected from the training data ensemble.

**Stage 1:** Fix a set \(\{\mathbf{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 hull3 of a set of points \(\{\mathbf{q}_{ij}\}_{j=1}^{N}\in\Omega\cup\partial\Omega\) local to \(\mathbf{p}_{i}\), an example of which is shown in Figure 5 (a). Specifically, for each \(i\), we assume that there exist points satisfying

Footnote 3: The convex hull of a set of points is the smallest convex subset containing all such points.

* \(\mathbf{q}_{ij}\in\Omega_{E}\cup\partial\Omega\), for each \(j=1,\ldots,N\);
* \(\mathbf{p}_{i}\in\mathrm{int}(Q_{i})\) where \(Q_{i}=\mathrm{conv}\{\mathbf{q}_{ij}:j=1,\ldots,N\}\);
* \(\mathbf{q}_{ij}\in\partial Q_{i}\), for each \(j=1,\ldots,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}{cc}e(\boldsymbol{u}_{ \boldsymbol{q}_{ij}}),&\quad\text{if }\boldsymbol{q}_{ij}\in\Omega;\\ e^{\partial}(\boldsymbol{b}_{\boldsymbol{q}_{ij}},\boldsymbol{\eta}_{ \boldsymbol{q}_{ij}}),&\quad\text{if }\boldsymbol{q}_{ij}\in\partial\Omega.\end{array} \right.,\qquad j=1,\ldots,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,\ldots,N\). This process is indicted in Figure 5 (b).

**Stage 3:** Given a generating matrix \(\mathcal{A}\in\mathbb{R}^{dr\times dr}\), define the associated linear

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\to\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}\).
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} \tag{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}\to 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], \tag{19}\]

where \(\alpha>0\) is a weighting parameter.

## 4 Numerical Implementation

In SS2 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.
### Interior generating functions

The aim is to create an ensemble of possible mappings \(e:L^{2}(E,\mathbb{R}^{n})\to\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}. \tag{20}\]

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

### 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})\to\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}. \tag{21}\]

Then, for any functions \(\mathbf{f}\in L^{2}(E_{\partial},\mathbb{R}^{n_{\partial}})\) and \(\mathbf{\eta}\in L^{2}(E_{\partial},\mathbb{R}^{d})\) we can form vectors \(\mathcal{D}^{\partial}\mathbf{f}\in\mathbb{R}^{n_{\partial}\times N_{E_{\partial}}}\) and \(\mathcal{D}^{\partial}\mathbf{\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}}}\to \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})\to\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.

\begin{table}
\begin{tabular}{|c|c|c|c c c|} \hline  & Interior Encoder & \multicolumn{3}{c|}{Boundary Encoder} \\ \hline \hline Semi-local & \(\begin{matrix}\epsilon:L^{2}(\Omega)&\to&C(\Omega)\\ \mathbf{u}&\overset{(\ref{eq:L^2})}{\mapsto}&e(\mathbf{u}_{\mathbf{x}})\end{matrix}\) & \(\begin{matrix}\epsilon^{\partial}:L^{2}(\partial\Omega)\times L^{2}(\partial \Omega)&\to&C(\partial\Omega)\\ (\mathbf{b},\mathbf{n})&\overset{(\ref{eq:L^1})}{\mapsto}&e^{\partial}(\mathbf{b}_{\mathbf{x }},\mathbf{n}_{\mathbf{x}})\end{matrix}\) \\ \hline GF & \(\begin{matrix}e:L^{2}(E)\to\mathbb{R}^{r}\\ e=\tilde{e}\circ\mathcal{D}\end{matrix}\) & \(\begin{matrix}e^{\partial}:L^{2}(E_{\partial})\times L^{2}(E_{\partial})\to \mathbb{R}^{r}\\ e^{\partial}=\tilde{e}^{\partial}\circ\mathcal{D}^{\partial}\end{matrix}\) \\ \hline Finite Dim. & \(\tilde{e}:\mathbb{R}^{n\times N_{E}}\to\mathbb{R}^{r}\) & \(\tilde{e}^{\partial}:\mathbb{R}^{(n_{\partial}+d)\times N_{E_{\partial}}}\to \mathbb{R}^{r}\) \\ \hline \end{tabular}
\end{table}
Table 1: Structure of the Interior and Boundary Encoders.
### Decoder Generating Functions

We aim to create functions \(d:\mathbb{R}^{r}\to 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}}\to 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}\to\mathbb{R}^{N_{D}}\) induces a decoder generating function via \(d:=\mathcal{I}\circ\tilde{d}\).

### 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=\operatorname{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\).

### 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}}\to\mathbb{R}^{r},\tilde{e}^{\partial}: \mathbb{R}^{(n_{\partial}+d)\times N_{E_{\partial}}}\to\mathbb{R}^{r},d: \mathbb{R}^{r}\to\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}}\to\mathbb{R}^{m_{O}}\) formed by repeated composition of a prescribed nonlinear _activation function_\(\phi:\mathbb{R}\to\mathbb{R}\) and a series of affine maps \(A_{i}:\boldsymbol{x}\mapsto W_{i}\boldsymbol{x}+\boldsymbol{b}_{i}\), for \(i=1,\ldots,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)=\operatorname{ReLU}(x)=\left\{\begin{array}{cc}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}}&\Psi( \mathcal{U},\mathcal{X})-\rho\log\left(\det\mathcal{A}\right)\\ &\Theta=(\Theta_{e},\Theta_{e^{\partial}},\Theta_{d}),\\ &\mathcal{A}\in\mathbb{R}^{dr\times dr},\\ & 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}}\\ &\rho>0,\alpha>0.\end{split} \tag{22}\]

The input and output dimensions of each neural network in (22) are prescribed by the encoder and decoder structure described in SS4.1, SS4.2, and SS4.3. The number of layers and neurons in each neural network is problem dependent and will be specified in SS5 for the particular numerical examples considered. The \(-\log\left(\det\mathcal{A}\right)\) term provides a _barrier_ in the sense that \(-\log\left(\det\mathcal{A}\right)\to\infty\) as \(\det A\to 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 SS4.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}/(\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.

### 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} \tag{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 \(\mathbf{z}=\mathbf{z}(x,y)\in\partial\Omega\), let \(\alpha(\mathbf{z})\) be the angle, measured anticlockwise, using a co-ordinate system with origin at the centre of the square domain \(\Omega\), namely

\[\alpha(\mathbf{z}(x,y))=\text{atan}2\left(x-0.5,y-0.5\right),\qquad\mathbf{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}(\mathbf{z})=\sum_{n=1}^{4}\frac{X_{n}\sin(n\alpha(\mathbf{z}))+Y_{n}\cos (n\alpha(\mathbf{z}))}{n},\qquad\mathbf{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(\mathbf{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(\mathbf{z})\) depend on the boundary values \(g(\mathbf{z})\) on the \(2m_{e}+1\) tile boundaries symmetrically surrounding \(\mathbf{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 SS4.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 SS2.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\{(0,0),(0,1),(1,0),(1,1)\}\]

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 \(\mathbf{\ell}:\Omega\to\mathbb{R}^{r}\) is then given by

\[D_{\mathcal{A}}\mathbf{\ell}=A_{11}\mathbf{l}_{xx}+2A_{12}\mathbf{l}_{xy}+A_{22}\mathbf{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

\[\mathbf{b}_{i}(\mathbf{z})=\left(g_{i}(\mathbf{z}),\mathbf{n}(\mathbf{z})\right),\qquad\mathbf{z}\in \partial\Omega, \tag{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

\[\mathbf{b}_{i}(\mathbf{z})=\left(g_{i}(\mathbf{z}),\frac{\partial u_{i}}{\partial\mathbf{n}}( \mathbf{z}),\mathbf{n}(\mathbf{z})\right),\qquad\mathbf{z}\in\partial\Omega, \tag{25}\]

The code available for generating results is available at: [https://github.com/jh6220/SINNs-for-boundary-observvtion-problems.git](https://github.com/jh6220/SINNs-for-boundary-observvtion-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}(\mathbf{b}_{i})\|_{L^{2}(\Omega,\mathbb{R})}^{2}\]

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.

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.0 1&-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\]

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \(r\) & 1 & 2 & 3 & 8 & 2 & 3 & 5 & 8 \\ \(N_{E}\) & 1 & 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 \\ \hline \(\mathcal{E}\) & \(4.54\times 10^{-3}\) & \(7.38\times 10^{-4}\) & \(5.3\times 10^{-4}\) & \(6.56\times 10^{-5}\) & \(2.18\times 10^{-2}\) & \(8.14\times 10^{-4}\) & \(6.74\times 10^{-4}\) & \(5.29\times 10^{-4}\) \\ \hline \end{tabular}
\end{table}
Table 2: Mean square error \(\mathcal{E}\) of SINNs using pure boundary data \(\mathbf{b}\in(u_{|\partial\Omega},\mathbf{n}_{|\partial\Omega}))\) for the nonlinear heat equation (23).
of the SINN is non-trivial, and the resulting error \(\mathcal{E}\) is an order of magnitude lower than that of the simplest model.

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,\mathbf{b})\), the SINN solution \(\mathcal{F}(\mathbf{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})\mathbf{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}(\mathbf{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

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)\mathbf{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
and the encoded data \(\epsilon^{\partial}\mathbf{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.

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\mathbf{\ell}}(\bar{\mathbf{\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.

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

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.
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\).

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}(\mathbf{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.

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 apply
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}(\mathbf{b})\); the lower plots (c–d) show errors \(u_{\text{error}}=u-\mathcal{F}(\mathbf{b})\).
ing 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 \((\mathbf{u}_{i},\mathbf{b}_{i})\). Consequently, the SINN approach appears robust to overfitting, even when employing only a relatively small training data ensemble.

### Steady laminar fluid flow

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

\[\mathbf{u}\cdot\nabla\mathbf{u}+\nabla p =\nu\Delta\mathbf{u}, \text{in}\ \Omega \tag{26}\] \[\nabla\cdot\mathbf{u} =0, \text{in}\ \Omega,\] \[\mathbf{u} =\mathbf{g}, \text{on}\ \partial\Omega.\]

Here, \(\mathbf{u}=(u_{x},u_{y})\) represents the velocity components 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}\mathbf{g}\cdot\mathbf{n}\,dS=\int_{\Omega}\nabla\cdot\mathbf{u}\,dV =0. \tag{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

\begin{table}
\begin{tabular}{c c c c} \hline \hline \(N_{\text{layers}}\) & 4 & 5 & 5 \\ \(N_{\text{nodes}}\) & 40 & 60 & 200 \\ \(N_{\text{parameters}}\) & \(5.53\times 10^{3}\) & \(1.55\times 10^{4}\) & \(1.64\times 10^{5}\) \\ \hline \(\mathcal{E}\) & \(1.03\times 10^{-3}\) & \(1.43\times 10^{-4}\) & \(2.34\times 10^{-5}\) \\ \hline \hline \end{tabular}
\end{table}
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)\).
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 SS5.1. Here, \(2\times 10^{3}\) random, sinusoidal, boundary velocity functions \(\mathbf{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 SS4.2.

The generating functions and training loop was implemented in an equivalent manner to the example in SS5.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 \(\mathbf{b}=(\mathbf{u}_{|\partial\Omega}=\mathbf{g},\frac{\partial\mathbf{u}}{\partial\mathbf{n}} _{|\partial\Omega},\mathbf{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 SS5.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 \((\mathbf{u}_{|_{\partial\Omega}},\mathbf{n})\) are available in model training although, for brevity, we do not discuss these results in detail here.

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}(\mathbf{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\).

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}(\mathbf{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

Figure 10: An indicative example of SINN performance and latent variables when solving the boundary value problem (26). The true solution, \(\mathbf{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\mathbf{u}\) and solved latent variables \((\mathcal{E}_{\mathcal{A}}\circ\epsilon^{\partial})(\mathbf{b})\) are shown in the left-hand and middle columns, respectively.
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.

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}}(\boldsymbol{\overline{ \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 \(\boldsymbol{\overline{\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 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.
Figure 12: Latent variable sensitivities \((\nabla\delta)(\bar{\mathbf{\ell}})\) at the mean ensemble latent variable value \(\bar{\mathbf{\ell}}\).
Discussion

The numerical examples discussed in SS5 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}})\to 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 SS5.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}\to\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 SS5.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.

### Proof of Lemma 1

_Regularity of \(\epsilon\mathbf{u}\):_ Given \(\mathbf{x},\mathbf{y}\in\Omega_{E}\), note that

\[|(\epsilon\mathbf{u})(\mathbf{x})-(\epsilon\mathbf{u})(\mathbf{y})|=|e(\mathbf{u_{x}})-e(\mathbf{u_{y }})|. \tag{28}\]

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

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

\[\sup_{x\in\Omega_{E}}\|\mathbf{u_{x}}\|_{L^{2}(E)}^{2}=\sup_{x\in\Omega_{E}}\int_ {E}|u(\mathbf{x}+\mathbf{y})|^{2}d\mathbf{y}\leq\|\mathbf{u}\|_{L^{2}(\Omega)}^{2}\]
Since \(e:L^{2}(E)\to\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\mathbf{u})(\mathbf{x})|=\sup_{x\in\Omega_{E}}|e(\mathbf{u}_{ \mathbf{x}})|_{2}<\infty.\]

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

### Proof of Lemma 2

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

\[D_{\mathbf{x}\setminus\mathbf{z}}=(D_{\mathbf{x}}\cap\Omega)\setminus D_{\mathbf{x}\mathbf{y}}, \quad D_{\mathbf{z}\setminus\mathbf{x}}=(D_{\mathbf{z}}\cap\Omega)\setminus D_{\mathbf{x}\mathbf{y}}\]

and set volumes by

\[c_{\mathbf{x}}=|D_{\mathbf{x}}\cap\Omega|,\quad c_{\mathbf{z}}=|D_{\mathbf{z}}\cap\Omega|.\]

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

\[|f(\mathbf{x})-f(\mathbf{z})| =\left|\frac{1}{c_{\mathbf{x}}}\int_{D_{\mathbf{x}}}g_{\mathbf{y}}(\mathbf{x}- \mathbf{y})d\mathbf{y}-\frac{1}{c_{\mathbf{z}}}\int_{D_{\mathbf{z}}}g_{\mathbf{y}}(\mathbf{z}-\mathbf{y}) d\mathbf{y}\right|\] \[\leq\underbrace{\frac{1}{c_{\mathbf{x}}}\int_{D_{\mathbf{x}\setminus\bm {z}}}|g_{\mathbf{y}}(\mathbf{x}-\mathbf{y})|d\mathbf{y}+\frac{1}{c_{\mathbf{z}}}\int_{D_{\mathbf{z} \setminus\mathbf{x}}}|g_{\mathbf{y}}(\mathbf{z}-\mathbf{y})|d\mathbf{y}}_{:=I_{1}}\] \[\quad+\underbrace{\frac{1}{c_{\mathbf{x}}}\int_{D_{\mathbf{x}\mathbf{z}}}|g_ {\mathbf{y}}(\mathbf{x}-\mathbf{y})-g_{\mathbf{y}}(\mathbf{z}-\mathbf{y})|\,d\mathbf{y}}_{:=I_{2}}\] \[\quad+\underbrace{\left|\frac{1}{c_{\mathbf{x}}}-\frac{1}{c_{\mathbf{z}} }\right|\int_{D_{\mathbf{x}\mathbf{z}}}|g_{\mathbf{y}}(\mathbf{z}-\mathbf{y})|d\mathbf{y}}_{:=I_{3}}\]

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

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

Then, since \(D_{\mathbf{x}\setminus\mathbf{z}},D_{\mathbf{z}\setminus\mathbf{x}}\to 0\) and \(c_{\mathbf{x}}-c_{\mathbf{z}}\to 0\) as \(\mathbf{x}\to\mathbf{z}\), it follows that there exists \(\delta_{1}>0\) such that

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

whenever \(|\mathbf{x}-\mathbf{y}|<\delta_{1}\).
Finally, since \(d\) is continuous and \(\boldsymbol{\ell}(\tilde{\Omega})\subset\mathbb{R}^{d}\) is compact, it follows that \(\mathcal{F}:=\{g_{\boldsymbol{y}}\}_{\boldsymbol{y}\in\tilde{\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). 
