# Enforcing exact physics in scientific machine learning: a data-driven
exterior calculus on graphs

Nathaniel Trask<EMAIL_ADDRESS>Andy Huang  Xiaozhe Hu  Center for
Computing Research, Sandia National Laboratories  Radiation and Electrical
Science, Sandia National Laboratories  Department of Mathematics, Tufts
University

###### Abstract

As traditional machine learning tools are increasingly applied to science and
engineering applications, physics-informed methods have emerged as effective
tools for endowing inferences with properties essential for physical
realizability. While promising, these methods generally enforce physics weakly
via penalization. To enforce physics strongly, we turn to the exterior
calculus framework underpinning combinatorial Hodge theory and physics-
compatible discretization of partial differential equations (PDEs).
Historically, these two fields have remained largely distinct, as graphs are
strictly topological objects lacking the metric information fundamental to PDE
discretization. We present an approach where this missing metric information
may be learned from data, using graphs as coarse-grained mesh surrogates that
inherit desirable conservation and exact sequence structure from the
combinatorial Hodge theory. The resulting data-driven exterior calculus (DDEC)
may be used to extract structure-preserving surrogate models with mathematical
guarantees of well-posedness. The approach admits a PDE-constrained
optimization training strategy which guarantees machine-learned models enforce
physics to machine precision, even for poorly trained models or small data
regimes. We provide analysis of the method for a class of models designed to
reproduce nonlinear perturbations of elliptic problems and provide examples of
learning $H(div)/H(curl)$ systems representative of subsurface flows and
electromagnetics.

##  1 Introduction

Network models of physical systems are ubiquitous throughout the sciences and
engineering. The electronic circuit models typically learned in first year
undergraduate physics may similarly be used to describe fluid, mechanical or
heat transfer systems in corresponding hydraulic circuit, mass-spring-damper,
or thermal circuit analogies, respectively [66, 67, 29, 41, 13, 60, 21, 51].
Such network models represent discrete representations of conservation laws
corresponding to classical control volume analysis and their physical
properties are intricately tied to the underlying graph topology [59, 7],
allowing modeling of a number of thermodynamic principles [19]. In
contemporary machine learning (ML), graph neural networks (GNNs) attach a more
”black-box” message passing model to a network, inferring “object-
relation”/causal relationships from unstructured data to obtain data
representations and model pair-wise interactions [32, 33, 34, 10, 20].
Generally however, applications of GNNs have not focused upon preserving
physical structures in their network topology, although some recent works have
considered how the interplay between microstructure and mechanics can be
encoded via graphs [63].

In this work we present a framework to extract efficient data-driven network
models which exactly preserve desirable mathematical structures of the
underlying physics. The framework we introduce is general, but we focus on a
particular application in which we assume access to high-fidelity data and aim
to learn a low dimensional network as an efficient structure-preserving
surrogate model. The resulting model bears similarities to reduced-order
models (ROMs), however the control volume analysis pursued in the current
approach may provide advantages in preserving structural properties that have
proven challenging for the variational derivation underpinning ROMs [53, 52].

Solutions of PDEs generally rely upon computational meshes partitioning space
into disjoint curvilinear cells. The shared topological structure offered by
both networks and meshes allow a unified analysis in terms of the exterior
calculus. For PDEs, the discrete exterior calculus (DEC) and related spatially
compatible discretizations encapsulate a range of so-called mimetic
discretizations providing discrete solutions mimicking physical properties of
the continuous problem [27, 4, 3, 16]. These methods generally preserve
conservation properties and spectral representations of operators, provide a
coordinate-free means of prescribing physics on manifolds, and allow handling
of the non-trivial null-spaces required in electromagnetics. In topological
data analysis, combinatorial Hodge theory has emerged as a tool for analyzing
flows on graphs along with their spectral and homological properties, e.g.,
[46, 37, 45, 8, 11, 44, 58, 62]. These techniques are supported by a graph
calculus providing generalizations of gradient, curl, and divergence operators
admitting interpretation as discrete exterior derivatives. As graphs are
purely topological however, this graph calculus lacks any metric information,
and the operators do not converge in a meaningful way to the familiar vector
calculus necessary to model physical systems.

The key observation of the current work is that one may exploit the shared
structure of graphs and meshes and use data to endow the graph calculus with
”missing” metric information, in the process learning a data-driven exterior
calculus (DDEC). Conservation laws are encoded onto the graph via trainable
div/curl operators satisfying a Stokes theorem, while ”black-box” DNNs may
parameterize fluxes. This guarantees that the mathematical structures related
to conservation and exact sequence properties are independent of the errors
incurred during training. Consequently, this allows us to prove a number of
properties of the resulting model, independent of the quantity of available
data or particular local minimizer found during training.

The incorporation of physical principles into scientific machine learning
tasks has surged in recent years, as it has become apparent that ”off-the-
shelf” ML tools often fail to provide robust predictions for science and
engineering applications [6]. Physics-informed approaches have achieved a
range of successes by introducing regularizers that penalize deviations from
physical properties [55, 54, 43]. While simple to implement and effective for
a range of problems, physics in these cases are only enforced to within
optimization error and may hold to a relatively coarse tolerance (for open
problems regarding their training, see e.g. [64, 65]). Some applications
however require constraints to hold to machine precision; e.g. in
electrodynamics it is critical that the solenoidal constraint hold to machine
precision to handle the involution condition [9, 26], and in forward modeling
contexts compatible/mimetic discretizations provide approximations
guaranteeing such properties hold by construction [68, 48, 17, 14]. The
current work provides an analogous means of designing machine learning
architectures which enforce physics by construction, therefore removing the
need for physics-informed regularizers.

###  1.1 Paper organization

We first recall necessary exterior calculus fundamentals before introducing
our data-driven exterior calculus. While the theory is abstract, we focus on
applying it to learn nonlinear perturbations of $div-grad$ and $div-curl$
model problems as canonical examples of physics requiring structure
preservation. We next provide numerical analysis, establishing conditions
under which the learned model has a unique solution. We establish that the
data-driven exterior calculus inherits the desirable properties of the graph
calculus and use them to analyze the well-posedness of a class of nonlinear
elliptic problems. A necessary implementation question is how to obtain a
graph to define the model upon. We provide a specific example considering
coarse-graining of a high-fidelity mesh, associated with either a finite
element simulation or histograms binning experimental data, and show how the
relevant commuting diagrams are preserved under coarsening. Finally, we
introduce a PDE-constrained optimizer to fit the model to data, allowing
enforcement of physics exactly via equality constraint. The numerical analysis
implies that the forward problem associated with the equality constraint is
always well-posed given mild conditions on the architecture. Finally, we
provide several numerical results demonstrating how one may learn efficient
physics-preserving surrogates from high-fidelity data.

##  2 The graph exterior calculus

We recall first the graph calculus, which serves as the foundation for DDEC.
Let $\mathbf{N}=\left\\{n_{i}\right\\}_{i=1}^{N_{N}}$ denote a set of nodes.
We embed $\mathbf{N}$ in $\mathbb{R}^{d}$ by associating with each node a
unique position $\mathbf{p}_{i}\subset\mathbf{R}^{d}$, $i\in 1,\dots,N_{N}$.
We define a $k$-clique as an ordered tuple consisting of $k$ nodes, i.e.
$t_{k}=[n_{1},...,n_{k}]$. A $k$-clique has positive orientation if
$\pi=\left\\{i_{1},...,i_{k}\right\\}$ is an even permutation of
$\left\\{1,...,k\right\\}$ and negative otherwise. Via the embedding of the
graph, we may associate with each $k$-clique the $(k-1)$-simplex defined as
the convex hull of the vertices $s_{k}=[\mathbf{p}_{1},...,\mathbf{p}_{k}]$,
for which the $k\leq d$ distinct points span a $k$-dimensional hyperplane. A
$k$-chain $c_{k}$ may then be defined as a linear combination of
(k+1)-cliques, and we denote the set of $k$-chains by $C_{k}$. One may
introduce a boundary operator $\partial_{k}:C_{k+1}\rightarrow C_{k}$ defined
via

$\partial_{k}[n_{1},...,n_{k}]=\sum_{i}^{k}(-1)^{i-1}[n_{1},...,n_{i-1},\widehat{n_{i}},n_{i+1},n_{k}],$
(1)

where $\hat{\cdot}$ denotes an omitted entry, and which satisfies the property
$\partial_{k-1}\circ\partial_{k}=0$. When the dependence upon $k$ is clear, we
will sometimes write the coboundary simply as $\partial$ \- we will adopt this
convention for similar operators throughout. With these definitions in hand we
may finally introduce the chain complex as the following exact sequence
pairing $k$-chains and boundary operators.

${C_{0}}$${C_{1}}$${\dots}$${C_{d}}$$\scriptstyle{\partial_{0}}$$\scriptstyle{\partial_{1}}$$\scriptstyle{\partial_{d-1}}$
(2)

with the standard convention that $\partial_{-1}$ maps $C_{0}$ to the empty
set.

Note that in this graph context, the specific realization of the chain complex
may be qualitatively different from the DEC setting. Traditionally in
compatible discretizations, the complex is realized by partitioning the domain
of interest into a collection of disjoint simplices to obtain $C_{d}$, and
then deriving lower dimensional mesh entities $C_{k}$, $k<d$ via the boundary
operator. In contrast, for the graph setting one may obtain overlapping
simplices, in the sense that given unique $t_{k},t_{k}^{\prime}\in C_{l}$,
$t_{k}\cap t_{k}^{\prime}\neq\emptyset$. We will later discuss details
regarding specific choice of chain complex, but for now keep the presentation
abstract.

We next associate real numbers with the graph entities constituting the chain
complex. For each set of chains $C_{k}$, we introduce the dual set of cochains
$C^{k}$ consisting of linear functionals acting on $C_{k}$. Given $\phi\in
C^{k}$, we denote the value associated with the $k$-chain
$t_{i_{1}i_{2}...i_{k}}$ via the shorthand
$\phi_{i_{1}i_{2}...i_{k}}:=\phi(t_{i_{1}i_{2}...i_{k}})$. Note that cochains
inherit the orientation of the underlying chains, e.g. $\phi_{ij}=-\phi_{ji}$
via the definition of $\pi$. Introducing the coboundary operator
$\delta_{k}:C^{k}\rightarrow C^{k+1}$, we next arrive at the following cochain
complex

${C^{0}}$${C^{1}}$${\dots}$${C^{d}}$$\scriptstyle{\delta_{0}}$$\scriptstyle{\delta_{1}}$$\scriptstyle{\delta_{d-1}}$
(3)

We will formally denote the pairing between boundary and coboundary operators
via the inner product

$\langle\phi,\partial_{k}t\rangle=\langle\delta_{k}\phi,t\rangle.$ (4)

In the traditional DEC setting, one would arrive at a definition of the
coboundary via the generalized Stokes theorem, defined by the dual pairing of
$\delta_{k}$ and $\partial_{k}$ via
$\int_{\partial\omega}w=\int_{\omega}\delta w$. In the graph setting, we
identify $\delta_{k}$ algebraically as the adjoint of the matrix representing
$\partial_{k}$. For example, this gives rise to the following combinatorial
gradient, assigning to the $2$-clique $t_{ij}$ the function

$\delta_{0}\phi_{ij}=\phi_{j}-\phi_{i}.$ (5)

Similarly we may obtain the combinatorial curl by assigning to the 3-clique
$t_{ijk}$ the function

$\delta_{1}\phi_{ijk}=\phi_{ij}+\phi_{jk}+\phi_{ki}.$ (6)

One may easily see that $\delta_{1}\circ\delta_{0}=0$. For the purposes of
this work, we will not require $\delta_{k}$ for $k>1$, however the definition
extends naturally to $1<k\leq d$, and one may show that
$\delta_{k}\circ\delta_{k-1}=0$. In this manner, the coboundary operator
inherits the exact sequence property of the boundary operator.

We next let $(\cdot,\cdot)_{k}$ denote an inner product mapping $C^{k}\times
C^{k}\rightarrow\mathbb{R}$. This inner product induces a codifferential
operator $\delta^{*}_{k}:C^{k+1}\rightarrow C^{k}$ via the pairing
$(v,\delta^{*}_{k}u)_{k}=(\delta_{k}v,u)_{k+1}$. In this manner, the careful
choice of $(\cdot,\cdot)_{k}$ will endow the codifferential with desirable
approximation properties, however we note that independent of the choice of
inner products the codifferential again inherits the exact sequence properties
of the coboundary operator so that $\delta^{*}_{k-1}\circ\delta^{*}_{k}=0$.
This follows trivially from the definition, so that for all $v\in C^{k-2}$ and
all $u\in C^{k}$

$(v,\delta^{*}_{k-1}\delta^{*}_{k}u)_{k-2}=(\delta_{k-1}v,\delta^{*}_{k}u)_{k-1}=(\delta_{k}\delta_{k-1}v,u)_{k}=0.$
(7)

Finally, we will contrast how the choice of inner product typically used in
the graph exterior calculus precludes the use of graph boundary/coboundary
operators in discretizing PDE. In the graph context, one selects as
$(\cdot,\cdot)_{k}$ the $\ell_{2}$ inner-product:

$(x,y)_{k}=\underset{{i\in\text{dim}(C^{k})}}{\sum}x_{i}y_{i}.$ (8)

And its induced norm is denoted by $\|\cdot\|_{k}$. In this case, the
codifferential $\delta^{*}_{k}$ may be identified as the transpose of the
matrix associated with the coboundary $\delta_{k}$. For the remainder of this
work, we will assume $\ell_{2}$ inner products in the definition of the
codifferential unless otherwise specified, so that $\delta_{0}^{*}$ and
$\delta_{1}^{*}$ correspond to the graph divergence and curl.

Besides the properties that $\delta_{k}\circ\delta_{k-1}=0$ and
$\delta_{k-1}^{*}\circ\delta_{k}^{*}=0$, the Hodge-Laplacian is defined as
follows

$L_{k}=\delta_{k-1}\delta_{k-1}^{*}+\delta_{k}^{*}\delta_{k}.$ (9)

From the combinatorial Hodge theorem, the combinatorial Laplacian decomposes
the space $C^{k}$ as follows

$C^{k}=\operatorname{im}(\delta_{k-1})\bigoplus\operatorname{ker}(L_{k})\bigoplus\operatorname{im}(\delta_{k}^{*}),$
(10)

which is the Hodge decomposition [28]. One direct consequence of the Hodge
decompostion (10) is the following Poincaré inequalities,

$\displaystyle\|z_{k}\|_{k}\leq c\|\delta_{k}z\|_{k+1},\quad
z_{k}\in\operatorname{im}(\delta_{k}^{*}),$ $\displaystyle\|z_{k}\|_{k}\leq
c^{*}\|\delta_{k-1}^{*}z_{k}\|_{k-1},\quad
z_{k}\in\operatorname{im}(\delta_{k-1}).$

Here $c$ and $c^{*}$ are Poincaré constants which depend upon the topological
structure of the graph. Another consequence of the Hodge decomposition (10) is
that the Hodge Laplacian $L_{k}$ is positive semidefinite and the dimension of
the null space is equal to the dimension of the homology
$\operatorname{ker}(\delta_{k})/\operatorname{im}(\delta_{k-1})$. This implies
the invertibility of the Hodge Laplacian $L_{k}$ on the orthogonal complement
of the null space.

In the next section, we will introduce a data-driven exterior calculus and
establish analogues of these properties.

##  3 The data-driven exterior calculus

We now parameterize these graph calculus operators in a manner which allows
recovery of traditional DEC schemes for PDE discretizations as a specific
instance. We consider the de Rham complex as a prototypical means of analyzing
structure-preserving differential operators in $\mathbb{R}^{d}$, and seek to
develop its discrete approximation on a graph. In general the construction
presented here may be used to approximate any exact sequence: we restrict our
presentation to the de Rham complex as a single example to support later
discretization of PDEs in $\mathbb{R}^{d}$. In three dimensions, the de Rham
complex for $\Omega\in\mathbb{R}^{3}$ is given by

${C^{\infty}(\Omega)}$${\left[C^{\infty}(\Omega)\right]^{3}}$${\left[C^{\infty}(\Omega)\right]^{3}}$${C^{\infty}(\Omega)}$$\scriptstyle{\operatorname{grad}}$$\scriptstyle{{\operatorname{div}}^{*}}$$\scriptstyle{\operatorname{curl}}$$\scriptstyle{{\operatorname{curl}}^{*}}$$\scriptstyle{\operatorname{div}}$$\scriptstyle{{\operatorname{grad}}^{*}}$
(11)

Here, coboundary operators $\operatorname{grad}$, $\operatorname{curl}$, and
$\operatorname{div}$ satisfying
$\operatorname{curl}\circ\operatorname{grad}=\operatorname{div}\circ\operatorname{curl}=0$;
codifferential operators $\operatorname{div}^{*}$, $\operatorname{curl}^{*}$,
and $\operatorname{grad}^{*}$ satisfying
$\operatorname{curl}^{*}\circ\operatorname{grad}^{*}=\operatorname{div}^{*}\circ\operatorname{curl}^{*}=0$.
To support later two-dimensional examples, we will also consider the
restriction to two dimensions $\Omega\in\mathbb{R}^{2}$.

${C^{\infty}(\Omega)}$${\left[C^{\infty}(\Omega)\right]^{2}}$${C^{\infty}(\Omega)}$$\scriptstyle{\operatorname{curl}}$$\scriptstyle{{\operatorname{curl}}^{*}}$$\scriptstyle{\operatorname{div}}$$\scriptstyle{{\operatorname{grad}}^{*}}$
(12)

For completeness, we recall the alternate complex in $\mathbb{R}^{2}$

$\leavevmode\hbox to236.59pt{\vbox
to26.92pt{\pgfpicture\makeatletter\hbox{\hskip
118.29692pt\lower-17.35913pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{}{}{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{\offinterlineskip{}{}{{{}}{{}}{{}}}{{{}}}{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-118.29692pt}{-9.56194pt}\pgfsys@invoke{
}\hbox{\vbox{\halign{\pgf@matrix@init@row\pgf@matrix@step@column{\pgf@matrix@startcell#\pgf@matrix@endcell}&#\pgf@matrix@padding&&\pgf@matrix@step@column{\pgf@matrix@startcell#\pgf@matrix@endcell}&#\pgf@matrix@padding\cr\hfil\hskip
53.04556pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-48.74002pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${C^{\infty}(\Omega)}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}}}&\hskip
53.04556pt\hfil&\hfil\hskip 46.71455pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-18.40904pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{${\left[C^{\infty}(\Omega)\right]^{2}}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}}}&\hskip
22.71458pt\hfil&\hfil\hskip 42.53677pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-14.23126pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${C^{\infty}(\Omega)}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip
18.5368pt\hfil\cr}}}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}}{{{{}}}{{}}{{}}{{}}}} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} {}{ {}{}{}}{}{ {}{}{}} {{{{{}}{
{}{}}{}{}{{}{}}}}}{}{{{{{}}{
{}{}}{}{}{{}{}}}}}{{}}{}{}{}{}{}{{{}{}}}{}{{}}{}{}{}{{{}{}}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.39998pt}\pgfsys@invoke{
}{{{}}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{1.72221pt}\pgfsys@invoke{
}}{}{}{}{}{{}}{}{}{{}}\pgfsys@moveto{-12.0058pt}{-7.06194pt}\pgfsys@lineto{11.19421pt}{-7.06194pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}}}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{11.3942pt}{-7.06194pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} {
}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-7.02109pt}{-3.34807pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{$\scriptstyle{\operatorname{grad}}$} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{ {}{}{}}{}{ {}{}{}}
{{{{{}}{ {}{}}{}{}{{}{}}}}}{}{{{{{}}{
{}{}}{}{}{{}{}}}}}{{}}{}{}{}{}{}{{{}{}}}{}{{}}{}{}{}{{{}{}}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.39998pt}\pgfsys@invoke{
}{{{}}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{-1.72221pt}\pgfsys@invoke{
}}{}{}{}{}{{}}{}{}{{}}\pgfsys@moveto{11.59418pt}{-7.06194pt}\pgfsys@lineto{-11.60583pt}{-7.06194pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}}}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{-1.0}{0.0}{0.0}{-1.0}{-11.80582pt}{-7.06194pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} {
}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-5.96971pt}{-15.20636pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{$\scriptstyle{{\operatorname{div}}^{*}}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{ {}{}{}}{}{ {}{}{}}
{{{{{}}{ {}{}}{}{}{{}{}}}}}{}{{{{{}}{
{}{}}{}{}{{}{}}}}}{{}}{}{}{}{}{}{{{}{}}}{}{{}}{}{}{}{{{}{}}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.39998pt}\pgfsys@invoke{
}{{{}}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{1.72221pt}\pgfsys@invoke{
}}{}{}{}{}{{}}{}{}{{}}\pgfsys@moveto{57.42334pt}{-7.06194pt}\pgfsys@lineto{80.62335pt}{-7.06194pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}}}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{80.82333pt}{-7.06194pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} {
}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{63.38026pt}{-4.70917pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{$\scriptstyle{\operatorname{curl}}$} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{ {}{}{}}{}{ {}{}{}}
{{{{{}}{ {}{}}{}{}{{}{}}}}}{}{{{{{}}{
{}{}}{}{}{{}{}}}}}{{}}{}{}{}{}{}{{{}{}}}{}{{}}{}{}{}{{{}{}}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.39998pt}\pgfsys@invoke{
}{{{}}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{-1.72221pt}\pgfsys@invoke{
}}{}{}{}{}{{}}{}{}{{}}\pgfsys@moveto{81.02332pt}{-7.06194pt}\pgfsys@lineto{57.8233pt}{-7.06194pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}}}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{-1.0}{0.0}{0.0}{-1.0}{57.62332pt}{-7.06194pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} {
}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{62.38026pt}{-15.20636pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{$\scriptstyle{{\operatorname{curl}}^{*}}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>(13)

In this work however, we consider only (12) as the natural complex to obtain
conservation properties.

###  3.1 Data-driven coboundaries and codifferentials

Consider the general case $(\Omega\in\mathbb{R}^{d})$, we define new
coboundary and codifferential operators by associating learnable metric
information with the graph operators as follows,

$\mathsf{d}_{k}:=\mathbf{B}_{k+1}\delta_{k}\mathbf{B}_{k}^{-1},\quad\text{and}\quad\mathsf{d}_{k}^{*}=\mathbf{D}_{k}^{-1}\delta_{k}^{*}\mathbf{D}_{k+1}.$
(14)

For example, in $\mathbb{R}^{3}$, we have

$GRAD=\mathbf{B}_{1}\delta_{0}\mathbf{B}_{0}^{-1}\hskip
28.45274ptCURL=\mathbf{B}_{2}\delta_{1}\mathbf{B}_{1}^{-1}\hskip
28.45274ptDIV=\mathbf{B}_{3}\delta_{2}\mathbf{B}_{2}^{-1}$ (15)
$DIV^{*}=\mathbf{D}_{0}^{-1}\delta^{*}_{0}\mathbf{D}_{1}\hskip
28.45274ptCURL^{*}=\mathbf{D}_{1}^{-1}\delta^{*}_{1}\mathbf{D}_{2}\hskip
28.45274ptGRAD^{*}=\mathbf{D}_{2}^{-1}\delta^{*}_{2}\mathbf{D}_{3}$ (16)

In $\mathbb{R}^{2}$, we have

$CURL=\mathbf{B}_{1}\delta_{0}\mathbf{B}_{0}^{-1},\hskip
28.45274ptDIV=\mathbf{B}_{2}\delta_{1}\mathbf{B}_{1}^{-1},$ (17)
$CURL^{*}=\mathbf{D}_{0}^{-1}\delta_{0}^{*}\mathbf{D}_{1},\hskip
28.45274ptGRAD^{*}=\mathbf{D}_{1}^{-1}\delta_{1}^{*}\mathbf{D}_{2},$ (18)

Here, we denote discrete graph operations in capital letters, and use lower
case for continuum counterparts. The $\mathbf{B}_{k}$ and $\mathbf{D}_{k}$
denote a diagonal tensor with positive entries weighting corresponding
elements of $C^{k}$. They naturally provide chain maps as follows,

${C^{0}}$${C^{1}}$${C^{2}}$${C^{3}}$${\cdots}$${C^{d}}$${C^{0}}$${C^{1}}$${C^{2}}$${C^{3}}$${\cdots}$${C^{d}}$${C^{0}}$${C^{1}}$${C^{2}}$${C^{3}}$${\cdots}$${C^{d}}$$\scriptstyle{\mathsf{d}_{0}}$$\scriptstyle{\mathsf{d}_{1}}$$\scriptstyle{\mathsf{d}_{2}}$$\scriptstyle{\mathsf{d}_{3}}$$\scriptstyle{\mathsf{d}_{d-1}}$$\scriptstyle{\delta_{0}}$$\scriptstyle{\mathbf{B}_{0}}$$\scriptstyle{\delta^{*}_{1}}$$\scriptstyle{\delta_{1}}$$\scriptstyle{\mathbf{B}_{1}}$$\scriptstyle{\delta^{*}_{2}}$$\scriptstyle{\delta_{2}}$$\scriptstyle{\mathbf{B}_{2}}$$\scriptstyle{\delta^{*}_{3}}$$\scriptstyle{\delta_{3}}$$\scriptstyle{\mathbf{B}_{3}}$$\scriptstyle{\delta^{*}_{d-1}}$$\scriptstyle{\delta_{d-1}}$$\scriptstyle{\delta^{*}_{d}}$$\scriptstyle{\mathbf{B}_{d}}$$\scriptstyle{\mathbf{D}_{0}}$$\scriptstyle{\mathbf{D}_{1}}$$\scriptstyle{\mathsf{d}_{0}^{*}}$$\scriptstyle{\mathbf{D}_{2}}$$\scriptstyle{\mathsf{d}_{1}^{*}}$$\scriptstyle{\mathbf{D}_{3}}$$\scriptstyle{\mathsf{d}_{2}^{*}}$$\scriptstyle{\mathsf{d}_{3}^{*}}$$\scriptstyle{\mathsf{d}_{d-1}^{*}}$$\scriptstyle{\mathbf{D}_{d}}$
(19)

Based on the definitions (15) and (16), it is easy to verify that the diagram
(19) is commutative, i.e.,

$\mathbf{B}_{k+1}\delta_{k}=\mathsf{d}_{k}\mathbf{B}_{k},\quad\text{and}\quad\mathbf{D}_{k}^{-1}\delta_{k}^{*}=\mathsf{d}_{k}^{*}\mathbf{D}_{k+1}^{-1}.$

In $\mathbb{R}^{3}$, we have,

$\mathbf{B}_{1}\delta_{0}=GRAD\mathbf{B}_{0},\quad\mathbf{B}_{2}\delta_{1}=CURL\mathbf{B}_{1},\quad\mathbf{B}_{3}\delta_{2}=DIV\mathbf{B}_{2},$

and

$\delta_{0}^{*}\mathbf{D}_{0}^{-1}=DIV^{*}\mathbf{D}_{1}^{-1},\quad\mathbf{D}_{1}^{-1}\delta_{1}^{*}=CURL^{*}\mathbf{D}_{2}^{-1},\quad\mathbf{D}_{2}^{-1}\delta_{2}^{*}=GRAD^{*}\mathbf{D}_{3}^{-1}.$

And in $\mathbb{R}^{2}$, similarly, we have

$\mathbf{B}_{1}\delta_{0}=CURL\mathbf{B}_{0},\quad\mathbf{B}_{2}\delta_{1}=DIV\mathbf{B}_{1},$

and

$\delta_{0}^{*}\mathbf{D}_{0}^{-1}=CURL^{*}\mathbf{D}_{1}^{-1},\quad\mathbf{D}_{1}^{-1}\delta_{1}^{*}=GRAD^{*}\mathbf{D}_{2}^{-1}.$

######  Theorem 3.1.

The discrete derivatives $\mathsf{d}_{k}$ in (14) form an exact sequence if
the simplicial complex is exact, and in particular
$\mathsf{d}_{k+1}\circ\mathsf{d}_{k}=0$.

###### Proof.

Because (19) is a commutative diagram, the chain maps take cycles to cycles
and boundaries to boundaries. Therefore, if the simplicial complex is exact,
then (15) forms an exact sequence. Moreover,
$\mathsf{d}_{k+1}\circ\mathsf{d}_{k}=0$ can be verified by the definitions
(14). ∎

######  Remark 3.1.

In $\mathbb{R}^{3}$, we have $CURL\circ GRAD=DIV\circ CURL=0$. And in
$\mathbb{R}^{2}$, we have $DIV\circ CURL=0$.

######  Theorem 3.2.

The discrete derivatives $\mathsf{d}_{k}^{*}$ in (14) form an exact sequence
if the simplicial complex is exact, and in particular
$\mathsf{d}_{k}^{*}\circ\mathsf{d}_{k+1}^{*}=0$.

###### Proof.

The conclusion follows from the fact that (19) is commutative and the
definitions (14). The proof is essentially the same as the proof of Theorem
3.1 ∎

######  Remark 3.2.

In $\mathbb{R}^{3}$, we have $DIV^{*}\circ CURL^{*}=CURL^{*}\circ GRAD^{*}=0$.
And in $\mathbb{R}^{2}$, we have $CURL^{*}\circ GRAD^{*}=0$

We will refer to this collection of operators as a data-driven exterior
calculus, with the understanding that the metric information encoded in
$\mathbf{B}_{k}$ and $\mathbf{D}_{k}$ will be learned from data. Note that in
the traditional low-order compatible/mimetic schemes, these metric tensors
contain geometric information related to the oriented measures of mesh
entities, such as cell volumes, face moments, etc. [1]. Following from the
exact sequence property, this exterior calculus structure inherits the
following other properties of the graph calculus.

Naturally, based on the Hodge Laplacians (9) and (19), we can define
generalized data-driven Hodge-Laplacians as follows,

$\Delta_{k}=\mathsf{d}_{k-1}\mathsf{d}_{k-1}^{*}+\mathsf{d}_{k}^{*}\mathsf{d}_{k}$

For example, for practical applications, we consider the following Hodge-
Laplacians in $\mathbb{R}^{2}$:

$\displaystyle\Delta_{1}$ $\displaystyle:=CURL\circ CURL^{*}+GRAD^{*}\circ
DIV$ $\displaystyle\Delta_{2}$ $\displaystyle:=DIV\circ GRAD^{*}.$

The Hodge decomposition (10) also can be generalized naturally. Here, we
choose $(\cdot,\cdot)_{k}:=(\cdot,\cdot)_{\mathbf{D}_{k}\mathbf{B}_{k}^{-1}}$
and denote its induced norm as $\|\cdot\|_{k}$.

######  Theorem 3.3 (Hodge Decomposition).

For $C^{k}$, the following decomposition holds

$C^{k}=\operatorname{im}(\mathsf{d}_{k-1}){\bigoplus}_{k}\operatorname{ker}(\Delta_{k}){\bigoplus}_{k}\operatorname{im}(\mathsf{d}_{k}^{*}),$
(20)

where ${\bigoplus}_{k}$ means the orthogonality with respect to the
$(\cdot,\cdot)_{k}$-inner product.

###### Proof.

Since (19) is a commutative diagram, following from the Hodge decomposition
(10), $C^{k}$ can be naturally decomposed into three parts,
$\operatorname{im}(\mathsf{d}_{k-1})$,
$\operatorname{im}(\mathsf{d}_{k}^{*})$, and $\operatorname{ker}(\Delta_{k})$.
Next we show this decomposition is orthogonal with respect to the
$(\cdot,\cdot)_{k}$-inner product, i.e.,
$(\cdot,\cdot)_{\mathbf{D}_{k}\mathbf{B}_{k}^{-1}}$-inner product. For
$\mathbf{u}_{k-1}\in C^{k-1}$ and $\mathbf{u}_{k+1}\in C^{k+1}$, we have

$\displaystyle\quad(\mathsf{d}_{k-1}\mathbf{u}_{k-1},\mathsf{d}_{k}^{*}\mathbf{u}_{k+1})_{k}$
$\displaystyle=(\mathbf{B}_{k}\delta_{k-1}\mathbf{B}_{k-1}^{-1}\mathbf{u}_{k-1},\mathbf{D}_{k}^{-1}\delta_{k}^{*}\mathbf{D}_{k+1}\mathbf{u}_{k+1})_{\mathbf{D}_{k}\mathbf{B}_{k}^{-1}}$
$\displaystyle=(\delta_{k-1}\mathbf{B}_{k-1}^{-1}\mathbf{u}_{k-1},\delta_{k}^{*}\mathbf{D}_{k+1}\mathbf{u}_{k+1})$
$\displaystyle=(\delta_{k}\delta_{k-1}\mathbf{B}_{k-1}^{-1}\mathbf{u}_{k-1},\mathbf{D}_{k+1}\mathbf{u}_{k+1})$
$\displaystyle=0$

For $\mathbf{h}_{k}\in\operatorname{ker}(\Delta_{k})$, we have
$\mathsf{d}_{k-1}^{*}\mathbf{h}_{k}=0$ and $\mathsf{d}_{k}\mathbf{h}_{k}=0$,
which implies $\delta_{k-1}^{*}\mathbf{D}_{k}\mathbf{h}_{k}=0$ and
$\delta_{k}\mathbf{B}_{k}^{-1}\mathbf{h}_{k}=0$. And then for
$\mathbf{u}_{k-1}\in C^{k-1}$,

$\displaystyle\quad(\mathsf{d}_{k-1}\mathbf{u}_{k-1},\mathbf{h}_{k})_{k}$
$\displaystyle=(\mathbf{B}_{k}\delta_{k-1}\mathbf{B}_{k-1}^{-1}\mathbf{u}_{k-1},\mathbf{h}_{k})_{\mathbf{D}_{k}\mathbf{B}_{k}^{-1}}$
$\displaystyle=(\mathbf{D}_{k}\delta_{k-1}\mathbf{B}_{k-1}^{-1}\mathbf{u}_{k-1},\mathbf{h}_{k})$
$\displaystyle=(\mathbf{B}_{k-1}^{-1}\mathbf{u}_{k-1},\delta_{k-1}^{*}\mathbf{D}_{k}\mathbf{h}_{k})$
$\displaystyle=0.$

On the other hand, for $\mathbf{u}_{k+1}\in C^{k+1}$, we have

$\displaystyle\quad(\mathsf{d}_{k}^{*}\mathbf{u}_{k+1},\mathbf{h}_{k})_{k}$
$\displaystyle=(\mathbf{D}_{k}^{-1}\delta_{k}^{*}\mathbf{D}_{k+1}\mathbf{u}_{k+1},\mathbf{h}_{k})_{\mathbf{D}_{k}\mathbf{B}_{k}^{-1}}$
$\displaystyle=(\mathbf{B}_{k}^{-1}\delta_{k}^{*}\mathbf{D}_{k+1}\mathbf{u}_{k+1},\mathbf{h}_{k})$
$\displaystyle=(\mathbf{D}_{k+1}\mathbf{u}_{k+1},\delta_{k}\mathbf{B}_{k}^{-1}\mathbf{h}_{k})$
$\displaystyle=0.$

Thus, the decomposition is orthogonal with respect to the
$(\cdot,\cdot)_{k}$-inner product, which completes the proof. ∎

For example, in $\mathbb{R}^{3}$, we have the following Hodge decomposition
when $k=1$,

$C^{1}=\operatorname{im}(GRAD){\bigoplus}_{k}\operatorname{ker}(\Delta_{1}){\bigoplus}_{k}\operatorname{im}(CURL^{*}),$

and when $k=2$

$C^{2}=\operatorname{im}(CURL){\bigoplus}_{k}\operatorname{ker}(\Delta_{2}){\bigoplus}_{k}\operatorname{im}(GRAD^{*}).$

In $\mathbb{R}^{2}$, we have the following Hodge decomposition when $k=1$

$C^{1}=\operatorname{im}(CURL){\bigoplus}_{k}\operatorname{ker}(\Delta_{1}){\bigoplus}_{k}\operatorname{im}(GRAD^{*}).$

Based on the Hodge decomposition, we have the following Poincaré inequality.

######  Theorem 3.4 (Poincaré inequality).

For each $k$, there exists a constant $c_{P,k}$ such that

$\|\mathbf{z}_{k}\|_{k}\leq
c_{P,k}\|\mathsf{d}_{k}\mathbf{z}_{k}\|_{k+1},\quad\mathbf{z}_{k}\in\operatorname{im}(\mathsf{d}_{k}^{*}),$

and another constant $c_{P,k}^{*}$ such that

$\|\mathbf{z}_{k}\|_{k}\leq
c^{*}_{P,k}\|\mathsf{d}^{*}_{k-1}\mathbf{z}_{k}\|_{k-1},\quad\mathbf{z}_{k}\in\operatorname{im}(\mathsf{d}_{k-1}).$

Thus, for $\mathbf{u}_{k}\in C^{k}$, we have

$\inf_{\mathbf{h}_{k}\in\operatorname{ker}(\Delta_{k})}\|\mathbf{u}_{k}-\mathbf{h}_{k}\|_{k}\leq
C\left(\|\mathsf{d}_{k}\mathbf{u}_{k}\|_{k+1}+\|\mathsf{d}^{*}_{k-1}\mathbf{u}_{k}\|_{k-1}\right),$

where constant $C>0$ only depends on $c_{P,k}$ and $c_{P,k}^{*}$.

###### Proof.

These inequalities are a direct consequence of the Hodge decomposition (20). ∎

Since we are considering matrix representation, by direct computation, we can
see that

$c_{P,k}=\lambda_{\min}(\mathsf{d}_{k}^{*}\mathsf{d}_{k})^{-1/2}.$

where $\lambda_{\min}$ denotes the smallest non-trivial eigenvalue. For
example, in $\mathbb{R}^{3}$, when $k=0$,

$c_{P,0}=\lambda_{\min}(\mathsf{d}_{0}^{*}\mathsf{d}_{0})^{-1/2}=\lambda_{\min}(\delta_{0}^{*}\mathbf{D}_{1}\mathbf{B}_{1}\delta_{0})^{-1/2}=\lambda_{\min}(DIV^{*}\circ
GRAD)^{-1/2}.$

Note that

$\min_{i}\\{(\mathbf{D}_{1}\mathbf{B}_{1})_{ii}\\}\lambda(\delta_{0}^{*}\delta_{0})\leq\lambda(\delta_{0}^{*}\mathbf{D}_{1}\mathbf{B}_{1}\delta_{0})\leq\max_{i}\\{(\mathbf{D}_{1}\mathbf{B}_{1})_{ii}\\}\lambda(\delta_{0}^{*}\delta_{0}).$

This implies

$\min_{i}\\{(\mathbf{D}_{1}\mathbf{B}_{1})_{ii}\\}\lambda_{\min}(\delta_{0}^{*}\delta_{0})\leq
c_{P,0}^{-2}\leq\max_{i}\\{(\mathbf{D}_{1}\mathbf{B}_{1})_{ii}\\}\lambda_{\min}(\delta_{0}^{*}\delta_{0}),$

which relates the Poincaré constant with
$\lambda_{\min}(\delta_{0}^{*}\delta_{0})$, also known as the Fielder value of
the graph Laplacian $L_{0}=\delta_{0}^{*}\delta_{0}$. Classical works provide
bounds upon the Fiedler eigenvalue in terms of the size and degree of a given
graph, see for example [22].

Another consequence of the Hodge decomposition (20) is the invertibility of
the Hodge Laplacian $\Delta_{k}$ once we account for its kernel.

######  Theorem 3.5 (Invertibility of Hodge Laplacian).

The $k^{th}$-order Hodge Laplacian $\Delta_{k}$ is positive-semidefinite, with
the dimension of its null-space equal to the dimension of the corresponding
homology
$H^{k}=\operatorname{ker}(\mathsf{d}_{k})/\operatorname{im}(\mathsf{d}_{k-1})$.

###### Proof.

For $\mathbf{u}_{k}\in C^{k}$, we have

$\displaystyle\quad(\Delta_{k}\mathbf{u}_{k},\mathbf{u}_{k})_{k}$
$\displaystyle=((\mathbf{B}_{k}\delta_{k-1}\mathbf{B}_{k-1}^{-1}\
\mathbf{D}_{k-1}^{-1}\delta_{k-1}^{*}\mathbf{D}_{k}+\mathbf{D}_{k}^{-1}\delta_{k}^{*}\mathbf{D}_{k+1}\
\mathbf{B}_{k+1}\delta_{k}\mathbf{B}_{k}^{-1})\mathbf{u}_{k},\mathbf{u}_{k})_{\mathbf{D}_{k}\mathbf{B}_{k}^{-1}}$
$\displaystyle=((\mathbf{D}_{k}\delta_{k-1}\mathbf{B}_{k-1}^{-1}\
\mathbf{D}_{k-1}^{-1}\delta_{k-1}^{*}\mathbf{D}_{k}+\mathbf{B}_{k}^{-1}\delta_{k}^{*}\mathbf{D}_{k+1}\
\mathbf{B}_{k+1}\delta_{k}\mathbf{B}_{k}^{-1})\mathbf{u}_{k},\mathbf{u}_{k})$
$\displaystyle=(\mathsf{d}_{k-1}^{*}\mathbf{u}_{k},\mathsf{d}_{k-1}^{*}\mathbf{u}_{k})_{k-1}+(\mathsf{d}_{k}\mathbf{u}_{k},\mathsf{d}_{k}\mathbf{u}_{k})_{k+1}$

which shows that $\Delta_{k}$ is positive-semidefinite. The second part
follows directly from the Hodge decomposition (20). ∎

Theorem 3.5 means that the following linear system of the Hodge Laplacian

$\Delta_{k}\mathbf{u}_{k}=\mathbf{f}_{k},$ (21)

is solvable as long as $\mathbf{f}_{k}\in C^{k}$ is orthogonal to
$\operatorname{ker}(\Delta_{k})$ with respect to the $(\cdot,\cdot)_{k}$-inner
product.

###  3.2 Nonlinear Perturbation of Hodge-Laplacian Problems

In many cases, it is helpful to consider the mixed form of the Hodge-Laplacian
problem (21) as it naturally provides connections to integral balance laws and
conservation principles [5, 2]. To this end, let us introduce a new variable
$\mathbf{w}_{k-1}:=\mathsf{d}_{k-1}^{*}\mathbf{u}_{k}$ and the mixed form of
(21) as follows,

$\displaystyle\mathbf{w}_{k-1}-\mathsf{d}_{k-1}^{*}\mathbf{u}_{k}$
$\displaystyle=0,$ (22)
$\displaystyle\mathsf{d}_{k-1}\mathbf{w}_{k-1}+\mathsf{d}_{k}^{*}\
\mathsf{d}_{k}\mathbf{u}_{k}$ $\displaystyle=\mathbf{f}_{k}.$ (23)

This class of problems describes several canonical second-order elliptic
operators; for example, in $\mathbb{R}^{2}$, for $k=2$ we obtain the Darcy
flow model

$\mathbf{F}+\kappa\nabla\phi=0$ | $\rightarrow$ | $\mathbf{w}_{1}-GRAD^{*}\mathbf{u}_{0}=0$
---|---|---
$\nabla\cdot\mathbf{F}=f$ |  | $DIV\mathbf{w}_{1}=\mathbf{f}_{0},$

and for $k=1$ we obtain the magnetostatics model, after applying a vector
potential for the magnetic field and applying a suitable gauge condition [15,
18].

$\nabla\times\mathbf{J}=\mathbf{f}$ | $\rightarrow$ |  $\mathbf{w}_{0}-CURL^{*}\mathbf{u}_{1}=0$,
---|---|---
$\nabla\cdot\mathbf{B}=0$ |  | $CURL\mathbf{w}_{0}+GRAD^{*}\circ DIV\mathbf{u}_{1}=\mathbf{f}_{1}.$
$\mathbf{J}=\mu\mathbf{B}$ |  |

While this model form is appropriate for learning, e.g. diffusion coefficients
corresponding to second-order elliptic problems, realistic problems require
accounting for nonlinearities. With this in mind, we introduce a nonlinear
perturbation of the fluxes while leaving the relevant conservation statements
untouched. This preserves the exterior calculus structure while incorporating
data into fluxes only, which are traditionally more empirical. Any
parameterization may be used for the nonlinearities, but we consider deep
neural networks. As a result, we obtain a nonlinear perturbation of a Hodge-
Laplacian problem in the mixed form as follows,

$\displaystyle\mathbf{w}_{k-1}=\mathsf{d}^{*}_{k-1}\mathbf{u}_{k}+\epsilon\mathcal{NN}(\mathsf{d}^{*}_{k-1}\mathbf{u}_{k};\xi),$
(24) $\displaystyle\mathsf{d}_{k-1}\mathbf{w}_{k-1}+\mathsf{d}_{k}^{*}\
\mathsf{d}_{k}\mathbf{u}_{k}=\mathbf{f}_{k}.$ (25)

The corresponding primal form is

$\Delta_{k}\mathbf{u}_{k}+\epsilon\mathsf{d}_{k-1}\circ\mathcal{NN}(\mathsf{d}^{*}_{k-1}\mathbf{u}_{k};\xi)=\mathbf{f}_{k}$
(26)

Later in this section, we will theoretically show that when $\epsilon>0$ is
sufficiently small, the nonlinear problem (24) and (25) remains well-posed.
First, let us look at some examples. In $\mathbb{R}^{2}$, when $k=2$, we have

$\displaystyle DIV\circ GRAD^{*}\mathbf{u}_{0}+\epsilon
DIV\circ\mathcal{NN}(GRAD^{*}\mathbf{u}_{0};\xi)=\mathbf{f}_{0}$
$\displaystyle=\Delta_{0}\mathbf{u}_{2}+\epsilon
DIV\circ\mathcal{NN}(GRAD^{*}\mathbf{u}_{0};\xi)=\mathbf{f}_{0}$

and when $k=1$, we have

$\displaystyle GRAD^{*}\circ DIV\mathbf{u}_{1}+CURL\circ
CURL^{*}\mathbf{u}_{1}+\epsilon
CURL\circ\mathcal{NN}(CURL^{*}\mathbf{u}_{1};\xi)=\mathbf{f}_{1},$
$\displaystyle=\Delta_{1}\mathbf{u}_{1}+\epsilon
CURL\circ\mathcal{NN}(CURL^{*}\mathbf{u}_{1};\xi)=\mathbf{f}_{1}.$

###  3.3 Well-posedness

Next we investigate the well-posedness of the perturbed Hodge-Laplacian
problem (24)-(25). We write the perturbed problem in the primal form, i.e.,

$\Delta_{k}\mathbf{u}_{k}+N(\mathbf{u}_{k})=\mathbf{f},$ (27)

where
$N(\mathbf{u}_{k}):=\epsilon\mathsf{d}_{k-1}\circ\mathcal{NN}(\mathsf{d}^{*}_{k-1}\mathbf{u}_{k};\xi)$.
Consider the space $\mathbb{V}=C^{k}\backslash\operatorname{ker}(\Delta_{k})$,
we introduce the following weak formulation of (27),

$a(\mathbf{u}_{k},\mathbf{v})+\epsilon(\mathcal{NN}(\mathsf{d}_{k-1}^{*}\mathbf{u}_{k}),\mathsf{d}_{k-1}^{*}\mathbf{v})_{k-1}=(\mathbf{f},\mathbf{v})_{k},\quad\mathbf{v}\in\mathbb{V},$
(28)

where

$a(\mathbf{u},\mathbf{v}):=(\mathsf{d}^{*}_{k-1}\mathbf{u},\mathsf{d}^{*}_{k-1}\mathbf{v})_{k-1}+(\mathsf{d}_{k}\mathbf{u},\mathsf{d}_{k}\mathbf{v})_{k+1},\quad\forall\,\mathbf{u},\,\mathbf{v}\in\mathbb{V},$

and its induced norm is
$\|\mathbf{u}\|_{a}:=\sqrt{a(\mathbf{u},\mathbf{u})}=\sqrt{\|\mathsf{d}_{k-1}^{*}\mathbf{u}\|_{k-1}+\|\mathsf{d}_{k}\mathbf{u}\|_{k+1}^{2}}$.
We assume Liptschitz continuity of the nonlinear perturbation, i.e. that there
exists a constant $L_{N}>0$, such that,

$\|\mathcal{NN}(\mathbf{v}_{k-1})-\mathcal{NN}(\mathbf{w}_{k-1})\|_{k-1}\leq
L_{N}\|\mathbf{v}_{k-1}-\mathbf{w}_{k-1}\|_{k-1},\quad\forall\mathbf{v}_{k-1},\,\mathbf{w}_{k-1}\in
C^{k-1}.$ (29)

In addition, we also assume that

$\mathcal{NN}(\mathbf{0})=\mathbf{0},$ (30)

which means that the nonlinear perturbation reduces to a linear gradient
closure in the limit as the gradient becomes small.

The main tool we use is the Leray-Schauder fixed point theorem [31] We define
$T:\mathbb{V}\mapsto\mathbb{V}$ such that for each $\mathbf{w}\in\mathbb{V}$,
$\widetilde{\mathbf{u}}:=T(\mathbf{w})\in\mathbb{V}$ is given as the solution
of the following linear problem,

$a(\widetilde{\mathbf{u}},\mathbf{v})+\epsilon(\mathcal{NN}(\mathsf{d}_{k-1}^{*}\mathbf{w}),\mathsf{d}_{k-1}^{*}\mathbf{v})_{k-1}=(\mathbf{f},\mathbf{v})_{k},\quad\forall\mathbf{v}\in\mathbb{V}.$
(31)

The map $T$ is clearly continuous and, therefore, compact in the finite
dimensional space $\mathbb{V}$. The sovability of (28) is thus equivalent to
the solvability of the equation $\mathbf{u}_{k}=T(\mathbf{u}_{k})$ in
$\mathbb{V}$, which is a fixed point problem.

If $\lambda>0$ and $\mathbf{w}$ satisfies $T(\mathbf{w})=\lambda\mathbf{w}$,
then

$\lambda
a(\mathbf{w},\mathbf{v})+\epsilon(\mathcal{NN}(\mathsf{d}_{k-1}^{*}\mathbf{w}),\mathsf{d}_{k-1}^{*}\mathbf{v})_{k-1}=(\mathbf{f},\mathbf{v})_{k},\quad\forall\mathbf{v}\in\mathbb{V}.$

By choosing $\mathbf{v}=\mathbf{w}$, we obtain

$\lambda\|\mathbf{w}\|_{a}^{2}\leq\|\mathbf{f}\|_{-a}\|\mathbf{w}\|_{a}+\epsilon\|\mathcal{NN}(\mathsf{d}_{k-1}^{*}\mathbf{w})\|_{k-1}\|\mathsf{d}_{k-1}^{*}\mathbf{w}\|_{k-1}$

By the Liptschitz continuity assumption (29) and assumption (30), we have

$\displaystyle\|\mathcal{NN}(\mathsf{d}_{k-1}^{*}\mathbf{w})\|_{k-1}=\|\mathcal{NN}(\mathsf{d}_{k-1}^{*}\mathbf{w})-\mathcal{NN}(\mathbf{0})\|_{k-1}\leq
L_{N}\|\mathsf{d}_{k-1}^{*}\mathbf{w}\|_{k-1}.$

Using the fact that
$\|\mathsf{d}_{k-1}^{*}\mathbf{w}\|_{k-1}\leq\|\mathbf{w}\|_{a}$, we have

$\lambda\|\mathbf{w}\|_{a}^{2}\leq\left(\|\mathbf{f}\|_{-a}+\epsilon
L_{N}\|\mathbf{w}\|_{a}\right)\|\mathbf{w}\|_{a},$

therefore,

$\lambda\leq\frac{\|\mathbf{f}\|_{-a}+\epsilon
L_{N}\|\mathbf{w}\|_{a}}{\|\mathbf{w}\|_{a}}.$

Thus, $\lambda<1$ holds true for any $\mathbf{w}$ being on the boundary of the
ball in $\mathbb{V}$ centered at the origin with radius

$\rho\geq\frac{\|\mathbf{f}\|_{-a}}{(1-\epsilon L_{N})}.$

Consequently, the Leray-Schauder fixed point theorem implies that the
nonlinear map $T$ has a fixed point in any ball centered at the origin with
radius $\rho\geq\frac{\|\mathbf{f}\|_{-a}}{(1-\epsilon L_{N})}$. This fixed
point is a solution of the equation (28).

######  Theorem 3.6.

Assume (29) and (30) hold. If $\epsilon L_{N}<1$, the equation (28) has an
unique solution $\mathbf{u}_{k}\in\mathbb{V}$ satisfies

$\|\mathbf{u}_{k}\|_{a}\leq\frac{\|\mathbf{f}\|_{-a}}{(1-\epsilon L_{N})}.$
(32)

###### Proof.

The existence has been discussed before the theorem. Let $\mathbf{u}_{k}$ be
the solution and $\mathbf{v}=\mathbf{u}_{k}$ in (28), we have

$\|\mathbf{u}_{k}\|_{a}^{2}\leq\left(\|\mathbf{f}\|_{-a}+\epsilon
L_{N}\|\mathbf{u}_{k}\|_{a}\right)\|\mathbf{u}_{k}\|_{a},$

which implies (32).

Next we prove uniqueness of (28) under the same assumptions. Let
$\mathbf{u}_{k}$ and $\bar{\mathbf{u}}_{k}$ be two solutions of (28). Denoting
$\mathbf{e}=\mathbf{u}_{k}-\bar{\mathbf{u}}_{k}$, we have,

$a(\mathbf{e},\mathbf{v})+\epsilon(\mathcal{NN}(\mathsf{d}_{k-1}^{*}\mathbf{u}_{k}),\mathsf{d}_{k-1}^{*}\mathbf{v})_{k-1}-\epsilon(\mathcal{NN}(\mathsf{d}_{k-1}^{*}\bar{\mathbf{u}}_{k}),\mathsf{d}_{k-1}^{*}\mathbf{v})_{k-1}=0.$

Letting $\mathbf{v}=\mathbf{e}$ and assuming Lipschitz continuity,

$\|\mathcal{NN}(\mathbf{v}_{k-1})-\mathcal{NN}(\mathbf{w}_{k-1})\|_{k-1}\leq
L_{N}\|\mathbf{v}_{k-1}-\mathbf{w}_{k-1}\|_{k-1},\quad\forall\mathbf{v}_{k-1},\,\mathbf{w}_{k-1}\in
C^{k-1},$ (33)

we arrive at,

$\displaystyle\|\mathbf{e}\|_{a}^{2}\leq\epsilon
L_{N}\|\mathsf{d}_{k-1}^{*}\mathbf{e}\|_{k-1}^{2}\leq\epsilon
L_{N}\|\mathbf{e}\|^{2}_{a}$

which implies the uniqueness of the solution since $\epsilon L_{N}<1$. This
completes the proof. ∎

We finally consider design of an architecture which satisfies this condition
to ensure the extracted model is solvable, considering multilayer perceptrons
as canonical architectures [56].

######  Remark 3.3.

Consider a $L$-layer neural net which has the following structure

$\mathcal{NN}(\mathbf{x}):=\psi_{L}\circ T_{L}\circ\cdots\circ\psi_{1}\circ
T_{1}(\mathbf{x})$

where
$T_{\ell}(\mathbf{x}_{\ell}):=M_{\ell}\mathbf{x}_{\ell}+\mathbf{b}_{\ell}$ is
an affine function and $\psi_{\ell}$ is a nonlinear activation function. If we
assume $\psi_{\ell}$ are $1$-Lipschitz nonlinear functions (e.g., ReLU, Leaky
ReLU, tanh, sigmoid) and define
$\mathcal{N}_{\ell}:=T_{\ell}\circ\psi_{\ell-1}\cdots\circ\psi_{1}\circ
T_{1}(\mathbf{x})$, then the Lipschitz constant $L_{N}$ (29) of $\mathcal{NN}$
can be estimated as follows,

$\displaystyle L_{N}$
$\displaystyle=\sup_{\mathbf{x}\in\mathbb{R}^{n}}\|\operatorname{diag}(\psi_{L}^{\prime}(\mathcal{N}_{L}))M_{L}\cdots
M_{2}\operatorname{diag}(\psi_{1}^{\prime}(\mathcal{N}_{1}))M_{1}\|_{k-1}$
$\displaystyle\leq\max_{i}\\{(\mathbf{D}_{k-1}\mathbf{B}_{k-1}^{-1})_{ii}\\}\sup_{\mathbf{x}\in\mathbb{R}^{n}}\|\operatorname{diag}(\psi_{L}^{\prime}(\mathcal{N}_{L}))M_{L}\cdots
M_{2}\operatorname{diag}(\psi_{1}^{\prime}(\mathcal{N}_{1}))M_{1}\|$
$\displaystyle\leq\max_{i}\\{(\mathbf{D}_{k-1}\mathbf{B}_{k-1}^{-1})_{ii}\\}\prod_{\ell=1}^{L}\|M_{\ell}\|$

Therefore, if we choose
$0<\epsilon<\left(\max_{i}\\{(\mathbf{D}_{k-1}\mathbf{B}_{k-1}^{-1})_{ii}\\}\prod_{\ell=1}^{L}\|M_{\ell}\|\right)^{-1}$,
then the assumption $\epsilon L_{n}<1$ would be satisfied. We also note that a
more accurate upper bound of the Lipschitz constant $L_{N}$ can be computed
numerically using the advanced algorithms developed in [57], providing a
tighter bound on $\epsilon$.

###  3.4 Construction of chain complex

The model introduced in the previous section assumes access to an underlying
graph to apply the DDEC to. Motivated by our surrogate modeling application,
we assume access to a very fine polygonal mesh and access to a high-fidelity
PDE solution defined as oriented moments of mesh entities (i.e. cell average
scalar potentials, face average fluxes, edge average circulations). We next
show how a coarsening of the underlying fine mesh preserves some structure
allowing a particularly simple implementation; essentially, given graph
grad/curl/div matrices on the fine mesh, we derive coarsening matrices that
encode proper orientations for the coarsened complex and the desired
$\mathsf{d}_{k}$ and $\mathsf{d}^{*}_{k}$ operators. We consider here graph-
cut coarsening available in packages such as METIS [39], which partitions the
domain $\Omega$ into $N_{c}$ disjoint volumetric subdomains (or $d$-cells),
and derive appropriately oriented lower degree $k$-cells. In Figure 1 we
provide a cartoon of the process. While the following presentation provides a
mathematical description of the process via a commutative diagram, it
practically will provide a simple implementation of the coarse coboundary
operators $\delta_{k}$ in terms of a few simple matrices, greatly simplifying
implementation.

######  Remark 3.4.

This is one particular construction appropriate for synthetic data which takes
advantage of available adjacency matrices of an underlying fine mesh. We
stress however that one may apply the calculus to any appropriately defined
graph. For example, in experimental contexts the bins associated with
histograms may be used instead, or the calculus may be applied to graphs with
no associated partition of space at all (see e.g. [38]). We further comment
that the graph-cut coarsening assumed here provides a quasi-uniform partition
of space that does not take advantage of the data; in another work we consider
spectral graph partitioning strategies to obtain data-driven partitions tuned
to give optimal representations of data.

$\longrightarrow$

$\longrightarrow$

Figure 1: The coarsening process takes a fine mesh discretizing $\Omega$ in
$d$-cells, partitions them into disjoint subdomains and uses the fine scale
mesh to derive coarsened $k$-cells and coboundary operators $\partial_{k}$
with consistent orientations. Values on the coarsened $k$-cells correspond to
sums of the constituent fine-grained moments with proper accounting of
orientation. Top: five $2$-cells coarsened into a pair of $2$-cells. Bottom: a
representative coarsening of a Cartesian mesh. Coarsened $0$-, $1$-, and
$2$-cells are each colored differently and are a formal linear sum of elements
of the fine mesh.

We assume the high-fidelity solution is associated with a $d$-dimensional
chain complex, with $k$-cells denoted by $F_{k}$ for $k<d$ and
$\partial_{k}^{fine}:F_{k+1}\rightarrow F_{k}$ as defined in Section 2, and
coboundary $\delta_{k}^{fine}:F^{k}\rightarrow F^{k+1}$ encoded via the
adjacency matrices generally available in mesh data structures. Here, we use
the symbol $F$ and superscript $fine$ to denote fine scale objects. Our goal
is to construct a coarse complex, denoted by the symbol $C$, that can be used
in the DDEC framework as we discussed in Section 3.1-3.3. The construction
contains two steps. First, we will coarsen the fine level set of the
$k$-chain, $F_{k}$, to obtain the coarse level set of the $k$-chain, $C_{k}$.
This can be done inductively by starting with coarsening the $d$-cells. Then
we will define the coarse boundary and coboundary operators via relating the
fine and coarse level complex properly.

#### Constructing $C_{k}$ and $C^{k}$.

To inductively define $C_{k}$ from $F_{k}$, we start with coarsening the
$d$-cells. Given a partitioning of the $N_{f}$ $d$-cells into $N_{c}$
disjoint, connected subsets, denoted by $F_{d}=\sqcup_{i=1}^{N_{c}}P_{i}$. we
define the space of coarse $d$-chain $C_{d}$ as a subspace of $F_{d}$ based on
the partition $\\{P_{i}\\}_{i=1}^{N_{c}}$ with the following natural inclusion

$\iota_{d}\in\mathbb{R}^{N_{f}\times
N_{c}},\quad(\iota_{d})_{ij}=\begin{cases}1,\qquad f_{i}\in P_{j},\\\ 0,\qquad
f_{i}\notin P_{j}.\end{cases}$

This has a dual $C^{d}$ identified as a subspace of $F^{d}$.

Now, consider the image of $C_{d}$ under $\partial_{d-1}^{fine}$, consisting
of the fine $d$-chains lying at partition interfaces. We define the element
$c_{ij}:=P_{i}\cap P_{j}$ (oriented consistently) for $i,j=1,2,\ldots N_{c}$
and their span as the the coarse $(d-1)$ chains $C_{d-1}$. Then we can proceed
inductively. With the chain and cochain spaces indexed by $d-k$, consider the
image of $C_{d-k}$ under the map $\partial_{d-k-1}^{fine}$. This consists of
the $(d-k-1)$-interfaces between the coarse $d-k$-chains. Define the element

$c_{i_{1}i_{2}\cdots i_{d-k-1}}=P_{i_{1}}\cap P_{i_{2}}\cap\cdots\cap
P_{i_{d-k-1}}$

The span of these elements is the space of coarse $(d-k-1)$-chains
$C_{d-k-1}$. As before, by identifying $C_{d-k-1}$ as a subset of $F_{d-k-1}$,
we can define the dual cochains $C^{d-k-1}$.

######  Remark 3.5.

It may be the case that $c_{i_{1}i_{2}\cdots i_{d-k-1}}$ is disjoint. As a
choice of implementation detail, each of the connected components can be taken
as a separate chain, or their disjoint union can be taken as a single chain.
What results from this choice is the treatment of parallel fluxes between
partitions as either: a sum of individual distinct fluxes between the
partitions, or the sum of those fluxes as a single effective quantity between
the partitions. We choose the former in our computational examples.

By construction, we have $C_{k}\subset F_{k}$, for $k=0,\cdots,d$. Therefore,
there is a natural inclusion $\iota_{k}:C_{k}\mapsto F_{k}$ whose entries are
$1$, $-1$, or $0$ depending on the partition and orientation. Furthermore, we
also have $C^{k}\subset F^{k}$ and the corresponding natural inclusions are
defined as $\iota^{k}:=\iota_{d-k}$, for $k=0,\cdots,d$.

#### Building boundary and coboundary operators.

In order to define the coarse boundary and coboundary operators that connects
$C_{k}$ and $C^{k}$, we first define the linear projections
$\pi_{k}:F_{k}\mapsto C_{k}$ as

$\pi_{k}=\left(\iota_{k}^{T}\iota_{k}\right)^{-1}\iota_{k}^{T},$

and $\pi^{k}:F^{k}\mapsto C^{k}$ as

$\pi^{k}=\left(\left(\iota^{k}\right)^{T}\iota^{k}\right)^{-1}\left(\iota^{k}\right)^{T}.$

Note that, $\pi^{k}=\pi_{d-k}$ since $\iota^{k}=\iota_{d-k}$.

With these in hand, we define the coarse boundary operators as follows,

$\displaystyle\partial_{k}:=$
$\displaystyle\pi_{k}\cdot\partial_{k}^{fine}\cdot\iota_{k+1},$

and the coarse coboundary operators

$\displaystyle\delta_{k}:=$
$\displaystyle\pi^{k+1}\cdot\delta_{k}^{fine}\cdot\iota^{k}.$

######  Remark 3.6.

The linear projections $\pi_{k}$ and $\pi^{k}$ are least-squares projections.
Since the natural inclusions $\iota_{k}$ and $\iota^{k}$ are constructed based
on partitions, $\iota_{k}^{T}\iota_{k}$ and $(\iota^{k})^{T}\iota^{k}$ are
diagonal matrices with diagonal entries equal to the number of $k$-chains in
the corresponding partition. Thus, in a practical implementation, we can
simply use $\pi_{k}=\iota_{k}^{T}$ and $\pi^{k}=(\iota^{k})^{T}$ and, based on
the definition (14), $(\iota_{k}^{T}\iota_{k})^{-1}$ and
$((\iota^{k})^{T}\iota^{k})^{-1}$ are implicitly absorbed in $\mathbf{B}_{k}$
and $\mathbf{D}_{k}$.

To summarize, in the case of $d=2$, all of the spaces discussed are related
through the following commutative diagram:

${0}$${F^{0}}$${F^{1}}$${F^{2}}$${0}$${0}$${C^{0}}$${C^{1}}$${C^{2}}$${0}$${0}$${F_{0}}$${F_{1}}$${F_{2}}$${0}$${0}$${C_{0}}$${C_{1}}$${C_{2}}$${0}$$\scriptstyle{\delta_{0}^{fine}}$$\scriptstyle{\pi^{0}}$$\scriptstyle{\delta_{1}^{fine}}$$\scriptstyle{\pi^{1}}$$\scriptstyle{\pi^{2}}$$\scriptstyle{\delta_{0}}$$\scriptstyle{\iota^{0}}$$\scriptstyle{\delta_{1}}$$\scriptstyle{\iota^{1}}$$\scriptstyle{\iota^{2}}$$\scriptstyle{\pi_{0}}$$\scriptstyle{\partial_{0}^{fine}}$$\scriptstyle{\pi_{1}}$$\scriptstyle{\pi_{2}}$$\scriptstyle{\partial_{1}^{fine}}$$\scriptstyle{\iota_{0}}$$\scriptstyle{\partial_{0}}$$\scriptstyle{\iota_{1}}$$\scriptstyle{\partial_{1}}$$\scriptstyle{\iota_{2}}$
(34)

The top-front row of this diagram corresponds to the middle row of the chain
complex (19). The DDEC complex results from perturbing these operators as in
(14).

##  4 PDE-constrained optimization

We finally turn toward the question of how to fit a model of the form (26) to
data. We assume access to data via the coarsening process of the previous
section, and denote by $\mathbf{u}_{data}$ a vector concatenating the
coarsened $\mathbf{w}_{k+1}$ and $\mathbf{u}_{k}$ degrees of freedom. We may
have only partial data: observations of possibly a single field
$\mathbf{w}_{k+1}$ or $\mathbf{u}_{k}$, or observations on a subset of the
chain complex. We then concisely express the boundary value problem in (26)
via the nonlinear operator $\mathcal{L}_{\xi}[\mathbf{u}]=\mathbf{f}$, where
we have lumped all model parameters into $\xi$ (i.e. $\mathbf{B}_{k}$,
$\mathbf{D}_{k}$, and neural network weights and biases), and refer to this as
the forward problem. We postpone a problem-specific discussion of how boundary
conditions will be imposed for the following section, and assume the forward
problem is prescribed such that BCs are imposed naturally. We seek a solution
to the following quadratic program with nonlinear equality constraints

$\displaystyle\underset{\xi}{\text{argmin}}||\mathbf{u}-\mathbf{u}_{data}||^{2}$
(35) $\displaystyle\text{such that }\mathcal{L}_{\xi}[\mathbf{u}]=\mathbf{f},$
(36)

where $||\cdot||$ denotes the $\ell_{2}$ norm.

To avoid handling the equality constraint, one may introduce a single scalar
penalty parameter $\lambda$ and use a gradient descent optimizer to solve

$\displaystyle\underset{\xi}{\text{argmin}}||\mathbf{u}-\mathbf{u}_{data}||^{2}+\lambda||\mathcal{L}_{\xi}[\mathbf{u}]-\mathbf{f}||^{2}.$
(37)

This approach resembles currently popular approaches such as physics-informed
neural networks and is simple to implement in machine learning libraries using
automatic differentiation to implement first-order optimization schemes.
However, $\lambda$ becomes a hyperparameter introducing well-known issues with
training, and ultimately results in $\mathcal{L}_{\xi}[\mathbf{u}]=\mathbf{f}$
holding only to within optimization error. We instead enforce the equality
constraint exactly, introducing Lagrange multipliers $\mathbf{\lambda}$ as
follows.

$\displaystyle\underset{\xi,\mathbf{u},\mathbf{\lambda}}{\text{argmin}}\,\mathbf{L}_{\xi}(\mathbf{u},\mathbf{\lambda})$
(38) $\displaystyle\text{such that }\mathcal{L}_{\xi}[\mathbf{u}]=\mathbf{f},$
(39) $\displaystyle\text{where
}\mathbf{L}_{\xi}(\mathbf{u},\mathbf{\lambda})=||\mathbf{u}-\mathbf{u}_{data}||^{2}+\mathbf{\lambda}^{\intercal}\left(\mathcal{L}_{\xi}[\mathbf{u}]-\mathbf{f}\right)$
(40)

Stationarity of the Karush-Kuhn Tucker conditions requires that the gradient
of $\mathbf{L}_{\xi}(\mathbf{u},\mathbf{\lambda})$ with respect to
$\mathbf{u}$, $\mathbf{\lambda}$ and $\xi$ be zero. This yields the following
set of three necessary conditions for a minimizer of (38), which we will
iteratively solve. If this fixed point iteration converges and all three are
satisfied, than one has obtained a minimizer.

  * 1.

Forward problem: The condition
$\partial_{\mathbf{\lambda}}\mathbf{L}_{\xi}(\mathbf{u},\mathbf{\lambda})=0$
requires that the forward problem is solved:
$\mathcal{L}_{\xi}[\mathbf{u}]=\mathbf{f}$. Assuming $\xi$ fixed, one may
solve with a Newton-Rhapson method, requiring calculation of the Jacobian of
$\mathcal{L}_{\xi}$, which we denote $J_{\xi,\bm{u}}$. Following the analysis
in Section 3, this problem is guaranteed to be solvable provided $\epsilon
L_{n}<1$.

  * 2.

Adjoint problem: Enforcing
$\partial_{\mathbf{u}}\mathbf{L}_{\xi}(\mathbf{u},\mathbf{\lambda})=0$
provides the linear adjoint problem
$J_{\xi,\bm{u}}^{\intercal}\mathbf{\lambda}=-2(\mathbf{u}-\mathbf{u}_{data})$
for the Lagrange multipliers. Having solved the forward problem in the
previous step, the Jacobian is already available.

  * 3.

Model update: The remaining condition
$\partial_{\xi}\mathbf{L}_{\xi}(\mathbf{u},\mathbf{\lambda})=0$ does not
readily admit solution with second-order optimizers due to the neural networks
embedded in the nonlinear perturbations of (27). It is well-known that neural
networks admit a complex optimization landscape due to linear dependence with
many suboptimal local minima. With this in mind, we apply a single step of a
first-order gradient optimizer instead, providing a small perturbation of the
model at each iteration.

This process ensures that the physics imposed by the carefully designed model
formed in the previous sections hold to machine precision at each iteration,
even in scenarios with limited training data. Asymptotically, the added
complexity compared to typical gradient-descent approaches is that of solving
the forward problem at each iteration, which is generally inexpensive for the
low-dimensional models under consideration and converges rapidly. The
remaining complexity lies in calculating the relevant derivatives for Newton,
which may be simply calculated with the same automatic differentiation used to
perform the gradient descent step.

Considering that the assumed model form is nonlinear, it is necessary to train
simultaneously on a variety of boundary conditions to learn the nonlinear
response across a range of conditions - we must therefore assimilate multiple
solutions to the nonlinear problem simultaneously. With that in mind, we
present in Algorithm 1 a batch-training strategy for handling training data
$\mathcal{D}=\left\\{\bm{u}_{data,i}\right\\}_{i=1}^{N_{batches}}$ consisting
of $N_{batches}$ solutions.

Data: Training solutions $\mathcal{D}$, parameter initialization $\xi_{0}$,
tolerance $\epsilon_{TOL}$

Result: $\xi$

for  _$e\in\\{1,...,\texttt{{epochs}}\\}$_ do

$\mathbf{u}\leftarrow\mathbf{0}$;

for  _$\mathbf{u}_{data,i}\in\mathcal{D}$_ do

Solve forward problem;

while  _$||\mathcal{L}_{\xi}[\mathbf{u}]-f|| >\epsilon_{TOL}$_ do

$J_{\xi,\bm{u}}\leftarrow\nabla_{\mathbf{u}}\mathcal{L}_{\xi}(\mathbf{u})$;

Solve $\Delta=-(J_{\xi,\bm{u}})^{-1}\mathcal{L}_{\xi}(\mathbf{u})$;

$\mathbf{u}\leftarrow\mathbf{u}+\Delta$;

end while

Solve adjoint problem;

$\lambda\leftarrow
2(J_{\xi,\bm{u}})^{-\intercal}\left(\mathbf{u}_{data,i}-\mathbf{u}\right)$;

Apply gradient descent update;

$\xi\leftarrow GD(\mathbf{L}_{\xi}(\mathbf{u},\mathbf{\lambda}),\xi)$;

end for

end for

Algorithm 1 Application of equality constrained optimizer batched over
multiple solutions $\mathbf{u}_{i,data}$. By $GD(l,p)$ we denote a gradient
descent update of the parameters $p$ to minimize the loss $l$. For this work
we apply the Adam optimizer [40].

##  5 Numerical results

Figure 2: We consider a cylindrical inclusion in a 2D domain for steady
diffusion problems in subsurface flow and magnetostatics (left), considering
diffusion coefficients which vary with the boundary conditions magnitude to
obtain a nonlinear response. The high-fidelity solution is coarsened,
obtaining a target solution on the corresponding $k-$chains (center). The
coarsened graph model is trained to reproduce moments on each subdomain
(right). Pictured here is the resulting magnetic field defined on coarsened
$0-$chains, with nearest neighbor interpolant of surrogate solution.

In the remainder we present results for two canonical $H(div)$ and $H(curl)$
problems from subsurface flow and magnetostatics. In both examples we consider
a cylindrical inclusion of radius $a$ embedded within the unit square centered
at the origin, whose material properties are prescribed by

$\mu_{\alpha}(\mathbf{x})=\begin{cases}\alpha,&\text{if }||\mathbf{x}||<a\\\
1,&\text{otherwise}\end{cases}$ (41)

Treatment of the material interface at $||\mathbf{x}||=a$ without interesting
spurious oscillations is a hallmark of mimetic discretizations and stems for
the exact treatment of interface conditions via Stokes theorem.

For both problems we will impose boundary conditions by identifying
appropriate cochains on the boundary, replacing their corresponding row of the
Jacobian matrix with a zero vector with one on the diagonal and setting the
desired value on the right hand side. For the introductory nature of this
paper this is sufficient, but we note that the imposition of boundary
conditions is a rich topic in the discrete exterior calculus literature, with
more complex applications requiring a deeper consideration of the interaction
between boundary conditions and the discrete exterior calculus [24, 42, 12].

###  5.1 Problem 1: Darcy

We consider as training data solutions to the system of equations

$\displaystyle\mathbf{F}+\mu_{\alpha}\nabla\phi=0$
$\displaystyle\qquad\bm{x}\in\Omega$ $\displaystyle\nabla\cdot\mathbf{F}=f$
$\displaystyle\qquad\bm{x}\in\Omega$ $\displaystyle\mathbf{F}\cdot\hat{n}=g$
$\displaystyle\qquad\bm{x}\in\partial\Omega$

obtained via the scheme in [50]. Boundary conditions are imposed upon the
1-cochains encoding the fluxes of $\mathbf{F}$ through subdomain boundaries,
and we stress that the resulting DEC method will guarantee that
$\int_{\Omega}\nabla\cdot\mathbf{F}dx=\int_{\partial\Omega}gdA$. For this
problem we consider as neural network a dense elu network [23] with two hidden
layers of width five, initialized with the He initializer [35]. We consider a
$50^{2}$ fine mesh to generate training data.

For this problem we will consider two scenarios.

Darcy problem 1 (D1): We take $f=0$ and $g=<1,0>\cdot<1,0>$ consistent with
applying a unit horizontal flux and study the effect of varying $\alpha$. For
a fixed $\alpha$ the PDE is linear and will not require a contribution from
the neural network, and we will use this to gauge the methods ability to
recover a PDE discretization consistent with a single solution (i.e.
$N_{data}=1$). Figure 3 provides a summary and discussion of results.

Figure 3: Solutions for problem D1. Top: RMS error vs number of subdomains in
coarsened mesh, using learning rate of $0.05,0.005,0.0005$ for
$3^{2},6,^{2},12^{2}$ subdomains, respectively, demonstrating ability to
capture training data to machine precision. Center: piecewise constant plots
of pressure for obtained model. Bottom: Comparison of pressure along $y=0.5$
for model solution (orange) and Delaunay interpolant (blue) to training
solution (dashed green) for increasing resolution, using coarsened subdomain
centroids to compute interpolant.

Darcy problem 2 (D2): We take $f=0$ and $g=<\alpha,0>\cdot<1,0>$ consistent
with applying a horizontal flux which scales with the diffusion parameter.
This corresponds to a material which becomes more conductive as the
current/flux increases. To capture the nonlinear behavior of this problem will
require the neural network to learn fluxes which depend upon the magnitude of
the potential. Figure 4 provides a summary and discussion of results.

Figure 4: Solutions for problem D2. A single nonlinear model is extracted
which may be solved for increasing $\alpha$. Top: Comparison of pressure along
$y=0.5$ for increasing $\alpha$. True solution given by solid line, while
learned surrogate is given by dashed line. Bottom: Convergence of nonlinear
model over three solutions during training.

###  5.2 Problem 2: Magnetostatics

We consider as training data solutions to the system of equations

$\displaystyle\nabla\times\mathbf{J}=\mathbf{f}$
$\displaystyle\qquad\bm{x}\in\Omega$ $\displaystyle\nabla\cdot\mathbf{B}=0$
$\displaystyle\qquad\bm{x}\in\Omega$
$\displaystyle\mathbf{B}=\mu_{\alpha}\mathbf{J}$
$\displaystyle\qquad\bm{x}\in\Omega$ $\displaystyle\mathbf{B}\cdot\hat{n}=g$
$\displaystyle\qquad\bm{x}\in\partial\Omega$

obtained via the scheme in [49]. For this 2D problem, the magnetic field may
be identified with a scalar ($\mathbf{B}=\left<0,0,B(x,y)\right>$). Boundary
conditions are imposed upon the 0-cochains encoding the fluxes of $\mathbf{B}$
through subdomain boundaries, and the 1-cochains encoding the magnetic
potential on the boundary are fixed to a value of zero. For this problem we
will consider only the nonlinear generalization of D2, taking $f=0$ and
$g=\alpha$ to obtain a nonlinear material whose permittivity increases with
the magnitude of applied magnetic field (Figure 5. For this problem, we take
as neural network a parametric ReLU activation and a single hidden layer of
width ten. In the same manner as the Darcy problem, the resulting model is
able to recover a range of $\alpha$’s and the resulting discontinuity in the
magnetic field.

Figure 5: Profile of magnetostatics solution along $y=0.5$ line for $5^{2}$
partitions. The trained model reproduces the coarsened moments of the problem
and avoid introduction of any spurious oscillations at the material interface.

##  6 Conclusions

We have presented a new data-driven exterior calculus which parameterizes the
classical graph calculus. The resulting framework allows training of discrete
exterior calculus operators which incorporates metric information from data
while preserving exact sequence structure and invertibility of Hodge
Laplacians - both of which are important for handling a range of physical
systems. Analysis of the resulting system shows that when nonlinear
perturbations are applied to a class of second-order elliptic problems the
system remains well-posed under general conditions, allowing the machine
learning of nonlinear elliptic systems. Numerical results demonstrate the
practical aspects of the approach.

For the sake of introductory exposition, we have restricted ourselves to
elementary elliptic problems in the current work. The DDEC framework however
may be applied to a broad range of more sophisticated problems; for example,
we are currently using it to discover surrogates for semiconductor physics
governed by nonlinear drift-diffusion equations (using [47]), and where the
resulting network model can be embedded within circuit simulators such as Xyce
[30, 36].

Abstractly, the DDEC framework provides a structure-preserving means of
parameterizing Dirichlet-to-Neumann maps [61] governing multiscale physics,
generalizing previous works restricted to resistor networks [25] while
supporting machine learning tasks. In future work we will provide details
regarding how this may be incorporated into a workflow to develop provably
stable multiscale models that preserve structure at both fine and coarse
scales.

## Acknowledgement

Sandia National Laboratories is a multimission laboratory managed and operated
by National Technology and Engineering Solutions of Sandia, LLC, a wholly
owned subsidiary of Honeywell International, Inc., for the U.S. Department of
Energy’s National Nuclear Security Administration under contract DE-NA0003525.
This paper describes objective technical results and analysis. Any subjective
views or opinions that might be expressed in the paper do not necessarily
represent the views of the U.S. Department of Energy or the United States
Government.

N. Trask has also been supported by the U.S. Department of Energy, Office of
Advanced Scientific Computing Research under the Collaboratory on Mathematics
and Physics-Informed Learning Machines for Multiscale and Multiphysics
Problems (PhILMs) project and the U.S. Department of Energy, Office of
Advanced Scientific Computing Research under the Early Career Research
Program. A. Huang has been supported under the Sandia National Laboratories
Laboratory Directed Research and Development (LDRD) program. The work of X. Hu
is partially supported by the National Science Foundation under grant
DMS-1812503 and CCF-1934553.

SAND Number: SAND2020-14261 O

## References

  * [1] James H. Adler, Casey Cavanaugh, Xiaozhe Hu, and Ludmil T. Zikatanov.  A finite-element framework for a mimetic finite-difference discretization of Maxwell’s equations.  arXiv:2012.03148 [cs, math], December 2020.
  * [2] Douglas Arnold, Richard Falk, and Ragnar Winther.  Finite element exterior calculus: From Hodge theory to numerical stability.  Bulletin of the American Mathematical Society, 47(2):281–354, 2010\.
  * [3] Douglas N Arnold.  Finite element exterior calculus.  SIAM, 2018.
  * [4] Douglas N Arnold, Pavel B Bochev, Richard B Lehoucq, Roy A Nicolaides, and Mikhail Shashkov.  Compatible spatial discretizations, volume 142.  Springer Science & Business Media, 2007.
  * [5] Douglas N. Arnold, Richard S. Falk, and Ragnar Winther.  Finite element exterior calculus, homological techniques, and applications.  Acta Numerica, 15:1–155, May 2006.
  * [6] Nathan Baker, Frank Alexander, Timo Bremer, Aric Hagberg, Yannis Kevrekidis, Habib Najm, Manish Parashar, Abani Patra, James Sethian, Stefan Wild, et al.  Workshop report on basic research needs for scientific machine learning: Core technologies for artificial intelligence.  Technical report, USDOE Office of Science (SC), Washington, DC (United States), 2019.
  * [7] Paul Bamberg and Shlomo Sternberg.  A Course in Mathematics for Students of Physics, volume 2.  1991\.
  * [8] Sergio Barbarossa and Stefania Sardellitti.  Topological Signal Processing over Simplicial Complexes.  arXiv:1907.11577 [eess], March 2020.
  * [9] Timothy Barth.  On the role of involutions in the discontinuous galerkin discretization of maxwell and magnetohydrodynamic systems.  In Compatible spatial discretizations, pages 69–88. Springer, 2006\.
  * [10] Peter Battaglia, Razvan Pascanu, Matthew Lai, Danilo Jimenez Rezende, et al.  Interaction networks for learning about objects, relations and physics.  Advances in neural information processing systems, 29:4502–4510, 2016.
  * [11] Federico Battiston, Giulia Cencetti, Iacopo Iacopini, Vito Latora, Maxime Lucas, Alice Patania, Jean-Gabriel Young, and Giovanni Petri.  Networks beyond pairwise interactions: Structure and dynamics.  Physics Reports, June 2020.
  * [12] René Beltman, MJH Anthonissen, and Barry Koren.  Conservative polytopal mimetic discretization of the incompressible navier–stokes equations.  Journal of Computational and Applied Mathematics, 340:443–473, 2018\.
  * [13] A Bloch.  Electromechanical analogies and their use for the analysis of mechanical and electromechanical systems.  Journal of the Institution of Electrical Engineers-Part I: General, 92(52):157–169, 1945.
  * [14] P Bochev and A Robinson.  Matching algorithms with physics: exact sequences of finite element spaces.  Collected Lectures on the Preservation of Stability Under Discretization, edited by D. Estep and S. Tavener, SIAM, Philadelphia, 2001.
  * [15] Pavel B Bochev, Jonathan J Hu, Christopher M Siefert, and Raymond S Tuminaro.  An algebraic multigrid approach based on a compatible gauge reformulation of maxwell’s equations.  SIAM Journal on Scientific Computing, 31(1):557–583, 2008.
  * [16] Pavel B Bochev and James M Hyman.  Principles of mimetic discretizations of differential operators.  In Compatible spatial discretizations, pages 89–119. Springer, 2006\.
  * [17] Alain Bossavit.  Computational electromagnetism: variational formulations, complementarity, edge elements.  Academic Press, 1998.
  * [18] Alain Bossavit.  ” stiff” problems in eddy-current theory and the regularization of maxwell’s equations.  IEEE transactions on magnetics, 37(5):3542–3545, 2001.
  * [19] Peter C Breedveld.  Multibond graph elements in physical systems theory.  Journal of the Franklin Institute, 319(1-2):1–36, 1985.
  * [20] Michael B Chang, Tomer Ullman, Antonio Torralba, and Joshua B Tenenbaum.  A compositional object-based approach to learning physical dynamics.  arXiv preprint arXiv:1612.00341, 2016.
  * [21] Qun Chen, Rong-Huan Fu, and Yun-Chao Xu.  Electrical circuit analogy for heat transfer analysis and optimization in heat exchanger networks.  Applied Energy, 139:81–92, 2015.
  * [22] Fan RK Chung and Fan Chung Graham.  Spectral graph theory.  Number 92. American Mathematical Soc., 1997.
  * [23] Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter.  Fast and accurate deep network learning by exponential linear units (elus).  arXiv preprint arXiv:1511.07289, 2015.
  * [24] Lorenzo Codecasa.  Refoundation of the cell method using augmented dual grids.  IEEE transactions on magnetics, 50(2):497–500, 2014.
  * [25] Edward B Curtis and James A Morrow.  The dirichlet to neumann map for a resistor network.  SIAM Journal on Applied Mathematics, 51(4):1011–1029, 1991.
  * [26] Constantine M Dafermos.  Quasilinear hyperbolic systems with involutions.  Archive for Rational Mechanics and Analysis, 94(4):373–389, 1986\.
  * [27] Mathieu Desbrun, Anil N. Hirani, Melvin Leok, and Jerrold E. Marsden.  Discrete exterior calculus.  2005\.
  * [28] Beno Eckmann.  Harmonische funktionen und randwertaufgaben in einem komplex.  Commentarii Mathematici Helvetici, 17(1):240–255, 1944.
  * [29] Floyd A Firestone.  A new analogy between mechanical and electrical systems.  The Journal of the Acoustical Society of America, 4(3):249–267, 1933.
  * [30] Xujiao Gao, Andy Huang, Nathaniel Trask, and Shahed Reza.  Physics-informed graph neural network for circuit compact model development.  In 2020 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD), pages 359–362. IEEE.
  * [31] David Gilbarg and Neil S Trudinger.  Elliptic partial differential equations of second order.  springer, 2015.
  * [32] Marco Gori, Gabriele Monfardini, and Franco Scarselli.  A new model for learning in graph domains.  In Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005., volume 2, pages 729–734. IEEE, 2005.
  * [33] Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  In Advances in neural information processing systems, pages 1024–1034, 2017.
  * [34] William L Hamilton, Rex Ying, and Jure Leskovec.  Representation learning on graphs: Methods and applications.  arXiv preprint arXiv:1709.05584, 2017.
  * [35] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In Proceedings of the IEEE international conference on computer vision, pages 1026–1034, 2015.
  * [36] S Hutchinson, E Keiter, R Hoekstra, H Watts, A Waters, T Russo, R Schells, S Wix, and C Bogdan.  The xyce™ parallel electronic simulator–an overview.  In Parallel Computing: Advances and Current Issues, pages 165–172. World Scientific, 2002.
  * [37] Xiaoye Jiang, Lek-Heng Lim, Yuan Yao, and Yinyu Ye.  Statistical ranking and combinatorial Hodge theory.  Mathematical Programming, 127(1):203–244, 2011.
  * [38] Xiaoye Jiang, Lek-Heng Lim, Yuan Yao, and Yinyu Ye.  Statistical ranking and combinatorial hodge theory.  Mathematical Programming, 127(1):203–244, 2011.
  * [39] George Karypis and Vipin Kumar.  Metis: A software package for partitioning unstructured graphs, partitioning meshes, and computing fill-reducing orderings of sparse matrices.  1997\.
  * [40] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [41] HE Koenig and WA Blackwell.  Linear graph theory-a fundamental engineering discipline.  IRE Transactions on Education, 3(2):42–49, 1960.
  * [42] Jasper Kreeft, Artur Palha, and Marc Gerritsma.  Mimetic framework on curvilinear quadrilaterals of arbitrary order.  arXiv preprint arXiv:1111.4304, 2011.
  * [43] Isaac E Lagaris, Aristidis Likas, and Dimitrios I Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE transactions on neural networks, 9(5):987–1000, 1998.
  * [44] Lek-Heng Lim.  Hodge Laplacians on Graphs.  SIAM Review, 62(3):685–715, January 2020.
  * [45] S. Maletić and M. Rajković.  Combinatorial Laplacian and entropy of simplicial complexes associated with complex networks.  The European Physical Journal Special Topics, 212(1):77–97, September 2012.
  * [46] Abubakr Muhammad and Magnus Egerstedt.  Control Using Higher Order Laplacians in Network Topologies.  page 15, 2006.
  * [47] Lawrence Musson, Xujiao Gao, Mihai Negoita, Andy Huang, and Gary L Hennigan.  Charon: A radiation aware massively parallel tcad modeling code.  Technical report, Sandia National Lab.(SNL-NM), Albuquerque, NM (United States), 2018.
  * [48] Jean-Claude Nédélec.  Mixed finite elements in $\mathbb{R}^{3}$.  Numerische Mathematik, 35(3):315–341, 1980.
  * [49] R. A. Nicolaides.  Direct Discretization of Planar Div-Curl Problems.  SIAM Journal on Numerical Analysis, 29(1):32–56, February 1992\.
  * [50] Roy A Nicolaides and Kathryn A Trapp.  Covolume discretization of differential forms.  In Compatible spatial discretizations, pages 161–171. Springer, 2006.
  * [51] Kwang W Oh, Kangsun Lee, Byungwook Ahn, and Edward P Furlani.  Design of pressure-driven microfluidic networks using electric circuit analogy.  Lab on a Chip, 12(3):515–545, 2012.
  * [52] Mario Ohlberger and Stephan Rave.  Reduced basis methods: Success, limitations and future challenges.  arXiv preprint arXiv:1511.02021, 2015.
  * [53] Alfio Quarteroni, Andrea Manzoni, and Federico Negri.  Reduced basis methods for partial differential equations: an introduction, volume 92.  Springer, 2015.
  * [54] Maziar Raissi, Paris Perdikaris, and George E 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.
  * [55] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations.  arXiv preprint arXiv:1711.10561, 2017.
  * [56] Frank Rosenblatt.  Principles of neurodynamics. perceptrons and the theory of brain mechanisms.  Technical report, Cornell Aeronautical Lab Inc Buffalo NY, 1961.
  * [57] Kevin Scaman and Aladin Virmaux.  Lipschitz regularity of deep neural networks: Analysis and efficient estimation.  arXiv:1805.10965 [cs, stat], October 2019.
  * [58] Michael T. Schaub, Austin R. Benson, Paul Horn, Gabor Lippner, and Ali Jadbabaie.  Random Walks on Simplicial Complexes and the Normalized Hodge 1-Laplacian.  SIAM Review, 62(2):353–391, January 2020.
  * [59] S. Smale.  On the mathematical foundations of electrical circuit theory.  J. Differential Geom., 7(1-2):193–210, 1972.
  * [60] Malcolm C Smith.  Synthesis of mechanical networks: the inerter.  IEEE Transactions on automatic control, 47(10):1648–1662, 2002\.
  * [61] John Sylvester and Gunther Uhlmann.  The dirichlet to neumann map and applications.  In Inverse problems in partial differential equations, volume 42, page 101. SIAM Publications, Philadelphia, 1990.
  * [62] Joaquín J. Torres and Ginestra Bianconi.  Simplicial complexes: Higher-order spectral dimension and dynamics.  Journal of Physics: Complexity, 1(1):015002, May 2020.
  * [63] Nikolaos Vlassis, Ran Ma, and WaiChing Sun.  Geometric deep learning for computational mechanics part i: Anisotropic hyperelasticity.  arXiv preprint arXiv:2001.04292, 2020.
  * [64] Sifan Wang, Yujun Teng, and Paris Perdikaris.  Understanding and mitigating gradient pathologies in physics-informed neural networks.  arXiv preprint arXiv:2001.04536, 2020.
  * [65] Sifan Wang, Xinling Yu, and Paris Perdikaris.  When and why pinns fail to train: A neural tangent kernel perspective.  arXiv preprint arXiv:2007.14527, 2020.
  * [66] Hermann Weyl.  Repartición de corriente en una red conductora.  Revista Matemática Hispano-Americana, 5:153–164, 1923.
  * [67] Herbert S. Wilf and Frank Hilary.  ”Mathematical aspects of electrical network analysis”, volume 3 of SIAM-AMS proceedings.  American Mathematical Society, 1971.
  * [68] Kane Yee.  Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media.  IEEE Transactions on antennas and propagation, 14(3):302–307, 1966\.