Greedy Fiedler Spectral Partitioning for Data-driven Discrete Exterior Calculus Andy Huang1 , Nathaniel Trask1 , Christopher Brissette 2 , Xiaozhe Hu 3 1 Sandia National Laboratories, NM 2 Rensselaer Polytechnic Institute, NY 3 Tufts University, MA ahuang@sandia.gov, natrask@sandia.gov, brissc@rpi.edu, xiaozhe.hu@tufts.edu Abstract graph exterior calculus which allows the application of pop- ular back propagation techniques from machine learning to The data-driven discrete exterior calculus (DDEC) structure provides a novel machine learning architecture for discover- learn coordinate-free vector calculus operators from data. ing structure-preserving models which govern data, allowing The resulting operators are compatible with a discrete ver- for example machine learning of reduced order models for sion of Stokes’ theorem, allowing the learning of physical complex continuum scale physical systems. In this work, we conservation laws which naturally preserve conservation prin- present a Greedy Fiedler Spectral (GFS) partitioning method ciples, provide solvability guarantees for second-order elliptic to obtain a chain complex structure to support DDEC mod- operators, and naturally handle diffusion problems. A key els, incorporating synthetic data obtained from high-fidelity ingredient for this calculus is a chain complex, which infor- solutions to partial differential equations. We provide justi- mally corresponds to a graph structure generalizing meshes fication for the effectiveness of the resulting chain complex typically used to partition space. We propose a greedy Fiedler and demonstrate its DDEC model trained for Darcy flow on a spectral (GFS) partitioning scheme which extracts this graph heterogeneous domain. structure from eigenspectrum information of a given data set. Data-driven Discrete Exterior Calculus A discrete exterior calculus structure (Desbrun et al. 2005) A motivating example: resistor networks provides a notion of differentiation beyond the differentiation of scalar-valued functions. The familiar div, grad, and curl Before introducing DDEC in the abstract, we provide a fa- operators from vector calculus in R3 satisfy the identities miliar example to motivate both the discrete exterior calculus div ◦ curl = 0 and curl ◦ grad = 0, and also appear in the and the chain complex structure on which it is defined. The integral Green’s, Stokes’, and divergence theorems. Although discrete exterior calculus is best exemplified by Kirchhoff’s back-propagation provides the differentiation of tensors by current law and Ohm’s law from circuit theory (see (Weyl regarding the components as scalar functions, preservation of 1923), (Smale 1972), and (Bamberg and Sternberg 1991)). these identities are critical to handling a number of physical As a concrete example, consider an electrical circuit on systems, particularly electromagnetic phenomena. the graph in Figure 1. Suppose nodal electric potentials are In the context of finite element approximation, compati- ble spatial discretizations construct approximations faithful to these underlying structures, and form the bedrock of dis- [2] cretization of multiphysics problems, particularly for second- order elliptic systems (Arnold et al. 2007; Bochev and Hyman [1] 2006; Arnold, Falk, and Winther 2006). We consider in this [0] work the design of machine learning architectures which preserve a similar structure. [3] Stokes’ theorem plays a central role for expressing conser- vation laws on a Riemannian manifold: given a differential Figure 1: A simple circuit on a graph with 4 nodes. form ω on a manifold M , we have Z Z dω = ω (1) known: the voltages at nodes [0], [1], [2], and [3] are 10V , 5V , Ω ∂Ω 3V , and 2V , respectively; we represent this fact by the vector The data-driven exterior calculus (DDEC) (introduced in V = [10, 5, 3, 2]T whose entry i corresponds to the voltage (Trask, Huang, and Hu 2020)) is a parameterization of the V ([i]) at node [i]. Further, suppose that the oriented current Copyright ©2021 for this paper by its authors. Use permitted under on the edge [0, 1] is 1A, on edge [1, 2] is 0.7A, and on edge Creative Commons License Attribution 4.0 International (CC BY [1, 3] is 0.3A; we order the edge labels lexicographically and 4.0) so represent these currents as I = [1.0, 0.7, 0.3]. The discrete exterior calculus boundary operator ∂ is The data-driven exterior calculus " −1 +1 # We gather the requisite definitions from DDEC to define ∂= −1 +1 (2) the chain complex and how it’s used, referring the inter- −1 +1 ested reader to (Trask, Huang, and Hu 2020) for details. Let NN N = {ni }i=1 denote a set of nodes. We embed N in Rd and the coboundary operator is defined as δ := ∂ T ; explicitly, by associating with each node a unique position pi ⊂ Rd , it maps the edge [i, j] to the linear combination of nodes i ∈ 1, . . . , NN . A k-clique is an ordered tuple of k nodes; it +[j] − [i] accounting for edge orientation. is assigned an orientation equal the sign of the permutation In this language of discrete exterior calculus, Kirchhoff’s which puts these k nodes in increasing order. We associate current law states that δI vanishes on interior nodes. Indeed, each k-clique with the (k − 1)-simplex spanned by its k node  −1   −1.0  positions. A k-chain ck is a formal R-linear combination of 1.0 " # +1 −1 −1  0.0  (k + 1)-cliques, and we denote the space of k-chains by Ck . δI =  +1  0.7 =  0.7  (3) The boundary operator ∂k : Ck+1 → Ck is defined via 0.3 +1 0.3 k X and the second entry vanishes since it is the discrete graph di- ∂k [n1 , ..., nk ] = (−1)i−1 [n1 , ..., nbi , . . . , nk ], (7) vergence of I at node [1], i.e., it is the sum of all edge currents i leaving node [1]. Ohm’s law is realized by the equation where ˆ· denotes an omitted entry. It satisfies the property ∂V = RI (4) ∂k−1 ◦ ∂k = 0. This forms an exact sequence for which R represents the resistances on the edges. R is a C0 ∂0 C1 ∂1 ... ∂d−1 Cd (8) matrix with diagonal entries R[0,1] , R[1,2] , and R[1,3] corre- sponding to the resistances at the respective edges. with the convention that ∂−1 maps C0 to the empty set. However, Ohm’s Law is more accurately expressed as The cochains C k consist of linear functionals acting on Ck . Given φ ∈ C k , we denote the value associated with the R−1 ∂V = I (5) k-chain ti1 i2 ...ik via the shorthand φi1 i2 ...ik := φ(ti1 i2 ...ik ). and the operator δ ∗ := R−1 ∂ is a metric-dependent gradient Note that the cochains inherit the orientation of the under- (as opposed to the purely combinatorial ∂). The discrete lying chains, e.g. φij = −φji via the definition of π. Intro- Stokes’ theorem states ducing the coboundary operator δk : C k → C k+1 , we next Z Z arrive at the following cochain complex δ∗ V = V (6) [0,1] ∂[0,1] δ δ δd−1 C0 0 C1 1 ... Cd (9) That this holds is the expression of Ohm’s law. For example, on the edge [0, 1], we see that the left-hand side the resulting coboundary operators result in the traditional Z Z graph div/grad/curl from the combinatorial Hodge theory [REF]. Duality between the combinatorial boundary and δ∗ V = R−1 I [0,1] [0,1] coboundary operators is expressed through the inner product = R[0,1] I([0, 1]) hφ, ∂k ti = hδk φ, ti. (10) is equal to the right-hand side The DDEC operators on the chain complex are parameterized Z Z metric perturbations of these, defined as V = V ∂0 [0,1] [1]−[0] dk := Bk+1 δk B−1 k , and d∗k = D−1 ∗ k δk Dk+1 (11) = V ([1]) − V ([0]) where Bk and Dk are diagonal matrices for each k. In this paper, we specialize to chain complexes whose since V ([1]) − V ([0]) = R[0,1] I([0, 1]). largest cells are 2-dimensional, where we introduce notation While simple, this resistor network demonstrates how con- alluding to the vector calculus operators: servation statements can be encoded onto a graph structure which exposes trainable metric information (the resistances) CU RLh = B1 δ0 B−1 0 , DIVh = B2 δ1 B−1 1 , (12) to parameterize diffusion processes. In this example, the circuit network is a 1-dimensional chain complex which sup- CU RL∗h = D−1 ∗ 0 δ0 D1 , GRADh∗ = D−1 ∗ 1 δ1 D2 , (13) ports conservation (encoded in Kirchoff’s law) and diffusion (encoded in Ohm’s law) expressed using graph divergence DDEC chain complex through mesh partitioning and gradient operators. More generally, the DDEC expresses Given a domain mesh M consisting of n-dimensional cells a class of such parameterized vector calculus models, pro- used in a finite element discretization of a physics PDE prob- vided a chain complex can be obtained. lem, it may be tempting to use M itself as a graph structure. The chain complex provided by the GFS partitioning pro- However, working with M directly is typically prohibitively posed in this paper provides in this example a circuit topology expensive since it is subject to the curse of dimensionality. for a given data set. In addition, with a solution to the PDE on M in hand, it is suitable to incorporate that physics information in obtaining tioning of the entire chain complex. We will see in Figure a coarse representation of M for use with a DDEC. 8 an example induced chain complex partitioning where the We now describe partitioning of the cells of M as a means coarse cells correspond to the color-coded regions. to obtain a chain complex which comes equipped with a DDEC structure. We consider the topological dual graph T Greedy Fiedler Spectral Partitioning to M . Its construction is as follows: for every n-cell of M , We aim to apply a spectral partitioning scheme on the nodes associate to it a node of T ; for every pair of n-cells of M of T . Spectral methods have been used for compression of sharing an (n − 1)-dimensional face, create an edge between information and partitioning in many domains. Examples their corresponding nodes of T ; by repeating this process, include JPEG image compression, manifold eigenfunction i.e., creating a k-cell for T if k + 1 corresponding n-cells parameterization (Jones, Maggioni, and Schul 2008), segmen- of M share an (n − k)-dimensional face, then T becomes tation of 3D models (Sharma et al. 2009), and partitioning the topological dual to the mesh M . This is a chain complex of finite element meshes for parallel processing (Kaveh and with combinatorial δk and δk∗ operators as described in the Davaran 1999). preceding section. On the one hand, this discrete exterior calculus structure Graph Laplacian and Fiedler partitioning on T itself can be perturbed into a DDEC by introducing the Bk and Dk metrization as in (11). However, T is typically To define spectral partitioning of T , we first define some some prohibitively large and subject to the curse of dimensionality. mathematical objects used for that purpose. Consider the Instead, by partitioning the nodes of T , a coarsened chain graph obtained from the nodes and edges of T . The Laplacian complex can be obtained which naturally admits a DDEC on T is defined as structure. As noticed in (Trask, Huang, and Hu 2020), nodal L = D − A, functions on the fine-scale chain complex are projected onto where D is the degree of the nodes and A is the (weighted) the coarsened chain complex through averaging over each adjacency matrix. For non-negative weights, it was first no- partition. Since DDEC training involves minimizing the error ticed in (Fiedler 1975) that the eigenvector (called the Fiedler through approximation by these coarse functions, we are eigenvector) of L corresponding to the second smallest eigen- motivated to obtain a coarse chain complex which has good value ψ can be used to partition T into two connected sub- function approximation quality. graphs, T+ and T− , for which ψ is positive on the nodes of We illustrate the induced chain complex coarsening by a T+ and negative on the nodes of T− , and whose edges of partitioning of the nodes of T in the case it is a 2-dimensional T not belonging to T+ or T− are exactly those on which ψ chain complex. In Figure 2, F0 , F1 , and F2 denote the 0-, 1-, takes a positive value on one terminal node and a negative and 2-chains on T , i.e., the nodes, edges, and cells; F0 , F 1 , value on the other (ignoring the special case when ψ may and F 2 denote the 0-, 1-, and 2-chains on T , i.e., the linear take a value of 0 on some nodes). In this way, the value 0 can functionals on the nodes, edges, and cells. be treated as a cut for this Fiedler eigenvector of the graph Laplacian and used for graph partitioning. Typically, a recursive hierarchical bi-partitioning is used δ0f ine δ1f ine to further partition a graph beyond the first cut. Each of 0 F0 F1 F2 0 ι0 ι1 ι2 the two partitioned subgraphs can be cut using their Fiedler π0 1 2 eigenvector, and so on. However, for the purposes of good δ0 π δ1 π 0 C0 C1 C2 0 approximation, since the cochains take the average value of the field on each partition, any partitioning of the mesh near small field variations is unnecessary. So, for a given nodal 0 F0 F1 F2 0 ι0 ∂0f ine ι1 ∂1f ine ι2 field V on T , we consider a different iterative approach: as π0 π1 π2 in adaptive quadrature, at each iteration, we choose to cut 0 C0 C1 C2 0 the partition with the worst approximation by a piecewise ∂0 ∂1 constant function (which is bounded a priori by the L2 norm of V on each partition). The GFS partitioning scheme is Figure 2: A diagram of domain partitioning for a 2- detailed in Algorithm 1. We use |∇V | for the node weights dimensional chain complex. This summarizes the relationship of T , and each partition weight is the sum of all nodal weights between the nodes, edges, and faces of a T and the nodes, within the partition. edges, and faces of a coarsened chain complex resulting from partitioning T . Implementation We implemented the GFS partitioning in P YTHON using With a partition π2 of the nodes of T , we agglomerate the the packages N UM P Y, S CI P Y, and N ETWORK X for graph fine-scale chains of F2 into coarse chains C2 . By applying manipulation. The pseudocode in Algorithm 1 can be im- the boundary operator to these chains, we establish the coarse plemented using NETWORKX in a straightforward manner - 1-chains C1 from F1 , etc. For more details of this partitioning calculuation of edge weights by partition, deletion of edges, process, we defer to (Trask, Huang, and Hu 2020). For the and obtaining the positive and negative connected compo- purposes of this paper, it suffices to produce a partitioning nents after a Fiedler eigenvector nodal cut are possible using of the highest dimensional cells because it induces a parti- NETWORKX. The Fiedler eigenvectors were calculated using Data: points X = {xn }, field values Vn = V (xn ) Result: partitioning Π = (P1 , . . . , PN parts ) such that FN parts i=1 Pi = X Construct graph T = (X, Xij = (xi , xj )) Assign edge weights |Vj − Vi | on Xij Initialize: Π = [G], Wparts = [WG ] for j ∈ {1, ..., Nparts} do P ←Gparts [−1] Compute Fiedler eigenvector vP of P for eij ∈ EP do if vP (i) · vP (j) < 0 then Remove edge eij from P end end Identify partition: P = P+ t P− Update: Gparts [−1] = P− Append: Gparts + = [P+ ] for G ∈ Gparts do Figure 3: GFS approximations of f (x) = tanh(x)+2 versus Calculate and update partition weights a uniform approximation with the same number of partitions. end The four plots display the cochain approximations of f de- Sort Gparts with respect to partition weights fined by the GFS approximation and the uniform approxi- end mation with the same number of partitions. The L2 errors Algorithm 1: Greedy Fiedler Spectral Partitioning algo- associated with each method are also shown. rithm for partitioning the domain of a field. networkx.fiedler_vector. We note that a more ef- ficient implementation can be made possible by evaluating subgraph Laplacians efficiently; since L = D − A, and A restricts to subgraphs simply by referincing the appropriate submatrix, only the degree matrix D needs to be updated properly between iterative cuts. One-dimensional function approximation example Figure 4: Error and partition weight versus number of parti- To better understand the GFS algorithm, we apply it to a tions, again comparing partitioning by GFS and a uniform simple one-dimensional function f (x) = tanh(x) + 2 on partitioning scheme on f (x) = tanh(x) + 2. The left plot x ∈ [−5, 5]. This function plotted in Figure 3 along with shows the maximum, mean, and minimum over all partition approximate functions given by uniform partitions and GFS weights. The right plot displays the resulting L2 error of GFS partitions for comparison (the average value of f over each approximation versus uniform approximation, as a function partition is used for the approximate functions). The fine- of the number of partitions. k scale graph is given by the nodes { 1000 |k = 0, 1, . . . , 1000} k k+1 and the edges are {[ 1000 , 1000 ]|k = 0, 1, . . . , 999}. Uniform partitioning here refers to cutting at the equally spaced points the GFS method performs better and has less L2 error than k n for k = 1, 2, . . . , n − 1 for an n-partition scheme. that of uniform partitioning. Intuitively, to minimize the error of a coarse cochain ap- proximation - which is piecewise constant on each partition - it is prudent for a partitioning algorithm to apply most re- Two-dimensional electrostatics/Darcy flow example finement between [−2, 2] where the derivative is great. Given the weights |Vj − Vi | in Algorithm 1 we can see that parti- Next, we consider an application of the GFS partitioning to tions within the interval [−2, 2] will carry the most partition data obtained from a textbook electrostatics problem. Con- weights, and thus the greedy iteration to partition by the sider an infinite dialectric cylinder in R3 co-axial with the heaviest weight will result in most approximation fidelity in z-axis, of radius b and relative capacivity K. Due to the in- this region. Indeed, we see this refinement of partitioning variance in the z-direction, the problem is analyzed on the over [−2, 2] in Figure 3. Furthermore, Figure 4 compares x − y plane. In polar coordinates, the background electric the GFS algorithm to a uniform partitioning scheme, and field is given by Vback = Einf r cos(θ). The induced electric demonstrates GFS yielding a lower L2 -error, as expected. As potential due to the charge on the cylinder which matches shown in Figure 4, while the rate of convergence is the same, the boundary conditions imposed by Vback on the cylinder surface interface is described analytically by    Einf − K−1 b2 cos(θ) for r > b K+1 r V (r, θ) =  Einf  (14)  K+1 − E inf r cos(θ) for r ≤ b Notably, a sharp interface is experienced by the gradient of the field at r = b as depicted n Figure 5. Figure 5: Electric potential V (left) and norm of the gradient of the electric potential ||∇V || (right) over a Cartesian mesh of [0, 1] × [0, 1]. We obtained point data by sampling the analytic expres- sion for the electric potential due to the dialectric cylinder. For the purpose of this example, we chose a simple Cartesian mesh of [0, 1] × [0, 1] with 70 cells in each dimension. Figure 6 displays the partitioning of the 70 × 70 cells, comparing the GFS partition to a uniform partitioning obtained by M ETIS. We see a notable clustering of partitions near r = b. Further- more, in Figure 7, we display the better approximation error of GFS compared to M ETIS. Using GFS partitioning for DDEC: Example Electrostatics/Darcy flow problem Finally, we display the utility of a GFS partitioning to obtain a chain complex supporting a DDEC. The DDEC will be trained to solve a system of PDEs which includes a conserva- tion law, resulting in a machine learned surrogate model. First, we apply GFS partitioning to obtain a coarse chain complex. In Figure 8, we display the GFS partitioning of the 70 × 70 unit square fine-scale cells into 16 coarse cells, and the resulting chain complex. Increased resolution is seen Figure 6: Juxtaposition of domain partitioning performed near the r = b interface. Nodes, edges, and faces shaded by by M ETIS (left column) versus GFS (right column), using different colors to distinguish the different components. electric field magnitude through a dialectric cylinder as node Next, on that GFS partitioning-obtained chain complex, weights. The rows correspond to 16, 32, 64, and 128 resulting we pose the electrostatics problem partitions, in increasing order from top to bottom. ~ D = −∇φ (15) ~ ∇·D = 0 (16) ~ φ,  are the electric displacement field, the electric where we have identified F as the electric displacement field where D, and the φ as the electric potential. The field κ was taken to potential, and the inhomogeneous capacivity. Mathematically be 10 inside the circle of radius 0.1 centered at (0.6, 0.57); equivalently, this system describes the Darcy flow of a fluid, the asymmetry of this configuration was chosen to break the where these variables correspond to the velocity vector field, symmetry of domain bisection. the pressure, and the permeability, respectively, and where Finally, we train this DDEC model on the boundary value the symbols ~u, p, and κ more commonly employed. In the data of the analytic solution (following the procedure in DDEC formulation, we have the system of equations (Trask, Huang, and Hu 2020) to learn F and φ). Figure 9 F + κ∇φ = 0 w1 − GRADh∗ u0 = 0 displays good agreement between the DDEC φ with the aver- −→ (17) age of the analytic solution over each partition. ∇·F=f DIVh w1 = f0 not necessarily represent the views of the U.S. Department of Energy or the United States Government. SAND Number: SAND2020-XXXX. The work of N. Trask is supported by the U.S. Depart- ment of Energy, Office of Advanced Scientific Comput- ing Research under the Collaboratory on Mathematics and Physics-Informed Learning Machines for Multiscale and Multiphysics Problems (PhILMs) project. N. Trask is sup- Figure 7: Approximation error comparison between piece- ported by the Department of Energy early career program. wise constant averaging from METIS and GFS. The work of A. Huang, N. Trask, and C. Brissette is supported by the Sandia National Laboratories Laboratory Directed Re- search and Development (LDRD) program. References Arnold, D. N.; Bochev, P. B.; Lehoucq, R. B.; Nicolaides, R. A.; and Shashkov, M. 2007. Compatible spatial discretiza- tions, volume 142. Springer Science & Business Media. Arnold, D. N.; Falk, R. S.; and Winther, R. 2006. Finite element exterior calculus, homological techniques, and appli- cations . Bamberg, P.; and Sternberg, S. 1991. A Course in Mathemat- ics for Students of Physics, volume 2. Figure 8: The DDEC complex obtained by GFS partitioning Bochev, P. B.; and Hyman, J. M. 2006. Principles of mimetic of the unit square mesh using the electric potential gradient. discretizations of differential operators. In Compatible spatial Partitioning results in coarse cells, edges, and nodes. discretizations, 89–119. Springer. Desbrun, M.; Hirani, A. N.; Leok, M.; and Marsden, J. E. 2005. Discrete Exterior Calculus . Conclusion and Further Work Fiedler, M. 1975. A property of eigenvectors of nonnega- In this article, we addressed the question of how to obtain tive symmetric matrices and its application to graph theory. a chain complex for use in the data-driven discrete exterior Czechoslovak Mathematical Journal 25: 619–633. calculus. We introduced a Greedy Fiedler Spectral partition- Jones, P. W.; Maggioni, M.; and Schul, R. 2008. Mani- ing scheme of a finite-element mesh to obtain a coarse chain fold parametrizations by eigenfunctions of the Laplacian and complex with favorable approximation of a fine-scale scalar heat kernels. Proceedings of the National Academy of Sci- field on the mesh. The question of how to obtain a suitable ences 105(6): 1803–1808. ISSN 0027-8424. doi:10.1073/ graph to support data-driven modeling is requisite for suc- pnas.0710175104. URL https://www.pnas.org/content/105/ cessful application of the DDEC. We demonstrate the utility 6/1803. of this calculus in a number of upcoming works. Kaveh, A.; and Davaran, A. 1999. Spectral bisection of adap- Acknowledgments. Sandia National Laboratories is a mul- tive finite element meshes for parallel processing. Computers timission laboratory managed and operated by National Tech- Structures 70(3): 315 – 323. ISSN 0045-7949. doi:https: nology and Engineering Solutions of Sandia, LLC, a wholly //doi.org/10.1016/S0045-7949(98)00170-9. URL http://www. owned subsidiary of Honeywell International Inc. for the sciencedirect.com/science/article/pii/S0045794998001709. U.S. Department of Energy’s National Nuclear Security Ad- Sharma, A.; Horaud, R.; Knossow, D.; and Von Lavante, E. ministration under contract DE-NA0003525. This paper de- 2009. Mesh Segmentation Using Laplacian Eigenvectors scribes objective technical results and analysis. Any subjec- and Gaussian Mixtures. In Souvenir, R., ed., AAAI Fall tive views or opinions that might be expressed in the paper do Symposium on Manifold Learning and its Applications, 50– 56. AAAI, Arlington, VA, United States: AAAI Press. URL https://hal.inria.fr/inria-00446990. Smale, S. 1972. On the mathematical foundations of elec- trical circuit theory. J. Differential Geom. 7(1-2): 193–210. doi:10.4310/jdg/1214430827. URL https://doi.org/10.4310/ jdg/1214430827. Trask, N.; Huang, A.; and Hu, X. 2020. Enforcing exact physics in scientific machine learning: a data-driven discrete exterior calculus on graphs. preprint arXiv:2012.11799. Figure 9: Trained DDEC model compared to analytic solution Weyl, H. 1923. Repartición de corriente en una red conduc- averaged over each partition (sliced through y = 0.53). tora. Revista Matemática Hispano-Americana 5: 153–164.