Reversible and irreversible bracket-based dynamics for deep graph neural networks
Authors: Anthony Gruber
Kookjin Lee
Nathaniel Trask
Venue: NeurIPS 2023 poster (2023)
Keywords: graph neural networks
structure preserving machine learning
neural ordinary differential equations
hamiltonian dynamics
metriplectic dynamics
TL;DR: Novel bracket-inspired GNN architectures are versatile and establish role of reversibility/irreversibility in stability of deep GNNs.

Abstract: Recent works have shown that physics-inspired architectures allow the training of deep graph neural networks (GNNs) without oversmoothing. The role of these physics is unclear, however, with successful examples of both reversible (e.g., Hamiltonian) and irreversible (e.g., diffusion) phenomena producing comparable results despite diametrically opposed mechanisms, and further complications arising due to empirical departures from mathematical theory. This work presents a series of novel GNN architectures based upon structure-preserving bracket-based dynamical systems, which are provably guaranteed to either conserve energy or generate positive dissipation with increasing depth.  It is shown that the theoretically principled framework employed here allows for inherently explainable constructions, which contextualize departures from theory in current architectures and better elucidate the roles of reversibility and irreversibility in network performance. Code is available at the Github repository \url{https://github.com/natrask/BracketGraphs}.

Reversible and irreversible bracket-based dynamics for deep graph neural networks
Anthony Gruber Center for Computing Research Sandia National Laboratories Albuquerque, NM. USA adgrube@sandia gov
Kookjin Lee
School of Computing and Augmented Intelligence Arizona State University
Tempe, AZ. USA kookjin in lee@asu. edu
Nathaniel Trask
School of Engineering and Applied Science University of Pennsylvania Philadelphia, PAA. USA ntrask@seas upenn edu
Abstract
Recent works have shown that physics-inspired architectures allow the train- ing of deep graph neural networks (GNNs) without oversmoothing. The role these physics is unclear, however, with successful examples of both re- versible (e.g., Hamiltonian) and irreversible (e.g., diffusion) phenomena pro- ducing comparable results despite diametrically opposed mechanisms, and fur- ther complications arising due empirical departures from mathematical theory. This work presents a series of novel GNN architectures based upon structure- preserving bracket-basee dynamical systems, which provably guaranteed either conserve energy or generate positive dissipation with increasing depth. It is shown that the theoretically principled framework employed here allows for inherently explainable constructions, which contextualize departures from theory in current architectures and better elucidate the roles of reversibility and irre- versibility in network performance. Code is available at the Github repository htwpsttub.coo//tttattttBBe/eeeee wwt.a//nnntaaatt//aaaeeerrrrrrr
'K. Lee acknowledges the support from the U.S. National Science Foundation under grant CNS2210137.
1 "Sandia National Laboratories a multimission laboratory managed and operated by National Technology 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 expressed the paper do not necessarily represent the views of the U.S. Department of Energy or United States Government. This article has been co-authored by an employee of National Technology & Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE). The employee all right, title and interest in and to the article and is solely responsible for its contents. The United States Government retains and publisher, by accepting the article for publication, acknowledges that the United States Government retains a non non-exclusive. paid-up, irrevocable, world-widd wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan https://www.ents https//wwwweeerry..ae wwwp://wwgg.ggo//wwoooooooooooo gow/downloads/dos dowaa///oo--uiii--ccc-sses
The work of N. Trask and Gruber supported by the U.S. Department Energy. Office of Advanced Computing Research under the "Scalable and Efficient Algorithms- Causal Reasoning, Operators and Graphs" (SEA-CROGS) project, the DoE Early Career Research Program, and the John von Neumann fellowship Sandia.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Introduction
Graph neural networks (GNNs) have emerged as a powerful learning paradigm able to treat un- structured data and extract subject-relation"caauss relationships while imparting inductive biases which preserve invariances through the underlying graph topology 2 4. This framework has proven effective for a wide range of both graph analytics and data-driven physics modeling problems. Despite successes, GNNs have generally struggle to achieve the improved performance with increasing depth typical of other architectures. Well-known pathologies, such as oversmoothing. oversquashing, bottlenecks, and aaploddnn/vanishiin gradients yield deep GNNs which are either unstable or lose performance as the number of layers increase
To combat this, number of works build architectures which mimic physical processes impart desirable numerical properties. For example, some works claim that posing message passing as either a diffusion process or reversible flow may promote stability or help retain information, respectively. These present opposite ends of spectrum between irreversible and reversible processes, which either dissipate or retain information. It unclear, however, what role (ir)reevrsibility plays [9. One could argue that dissipation entropically destroys information and could promote oversmoothing, so should be avoided. Alternatively, in dynamical systems theory, dissipation crucial to realize a low-dimensional attractor, and thus dissipation may play an important role in realizing dimensionality reduction. Moreover, recent work has shown that dissipative phenomena can actually sharpen information well as smooth it [0O), although this is not often noticed in practice since typical empirical tricks (batch norm, etc.) lead to a departure from the governing mathematical theory.
In physics, Poisson brackets and their metriplectic/port moott-Hamitooian generalization to dissipative systems provide an abstract framework for studying conservation and entropy production in dynamical systems. In this work, we construct four novel architectures which span the (ir)reversibility spectrum, using geometric brackets as means of parameterizing dynamics abstractly without empirically assuming a physical model. This relies on an application of the data-driven exterior calculus (DDEC) [1, which allows a reinterpretation of the message-passing and aggregation of graph attention networks 12] as the fluxes and conservation balances of physics simulators 13), providing simple but powerful framework for mathematical analysis. this context, we recast graph attention as an inner-product graph features, inducing graph derivative muulldigg-blokss which may be used to build geometric brackets. In the process, we generalize classical graph attention [12] to higher-order clique cochains (e.g.. labels edges and loops). The four architectures proposed here scale with identical complexity to classical graph attention networks, and possess desirable properties that have proven elusive in current architectures. On the reversible and irreversible end of the spectrum we have Hamiltonian and Gradient networks. In the middle of the spectrum, Double Bracket and Metriplectic architectures combine both reversibility and irreversibility, dissipating energy to either the environment or an entropic variable, respectively, in a manner consistent with the second law of thermodynamics. We summarize these brackets in Table I providing diagram of their architecture in Figure
Primary contributions:
Theoretical analysis of GAT in terms of exterior calculus. Using DDEC we establish a unified framework for construction and analysis of pressage-passing graph attention networks, and provide an extensive introductory primer to the theory in the appendices. In this setting, we show that with our modified attention mechanism, GATs amount to a diffusion process for special choice of activation and weights.
Generalized attention mechanism. Within this framework, we obtain a natural and flexible extension of graph attention from nodal features higher order cliques (e.g. edge features). We show attention and introduce a novel and
must have symmetric numerator to be formally srructue--reesrvvig flexible graph attention mechanism parameterized in terms of learnable inner products on nodes and edges.
Novel prructure-presevingg extensions. We develop four GNN architectures based upon bracket- based dynamical systems. In the metriplectic case, we obtain the first architecture with linear complexity in the size of the graph while previous works are O(N$).
Unified evaluation of dissipation. We use these architectures to systematically evaluate the role of (ir)reveriiiility in the performance of deep GNNs. We observe best performance for partially
2

---

Formalism Hamiltonian
Equation
Requirements
Completeness complete
Character conservative
x = fx, E}
L = -L,
Jacobi's identity
Gradient
X xx Ej
M' M
incomplete incomplete complete
totally dissipative
Double Bracket
* = xx E) + [xx, EJ)
= L
partially dissipative
Metriplectic
& - [X, E) + []
L = L.M = M.
partially dissipative
LVS = MVF=
Table I: The abstract bracket formulations employed in this work. Here x represents a state variable, while E E S S(x) are energy and entropy functions. "Conservative" indicates purely reversible motion, "totally dissipative" indicates purely irreversible motion, and "partially dissipative" indicates motion which either dissipates E' (in the double bracket case) or generates S (in the metriplectic case).
dissipative systems, indicating that a combination both reversibility and irreversibility are important. Pure diffusion is the least performant across all benchmarks. For physics-based problems including optimal control, there a distinct improvement. All models provide near state-of-the-art performance and marked improvements over black-box GAT/NODE networks.
2 Previous works
Neural ODEs: Many works use neural networks to fit dynamics of the form à = fæx 0) to time series data. Model calibration (e.g., UDE [4). dictionary-basee learning (e.g.. SINDy [5]), and neural ordinary differential equations (e.g., NODE [6]) pose a spectrum of inductive biases requiring progressively less domain expertise. Structure preservvttonssstioo provides a means of obtaining stable training without requiring domain knowledge, ideally achieving the flexibility of NODE with the robustness of UDE/SINDY. The current work learns dynamics on graph while using a modern NODE library to exploit the improved accuracy of high-order integrators [17 19)
Structure-preservingg dense networks: For dense networks,
relatively straightforward to parameterize reversible dynamics, see for example: Hamiltonian neural networks [20 22 23),
Hamiltonian generative networks (24), Hamiltonian with Control (SymODEN) 25 Deep Lagrangian networks [26] and Lagrangian neural networks
(27). Structure-preservinn extensions to dissipative systems are more challenging, particularly for metriplectic dynamics 28 which require a delicate degeneracy condition to preserve discrete notions of the first and second laws of thermodynamics. For
dense networks such constructions are intensive, suffering from O(N*) complexity the number of features [29 30 31,
In the graph setting we avoid this and achieve linear complexity by exploiting exact sequence structure, Alternative dissipative frameworks include Dissipative SymODEN [32 and Hort-Hamiltonian [33] We choose to focus metriplectic parameterizations due to their broad potential impact data-driven physics modeling, and ability to naturally treat fluctuations in multiscale systems 34].
Physics-ifformed vs struucture-preserving:: (Physics-informee learning imposes physics by penalty, adding regularizer corresponding to physics residual. The technique is simple to implement and has been successfully applied to solve a range of PDEs 335, discover data-driven models complement
first-principles simulators [36] 37 38) learn metriplectic dynamics (39], and perform uncertainty
quantification (40 41). Penalization poses a multiobjective optimization problem, however, with
parameters weighting competing objectives inducing pathologies during training, often resulting in physics being imposed to a coarse tolerance and qualitatively poor predictions [42 43 In contrast, prructure-preserviie architectures exactly impose physics by construction via carefully designed networks. Several works have shown that penalty-based approaches suffer in comparison, with structure-pressevvation providing improved long term stability, extrapolation and physical realizability,
Structure-preservving graph networks: Several works use discretizations of specific PDEs to combat oversmoothing or exploding/vanishiin gradients, e.g. telegraph equations 44 or various reaction-diffusion systems [45). Several works develop Hamiltonian flows on graphs [4 47.. For metriplectic dynamics, [88 poses a penalty based formulation on graphs. We particularly focus on
GRAND, which poses graph learning as a diffusive process [49], using a similar exterior calculus framework and interpreting attention as diffusion coefficient. We show in Appendix A.5 that their

---

Graph derivalives
Inpub/Cupput
da: dl
Attention
Metric acjoints
Bracket matrices
a(ta)
Encoder
#(ts)
d #(fa.))
Decoder
>n((te)))
Ao A
dyj
EnergyEEreoyy
Gradients
((ts)
Intogrator
S,S
Figure I: A diagrammatic illustration of the bracket-based architectures introduced in Section
analysis fails to account for the asymmetry in the attention mechanism, leading to a departure from the governing theory. To account for this, introduce modified attention mechanism which retains interpretation as part of diffusion PDE. In this purely irreversible case, adherence the theory provides improved results, GRAND's success is driven by something other than structure-preservatie preservation.
of interest whether
3 Theory and fundamentals
Here we introduce the two essential ingredients to our approach: bracket-based dynamical systems for neural differential equations, and the data-driven exterior calculus which enables their construction.
thorough introduction to this material provided in Appendices A.I A.3 and A.2
Bracket-based dynamics: Originally introduced as an extension of Hamilooiann/Larranne namies to include dissipation [50]. bracket formulations are used inform a dynamical system with certain structural properties, e.g., time-reversibility. invariant differential forms, or property preservation. Even without dissipation, bracket formulations may compactly describe dynamics while preserving core mathematical properties, making them ideal for designing neural architectures.
dy-
Bracket formulations are usually specified via some combination of reversible brackets {F,G) (VF, LVG) and irreversible brackets [F, G| F. MVG). ([F,G)} (VF, for poten- tially state-dependent operators L = and M' M. The particular brackets which are used in the present network architectures summarized in Table [ Note that complete systems the dynamical extensions isolated thermodynamical systems: they conserve energy and produce entropy, with nothing lost the ambient environment. Conversely, incomplete systems do not account for any lost energy: they only require that it vanish in a prescribed way. The choice of completeness an application-deeendent modeling assumption.
Exterior calculus: In the combinatorial Hodge theory [51]. an oriented graph g = (V,E) carries sets of k-cllqqes, denoted Gk. which are collections of ordered subgraphs generated by (k: 1) nodes. This induces natural exterior derivative operators de Q2 + Q4++1, acting on the spaces functions on Ge. which the signed incidence matrices between k-cliques and (k + 1)-cliques. An explicit representation these derivatives is given in Appendix [A.] from which is easy to check the exact sequence property dx+1 d, for any k. This yields discrete de Rham complex (Figure 2). Moreover, given a choice of inner product (say, (2) on Sr. there is an obvious dual de Rham complex which comes directly from adjointness. In particular, can define dual derivatives
the graph G
di P2+1 + 24 via the equality
(defsg)<-1= (f, dig)
from which nontrivial results such as the Hodge decomposition. Poincaré inequality, and coerciv- ity/nnvrtiibility of the Hodge Laplacian A di dx-1d follow (see e.g. (II). Using the derivatives dy. di. possible to build compatible discretizations of PDEs on G which are guaranteed preserve exactness properties such as, e.g., d, do curl grad 0.
The choice of inner product ( thus induces a definition of the dual derivatives dz. In the graph setting 52, one typically selects the e inner product, obtaining the adjoints of the signed incidence matrices as di dI. By instead working with the modified inner product (v, (v.w) v A for machine-learnate Ag. obtain di A dī Ak+1 (see Appendix A.3 This parameterization

---

do
d1
d2
dx 1
20
Ω
O
d
d
d2
d k-1
Ao
A1
A2
Ak
2
Ω
d
d
di
d
Figure A commutative diagram illustrating the relationship between the graph derivatives dæ, their e adjoints df and the learnable adjoints d These operators form de Rham complex due to the exact sequence property di+1 o d, = dd d d; o O. We show that the learnable Ag may encode attention mechanisms, without impacting the preservation of exact sequence structure.
inherits the exact sequence property from the graph topology encoded in dx while allowing for incorporation of geometric information from data. This leads directly to the following result, which holds for any (potentially ceature-dependent symmetric positive definite matrix Ag.
Theorem 3.1. The dual derivatives di. Sk+1 + Sy adjoint to dg 94 + SK+1 with respect to the learnable inner products A S2 S2g. satisfy exact sequence property.
Proof.
di di A df AAA, 'fA A (dadk-1)) Au+1 =0.
As will be shown in Section 4 by encoding graph attention into the As, may exploit the exact sequence property to obtain symmetric positive definite diffusion operators, as well as conduct the cancellations necessary enforce degeneracy conditions necessary for metriplectic dynamics.
For thorough review DDEC, we direct readers to Appendix A. and For exterior calculus in
topological data analysis see [22] and an overview in the context of PDEs see 53 54
4 Structure--resseringg bracket parameterizations
We next summarize properties of the bracket dynamics introduced in Section and displayed in Table I postponing details and rigorous discussion Appendices A.3 denote node-edge feature pairs, the following operators will be used to generate our brackets.
A.3 and A.6 Letting = (q. p)
do
A 0
0
di
O di + dodg)
= -
0 Aidid1A
do
A
As mentioned before, the inner products Ag, A1, A2 on Q4 which induce the dual derivatives di di. are chosen in such a way that their combination generalizes graph attention mechanism. The precise details of this construction are given below, and its relationship
the standard GAT network from
12] is shown in Appendix A.5 Notice that L = -L, L, while G = G, M M are positive semi- definite with respect to the block block-diagonal inner product (--)) defined by A diag (Ag. A1) (details are provided in Appendix A.6). Therefore, L generates purely reversible (Hamiltonian) dynamics and
G.M generate irreversible (dissipative) ones. Additionally, note that nate-dependence in L, M, G enters only through the adjoint differential operators, meaning that any structural properties induced by the topology of the graph G (such as the exact sequence property mentioned in Theorem 3.1 are automatically preserved.
Remark 4.1. Strictly speaking, L is guaranteed to be truly Hamiltonian system only when dă is state-independent, since may otherwise fail to satisfy Jacobi's 's identity. On the other hand, energy conservation is always guaranteed due to the fact that L is skew-adjoint.
In addition to the bracket matrices M, G. it is necessary to have access to energy and entropy functions S S and their associated functional derivatives with respect to the inner product on Po @ Ω defined by A. For the Hamiltonian, gradient, and double brackets, chosen simply as the "total kinetic energy"
E(q.p) p) - ((a1 +lp1²)
1 14:
Σ
Pa
5

---

whose A- A-gradient (computed in Appendix A.6 is just () = ( AT metriplectic bracket uses parameterizations is deferred to later in this Section.
Since the which are more involved, discussion this case
Attention learnable inner product: Before describing the dynamics, it remains diseuss how the matrices A; 0 < i S. 2, are computed in practice, and how they relate to the idea of graph attention. Recall that if ny > 0 denotes the nodal feature dimension, graph attention mechanism takes the form a(q; qi = (āzz) -; f (ā:j) for some differentiable pre-attention function à X ny - R (e.g., for scaled dot product 55 one typically represents a(q, 4j) a softmax, so that f exp(g)). This suggests decomposition a(q: qj) A, Ai where Ag = (d0.s) is diagonal on nodes and A = (a1,5)) is diagonal on edges,
d0,3
(ā, JEN())
Σ
d1,5
(q: 9,)).
Treating the numerator and denominator of the standard attention mechanism separately in Ao, A1 allows for flexible and theoretically sound incorporation of graph attention directly into the adjoint differential operators on G. In particular, if A is symmetric with respect to edge-orientation and is an edge feature which antisymmetric, follows that
(di), = (A, de A.P))
Σ a (q: 9j) P5s JEN()
which is just graph attention combined with edge aggregation. This makes it possible give the following informal statement regarding graph attention networks which is explained and proven in Appendix A..5
Remark 4.2. The GAT layer from 12 is almost the forward Euler discretization of metric heat equation.
The "almost" appearing here has to do with the fact that (1) the attentional numerator f (ā(qi, q:) is generally asymmetric in i, and is therefore symmetrized by the divergence operator dj (2) the activation function between layers is not included, and (3) learnable weight matrices W GAT are set to the identity.
Remark 4.3. The interprelation of graph attention combination of learnable inner products admits direct generalization to higher-order cliques, which discussed in Appendix A.4
Hamiltonian case: A purely conservative system is generated by solving x V(x))
di
A 0
q
(
A A p
do
A-
-1 A
This is noncanonical Hamiltonian system which generates purely reversible flow. In particular, it can be shown that
(x) = (x, VE(x)) (L(x)VE(x), VE(x)) (())E))
so that energy is conserved due to the skew-adjointness of L.
Gradient case: On the opposite end of the spectrum generalized gradient flows, which are totally dissipative. Consider solving x -G(x)V or
- Ao 0 A 0
A 0
9 P 2))((1)) A
= A A 1.Ai1p)
This system metric diffusion process nodes and edges separately. Moreover, it corresponds a generalized gradient flow, since
É(x) (x. Ex))) = (G(x)V VE(x)) = =VIV((x 00
due to the self-adjoint and positive semi-definite nature G.
Remark 4.4. The architecture in GRAND [49) is almost a gradient flow, however the pre-attention mechanism lacks the requisite symmetry to formally induce a valid inner product.

---

Double bracket case: Another useful formulation for incomplete systems is the so-called double- bracket formalism. Consider solving LVE LLVE. or
0 - do
00)
An A
(-6 do
Ao 9 A
d - 9 dadÿ A
A
- -dog
A
This provides a dissipative relationship which preserves the Casimirs of the Poisson bracket generated by L, since LVC implies L'VC 0. particular, follows that
(xx) = (x, V FE()) = (L(x) + V = [x) S O,
since L is skew-adjoint and therefore L is self- adjont.
Remark 4.5. It is interesting to note that the matrix L is essentially Dirac operator (square root of the Hodge Laplacian A = (d + d²)²) restricted cliques of degree most 1. However, here L = -A so that L is in some sense "pure imaginary".
Metriplectic case: Metriplectic systems expressible X LVE MVS where S energy resp. entropy functions which satisfy the degeneracy conditions LVS = MVE = 0. One way of setting this up in the present case is define the energy and entropy functions
E(q. p) fE (s(q)) 9E (s (dod)p)),
S(q. p) = 9s (s (d)d)p)),
where s is sum aggregation over nodes resp. edges, fE R + R acts on node features, and IE. 99 R + R act on edge features. Denoting the "all ones" vector (of variable length) by 1, it shown in Appendix A.6 that the A A-gradients of energy and entropy can be computed as
AT VE(x)= A
SE (h(q)) Vge (h (dodf
VS(x) (A- df (A-'df41 Vgs (h (d) (dTddp)))
dodf do-
Similarly, it is shown in Appendix A.6 that the degeneracy conditions LVS = MV E 0 are satisfied by construction. Therefore, the governing dynamical system becomes
( (doA5'1e fi f + Aididid di 1 & Vgs (d(dip)))
-A & Vye (s (dodp))
LVE MVS
With this, it follows that the system obeys version of the first and second laws of thermodynamics,
(x) (x. VF(()) (LVF(x), VE(x)) (MVS(x), ( (vS(x), MVF(()) 0,
$(x) (xx VS(x)) = (LVE(x), VS())) + (M) V + ( 2
Remark 4.6. As seen in the increased complexity of this formulation, enforcing the degeneracy conditions necessary for metriplectic structure is nontrivial. This is accomplished presently via an application of the exact sequence property in Theorem 3.1 which we derive in Appendix AA66
Remark 4.7. It is worth mentioning that, similar to the other architectures presented in this Section, the metriplectic network proposed here exhibits linear O(N) scaling the graph size. This is in notable contrast to 29 371 which scale as O(N3).
5 Experiments
This section reports results on experiments designed to probe the influence bracket structure on trajectory prediction and nodal feature classification. Additional experimental details can be found in Appendix B each Table, orange indicates the best result by our models, and blue indicates the best of those compared. We consider both physical systems, where the role of structure preservation is explicit, as well as graph-analytic problems.
5.1 Damped double pendulum
As first experiment, consider applying one of these architectures to reproduce the trajectory of double pendulum with damping force proportional angular momenta of the pendulum masses (see Appendix B.] for details).

---

Double pendulum
MAE q
MAE P
Total MAE
NODE
0,0240 = 0.015
0.0299 + 0.0091
0.0269 = 0.012
NODE+AE
0.0532 + 0.029
0.0671 + 0.043
0.0602 + 0.035
Hamiltonian
0.00368 + 0.0015
0.00402 = 0.0015
0.00369 + 0.0013
Gradient 0.00762 0.0023
0.0339 + 0.012
0.0208 = 0.0067
Double Bracket
0.00584 = 0.0013
0.0183 + 0.0071
0.0120 + 0.0037
Metriplectic
0.00364 = 0.00064
0.00553 = 0.00029
0.00459 0.00020
Table Mean absolute errors (MAEs) of the network predictions in the damped double pendulum case, reported as atgee over runs.
Since this system is metriplectic when expressed in rosition-momeenuu-entrooo coordinates (c.f. (56]), it is useful to see if any of the brackets from Section Joctionn4caan adequately capture these dynamics without an entropic variable. The results of applying the architectures of Section 4 to reproduce trajectory of five periods are displayed in Table 2 alongside comparisons with a black-box NODE network and latent NODE with feature encoder/decoder, While each network capable of producing a small mean absolute eπтоr, it is clear that the metriplectic and Hamiltonian networks produce the most accurate trajectories. is remarkable both that the Hamiltonian bracket does so well here and that the gradient bracket does so poorly, being that the damped double pendulum system is quite dissipative. On the other hand, it unlikely to be only the feature encoder/decoder leading to good performance here, both the NODE and NODE+AE architectures perform worse on this task by about one order of magnitude.
5.2 MuJoCo Dynamics
Next we test the proposed models on more complex physical systems that are generated by the Multi-Joint dynamics with Contact (MuJoCo) physics simulator [57]. We consider the modified versions of Open AI Gym environments [23]: HalfCheetah, Hopper, and Swimmer.
We represent an object in an environment as a fully-connected graph, where a node corresponds to a position of body part or an
body part of the object and, thus, the nodal feature qi corresponds angle of joint As the edge features, a pair of nodal velocities Pa = (Uure(a): Vas((a)) are provided, where (a)((a) and Vas((a)) denote velocities of the source and destination nodes connected the edge.
Since the MuloCo environments contain an actor applying controls, additional control input is accounted for with an additive forcing term which is parameterized by multi-layer perceptron and introduced into the bracket-based dynamics models. See Appendix [B.2 for additional experimental details. The problem therefore consists of finding optimal control MLP, and evaluate the improvement which comes from representing the physics surrogate with bracket dynamics over NODE.
All models are trained via minimizing the MSE between the predicted positions q and the ground truth positions q and are tested the test set measured in the relative le norm, ||q - alle 2// |||p////ll2. Similar to the double pendulum experiments, all models are able to produce accurate predictions with around or less than 10% errors. While the gradient bracket makes little to no improvements over NODEs, the Hamiltonian, double, and metriplectic brackets produce more accurate predictions. Interestingly, the Hamiltonian bracket performs the best in this case as well, meaning that any dissipation present is effectively compensated for by the autoencoder which transforms the features.
unseen test set. Table 3 reports the eπrоrs of network predictions on
5.3 Node classification
Moving beyond physics- based examples, it remains to see how saacket-based architectures perform "black-box" node classification problems. Table and present results on common benchmark problems including the citation networks Cora [38 Citeseer [59]. and Pubmed [60], well as the coauthor graph, CoauthorCS [6]], and the Amazon co-purchasing graphs, Computer and Photo For comparison, we report results on standard GAT [2]. neural graph differential equation
62
Results of experiment with alternative embedding (i.e., q (gi, v)) are reported in Appendix B.2.2

---

Dataset
HalfCheetah
Hopper
Swimmer
NODE+AE
0.106 + 0.0011
0.0780 = 0.0021
0.0297 + 0.0036
Hamiltonian
0.0566 0.013
0.0279 + 0.0019
0.0122 0.00044
Gradient
0.105 + 0.0076
0.0848 + 0.0011
0.0290 + 0.0011
Double Bracket
0.0621 0.0096
0.0297 = 0.0048
0.0128 0.00070
Metriplectic
0.105 = 0.0091
0.0398 + 0.0057
0.0179 + 0.00059
Table 3: Relative error of network predictions for the MuloCo environment on the test set, reported as stdev over 4 runs.
Planetoid splits
CORA
CiteSeer
PubMed
GAT
82.8 0.5
79.0 = 0.5 79.9 0.3 72.5 = 0.5
69.5 = 0.9
GDE
83.8 0.5
GRAND-nl
83.6 0.5
70.8 1.1
79.7 + 0.3
Hamiltonian
7
73.0 1.2
78.5 + 0.3
Gradient
7
71.8 = 1.4
78.6 + 0.7
Double Bracket
62.6=
74.2 1.4
79.6 + 0.6
Metriplectic
1.0
60.5 = 1.1
69.8 + 0.7
Table 4: Test accuracy and standard deviations (averaged over 20 randomly initialized runs) using the original Planetoid train-valid-test splits. Comparisons use the numbers reported in
architecture (GDE) [63]. and the nonlinear GRAND architecture (GRAND-nI) from 44 which is closest ours. Since our experimental setting similar to that of (49, the numbers reported for GAT, GDE, and GRAND-nl are taken directly from this paper. Note that, despite the similar O(N) scaling the metriplectic architecture, the high dimension of the node and edge features on the latter three datasets led to trainable E, S' functions which exhausted the memory on our available machines, and therefore results are not reported for these cases. A full description of experimental details is provided in Appendix B.3
Remark 5.1. To highlight the effect of bracket structure on network performance, only minimal modifications are employed during network training. In particular, we do not include any additional regularization, positional encoding, graph rewiring. extraction of connected components, extra terms on the right-hand side, or early stopping. While it is likely that better classification performance could be achieved with some of these modifications included, it becomes very difficult to isolate the effect of seuuttueepreservatiinn A complete list of tunable hyperparameters is given in Appendix B.3
The results show different behavior produced by each bracket architecture. It is empirically clear that there is some value in full or partial reversibility, since the Hamiltonian and double bracket architectures both perform better than the corresponding gradient architecture on datasets such as Computer and Photo. Moreover, it appears that the partially reversible double bracket performs the best of the bracket architectures every case, which consistent with the idea that both reversible and irreversible dynamics are critical for capturing the behavior of general dynamical systems. Interestingly, the metriplectic bracket performs worse on these tasks by large margin. We conjecture this architecture may be harder to train for larger problems despite its O(N) complexity in the graph size, suggesting that more sophisticated training strategies may be required for large problems.
6 Conclusion
This work presents unified theoretical framework for analysis and construction of graph attention networks. The exact sequence property of graph derivatives and coercivity of Hodge Laplacians which follow from the theory allow the construction of four structure-pres preserrip brackets, which we use to evaluate the role of irreversibility in both data-driven physics simulators and graph analytics problems. contexts, the pure diffusion bracket performed most poorly, with mixed results between purely reversible and partially dissipative brackets.

---

Random splits
CORA
CiteSeer
PubMed
Coauthor CS
Computer
Photo
GAT
78.7 2.3 1.94 1.1 2.2
1.3
90.5 0.6
78.0 19.0
85.7 20.3 2.0
GDE
91.6 3.7
82.9
GRAND-nl
8.6
1.0
77.5 1.8
92.4
82.4
92.4
Hamiltonian
2.12
9.2 1.1 1.9
84.0 1.0
9.2
Gradient $1.3122
72.1 1.7
77.2 2.1
92.2
78.1 1.2
88.2 0.6
Double Bracket
1.1
71.2
+: 2.5
78.2 + 2.0
92.5 + 0.2
8488 0.5
92.1 0.3
Metriplectic 59.6 99.6:2.00 2.0
63.1 + 2.4
69.8 + 2.1
Table 5: Test accuracy and standard deviations averaged over 20 runs with random 80/10/10 train/val/test splits. Comparisons use the numbers reported in [49].
The linear scaling achieved by the metriplectic brackets has potential major impact for data-driven physics modeling. Metriplectic systems emerge naturally when coarse gratsiin multiscale systems. With increasing interest in using ML to construct digital twins, fast data-driven surrogates for complex mulli-physics acting over multiple scales will become crucial. In this setting the stability encoded by metriplectic dynamics translates to robust surrogates, with linear complexity suggesting the possibility of scaling up to millions of degrees of freedom.
Limitations: All analysis holds under the assumption of modified attention mechanisms which allow interpretation of GAT networks as diffusion processes; readers should take care that the analysis is for non-standard attention. Secondly, for all brackets we did not introduce empirical modifications (e.g. regularization, forcing, etc) to optimize performance so that we could study the role (ir)reversibility isolation. With this in mind, one may be able to add "tricks" to e.g. obtain a diffusion architecture which outperforms those presented here. Finally, note that the use of a feature autoencoder the bracket architectures means that structure is enforced in the transformed space. This allows for applicability to more general systems, and can be easily removed when appropriate features are known.
Broader impacts: The work performed here is strictly foundational mathematics and intended to improve the performance of GNNs in the context of graph analysis and data-driven physics modeling. Subsequent application of the theory may have societal impact, but the current work anticipated to improve the performance of machine learning in graph settings only at a foundational level.
References
[1] Marco Gori, Gabriele Monfardini, and Franco Scarselli. new model for learning graph domains. In Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005., volume 2, pages 729-734. IEEE, 2005.
[2] Jie Zhou, Ganqu Cui, Shengding Hu, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang. Changeheng Li, and Maosong Sun. Graph neural networks: A review of methods and applications. Al open, 11, 2020.
[3] Will Hamilton, Zhitao Ying, and Jure Leskovec. Inductive representation learning on large graphs. Advances in neural information processing systems, 30, 2017.
[4] William L Hamilton, Rex Ying. and Jure Leskovec. Representation learning graphs: Methods and applications. arXiv preprint (oXiv:7009005888 2017.
[5] Deli Chen, Yankai Lin, Wei Li, Peng Li, Jie Zhou, and Xu Sun. Mcasuring and relieving the over-smoothing problem for graph neural networks from the topological view. In Proceedings of the AAAI conference on artificial intelligence, volume 34, pages 3438-445,, 2020.
[6] Michael M Bronstein, Joan Bruna, Taco Cohen, and Petar Veličković. Geometric deep learning: Grids, groups, graphs, geodesics, and gauges. arXiv preprint arXiv:2 13478, 2021.
[7] Kaixiong Zhou, Xiao Huang, Yuening Li, Daochen Zha, Rui Chen, and Xia Hu. Towards deeper graph neural networks with differentiable group normalization. Advances in neural information processing systems, 33:4917-4928 2020.
10

---

[8] Chen Cai and Yusu Wang. A note over-smoothing for graph neural networks, arXiv preprint Ttti:2000.133116 2020.
[9] Yifei Wang. Yisen Wang, Jiansheng Yang, and Zhouchen Lin. Dissecting the diffusion process in linear graph convolutional networks. Advances in Neural Information Processing Systems, 34:5758- 5769, 2021.
[10] Francesco Giovanni, James Rowbottom, Benjamin Chamberlain, Thomas Markovich, and Michael M.
Bronstein. Understanding convolution on graphs via energies, 2023.
[11] Nathaniel Trask, Andy Huang, and Xiaozhe Hu. Enforeing exact physics scientific machine learning: a data-driven exterior calculus on graphs. arXiv preprint arxiv: 1799,, 2020.
[12] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio. Graph attention networks. arXiv preprint
2017.
[13] Khemraj Shukla, Mengjia Xu, Nathaniel Trask, and George E Karniadakis. Scalable algorithms for physics-informed neural and graph networks. Data-Centric Engineering, 3:e24, 2022.
[14] Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman. Universal differential equations for scientific machine learning. arXiv preprint Ttt::200110438533 2020.
[15] Steven Brunton, Joshua L Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences. 13((15)33932-39) 3937, 2016.
[16] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 6572-6583, 2018
[17] Michael Poli, Stefano Massaroli, Atsushi Yamashita, Hajime Asama, Jinkyoo Park, and Stefano Ermon. Torchdyn: Implicit models and neural numerical methods pytorch.
[18] Louis-Pascal Xhonneux. Meng Qu, and Jian Tang. Continuous graph neural networks. In International Conference on Machine Learning, pages 10432- 10441. PMLR, 2020.
[19] Fangda Gu, Heng Chang, Wenwu 7hu, Somaych Sojoudi, and Laurent EI Ghaoui. Implicit graph neural networks. Advances in Neural Information Processing Systems, 33:1 1194-11995, 2020.
[20] Samuel Greydamus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
[21] Mare Finzi, Ke Alexander Wang, and Andrew G Wilson. Simplifying hamiltonian and lagrangian neural networks via explicit constraints. Advances neural information processing systems, 33:13880-13889, 2020.
[22] Renyi Chen and Molei Tao. Data-driven prediction of general hamiltonian dynamics via learning exactly- symplectic maps. In International Conference on Machine Learning, pages 1717-1727. PMLR, 2021.
[23] Nate Gruver, Mare Anton Finzi, Samuel Don Stanton, and Andrew Gordon Wilson. Deconstructing the inductive biases of hamiltonian neural networks. In International Conference Learning Reprexentations.
[24] Peter Toth, Danilo Rezende, Andrew Jaegle, Sébastien Racanière, Aleksandar Botev, and Irina Higgins. Hamiltonian generative networks. In International Conference Learning Representations, 2019.
[25] Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty. Symplectic ODE-Net: Learning Hamiltonian dynamies with control. In International Conference on Learning Representations.
[26] Michael Lutter, Christian Ritter, and Jan Peters. Deep lagrangian networks: Using physics as model prior for deep learning. International Conference on Learning Representations, 2018
[27] Miles Cranmer, Sam Greydanus. Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. La- grangian neural networks. In ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations, 2020.
[28] Partha Guha. Metriplectic structure, leibniz dynamics and dissipative systems. Journal of Mathematical Analysis and Applications, 326(1):121-15)) 2007.

---

[29] Kookjin Lee, Nathaniel Trask, and Panos Stinis. Machine learning structure preserving brackets for forecasting inreversible processes. Advances in Neural Information Processing Systems, 34:5696-5707. 2021.
[30] Kookjin Lee, Nathaniel Trask, and Panos Stinis. Structure-preserving sparse identification nonlinear dynamies for data-driven modeling. In Mathematical and Scientific Machine Learning pages 65-80.
PMLR, 2022.
[31] Zhen Zhang, Yeonjong Shin, and George Em Karniadakis. Gfinns: Generie formalism informed neural networks for deterministic and stochastic dynamical systems. Philosophical Transactions of the Royal Society A, 3800222)) 2022.
[32] Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty. Dissipative symoden: Encoding hamiltonian dynamics with dissipation and control into deep learning. In ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations.
[33] Shaan A Desai, Marios Mattheakis, David Sondak, Pavlos Protopapas, and Stephen J Roberts. Port- hamiltonian neural networks for learning explicit time-dependent dynamical systems. Physical Review 1043))034312, 2021.
[34] Miroslav Grmela. Generic guide to the multiscale dynamics and thermodynamics. Journal of Physics Communications, 2(3):032001, 2018.
[35] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed 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.
[36] Maziar Raissi and George Em Karniadakis. Hidden physics models: Machine learning of nonlinear partial differential equations. Journal of Computational Physics, 357:125- 14115-1112221 2018.
[37] Filippo Masi and Ioannis Stefanou. Multiscale modeling of inelastic materials with thermodynanics--base based artificial neural networks (lann). Computer Methods in Applied Mechanics and Engineering, 398: 15190. 2022.
[38] Ravi G Patel, Indu Manickam, Nathaniel A Trask, Mitchell A Wood, Myoungkyu Lee, Ignacio Tomas, and Eric C Cyr. Thermodynamically consistent physics-informed neural networks for hyperbolic systems. Journal of Computational Physics, 449:110754 2022.
[39] Quercus Heránndzz, Alberto Badías, David González, Francisco Chinesta, and Elías Cueto. Structure- preserving neural networks. Journal of Computational Physics, 426: 109950, 2021.
[40] Yibo Yang and Paris Perdikaris. Adversarial uncertainty quantification in physics- informed neural networks. Journal of Computational Physics, 394:136-152, 2019.
[41] Dongkun Zhang. Lu, Ling Guo, and George Em Karniadakis. Quantifying total uncertainty physics- informed neural networks for solving forward and inverse stochastic problems, Journal of Computational Physics, 397:108850, 2019.
[42] Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow pathologies in mhysics-informer neural networks. SIAM Journal on Scientific Computing, 43(5):A3055-A4 A3081, 2021.
[43] Sifan Wang. Xinling Yu, and Paris Perdikaris. When and why pinns fail to train: A neural tangent kernel perspective. Journal of Computational Physics, 449:110768 2022.
[44] T Konstantin Rusch, Ben Chamberlain, James Rowbottom, Siddhartha Mishra, and Michael Bronstein. Graph-coupled oscillator networks. In International Conference on Machine Learning, pages 18888-18909. PMLR, 2022.
[45] Jeongwhan Choi, Scoyoung Hong, Noseong Park, and Sung-Bae Cho. Gread: Graph neural reaction- diffusion equations. arXiv preprint arXiv:2211. 14208, 2022.
[46] Alvaro GonzazzZ Victor Bapst, Kyle Cranmer, and Peter Battaglia. Hamiltonian graph networks with ode integrators. arXiv preprint arXiv:1909 12790, 2019.
[47] Suresh Bishnoi, Ravinder Bhattoo, Jayadeva Jayadeva, Sayan Ranu, and NM Anoop Krishnan. Enhancing the inductive biases of graph neural ode for modeling physical systems. In The Eleventh International Conference on Learning Representations.
12

---

[48] Quercus Hernandez. Alberto Badias, Francisco Chinesta, and Elias Cueto. Thermeolyammis-inforrreed graph neural networks. IEEE Trunsactions Artificial Intelligence, 1(01):1-1,200 2022.
[49] Ben Chamberlain, James Rowbottom, Maria I Gorinova, Michael Bronstein, Stefan Webb, and Emanuele Rossi. Grand: Graph neural diffusion. In International Conference on Machine Learning, pages 1407-1418. PMI.R, 2021.
[50] PJ Morrison. Thoughts on brackets and dissipation: old and new. In Journal of Physics: Conference Series, volume 169, page 012006. IOP Publishing, 2009.
[51] Oliver Knill. The dirac operator of graph, 2013.
[52] Xiaoye Jiang, Hen- Lim, Yuan Yao, and Yinyu Ye. Statistical ranking and combinatorial hodge theory. Mathematical Progranming, 27711:203-244- 2011.
[53] Pavel B Bochev and James M Hyman. Principles of mimetic discretizations of differential operators. In Compatible spatial discretizations, pages 199. Springer, 2006.
[54] Douglas N Amold, Finite element exterior calculus. SIAM, 2018.
[55] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, E. ukasz Kaiser, and Illia Polosukhin, Attention is all you need. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors,Avaancee in Neural Information Processing Systems, volume 30. Curran Associates, Inc.. 2017.
[56] Ignacio Romero.
Thernodynamically consistent time- step--ttpping algorithms for non non-lineer thermomechanical systems. International Journal for Numerical Methods in Engineering, 79(6):706-734) 2023/05/14 2009,
[57] Emanuel Todorov, Tom Erez, and Yuval Tassa, Mujoco: physics engine for model-based control. In 2012 IEEF/RSJ international conference on intelligent robots and systems, pages 5026-5033 IEEE, 2012.
[58] Andrew Kachites McCallum, Kamal Nigam, Jason Rennie, and Kristie Seymore. Automating the construc- tion of internet portals with machine learning. Information Retrieval, 3:127-163, 2000.
[59] Prithviraj Sen, Galileo Namata, Mustafa Bilgic, Lise Getoor, Brian Galligher, and Tina Eliassi-Rad. Collective classification in network data. Al magazine, 29(3):93-99) 93, 2008
[60] Galileo Namata, Ben London, Lise Getoor, Bert Huang, and U Edu. Query-driven active surveying for collective classification. In 10th international workshop on mining and learning with graphs, volume 8, page 1, 2012.
[61] Oleksandr Shchur, Maximilian Mumme, Aleksandar Bojchevski, and Stephan Günnemann. Pitfalls of graph neural network evaluation. ArXiv, abs//s11l.0000 1.05868, 2018
[62] Julian McAuley, Christopher Targett, Qinfeng Shi, and Anton Van Den Hengel. Image-based recommen- dations on styles and substitutes. In Proceedings of the 38th international ACM SIGIR conference on research and development in information retrieval, pages 43 52, 2015.
[63] Michael Poli, Stefano Massaroli, Junyoung Park, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park. Graph neural ordinary differential equations. arXiv preprint arXiv:911..05322 2019.
[64] Xiaoye Jiang, Lek- Lim, Yuan Yao, and Yinyu Ye. Statistical ranking and combinatorial hodge theory. Mathematical Programming, 127(1):203-22 nov 2010.
[65] Anthony Bloch, P. Krishnaprasad, Jerrlld E. Marsden, and Tudor Ratiu. The euler-poincaré equations and double bracket dissipation. Communications in Mathematical Physics, 175(1):1 42, 1996.
[66] Darryl Holm, Jerrold E Marsden, and Tudor Ratiu. The coier-poincaré equations and semidirect products with applications to continuum thcorics. Advances in Mathematics, 137(1):1 -81, 1998.
[67] Hans Christian Octtinger. Irreversible dynamics, onsager-casimir symmetry, and an application to turbu- lence. Physical Review F, 9044)10 2014.
[68] Anthony Gruber, Max Gunzburger, Lili Ju, and Zhu Wang. Energetically consistent model reduction for metriplectic systems. Computer Methods in Applied Mechanics and Engineering, 404:115709, 2023.
[69] Robert J. Renka. A simple explanation of the sobolev gradient method. 2006.
[70] Chris Yu, Henrik Schumacher, and Keenan Crane. Repulsive curves. ACM Trans, Graph., 40(2), may 2021.
13

---

[71] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chaman, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances neural information processing systems, 32, 2019.
[72] Joe Chen. Chaos from simplicity: introduction the double pendulum. Technical report, University of Canterbury, 2008.
[73] Diederik P Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint arXiv: 2.600, 2014.
[74] Romeo Ortega, Arjan Van Der Schaft, Bernhard Maschke, and Gerardo Escobar. Interconnection and damping assignment passivity-based control port-controlled hamiltonian systems. Automatica, 38(4):585- 596, 2002.
[75] Greg Brockman, Vicki Cheung. Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojciech Zaremba. OpenAl Gym. arXiv preprint arXiv:1606.01500, 2016.
[76] Lukas Biewald. Experiment tracking with weights and biases, 2020. Software available from wandb.com.
14

---

Glossary of Notation and Symbols
The next list describes several symbols that will be later used within the body of the document
(5+))
(Irreversible) double bracket on functions with generator L
[1] Degenerate (irreversible) metric bracket on functions with generator M M
(s-) Poisson (reversible) bracket on functions with generator L = L
(i) Neighbors of node E V, neighbors node E V including
[S]
Indicator function of the statement S
5f.
V
Adjoint of df with respect to (., 1) adjoint of df with respect to (s.))
Ax
Hodge Laplacian dæ dz
5i
Kronecker delta
Derivative of f with respect to time
(5))
Learnable metric inner product on k-cliques with matrix representation A
Euclidean @2 inner product on k-cliques
G, V
Oriented graph, set of nodes, set of edges
Gk. S2
Set of k-cliques, vector space of real-valued functions on k-cliques
d, d2
Exterior derivative operator on functions, exterior derivative operator on k-cliques
Adjoint of da with respect (A)) adjoint of de with respect to (:)
A
Mathematical foundations
This Appendix provides the following: (1) introduction the ideas of graph exterior calculus, A.1 and bracket-based dynamical systems, A.2 necessary for understanding the results in the body, (2) additional explanation regarding adjoints with respect generic inner products and associated computations, A.3 (3) a mechanism for higher-order attention expressed in terms of learnable inner products, A.4 (4) discussion of GATs in the context of exterior calculus, A.5 and (5) proofs which
are deferred from Section 4 A.6
A.1
Graph exterior calculus
Here some basic notions from the graph exterior calculus recalled. More details be found in, e.g., M 51 64
Figure 3: A toy graph with six O-cliques (nodes), six I-cliques (edges), and one 2-clique.
As mentioned in Section 3 oriented graph G = (V, E} carries sets of k-cliques, denoted Gk. which are collections of ordered subgraphs generated by (A: + 1) nodes. For example, the graph in Figure contains six O-cliques (nodes), six 1-cliques (edges), and one 2-clique. A notion of combinatorial derivative is then given by the signed incidence matrices dg Ω4 - Ω operating on the space Og of differentiable functions on k-cliques, whose entries (dy) are or -1 if the jth k-clique is
15

---

incident on the ith (k + 1) cllique, and zero otherwise. For the example in Figure these are:
1
0
0
1
1
0
1
0
1
do 0
1 = 0
0
0
-1
1
1 0
0
-1
0
0
-1
Remark A.1. While the one-hop neighborhood of node i in G. denoted (i), does not include node itself. many machine learning algorithms employ the extended neighborhood ( N(i) U {i}. Since this is equivalent to considering the one-hop neighborhood of node in the self-looped graph G. this modification does not change the analysis of functions graphs.
It can be shown that the action of these matrices be conveniently expressed in terms of totally antisymmetric functions f E Ωg. via the expression
A+1
Σ C-1 (i..... Ey.... (nn1), (x).,.........,.., 1=0
(dxf) (1o. 21, ix+1)
where (io...., 11)) denotes a (k I--liqque of vertices E V. As convenient shorthand, we often write subseripts, e.g.. (d&f (orl)iiiin..... i11 instead of explicit function arguments. Using [S] to denote the indicator function the statement S, it is straightforward to check that o = 0,
k+1
(dada- 1. Diom.nik
(-1)² (dx-1 -oonsspeee.........e.
==
k+1
<-](- S Ron.neenddynnae
k+1
11 51 Gy
1=0
11)
(-1jst) S San 0,
since labels L. This shows that the k-cliques on & form a de Rham complex [53]: collection of function
and the final follows from swapping the
spaces Qg equipped with mappings d₃ satisfying Im dx-1 Ker dç as shown in Figure When
d2
dx
Ω
ΩK
Figure 4: Illustration of the de Rham complex G induced the combinatorial derivatives, where K is the maximal clique degree.
K 3, this is precisely the graph calculus analogue of the de Rham complex on R2 formed by the Sobolev spaces H² (curl), (div), 1,2 which satisfies div curl curl grad
While the construction of the graph derivatives and their associated de Rham complex is purely topological, building elliptic differential operators such as the Laplacian relies dual de Rham complex, which is specified by an inner product on 2z. In the case of E2, this leads to dual derivatives which are the matrix transposes of the dy having the following explicit expression.
Proposition A.1. The dual derivatives df. Ω4+1 + given by
S2g adjoint to de through the inner product
k+1
-
1
(d]f) (io. i1, ix)
..... [G i++1)),
k+2
where (i,, ix+1] iet1, ig. "........ indicates cyclic permutation forward by one index.
16

---

Proof. This is a direct calculation using the representation of dy in terms of antisymmetric functions. More precisely, let an empty sum Σ denote summation over all unspecified indices. Then, for any g E 9k:
(dxf.9)
Σ
(dkffio... Pig
f E0x+1
R+1
1
1) G.odee Bou
( + 2)
k+1
1
Σ fio.ni
ΣΣ
11 91o...(65mm(i11))))
(k 2)
Pe-15=0
k+1
Σ ftande
ΣΣ 1E41
- 1) gig (....I1g..dettr))
k
Σ
Sig (119) = (1.dgg)
which establishes the result.
Proposition A. is perhaps best illustrated with concrete example. Consider the graph gradient, defined for edge a = (i.j) as (duf) = (dof).) = -j - fi. Notice that this object is antisymmetric with respect to edge orientation, and measures outflow of information from source to target nodes. From this, it is easy to compute the -addooint of do. known as the graph divergence, via
(dof.g) = Σ (fr - f1)gs Σ Σ
gisfs - gis
a=(i
(j>i)EN()
1 4 Si fi 9is) (f,
2
JEN(()
where we have re-indexed under the double sum, used that i E (j) if and only if E (i)), and used that there are no self-edges in E. Therefore, it follows that the graph divergence at node is given by
1 Σ 9jz 9bj.
(d5.);; Σ9-a -Ja
--a 03i
jEN(s)
which reduces to the common form (dhg); = E, Dij and only if the edge feature Pij is antisymmetric.
Remark A.2. When the inner product on edges & is not I.2, but defined in terms of nonnegative, intentationiivariinn and (edge-wise) diagonal weight matrix W = (Wig), similar computation shows that the divergence becomes
1 (G
Σ
Wiy (fzr fiz).
2
jEN(()
The more general case of arbitrary inner products on V. & is discussed in section A.3
The differential operators df. induce a dual de Rham complex since di_,dī = (dxdx 0, which well as
enables both the construction of Laplace operators on k-cliques, A d. + d - 1-1 Theorem 3.3]. the celebrated Hodge decomposition theorem, stated below. For a proof, see, e.g., [I]]
Theorem A.3. (Hodge Decomposition Theorem) The de Rham c the following direct sum decomposition of the function space Q2.
by de, dj. induce
Ωr = Im d e Ker A @ Im di
In the case where the dual derivatives dž are adjoint with respect to a learnable inner product which does not depend on graph features, the conclusion of Theorem A.3 continues to hold, leading to an interesting rell-posedness result proved in [T] involving nonlinear perturbations of Hodge-Laplace problem in mixed form.
17

---

Theorem A.4. ([I Theorem 3.6]) Suppose f E s2g, and g (x: E) is a neural network with parame- ters & which is Lipschitz continuous and satisfies g(0) 0. Then, the problem
Wk k = d u + eg 1 u):
f. = 1 Wk + dg dg uk:
has unique solution on 2 Ax.
This result shows that initial-value problems involving the Hodge-Laplace are stable under nonlinear perturbations. Moreover, when Ag is the Hodge Laplacian on nodes, there is useful connection between Ao and the degree and adjacency matrices of the graph G. Recall that the degree matrix
D (dey) is diagonal with entries du = jeevio 1, while the adjacency matrix = (ais) satisfies
a. when E (i)) and a:j 0 otherwise. Proposition A.2. The combinatorial Laplacian on V, denoted Do = dodo. satisfies Do D A.
Proof. Notice that
(dīdo))e = Σ (do) (do) = =1==1] Σ +b OEE
1i Σ [i 5] dis - * 11 alj dis aij =D
(do) (do)
aEE
a=((i))
where we used that D diagonal, is diagonal-fres free, and (do).: (do) whenever a = (1, is an edge in E, since one of (do).. (do).. is I and the other is -1.
A.2 Bracket-based dynamical systems
Here we mention some additional facts regarding bracket-based dynamical systems. More information
can be found in, e.g., (50 65 66 677.
As mentioned before, the goal of bracket formalisms is to extend the Hamiltonian formalism to systems with dissipation. To understand where this originates, consider an action functional A(q) =
Sa L (q. 9) dt on the space of curves (t),, defined in terms of Lagrangian L on the tangent bundle to some Riemannian manifold. Using Lg. Lg to denote partial derivatives with respect to the subscripted variable, it is straightforward to show that, for any compactly supported variation åg of q. we have
dA(g))g - 1." dL (g.
[ LL,0q + LGōq = [((L, (Lg- 0,L) 64.
-1."
where the final equality follows from integration-by-pan derivatives commute in this setting. It follows that A is stationary (i.e., dA = for all variations only when Lg = Lg- These are the classical Buler-Lagrange equations which are (under some regularity conditions) transformed Hamiltonian form via a Legendre transformation,
and the fact that variational and temporal
H (q. p) sup ((p. L(g.
which defines the Hamiltonian functional H. on phase space, and yields the conjugate momentum vector p = Lg. Substituting (p. 4) H into the previously derived Euler-Lagrange equations leads immediately to Hamilton's equations for the state x (q T)_
-
!)
Hg Hp
x
JVH,
which are an equivalent description of the system in question in terms of the anti-involution J and the functional gradient V..
An advantage of the Hamiltonian description is its compact bracket-based formulation, x = JVH = {x, H). which requires only the specification of antisymmetric Poisson bracket {--} and Hamiltonian functional H. Besides admitting a direct generalization to more complex systems such as Korteweg-de Vries or incompressible Euler, where the involved bracket -taeedeepnn,nnt, this formulation makes the energy conservation property of the system obvious. In particular, it follows immediately from the antisymmetry of (--)} that
H (x, VH) (H, H} 0,
18

---

while it is more difficult to see immediately that the Euler-Lagrange system obeys this same property. The utility and ease-of-use of bracket formulations what inspired their extension to other systems of interest which do not conserve energy. On the opposite end of this spectrum are the generalized gradient flows, which can be written in terms of bracket which is purely dissipative. An example of this is heat flow q = Aq = [g. D), which is the I,2 -.A-gradient flow of Dirichlet energy D(q)= (1/2) Sa 1a'1 di (c.f. Appendix A.3). In this case, the functional gradient VD = - the negative of the usual Laplace operator, so that the positive-definite bracket [:: is generated by the identity operator M = id. It is interesting to note that the same system could be expressed using the usual
kinetic energy E (1/2) fa |g dl instead, provided that the corresponding bracket is generated by M = A. This is a good illustration of the flcxibility afforded by bracket-based dynamical systems.
Since physical systems are not always purely reversible or irreversible, other useful bracket formalisms have been introduced to capture dynamics which are mix of these two. The double bracket x = (x, E + {xx, E}] = LVE + L' is nice extension of the Hamiltonian bracket particularly because it is Casimir preserving, i.e., those quantities which annihilate the Poisson bracket {-, also annihilate the double bracket. This allows for the incorporation of dissipative phenomena into idealized Hamiltonian systems without affecting desirable properties such mass conservation, and has been used to model, e.g., the Landau-I .ifschitz dissipative mechanism, as well as a mechanism for fluids where energy decays but entrophy is preserved (see [65] for additional discussion). A complementary but alternative point of view is taken by the metriplectic bracket formalism, which requires that any dissipation generated by the system is accounted for within the system itself through the generation of entropy. In the metriplectic formalism, the equations of motion are fx, E) ( S] = LVE MVS, along with important and nontrivial compatibility conditions LVS = MVE = O, also called degeneracy conditions, which ensure that the reversible and irreversible mechanisms do not consacontame As shown in the body of the paper, this guarantees that metriplectic systems obey form of the first and second thermodynamical laws. Practically, the degeneracy conditions enforce a good deal of structure on the operators L, MM which has been exploited to generate surrogate models [29] 68 31 In particular, it can be shown that the reversible and irreversible brackets can be parameterized in terms of a totally antisymmetric order-3 tensor & (Eijk) and partially symmetric order-4 tensor Ç = (S through the relations (Einstein summation assumed)
(A, B} = Eiيk 0,AO,
[A, B] = c skزل 0, 0,A0
Moreover, using the symmetries of Ç. it follows (see 67 that this tensor decomposes into the product Sik.jt A DmnAAj Ai of a symmetric matrix and an order-3 3 tensor A which is skew-symmetric its lower indices. Thus, by applying symmetry relationships, it is easy to check that (} = F =
Remark A.5. In (29), trainable 4- and 3- tensors glik and Sik are constructed to achieve the degeneracy conditions, mandating costly O(N3) computational complexity. In the current work overcome this by instead achieving degeneracy through the exact sequence property.
A.3 Adjoints and gradients
Beyond the basic caleulus operations discussed in section A. which depend only on graph topology, the network architectures discussed in the body also make extensive use of learnable metric infor- mation coming from the nodal features. To understand this, it is useful to recall some information about general inner products and the derivative operators that they induce. First, recall that the usual e² inner product on node features a, b E R/M, (a, b) aTb, is (in this context) a discretization of the standard T.2 inner product fy ab du which aggregates information from across the vertex set V. While this construction is clearly dependent only on the graph structure (i.e., topology), any symmetric positive definite (SPD) matrix An
24 > Ωg also defines an inner product on functions
E Go through the equality
(a. b), = (a, Abb) = - A,b.
which gives a different way of measuring the distance between a and b. The advantage of this construction is that An can be chosen in way that incorporates geometric information which implicitly regularizes systems obcying variational principle. This follows from the following
19

---

intuitive fact: the Taylor series of a function does not change, regardless of the inner product on its domain. For any differentiable function(al) E Sg + R, using d denote the exterior derivative,
this means that the following equality holds
dE(a)b lim E(a cb) E(a) (óE(a), b) (VE(a),
c-
E
where SE denotes the e -graddeent of E and VE denotes Ao-gradient, i.e., its gradient with respect to the derivative operator induced by the inner product involving Ao. From this, it is clear that & = VE,, so that the Ao-gradient is just an anisotropic rescaling of the l version. The advantage of working with V over å in the present case of graph networks is that Ao can be learned based on the features of the graph. This means that learnable feature information (i.e., graph attention) can be directly incorporated into the differential operators governing our bracket-basse dynamical systems by construction.
The prototypical example of where this technique is useful is seen in the gradient flow of Dirichlet energy. Recall that the Dirichlet energy of a differentiable function D(u)= = (1/2) f /uul" dye, where now denotes the usual L --gradient of the function u on IR".. Using suvegratiionbyy-artt it is easy to see that dD(u = fffii vAu for any test function e with compact support, implying that the L'rggraieet of Dirichlet energy: the motion which decreases the quantity D(u) the fastest as measured by the L norm. It can be shown that high high-frequency modes decay quickly under this flow, while low-frequency information takes much longer to dissipate. On the other hand, we could alternatively run the H' -rradddientt flow of DD which motion of fastest decrease with respect to the H inner product (u, - (Vu, Vr) dj. This motion is prescribed in terms of the H -gradiet of D, which by the discussion above with Ao - A is easily seen to be the identity. This means that the H -gradient flow is given by ú = -, which retains the minimizers of the flo but with quite different intermediate character, since it functions by simultaneously flattening all spatial frequencies. The process of preconditioning a gradient flow by matching derivatives is known as a Sobolev gradient method (c.f. (69). and these methods often exhibit faster convergence and better numerical behavior than their 7,2 counterparts 77).
R + R is given by
A and ū Au is the Li-gradient flow of
Returning to the graph setting, our learnable matrices A& on k-cliques will lead to inner products (, on functions in S2g, and this will induce dual derivatives described in Appendix A.I However, in this case we will not have dj = dg. but instead the expression given by the following result: Proposition A.3. The Ag-adjoints di to the graph derivative operators dę given by di Ao di Ak11- Similarly. for any linear operator B 2g - S4 the Ax-adjoint B A, A A.
Proof. Let q. p denote vectors of k-clique resp. (k + I)-clique features. follows that (d=q. = (d. A = (q. 1p)) (q. A,d (q. (g.))) Therefore, we see that A A and hence di A5'd[AA Similarly, if q, q' denote vectors of k-clique features, it follows from the sesefffadjjiinte of A that
(q.Ba') j- = AB. = (BT A.q.q') = (A (A AT Ard') = (B'q.q)).
establishing that B A A Aç.
Remark A.6. It is common in graph theory to encounter the case where a, > are nodal weights and Wij > 0 are edge weights. These nothing more than the (diagonal) inner products As, Ai in disguise, and so Proposition A.3 immediately yields the familiar formula for the induced divergence
(dhp), =
Σ J(i-))EE
Wij (p,, -P;))
Note that all of these notions extend to the case of block inner products in the obvious way. For example, if q, p are node resp. edge features, it follows that A diag ( (Ap, A1) ) is an inner product on node-edge feature pairs, and the adjoints of node-edge operators with respect to A are computed as according to Proposition A.3
Remark A.7. For convenience, this work restricts to diagonal matrices Ag, A1. However, note that mairix which is diagonal in "edge space" G2 is generally full in nodal representation. This is because an (undirected) edge is uniquely specified by the two nodes which it connects, meaning that purely local quantity on edges is necessarily nonlocal nodes.
20

---

A.4 Higher order attention
As mentioned in the body, when f exp and ((),9)) (1/d) (Waq, Woq),, defining the learablle inner products Au = (0u,si) A1 (01.1)) as
Σ f(ā (q:
0 f (q:9)))
a0,i
JEN(1)
recovers scaled dot product attention as A A1.
Remark A.8. Technically, A1 an inner product only with respect to a predefined ordering of the edges (i. since do not require A₁ be orientation-invarinnt. On the other hand, is both unnecessary and distracting to enforce symmetry on A₁ in this context, since any necessary symmetrizzaioo will be handled automatically by the differential operator dz-
Similarly, other common attention mechanisms are produced by modifying the pre-attention function à. While A, 'A never appears the brackets of Section letting a (i.j) denote a global edge with endpoints j,j, it is straightforward to calculate the divergence of an antisymmetric edge feature p node
(hp),; = (Ar (.5'd[fAP) A1)) = 00,5 Σ (dh) (A1P)a
Σ (AiP) - (AiP), =
Σ
01.fi +
As.if
Pij.
00,5
a31
JEV(()
This shows that b(q, 9,) (a1.3) + a1.fs) /u0.:: appears under the divergence di A5' A1 which is the usual graph attention up to symmetrization in A1.
Remark A.9. While A1 is diagonal on global edges a = (i,j), it appears sparse nondiagonal in its nodal representation. Similarly, any diagonal extension A2 2-cliques will appear as sparse 3-tensor A2 = (02,ijk) when specified by its nodes.
This inspires straightforward extension of graph attention to higher-ord order cliques. In particular, de- note by K the highest degree of clique under consideration, and define A = (O -xf11.12....kk by
0K 0K = f (W (qi qi 9ix)). ))
where W R is learnable K K-tensor. Then, for any 0 l S K - define Aş (Ar.iria...en)) by
Σ SK KK...iK-k--
OK UK-1-19....
This recovers the matrices Ag, A₁ from before when K = 2, and otherwise extends the same core idea to higher-order cliques. It's attractive that the attention mechanism captured by di remains asymmetric, meaning that the attention of any one node to the others a k-clique need not equal the attention of the others that particular node.
Remark A.10. A more obvious but less expressive option for higher-order attention is to let
0.K i.1g-.i
Eix aK 1.4
for any 0 k However, application of the combinatorial codifferential dT_,, appearing in di 1 will necessarily symmetrize this quantity, that the asymmetry behind the attention mechanism lost in this formulation.
To illustrate how this works more concretely, consider the extension K = 3 to 2-cliques, and let (i..j) nN(j). We have the tensors Aq (a2,izk), A (a1,i)). and Ao (ao.i) defined by
02,igk f (W (qi, qi qx)).
01.1
02.ig
Σ
02.ikk.
kEA(())
JEN(()
REA
This provides way for (features on) 3-node subgraphs of G attend to each other, and can be similarly built-in to the differential operator di A d Ag.
21

---

A.5 Exterior calculus interpretation of GATs
Let (i() denote the one-hop neighborhood of node i, and let let (i) uit). Recall the standard (single-headed) graph attention network (GAT) described in [12). described layer-wise as
(
4 k+1 -
Σ a (ql, q²) W' (EX
where a is an element-wise nonlinearity, w is layer-dependent embedding matrix, and a (qis qj) denotes the attention node i pays to node j. Traditionally, the attention mechanism is computed through
aa(g).
(q, 9)) Softmax, ā (q. 4))
di
where the pre-attention coefficients (qe, 4j) and nodal weights or defined as
(q; q,) LeakyReT.U (aT (WTq: WTq)))
0= Σ (وq..)))
JEN()
However, the exponentials the outer Softmax are often replaced with other nonlinear functions, c.g. Squareplus, and the pre-attention coefficients à appear as variable (but learnable) functions of the nodal features. First, notice that (I) the attention coefficients (qi, qj) depend on the node features and not simply the topology of the graph, and (2) the attention coefficients are not symmetric, reflecting the fact that the attention paid by node node need not equal the attention paid by node to node i. A direct consequence of this is that GAT's are not purely diffusive under any circumstances, since it was shown in Appendix A.1 that the combinatorial divergence df will antisymmetrize the edge features it acts on. In particular, it is clear that the product (q: qi) (q: qi) is asymmetric in i under the standard attention mechanism, since even the pre-attention coefficients ā (q: qi) are not symmetric, meaning that there will be two distinct terms after application of the divergence. More precisely, there the following subtle result.
Proposition A.4. Let q E xnnv denote an array of nodal features. The expression
Σ a (q., q,) (q:-)()) -q),
JEN(
where A Ai is not the action of a Laplace operator whenever A is not symmetric.
Proof. From Appendix A.3 we know that any Laplace operator on nodes is expressible as di A df A2 for some positive definite Ag, A1. So, we compute the action of the Laplacian at node
i
(Aog), = (didoa), = (Anid AgAA = -0 Σ (df) (A)d)))
a
Σ (A1doq) (A1dog)a 1
anji + aisi (qi (q)-));
0
2
a3i
jEA(()
= Σ a (qi qi) (qi - -i),
JEN())
which shows that a (q, 9;) = (1/2) (01,ft + a1.15) /a0.5 must have symmetric numerator.
While this result shows that GATs (and their derivatives, e.g., GRAND) not purely diffusive, it also shows that it is possible to get close to GAT (at least syntactically) with a learnable diffusion mechanism. In fact, setting
W I in (I) yields precisely single-step diffusion equation provided that a (gf of) is right-stochastic (i.e., Σ, (q. q:) = 1;) and built as dictated by Proposition A.4
Theorem A.11. The GAT layer ① a single-step diffusion equation provided that W I. and the attention mechanism a (q. q,) (1/2) (a1 + a1,4)) / /ao,si is right-stochastic.
22

---

Proof. First, notice that the Laplacian with respect to an edge set which contains self-loops is computable via
(Aoq):= =- E a (qi, qi) (qy qi) = qi - Σ a (q. q,) qi
jeN())
jEN(()
Therefore, taking single step of heat flow q -oq with forward Euler discretization and time step r is equivalent to
q k+1 q. - (Aoq")
Σ a (q. 9.)
jEN(()
which is just a modified and non-activated GAT layer with W I and attention mechanism a.
Remark A.12. Since Softmax and its variants right-stochastic, Theorem A.1l is what establishes equivalence between the non-divee con-divergence equation
di =
Σ a (q: (4) -g), (EN
and the standard GAT layer seen in, e.g., 149), when a(q; qi) is the usual attention mechanism.
Remark A.13. In the literature, there an important (and often overlooked) distinction between the positive graph/Hodge Laplacian Ag and the negalive "geometer's 's Laplacian" A which is worth noting here. Particularly, we have from integration perratiio----aaee pyr-parts that the gradient V = dg is L2 -ajoiint to minus the divergence - dh. so that the two Laplace operators Do dīdo and A V V differ by a sign. This is why the same l2 -raaddiennt flow of Dirichlet energy can be equivalently expressed as q Aq -Anq but not by, e.g.. q Dnq.
This shows that, while they are not equivalent, there is close relationship between attention and diffusion mechanisms on graphs. The closest analogue to the standard attention expressible this format is perhaps the choice a = f (ā (qi, 4;)), ao,i E) a1.ij, discussed in Section and Appendix A.4 where f is any scalar-valued positive function. For example, when f(x) e.. .1 follows that
(Ang),
1
Σ jEN(()
((qq.a)) + (()))
(q: -q)
di
1
Σ JEN(()
a aq:
cala ad)
(q; -q),
2
which leads the standard GAT propagation mechanism plus an the attention a is not symmetric.
arising from the fact that
Remark A.14. Practically, GATs and their variants rypically make use of multi-head attention, defined in terms of an attention mechanism which is averaged over some number [h| of independent "heads",
a (q:19))
Σ ad (q:,9,)
T
which are distinct only in their learnable parameters. While the results of this section were presented in terms of hh] = 1, the reader can check that multiple attention heads can be used in this framework provided is the pre-attention à that is averaged instead.
A.6 Bracket derivations and properties
Here the architectures in the body are derived greater detail. First, will be shown that L* -L, G G. and M M, as required for structure preeturr-peeeevattion
Proposition A.5. For L, G, M defined in Section f have = -L, G* G. and M M.
23

---

Proof. First, denoting A diag (Ag. A1), it was shown in section A.3 that B" = AB' for
any linear operator B of appropriate dimensions. So, applying this to L, it follows that
A
1
0 dg
db
Ao
0
L
0
A
1
A1
A
dg
A1
0
di
-L.
A
(g))
Ao
0
dp
Similarly, it follows that
A -1
0
do
0
Ao
0
G
0
0
Ai
did,)
A,
A 0
dg
(d))
T
Ao
0
A )
doda
0
= G.
A df (d))
did1
A
0
0
Ao
M
AT
0
A
0
A1
M,
0
o (d)) A?
di
where the second-to-las equality used that AA II
Remark A.15. Note that the choice of zero blocks in L, G is sufficient but not necessary for these adjoinmess relationships to hold. For example, one could alternatively choose the diagonal blocks of L to contain terms like B B for appropriate message passiigg network B.
Next, we compute the gradients of energy and entropy with respect to (:,-). Proposition A.6. The A- A-gradient of the energy
E(q.p) 2 (lal + +lpl²)
1
Σ
2
Σ
Ppal
2
2
LEV
AEE
satisfies
VE(q.p)
A
9
A 9 A
0
AT
p
Moreover. given the energy and entropy defined as
Elq. p) = fE (s(q)) + 98 (s (dodfp))
S p) = gs (s (d]dip)).
where fE : RV + R acts on node features, 9E: 9s R" + R act on features, and s denotes sum aggregation over nodes or edges, the A A-gradients are
A A & fE (s(q)) VE(q.P) = (A (00;
dodj 1 & Vgz (s (dod
V((q.))) (A d @ Vgs (s (d d (d)))))
Proof. Since the theory of A-gradients in section A.3 establishes that VE = A SE, is only (q p) it follows for the first definition
necessary to compute the L'-grrdients. First, letting of energy that
dE(x)= = Σ (q, dq:) + Σ (p.. = (q. dq) + (p. dp) = (x, dx),
iEV
aEE
showing that SE p) (q T)_ as desired. Moving to the metriplectic definitions, since each term of E, has the same functional form, suffices compute the gradient of f (s (Bq)) for some function f m
RI - RIV, To that end, adopting the Einstein summation convention where repeated indices appearing up-and-down in an expression are implicitly summed, if S a, b < ny and S i,j I [V], we have
and matrix
(s(())) = Σ -qp = Σ
ddq = 1d = (1 dq = -V (s(g)) : dq.
KE
$E V
24

---

implying that (s(g) = 1 @ I. Continuing, it follows that
d( oss Bq) f' (s (Bq))a s' (Bq); Bd dgj (s (Bq)), ea ( +VV (Bq)) & BT1, dq) (V (f o bq),da)
1 dgf
showing that V (fSon o B) decomposes into an outer product across modalities. Applying this formula to the each term of E, S then yields the 12=gradients,
o dudh x Vge (s (dodfp)) ) 8S(g. (1 (
1 & fE (s(q))
EE p)
from which the desired A A-gradients follow directly.
Finally, can show that the degeneracy conditions for metriplectic structure are satisfied by the network in Section
Theorem A.16. The degeneracy conditions LVS MVE = 0 are satisfied by the metriplectic A.6
bracket in Section
Proof. This is a direct calculation using Theorem 3.1 and Proposition A.6
In particular, it follows
that
-
Aj Ad'dddis (A Vgs (s (dfdip))
LVS
do
(-A AJ (dido) d11 & Vgs (8 (d ((Hd))
0
0 0 AididiAi)
A 8 VJz (s(g))
MVE
A Aj
dodj
@ Vge (s (dodfp))
A A (di dif @ Vg (dodfp))
since dj do 0 as consequence of the graph calculus. These calculations establish the validity of the energy conservation and entropy generation properties seen previously in the manuscript body. O
Remark A.17. Clearly, this is not the only possible metriplectic formulation for GNNs. On the other hand, this choice is in some sense maximally general with respect to the chosen operators L, G. since only constants are the kernel of do (hence there is no reason to include a nodal term S), and only elements in the image of df (which do not exist in our setting) are guaranteed to be in the kernel of dī for any graph. Therefore, M is chosen to be essentially without the Do term, whose kernel is praph-dependent and hence difficult to design.
B Experimental details and more results
This Appendix provides details regarding the experiments in Section 5 as well any additional information necessary for reproducing them. We implement the proposed algorithms with PYTHON and PYTORCH [7I] that supports CUDA. The experiments are conducted on systems that are equipped with NVIDIA RTX A100 and V 100 GPUs. For NODEs capabilities, we use the TORCHDIFFEQ library 16
B.1 Damped double pendulum
The governing equations for the damped double pendulum can written in terms of four coupled first-order ODEs for the angles that the two pendula make with the vertical axis 02, 0g and their associated angular momenta W1, W2 (see (22).
0, = Wii i
(2)
W1 m sin (2A0) + 2mzlzw3 sin (AA) + 2gm2 cos 0g sin AO 2gm; sin 01 + 7/1
(3)
21, (mi + m2 sin? A0) 42 malaw? sin (2A0) + (m m2) lww sin A0 2g (mi + m2) cos 0 sin A0 + 72
(4)
21z (mi + m2 sin? A0)
25

---

where mi, m2, 11,12 l2 are the masses resp. lengths of the pendula, A 81 82 is the (signed) difference in vertical angle, g is the acceleration due to gravity, and
71 - 2 01 2k202 co8 A0.
Y/2 = 2k 01 cos Ao (m) + m2) k202,
m2
for damping constants k k2-
Dataset. A trajectory of the damped double pendulum by solving an initial value problem associated
with the ODE The initial condition used is (1.0, m/2, 0.0, 0.0), and the parameters m7 m2
1 = 1,4 = T, 12 = 0.9, k = l2 = 0.1. For time integrator, we use the TorchDiffeq library [16] 5 (DOPRI5) the numerical solver. The total simulation time is 50 (long
with Dormand- enough for significant dissipation to occur), and solution snapshots are collected at 500 evenly-spaced temporal indices.
To simulate the practical case where only positional data for the system is available, the double pendulum solution is integrated to time T = 50 (long enough for significant dissipation to occur) and proccesod once the angles and angular momenta are determined from the equations above, yielding the (Ir coordinaass of each pendulum mass at intervals of 0.1s. This is accomplished using the relationships
T = 1 sin 0
y1 - 11 cos 01
T2 = xi + l2 sin 02 1, sin 01 + 12 sin 02
1/2 1/1 12 cos 02 -1 cos 0 12 cos 02.
The double pendulum is then treated as a fully connected three-node graph with positional coordinates q, = (Ii y4) as nodal features, and relative velocities Pa = (doq)a as edge features. Note that the positional coordinates (zo, yo) = (0, 0) of the anchor node are held constant during training. To allow for the necessary flexibility of coordinate changes, each architecture from Section 4 makes use of a message-passee passing feature encoder before time integration, acting on node features and edge features separately, with corresponding decoder returning the original features after time integration.
To elicit a fair comparison, both the NODE and NODE+AE architectures are chosen to contain comparable numbers of parameters the bracket architectures (~~ 30k), and all networks are trained for 100,000 epochs. For each network, the configuration of weights producing the lowest overall error during training is used for prediction.
Hyperparameters. The networks are trained to reconstruct the node/edge features in mean absolute error (MAE) using the Adam optimizer 173) The NODEs and metriplectic bracket use an initial learning rate of 10-4, while the other models use an initial learning rate of 10-3, The width of the hidden layers in the message passing encoder/decoder 64, and the number of hidden features for nodes/edges is 32. The time integrator used is simple forward Euler.
Network architectures. The message passing encoders/decoders are 3-layer MLPs mapping, in the node case, nodal features and their graph derivatives, and in the edge case, edge features and their graph coderivatives, mechanism used in the learnable coderivatives is scaled dot product. The metriplectic network uses 2-layer MLPs JF, gF; Is with scalar output and hidden width 64. For the basic NODE, node and edge features are concatenated, flattened, and passed through 4-layer fully connected network of width 128 in each hidden layer, before being reshaped at the end. The NODE+AE architecture uses a 3-layer fully connected network which operates on the concatenated and flattened latent embedding of size 32 6 192, with constant width throughout all layers.
hidden representation. For the bracket architectures, the attention
B.2 Mujoco
We represent an object as a fully-connected graph, where node corresponds to a body part of the object and, thus, the nodal feature corresponds position of body part or joint. To learn the dynamics of object, we again follow the encoder--decoder-type architecture considered in the
26

---

(a) q
(b)
(c) E and
Figure 5: [Double pendulum] Trajectories of q and p: ground-truth (solid lines) and predictions of the metriplectic bracket model (dashed lines). The evolution of the energy E and the entropy S over the simulation time. Note that slight fluctuations appear in due the fact that forward Euler not a symplectic integrator.
double-pentn ceabbll-enuulm experiments. First we employ node-wise linear layer to embed the nodal feature into node-wise hidden representations (i.e., the nodal feature Qi corresponds to position of body part an angle of a joint.). As an alternative encoding scheme for the nodal feature, in addition to the position or the angle, nodal velocities are considered additional nodal features, i.e., qi = (gi, U.). The experimental results of the alternative scheme is represented in the following section (ection[.2..2
The proposed dynamies models also require edge features (e.g., edge velocity), which are not presented in the dataset. Thus, to extract a hidden representation for an edge, we employ a linear layer, which takes velocities of the source and destination nodes of the edge as an input and outputs edge-wise hidden representations, i.e., the edge feature correspond to pair of nodal velocities Pa (P)re(a) F) where (n) and Vasa(a)) denote velocities of the source and destination nodes connected to the edge.
The MuloCo trajectories
generated in the presence of an actor applying controls. To handle the changes in dynamics due to the control input, we introduce an additive forcing term, parameterized by MLP, to the dynamics models, which a similar approach considered in dissipative SymODEN [32] In dissipative SymODEN, the forcing term is designed to affect only the change of the generalized momenta (also known the port-Hamiltonian dynamics (74)). As opposed this approach, our proposed forcing term affects the evolution of both the generalized coordinates that are defined in the latent space. Once the latent states are computed at specified time indices, a node-wise linear decoder is applied to reconstruct the position of body parts of the object. Then the models are trained based on the data matching loss measured in mean suanrquuar errors between the reconstructed and the ground-truth positions.
B.2.1 Experiment details
We largely follow the experimental settings considered in 23
Dataset. As elaborated in 23 the standard Open AI Gym 75 environments preprocess observa- tions in ad-hoe ways, e.g.. Hopper clips the velocity observations [-10, 10]d. Thus, the authors in (23] modified the environments to simply return the position and the velocity (g. as the observation and we use the same dataset, which is made publicly available by the authors. The dataset consists of training and test data, which are constructed by randomly splitting the episodes in the replay buffer into training and test data. Training and test data consist of ~40K and ~300 or ~ 85 trajectories, respectively. For both training and test data, include the first 20 measurements (i.c., 19 transitions) in each trajectory.
Hyperparameters. For training, we use the Adam optimizer 73 with the initial learning rate 5e-3 and weight decay 1e-4. With the batch size of 200 trajectories, we train the models for 256 epochs. We also employ a cosine annealing lcarning rate scheduler with the minimum learning rate 1e-6. For time integrator, we use the Torchdiffeq library with the Euler method.
Network architectures. The encoder and decoder networks are parameterized as linear layer and the dimension of the hidden representations set to 80. For attention, the scaled dot-product attention is used with heads and the embedding dimension is set to 16. The MLP for handling the
27

---

forcing tern consists of three fully-connected layers (i.e., input, output layers and one hidden layer with 128 neurons). The MLP used for parameterizing the "black-box" NODEs also consists of three fully-connected layers with 128 neurons in each layer.
B.2.2 Additional results
Figure reportt the loss trajectories for all considered dynamies models. For the given number of maximum epochs (i.e., 256), the Hamiltonian and double bracket models tend to reach much lower training losses (an order of magnitude smaller) errors than the NODE and Gradient models do. The metriplectic model produces smaller training losses compared to the NODE and gradient models after a certain number of epochs (e.g., 100 epochs for Hopper).
Hamiltonian
Gradient
Double
Metripletic
101
trai
13a1
10-2
spoch
(a) HalfCheetah
(b) Hopper
(c) Swimmer
Figure [Mujoco] Train MSE over epoch for all considered dynamics models. For the nodal feature, only the position or the angle of the body part/joint is considered.
In the set of experiments, we provide only positions/angles of body parts nodal features, but also velocities of the body parts as nodal features (i.e., q = (gs= v))). Table 6 reports the relative errors measured in L2-norm; again, the Hamiltonian, double bracket, and metriplectic outperform other dynamics models. In particular, the metriplectic bracket produces the most accurate predictions in the Hopper and Swimmer environments. Figure repoorts the loss trajectories for all considered models. Similar to the previous experiments with the position the Hamiltonian and Double bracket produces the lower training losses than the NODE and Gradient models do. For the Hopper and Swimmer environments, however, among all considered models, the metriplectic model produces the lowest training MSEs after 256 training epochs.
the only nodal feature,
Dataset
HalfCheetah
Hopper
Swimmer
NODE+AE
0.0848 + 0.0011
0.0421 + 0.0041
0.0135 +:2 0.00082
Hamiltonian
0.0403 0.0052
0.0294 + 0.0028
0.0120 0.00022
Gradient
0.0846 + 0.00358
0.0490 + 0.0013
0.0158 + 0,00030
Double Bracket
0.0653 = 0.010
0.0274 + 0.00090
0.0120 + 0.00060
Metriplectic
0.0757 + 0.0021
0.0269 + 0.00035
0.0114 + 0.00067
Table 6: Relative enrors of the network predictions of the MuJoCo environments on the test set,
reported as stdev over runs.
B.3
Node classification
To facilitate comparison with previous work, we follow the experimental methodology of 61
28

---

NODE
Hamiltonian
Gradient
Double
Metripletic
101
0331 V62
107
10m1
10"
opoch
aooch
(a) HalfCheetah
(b) Hopper
(c) Swimmer
Figure 7: [Mujoco] Train MSE over epoch for all considered dynamics models. For the nodal feature, along with the position or the angle of the body part/joint, the node velocities are also considered.
B.3.1 Experiment details
Datasets. We consider the three well-known citation networks, Cora 588 Citeseer 55 and
Pubmed [60|; the proposed models are tested on the datasets with the original fixed Planetoid traing/test splits, as well as random train/test splits. In addition, we also consider the coauthor graph, CoauthorCSS [6]] and the Amazon co-purchasing graphs, Computer and Photo 62 Table provides some basic statistics about each dataset.
Dataset
Cora
Citeseer
PubMed
CoauthorCS
Computer
Photo
Classes
7
3
10
Features
1433
3703
500
6805
767
745
Nodes
2485
2120
19717
18333
13381
7487
Edges
5069
3679
44324
81894
245778
119043
Table 7: Dataset statistics.
Hyperparameters The bracket architectures employed for this task are identical to those in Section 4 except for that the right-hand side of the Hamiltonian, gradient, and double bracket networks is scaled by learnable parameter Sigmoid (a) > O, and the matrix A, = I used as the inner product on 3-cliques. It is easy to verify that this does not affect the structural properties or conservation character the networks. Nodal features q: are specified by the datasets, and edge features Р., (dog)., are taken as the combinatorial gradient of the nodal features. In order to determine good hyperparameter configurations for each bracket, Bayesian search is conducted using Weights and Biases (76] for each bracket and each dataset using a random 80/10/10 train/valid/test split with random seed 123. The number of runs per bracket was 500 for CORA, CiteSeer, and PubMed, and 250 for CoauthorCS, Computer, and Photo. The hyperparameter configurations leading to the best validation accuracy are used when carrying out the experiments Table 4 and Table
Specifically, the hyperparameters that are optimized are as follows: initial learning rate (from 0.0005 0.05), number of training epochs (from 25 to 150), method of integration (rk4 or dopri5), integration time (from 1 to 5), latent dimension (from 10 to 150 in increments of 10), pre-attention mechanism (see below), positive function f (either exp or Squareplus), number of pre-attention heads (from to 15, c.f. Remark 1.14 attention embedding dimension (from 1x 15x the number of heads), weight decay rate (from 0.05), dropout/input dropout rates (from to 0.8), and the MI.P activation function for the metriplectic bracket (cither relu, tanh, squareplus). The pre- pre-attention is chosen
29

---

from one of four choices, defined as follows:
(Wqq.) ) Woqi
a (qis qi)
scaled dot product,
d
(W ) Woq
a (qi q,)
cosine similarity,
W xq Wooi
(W kq: Wxq) (Woq: - Woq))
ä (q: .)) =
Pearson correlation
WKq: Wxq
Waq: Woq
Wu Wou,
W KX, W "oxx1""
à (q: 9,) (o. exp
exp
exponential kernel
2
2
Network architectures. The architectures used for this experiment follow that of GRAND [49]. consisting of the learnable affine encoder/decoder networks Φ, 1 and learnable bracket- dynam- ics in the latent space. However, recall that the bracket-based dynamics require edge features, which are manufactured as Pa = (dog).- In summary, the inference procedure as follows:
q(0) = p(q)
(nodal feature encoding);
p(0) doq(0)
(edge feature manufacturing),
//" (q, (q. (latent dynamics)
= (q(0). +
q (q(T)),
(nodal feature decoding)
y = c(q).
(class prediction)
Training is accomplished using the standard cross entropy
H (t,y)
Σ t log
where t; is the one-hot truth vector corresponding to the 7th node. In the case of the metriplectic network, the networks J', gF; Is are 2-layer MLPs with hidden dimension equal to the latent feature dimension and output dimension 1.
B.3.2 Additional depth study
Here we report the results of a depth study on Cora with the Planetoid train/val/test split. Table shows the train/test accuracy of the different bracket architectures under two different increased depth conditions, labeled Task and Task respectively. Task 1 refers to fixing the integration step-size at variable final time T', while Task 2 instead fixes the final time T to the
At = 1 and integrating value identified by the hyperparameter search (see Table and instead varies the step-size At. Notice that both tasks involve repeatedly composing the trained network and hence simulate increasing depth, so that any negative effects of network construction such as oversmoothing, oversquashing, or vanishing/explodien gradients should appear in both cases. For more information, Table 10 provides a runtime comparison corresonding to the depth studies in Table
Observe that every bracket-based architecture exhibits very stable performance in Task 2, where the final time is held fixed while the depth is increased. This suggests that our proposed networks are dynamically stable and effectively mitigate negative effects like oversmoothing which are brought by repeated composition and often seen in more standard GNNs. Interestingly, despite success on Task 2, only the gradient and metriplectic architectures perfectly maintain improve their performance during the more adversarial Task I where the final time increased with fixed step-size. This suggests that, without strong diffusion, the advection experienced during conservative dynamics has the potential to radically change label classification over time, as information is moved through the feature domain in loss-less fashion.
Remark B.1. It is interesting that the architecture most known for oversmoothing (ie., grudient) exhibits the most improved classification performance with increasing depth on Task This is perhaps due to the fact that the gradient system decouples over nodes and edges, while the others
30

---

do not, meaning that the gradient network does not have the added challenge of learning a useful association between the manufactured edge feature information and the nodal labels. It remains to be seen if purely node-based bracket dynamics exhibit the same characteristics as the node-edge formulations presented here.
CORA networks
Trainable Parameters
Integration Time
Hamiltonian
60723
1.49625
Gradient
160772
14.82404
Double Bracket
30718
5.36151
Metriplectic 104088
7.53107
Table 8: The integration time and number of trainable parameters corresponding to the best networks trained on CORA. Note that integration time can be considered as surrogate for depth, since temporal step-size of each network is fixed to 1.
Depth study (acc)
layer
2 layers
4 layers
8 layers
layers
layers
layers
Hamiltonian
75.0 + 1.4
77.9 1.0
62.000.6
38.8 1.0
32.0 0.8
0.5
17.1
Gradient
69.6 1.1
11 1.
7.0
8.2
82.0 1.
Double Bracket
1.1
8.2
13.9
1.
1.3
1.
1.. 61.2
Metriplectic
61.0 1.6
6..4
0.8
1.3
61.6 1.0
60.6
Hamiltonian
77.2
77.5
12
7.7 0.6
7..4 1.
Gradient
1..9
79.9
79.6
7..6+1.1
80.41.11
0.9
79.8 + 1.1
Double Bracket
84.2
84.5
84.1
84.2
84.1
0.8
83.8
84.2 + 1.0
Metriplectic 6..71144
0.9
61.0
1.
1.
61.6 1.4
61.9 1.2
Table 9: Accuracy results depth study on CORA. Test accuracies reported as m over 10 runs with random train/valid/test splits. Top/bottom groups correspond tasks I and 2 of the study, respectively, where Task uses a fixed step-size while Task 2 uses a fixed integration domain.
Depth study (time)
layer
2 layers
4 layers
layers
16 layers
32 layers
64 lavers
Hamiltonian
0.008 + 0.003
0.060
0.006
0.106 0.011
0.191 0.005
0.308
0.035
0.497
0.045
0.874 0.035
Gradient
0.053 1 0.005
0.091
0.011
0.159 1: 0.015
0.273 0.023
0.470 1: 0.031
0.809
0.037
1.44 1.4.1: 0.038
Double Bracket
0.068 1 0.007
0.125
a 0.010
0.232
1: 0.023
0,009
0.770
0.688
1366 1: 0.098
2.52 0.100
Metriplectic
0.161 = 0.010
0.305
= 0.011
0.574 + 0.014
1.10 0.012
2.13 0.010
111 + 0.075
7,86 0.129
Hamilfonian
0.052 # 0.009
0.067
0.013
0.114 + 0.021
0.200
0.350 = 0.051
0.613 = 0.056
1.7
Gradient
0.365 0.020
0.012
1.20
2.26
0.058
0.009
8.36 +
19.4 0.152
Double Bracket
0.241
+: 0.036
0.304 0.046
0.731 +: 0.057
1.44 1: 0.132
2.98 0.286
5.42 4: 0.609
9.44 0.487
Metriplectic
2.06 0.097
Table 10: Runtime results of a depth study on CORA. Wall clock times reported as mean=stdev over 10 runs with random train/valid/test splits. Top/bottom groups correspond tasks I and 2 of the study, respectively, where Task uses a fixed step-size while Task 2 uses a fixed integration domain.
31


Reviews:
[{"review_id": "Mf7bbxGzTQ", "reviewer": "Reviewer_Ui2h", "summary": "This paper proposes a bracket-based dynamical system framework to design structure-preserving graph neural networks (GNNs). Specifically, the authors leverage four formalisms: Hamiltonian, Gradient, Double-Bracket, and Metriplectic, that model physical systems with different completeness and dissipation characteristics, and implement them via discretizations from exterior calculus. The authors parameterize GNNs with a chosen dynamic by specifying the bracket matrices and learning the energy/entropy functions via graph attention mechanisms. Experiments on real-world graph inference tasks show that structure-preserving GNNs can perform better than those that are not strictly structure-preserving, and the best dynamics depend on the downstream task. ", "questions": "1. What's the implications for practitioners? What applications would we expect imposing certain bracket structure outperforms the unconstrained ones? What structural requirement one should choose for different applications (e.g., add a column in Table 1 giving example applications)? The current set of experiments reads a bit confusing to me. For example, experiment 2 considers the node classification task, which requires some sense of clustering / smoothing, but not over-smoothing. So the findings in Table 4 showing that Double Bracket performs the best aligns with our expectations? But in Table 8, the Double Bracket formalism cannot be computed for deeper architectures - Is it merely a computational constraint? More explicit remarks and discussions will be great!\n\n2. The role of (learnable inner product) $A$: the authors connect this with graph attention mechanism; but a priori one can also choose other forms of $A$? Have the authors compare some simple baselines, say fixed $A$ or even just set $A$ as identity? This could also help address the memory issue in experiments (Sec 5.3) where the authors cannot compute $A$ (and thus $E, S$) due to the high-dimensionality of the node features?\n\n3. The authors discussed higher-order attention mechanism based on higher-order cliques (Appendix 4). Will this be useful in practice or is this a pure theoretical discussion? If the former, can the authors provide experimental demonstration? If the latter, the authors might consider removing this from the primary contributions since this does not seem like the focus of the paper?", "limitations": "This work does not have negative societal impact. The limitations are discussed in Section 6.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "Strengths:\n\n1. It is novel to use bracket-based parameterizations to design physics-inspired GNNs.\n\n2. It is fruitful to analyze existing GNNs (e.g., GAT, GRAND) under the bracket-based dynamical framework.", "weaknesses": "Weaknesses:\n\n1. Since the paper is intended for general audience in graph machine learning, it can be written in a clearer manner with more motivating examples and targeting applications. For starters, include a notation section.\n\n2. One of the key motivations of the paper is to improve performance of deep GNNs, yet all the experiments in the main manuscripts do not investigate explicitly the role of depth - are the GNNs used in these experiments actually deep? The more related experiments are described and discussed in Appendix B.3.2 as ablation study, but results (Table 8) show that for node classification task, only Gradient-based GNNs perform well when the number of layers increases, whereas other dynamics exhibit (significant) performance drop. If the paper intends to mitigate over-smoothing, the authors should provide more compelling experiments. On the other hand, if the paper simply intends to describe a physics-inspired framework to parameterize GNNs, then the authors should re-write the discussion of deep GNNs.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review.  We are happy to hear that you enjoyed our bracket-based parameterizations and that you appreciate the analysis of existing GNN architectures in terms of this framework.  Your weakness related to the clarity of the exposition has been well received: we plan to include a section at the beginning of Appendix A which summarizes our notation, and we will also add some discussion of additional examples and applications in Appendices A.1 and A.2, along with references where the interested reader can find more information.  Additionally, we have gone through the body of the paper again and made sure that the takeaways of our work are clear even if the analysis cannot be precisely followed.  If you have other suggestions, we are happy to incorporate them.\n\nTo address your weakness related to the role of depth in our networks, we have modified the depth study as described in the global reviewer response, so that the role of depth in our architectures is more apparent.  Particularly, Tables 9-11 provide additional information related to the role of depth in the trained CORA models, and it can be seen from Table 10 that the performance of all bracket-based architectures does not depend critically on the depth in terms of number of layers, maintaining itself very well when the depth is increased while the interval of integration is held constant.\n\nTo address your specific questions, please see below:\n\n1) *What's the implications for practitioners? What applications would we expect imposing certain bracket structure outperforms the unconstrained ones? What structural requirement one should choose for different applications (e.g., add a column in Table 1 giving example applications)? The current set of experiments reads a bit confusing to me. For example, experiment 2 considers the node classification task, which requires some sense of clustering / smoothing, but not over-smoothing. So the findings in Table 4 showing that Double Bracket performs the best aligns with our expectations? But in Table 8, the Double Bracket formalism cannot be computed for deeper architectures - Is it merely a computational constraint? More explicit remarks and discussions will be great!*\n\nFirst, please note that the mentioned deficiency in the Double Bracket architecture was artificial and has been addressed since our submission (please see the global reviewer response for details).  As to the impact of our work for practitioners, this is a complex question.  We notice that, in almost every case, a combination of reversible and irreversible dynamics is necessary for optimal performance, which agrees with the intuition that network dynamics should simplify feature information without homogenizing it.  As for the question of when to choose structure-preserving networks over unconstrained ones: we chose not to focus on this at the present time, as there are many papers in the literature (e.g., [6, 19-26]) which have established the benefits of structure-preservation in machine learning – including dynamical stability, improved performance with depth, and relative insensitivity to initialization – on several different classes of problems.  Our work adds to this literature by showing that, perhaps counterintuitively, comparable performance can also be obtained on node classification experiments with architectures based on several different physical mechanisms.\n\n2) *The role of (learnable inner product) A: the authors connect this with graph attention mechanism; but a priori one can also choose other forms of A? Have the authors compare some simple baselines, say fixed A or even just set A as identity? This could also help address the memory issue in experiments (Sec 5.3) where the authors cannot compute A (and thus E,S) due to the high-dimensionality of the node features?*\n\nIt is certainly possible to choose other forms for the learnable matrix A which do not connect to the standard idea of graph attention: the only fundamental requirement A should have is positive definiteness to induce a valid norm.  For example, one could choose learnable diagonal matrices as in [10], or parameterize A in terms of learnable Cholesky factors.  We have indeed experimented with setting A to identity, but find that including an attentional A improves performance, likely for the same reasons that it is advantageous to incorporate graph attention in standard GNNs.  We have also experimented with setting the A matrices to update only at the start of every epoch, so that they are constant during the integration (this is referred to as “constant attention” in [48]).  We find that this choice makes very little difference in our experiments, so we have not included a comparative study to this effect.\n\n3) *The authors discussed higher-order attention mechanism based on higher-order cliques (Appendix 4). Will this be useful in practice or is this a pure theoretical discussion? If the former, can the authors provide experimental demonstration? If the latter, the authors might consider removing this from the primary contributions since this does not seem like the focus of the paper?*\n\nWhile the discussion in the paper is purely theoretical, we expect higher-order attention to be useful on datasets with, e.g., explicit edge features or relatively dense connectivity.  We had difficulty identifying an appropriate edge feature dataset when we were designing the experiments for this work.  However, since the time of submission, we have improved the implementation of our attentional inner products, including the 2-clique attention which appears in the edge equation of the Gradient bracket.  These changes do appear to yield a modest performance increase, and this code will be released following the review process.\n"}, {"review_id": "KKWRklhlqv", "reviewer": "Reviewer_y9dF", "summary": "The work proposes structure-preserving bracket-based dynamical systems to learn physical systems using GNNs. The authors proposed four formulations depending on completeness and character (conservative or dissipative). The models are demonstrated to be effective in physical system and node classification tasks.", "questions": "* Does the method work well for PDEs, e.g., heat and Navier–Stokes equations?", "limitations": "The authors adequately addressed the limitations.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The method is constructed based on solid mathematics, and incorporating attention as a learnable inner product seems quite natural.\n* The authors suggested various formalisms depending on completeness and character under the notion of the abstract bracket formulations, which offers a general viewpoint for physical systems.\n* The method is quite adaptive because it is demonstrated to be effective in various tasks, including physical system and node classification.", "weaknesses": "* It is unclear to the reviewer where is the authors' contribution and where it is from prior works, particularly the mathematical part. The presentation could be more clear to distinguish existing work and the contribution of the work.\n* The authors claim, \"We use these architectures to systematically evaluate the role of (ir)reversibility in the performance of deep GNNs.\" However, the reviewer could not catch \"the role of (ir)reversibility\" from the experiments. The authors could elaborate on this point in the experiments part.\n* Despite the broad range of experiments considered, the performance of the models is moderate. Therefore, it is unclear in which condition the proposed models are effective (see the following weakness for more details).\n* In Section 5.1, the Hamiltonian model works better than the gradient model on the dissipative phenomena. The author could explain why the conservative works better here as it could be key to understanding \"the role of (ir)reversibility.\"\n* In addition, there are no experiments where the gradient model works best. With these results, it is unclear when the gradient model is effective. The authors could add experiments where the gradient works fine, e.g., incomplete and totally dissipative system.\n* In Section 5.3, the authors reported that the proposed model consumed much memory. In this case, computation time for prediction could also be reported because there are typically tradeoffs between speed and accuracy.\n\nMinor points:\n* Pointers to \"Table 3\" (49, 124, and 159) might be \"Table 1.\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review. We are glad that you appreciate the idea of attention as a learnable inner product, and the versatility of our bracket-based architectures in capturing physical principles and performing a variety of tasks.\n\nWe believe the primary weakness mentioned in your review – that it is not obvious where our mathematical contributions are relative to related works – can be cleared up easily. First, the work on data-driven exterior calculus in [10] deals fundamentally with identifying elliptic PDEs for data-driven physics modeling, and has nothing to do with dynamical systems, conservation properties, or bracket formalisms. Moreover, [10] presents no analysis related to graph attention or existing GNN architectures in terms of learnable inner products. So, while it is true that both [10] and our paper rely on the graph exterior calculus, their goals and contributions are essentially disjoint. Similarly, the work in [28] has a fundamentally different focus: they consider one specific type of bracket-based dynamical system, the metriplectic system, with the goal of learning physics outside the graph setting. Moreover, that work makes no use of the graph calculus or learnable inner products, meaning that [28] also has no direct connection to graph attention or existing GNN architectures.\n\nResponses to your remaining concerns are below. Please note truncation due to the character limit.\n\n*The authors claim, \"We use these architectures to…*\n\nYou are right that it is not easy to precisely characterize when one should use a reversible GNN architecture over an irreversible one. Our experiments show that both types of architectures can achieve competitive performance, though it seems that combining both conservative and dissipative dynamics is often the best choice. This suggests that the precise role of (ir)reversibility is context-dependent:  in physics-based modeling there are obvious reasons to use a network which is invariant-preserving, but enforcing physics may not be as critical for graph classification. Regarding depth, please see the attached Tables 9-11 and the global reviewer response.\n\n*Despite the broad range of experiments…*\n\nThis is fair criticism, as our architectures only achieve the best performance in some of the experiments discussed. However, note that the goal of this work was not to develop an architecture which achieves state-of-the-art performance: this would require significantly more hyperparameter tuning, as well as additional ad hoc modifications in terms of feature encoding, graph rewiring, etc. Instead, we aim to establish a unified framework for understanding graph networks, in the process bringing some clarity to the roles of reversibility and irreversibility in modern GNN architectures using benchmarks from both graph analytics and physics-based modeling. By connecting ideas such as gradient flows and graph attention to bracket-based dynamical systems, we have shown that architectures such as GAT and GRAND are not conservative or totally dissipative, and that this has non-obvious consequences on network performance.\n\n*In Section 5.1, the Hamiltonian model works better…*\n\nThis is a good point. The key to understanding this is to note that, in order to achieve coordinate-independence, our brackets are wrapped in a feature encoder/decoder.  This is necessary for learning bracket-based dynamical systems from arbitrary data: it is highly unlikely that a given set of data will match the form of our structure-preserving brackets in any arbitrary set of user-collected features, even if this data is drawn from a system which is physics-constrained. So, the autoencoder provides additional flexibility to our architectures, but comes with a challenge: the bracket structure we prescribe is enforced only in the latent space. Section 5.1 shows that it is possible for machine learned conservative dynamics to still produce an effective model for the dissipative double pendulum system when combined with a learnable encoding scheme.  \n\nIn simple problems, this issue is circumvented by using data with specially designed features which match a known bracket structure, e.g., position and momentum for canonical Hamiltonian systems. This removes the need for an autoencoder, but also limits the ability of bracket-based architectures to discover unknown physics or apply to less physical problems such as graph classification. For this reason, we have chosen to focus on the more general case of arbitrary features.\n\n*In addition, there are no experiments…*\n\nYou are right that our experiments do not show a case where the gradient system performs the best. On the other hand, the original depth study in Table 8 shows that the gradient system appears to be the highly stable on node classification problems with respect to perturbations in the integration domain. See, e.g., Remark B.1.\n\n*In Section 5.3, the authors reported…*\n\nThis is a good suggestion. We have included some computational timings on CORA in Table 11 as part of the new depth study. Note that the memory issues which were observed have been mitigated: they no longer appear in the new version of the code, though the metriplectic bracket is still prohibitively expensive on the large node classification problems despite its linear scaling in the graph size, due to its reliance on automatic differentiation during the forward pass. We anticipate that the metriplectic model will still have substantial impact in the context of data-driven physics models: in these settings feature dimensions are low. For evidence in this direction, please refer to the results in Tables 2 and 3.\n\n*Does the method work well for PDEs…?*\n\nMetriplectic dynamics have tremendous potential for multiphysics modeling. Reduced-order physics models obey a fluctuation-dissipation theorem, and the framework introduced here provides a means of learning non-equilibrium statistical mechanics. We are currently investigating this in other projects.\n"}, {"review_id": "XAbmK2HS0q", "reviewer": "Reviewer_7QNL", "summary": "This work provides a comprehensive overview of graph attention networks (GATs), shedding light on their fundamental concepts and principles. Additionally, the authors introduce a set of novel GNN architectures that leverage structure-preserving bracket-based dynamical systems. By incorporating these systems into GNNs, the proposed architectures aim to enhance the modeling capabilities and performance of graph-based machine learning tasks. Overall, the paper presents both a broad perspective on GATs and innovative approaches to further advance the field.", "questions": "I'd like to know whether your analysis can explain why \"Gradient“ model performs worse on heterophily graphs like Computer and Photo. Which physics system is suitable for heterophily graphs and which is suitable for homophily graphs?", "limitations": "They have already listed their limitations.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "They give a detailed theoretical analysis of graph attention networks with exterior calculus techniques. They also provide a generalized attention mechanism and propose several novel structure-preserving GNNs with satisfying performance.", "weaknesses": "The table writing is not clear. Sometimes orange result is better than blue result while sometime blue is better.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review.  We are glad that you appreciate our analysis of graph attention networks in the context of the graph exterior calculus, as well as our interpretation of the attention mechanism as a learnable inner product on features.  Please see below for responses to your specific concerns:\n\n*The table writing is not clear. Sometimes orange result is better than blue result while sometime blue is better.*\n\nIt is true that blue and orange do not always mean “best” and “second best” in our Tables.  In particular, we have used orange to represent the best result of the bracket-based architectures, and blue to represent the best result used for comparison.  Although this has been mentioned at the beginning of Section 5, we will reiterate this in the captions of our Tables since it was confusing.\n\n*I'd like to know whether your analysis can explain why \"Gradient“ model performs worse on heterophily graphs like Computer and Photo. Which physics system is suitable for heterophily graphs and which is suitable for homophily graphs?*\n\nThis is a good question.  In general, there is not an easy way to determine whether a particular structure-preserving architecture will be suitable for a particular degree of homophily or heterophily, as it has been shown (see, e.g., [*] below) that the ability of a network to distinguish between neighboring nodes relies crucially on the spectra of the learned weights.  On the other hand, it can be informally expected that architectures which are fully or partially conservative – such as the Hamiltonian, Double Bracket, and Metriplectic systems – will have an easier time learning on highly heterophilic data.  Indeed, this is what we observe on Computer and Photo in Table 5, where we see that the double bracket architecture is most performant.  This is because the totally dissipative Gradient architecture, along with many other architectures based on gradient flows, has a strong tendency to assign similar predictions to nodes in the same local neighborhood, which is colloquially known as “oversmoothing”.  While it is possible in theory for a totally dissipative architecture to correctly classify highly heterophilic data, the connection of dissipation/diffusion to graph convolution appears to make this difficult in practice.\n\n[*]  Francesco Di Giovanni, James Rowbottom, Benjamin P. Chamberlain, Thomas Markovich, and Michael M. Bronstein, “Understanding convolution on graphs via energies”, *arXiv:2206.10991*, 2023.\n"}, {"review_id": "vCjAN8NSZP", "reviewer": "Reviewer_dWFZ", "summary": "This paper provides a unified framework inspired the bracket-based dynamical system to analysis the oversmoothing problem in GNN. The past work may leverage the opposite physics concept such as the reversible processes, irreversible process and therefore it is not clear how such concept help to design GNNs, while the framework in this paper gives a deeper understanding. Leveraging the data-driven exterior calculus, this paper constructs four novel architectures which span the both reversibility and irreversibility spectrum using geometric brackets as a means of parameterizing dynamics abstractly without empirically assuming a physical model. Interestingly, it can reinterpretate the message-passing  and GAT as the fluxes and conservation balances of physics simulator.  It also generalize the attention mechanism which extends the nodal feature to higher order cliques and provide a unified evaluation of dissipation.  Emperically, the author compare the architectures derived from this framework with classical GNN,e.g., GAT, GDE over several benchmarks.", "questions": "I am curious about the precise computational complexity of the architectures proposed by the author, specifically in the context of N nodes or N neighbors. It would be beneficial to have concrete values rather than just the order of complexity. Additionally, it would be valuable to understand the advantages of these architectures in comparison to traditional graph neural networks (GNNs).\n\nRegarding the node classification experiment, the algorithm's performance is generally comparable to the baseline, although it does appear weaker in some cases. However, it showcases excellent performance in physical systems like the double pendulum and MuJoCo dynamics. Could this be attributed to the notion that physics-inspired frameworks are inherently more suitable for describing and modeling physical systems rather than social networks?\n", "limitations": "Yes", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. This paper is well-written and offers valuable insights into the construction of graph neural networks. As far as I am aware, the idea presented in this paper is innovative. However, since I lack a background in physics, I would appreciate hearing the suggestions of other reviewers who are familiar with bracket-based dynamical systems.\n\n2. The experimental result in damped double pendulum and MuJoCO dynamics looks good.\n\n3. The author made the code available, promoting reproducibility.", "weaknesses": "1. This paper may be a little hard to understand for the reader without physics background.\n\n2. It seems that the algorithm proposed by the author does not outperform the baselines in node classification problem.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review.  We are happy that you have enjoyed our paper and that you appreciate the results of the physics-based simulations.  Below are responses to your specific concerns in the “Weaknesses” section:\n\n1. *This paper may be a little hard to understand for the reader without physics background.*\n\nThis is a fair point, as the physics underlying our architectures is reasonably involved.  Unfortunately, we are limited by space constraints in the main manuscript and cannot give a comprehensive explanation whenever concepts are introduced.  However, note that we have included the Appendices A.1 and A.2 which provide an introductory primer to the graph exterior calculus and bracket-based dynamical systems.  Additionally, we also plan to add some more motivating examples and literature references to these Appendices, as well as go carefully through the body to ensure that the meaning of our results is preserved even if the analysis cannot be precisely followed.  If you have additional suggestions regarding how to improve the readability of our paper, we are happy to incorporate them.\n\n2. *It seems that the algorithm proposed by the author does not outperform the baselines in node classification problem.*\n\nThis is true in most cases.  However, note that the goal of this work was not to develop an architecture which achieves state-of-the-art performance: this would require more significantly more hyperparameter tuning than we have done here, as well as additional empirical modifications in terms of feature encoding, graph rewiring, etc.  Instead, we aim to bring some clarity to the roles of reversibility and irreversibility in modern GNN architectures using common benchmarks from both graph analytics and physics-based modeling.  By connecting ideas such as gradient flows and graph attention to bracket-based dynamical systems, we have shown that architectures such as GAT and GRAND are not conservative or totally dissipative, and that this has non-obvious consequences on network performance.  \n\nRegarding your questions, please see the following responses:\n\n*I am curious about the precise computational complexity of the architectures proposed by the author, specifically in the context of N nodes or N neighbors. It would be beneficial to have concrete values rather than just the order of complexity. Additionally, it would be valuable to understand the advantages of these architectures in comparison to traditional graph neural networks (GNNs).*\n\nThe primary advantages of the architectures presented here over traditional GNNs are (1) explicit and automatic preservation of desirable physical properties such as the first and second laws of thermodynamics, and (2) straightforward incorporation of the graph attention mechanism directly into the differential operators governing network evolution.  By posing attentional GNNs as bracket-based dynamical systems, we make these networks more interpretable, conceptually simpler, and guaranteed to satisfy useful properties (e.g., dynamical stability) which follow from established theory.   This leads to better analysis of network properties and performance, as well as architectures which provably respect physical laws often found in training data.\n\nTo address the question of computational complexity:  note that our architectures are predicated on the graph exterior calculus, which encodes derivatives in terms of signed incidence matrices d_0, d_1, etc.  These derivatives are always local on k-cliques, and our attention matrices A_0, A_1, etc. are always diagonal.  This means that operations on graph entities are always sparse, and the complexity will therefore be linear in the graph size.  On the other hand, it will also depend on the number of matrix multiplies called for by the evolution equations, the dimension of the features, and the time integrator itself, along with any need to, e.g., compute automatic derivatives during the forward pass.  Based on this information, it can be seen that the Hamiltonian bracket is relatively cheap, and the Metriplectic bracket is the most expensive.  A rough estimate of this complexity in terms of numerical runtime on CORA is given in the new Table 11.\n\n*Regarding the node classification experiment, the algorithm's performance is generally comparable to the baseline, although it does appear weaker in some cases. However, it showcases excellent performance in physical systems like the double pendulum and MuJoCo dynamics. Could this be attributed to the notion that physics-inspired frameworks are inherently more suitable for describing and modeling physical systems rather than social networks?*\n\nThere may be some truth to this.  You are right that we do notice better results from structure-preservation on physics-based examples where there is a clear notion of energy/entropy, and the goal is to reproduce the trajectories of dynamical systems.  Conversely, it appears that graph classification networks can achieve comparable performance with many different types of structure-preserving and structure-agnostic architectures, perhaps indicating that the dynamics of the internal state in these networks are not critical to understanding the underlying classification problem.  With that said, we also notice remarkable stability of our structure-preserving architectures with increasing depth (see the attached Table 10), which is a notable advantage over standard GNN architectures even in this case.\n"}]

Author rebuttal:
Thank you all for your insightful comments and useful feedback on our work.  We have heard the shared criticism that (1) parts of our manuscript could be difficult to understand for a general machine learning audience, and (2) that the role of depth in our architectures could be made clearer in the text.   We are taking steps to address these concerns.  

Unfortunately, some technical complexity is unavoidable due to the mathematical machinery employed.  However, there are some additional things that we can do to ease the burden on our readers.  We have gone through the main manuscript and ensured that the practical implications of our work are clear, even if the analysis is difficult to follow.  Due to space constraints, it is difficult to include any more introductory details in the main body of the manuscript.  However, note that we have already written introductions to the graph exterior calculus in Appendix A.1 and to bracket-based dynamical systems in Appendix A.2, with the goal of making our manuscript accessible to the general machine learning community.  Following one reviewer’s suggestion, we have added a glossary of terms for ease of reference. We also plan to add more examples and common applications to these Appendices as per your suggestions.

Regarding the role of depth: we find that there is not an obvious relationship between depth and the performance of structure-preserving network architectures.  On the other hand, our experiments show that the bracket-based architectures presented here **can maintain their performance up to (and potentially beyond) 64 layers**, as demonstrated clearly in the new Table 10, while standard GNNs, including GATs, are known to experience significant degradation after ~4 layers.   In particular, many of the networks trained in Section 5 use the equivalent of 5-15 layers to reach their optimal predictions (see, e.g., the new Table 9).  Note that these new Tables correspond to a redesigned depth experiment which clarifies the influence that depth has on our architectures.  This new study, which has been designed to supplement Table 8 in Appendix B.3.2, takes our bracket-based networks trained on CORA with Planetoid splits, but now progressively reduces the step-size of their 4th-order Runge-Kutta time integration from 1 to 1/64 and tests their predictive accuracy on random splits.  This experiment tests the effect of increasing the depth of a network while holding its optimized domain of integration fixed and measures dynamical stability with increasing depth.  In the numerical analysis community, it is well known that the temporal step-size must be sufficiently small to achieve a stable result (the so called CFL conditions), which may not be respected by the experiment in Table 8 which fixes the temporal step-size to 1 and performs integration over progressively longer domains. The previous experiment (Table 8) was originally taken to match the methodology used in the GRAND paper, but we believe that Table 10 provides a more meaningful comparison.

In addition to these changes, we would like to reiterate that a primary contribution of our work is the principled reframing of graph neural networks in terms of bracket-based dynamical systems, which provides a rigorous way to both construct new network architectures as well as contextualize current ones.  Besides simplifying the analysis of existing architectures such as GAT and GRAND, this framework has great potential for data-driven physics applications, since reducing the order of a physics model naturally introduces dissipation in the form of information loss even when the full model is conservative.  Some evidence for the utility of the bracket formalism on physical problems is presented in Tables 2,3, and 6, and we are actively seeing success with this approach in other projects as well.

Finally, we would like to call your attention to the difference in results between the new depth study in Table 10 and Table 8 of the original submission.  This is because, in addition to reflecting a slightly modified experiment, the new results were obtained following a complete refactor of the code. Our original code was written around the publicly available GRAND codebase, but we had concerns related to memory usage and scalability of the GRAND code for deep networks (see missing data in tables 8 and 5). After reimplementing from scratch, we were able to remove the many hidden features of the GRAND code and avoid these issues. We are repeating all experiments with the new code but have completed preliminary tests that demonstrate the issues have been resolved, and this code will be made available following the review process. 


Meta-review: This paper proposes novel GNN architectures based on structure-preserving bracket-based dynamical systems. The reviewers appreciate the novel and solid contributions, and are generally positive about the paper. Most questions from the initial reviews were satisfactorily addressed. The authors committed to making the paper more accessible to a broader audience. We suggest the authors to incorporate the reviewers' suggestions in the final version.
Decision: Accept (poster)