Learning Interface Conditions in Domain Decomposition Solvers
Authors: Ali Taghibakhshi
Nicolas Nytko
Tareq Uz Zaman
Scott MacLachlan
Luke Olson
Matthew West
Venue: NeurIPS 2022 Accept (2022)
Keywords: Domain Decomposition
Optimized Schwarz Methods
Graph Neural Networks
Unsupervised Learning

Abstract: Domain decomposition methods are widely used and effective in the approximation of solutions to partial differential equations.  Yet the \textit{optimal} construction of these methods requires tedious analysis and is often available only in simplified, structured-grid settings, limiting their use for more complex problems. In this work, we generalize optimized Schwarz domain decomposition methods to unstructured-grid problems, using Graph Convolutional Neural Networks (GCNNs) and unsupervised learning to learn optimal modifications at subdomain interfaces. A key ingredient in our approach is an improved loss function, enabling effective training on relatively small problems, but robust performance on arbitrarily large problems, with computational cost linear in problem size. The performance of the learned linear solvers is compared with both classical and optimized domain decomposition algorithms, for both structured- and unstructured-grid problems.


Learning Interface Conditions in Domain Decomposition Solvers
Ali Taghibakhshi
Nicolas Nytko
Mechanical Science and Engineering University of Illinois Urbana- Champaagg Urbana, IL 61801, USA alit2@illinois. edu
Computer Science
University of Illinois Urbana-Chaman Chamaa-Champain
Urbana, IL 61801, USA myykoo2illinmg
Tareq Zaman
Scott MacLachlan
Scientific Computing Program Memorial University of Newfoundland and Labrador
Mathematics and Statistics
Memorial University of Newfoundland and Labrador
St. John's, s, NL, Canada tzaman@mun.co
St. John't 's NL, Canada smaclachl an@mun
Luke Olson
Matthew West
Computer Science
Mechanical Science and Engineering University of Illinois Urbana-CHie Chammapan
University of Illinois Urbana- Chaamaaig Urbana, IL 61801, USA Iukeo@illinois edu
Urbana, IL 61801, USA mwest@illinois. edu
Abstract
Domain decomposition methods widely used and effective in the approxima- tion solutions to partial differential equations. Yet the optimal construction of these methods requires tedious analysis and is often available only in simpli- fied, structured-gret grid settings, limiting their use for more complex problems. this work, generalize optimized Schwarz domain decomposition methods unstructured-gnic problems, using Graph Convolutional Neural Networks (GCNNs) and unsupervised learning learn optimal modifications subdomain interfaces. A key ingredient in our approach an improved loss function, enabling effective training relatively small problems, but robust performance arbitrarily large problems, with computational cost linear in problem size. The performance of the learned linear solvers compared with both classical and optimized domain decomposition algorithms, for both structured- and unstructured-grit problems.
1 Introduction
Domain decomposition methods (DDMs) 2 3 are highly effective in solving the linear and nonlinear systems of equations that arise from the numerical approximation of solutions to partial differential equations (PDEs). While most effective on elliptic boundary-value problems, DDMs can also be applied to nonlinear problems, either using their nonlinear variants, successively solving linearizations. Time-dependent problems are normally solved by using a time stepping algorithm in the time domain for implicit methods, which require the solution of a spatial problem for each time step. Of these methods, Schwarz methods are particularly popular given their relative simplicity and ease of parallelization. The common theme break the global problem into subproblems, derived either by projection or by discretizing the same PDE over a physical subdomain, and to use solutions
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

on the subdomains as a preconditioner for the global discretization. Classical Schwarz methods generally make use of Dirichlet or Neumann boundary conditions for these subdomain problems, while Optimized Schwarz Methods (OSMs) aim to improve the convergence of the algorithm by using more general interface conditions 4 Notably, [5] demonstrates that optimal, but non-local, interface conditions exist for more general decompositions.
Much of the OSM literature considers only one-level additive Schwarz methods, although multilevel extensions do exist. For one-level methods (i.e., domain decomposition approaches without "coarse grid"), restricted additive Schwarz (RAS) approaches [6] are arguably the most common; optimized restricted additive Schwarz (ORAS) methods are considered in [Z]. The OSM idea has also been extended to asynchronous Schwarz methods (8) where the computations on each subdomain are done using the newest information available in parallel computing environment without synchronizing the solves on each subdomain.
With a recent focus on machine learning (ML) techniques for solving PDE systems 9 T0], there is also effort to apply learning-based methods improve the performance of iterative solvers for PDEs, including DDM and algebraic multigrid (AMG) methods, Within AMG methods, ML techniques have been applied learning interpolation operators ( and to coarse-grid selection in reduction- based AMG 13]. Of particular note here is the loss function employed in I1 12]. where they use unsupervised learning to train graph neural network minimize the Frobenius norm of the error- proropropagtttio operator of their iterative method. Within DDM, significant effort has been invested in combining ML techniques with DDM, as in 14]. where two main families of approaches are given: 1) using MI. within classical DDM methods to improve convergence, and 2) using deep neural networks, such as Physics Informed Neural Networks (PINNs) 9. as a discretization module and solver for DDM problems. In [15] fully connected neural network is used to predict the geometric locations of constraints for coarse space enhancements in an adaptive Finite Element Tearing and Interconnecting- -Dal Primal (FETI-DP) method. Using the continuous formulation of DDM, the so-called D3M [6] uses a variational deep learning solver, implementing local neural networks on physical subdomains in a parallel fashion. Likewise, Deep-DDM [1 utilizes PINNs to discretize
and solve DDM problems, with coarse space corrections 18] being used to improve scalability.
In this paper, we advance DDM-based solvers by developing a framework for learning optimized Schwarz preconditioners. key aspect of this reconsidering the loss function to use a more effective
relaxation of the ideal objective function than Frobenius norm minimization [1] 12 Moreover, the approach introduced here offers an opportunity to reconsider existing limitations of optimized Schwarz methods, where optimal parameter choice is based on Fourier analysis and requires a highly regular subdomain structure, such as in the classical cases of square domains split into two equal subdomains or into regular grids of square subdomains. Our framework learns the optimized Schwarz parameters via training on small problem sizes, in a way that generalizes effectively to large problems, and in way that allows us to consider both irregular domains and unstructured grids, with no extraordinary limitations on subdomain shape. Furthermore, the evaluation time of our algorithm scales linearly with problem size. This allows significant freedom in defining optimized Schwarz methods, in comparison to classical approaches, allowing us to explore the potential benefits of optimized Schwarz over classical (restricted) additive methods unstructured grids for the first time.
2 Background
Let Z C R2 be an open set, and consider the positive-defintt Helmholtz problem
Lu = (1) A)u = in S2,
(1)
with inhomogeneous Dirichlet conditions imposed the boundary Un In (T. the parameter 1) 0 represents a shift in the Helmholtz problem. In the numerical results below, we consider both finite-difference discretizations of ( on regular grids, well as piecewise linear finite-element (FE) discretizations on arbitrary triangulations. In both cases, we denote the set of degrees of freedom as D, and note that these in one-to-one correspondence with the nodes of the underlying mesh. Consider a decomposition of D into non-overlapping subdomains D i E {1, 2,,,, S} such that each node is contained within exactly one subdomain D.. yielding UD! = D. In this subdomain notation, the superscript denotes the amount of overlap in the subdomains, which is zero for the non-overlapping subdomains that we first consider. Let R. be the restriction operator onto the set of degrees of freedom (DoFs) in D.. and let (RQ)² T be the corresponding extension operator from DM

---

D2
D
Optimized Schwarz DoFs
Figure 1: Two subdomains with overlap 5 = on 58-node unstructured grid. The blue and orange shading denotes the original non-overlapping partitions (Df and D9), while the blue and orange dashed outlines show the overlapping subdomains (D) and D2). Nodes belonging to only one subdomain are marked with a solid circle, while those white belong to both subdomains. The connections in the boundary matrices L1 and L2 are denoted by the edges shaded in blue (for L1) and orange (for L2).
into set D. Then, FE discretization of the Helmholtz problem leads to linear system of the form Ax = where is the global stiffness matrix and A R (R)) T is the subdomain stiffness matrix for D. We note that alternate definitions to the Galerkin projection for Aه are possible, and commonly considered in optimized Schwarz settings (as noted below).
In the case of restricted additive Schwarz (RAS) 16 the subdomains are extended to allow for overlap: nodes near the "boundary" of their subdomain are potentially included in two or more subdomains. We denote the amount of overlap by e N, defining subdomains D recursively, by D = DØ-1 U {j dxj + and D for 0. The conventional RAS preconditioner
defined as
S MRAS
((E))
where R2: is the standard restriction operator subdomain 2 fe. is modified restriction operator from D to the DoFs in D that takes nonzero values only for DoFs in D and A = (R.) AR Figure [ shows example unstructured grid with two subdomains and overlap 8 = 1.
In an optimized Schwarz setting, we modify the subdomain systems, A². Rather than using Galerkin projection onto D. rediscretize [[] over the subdomain of Q corresponding these DoFs, imposing Robin-type boundary condition on the boundary of the subdomain. We define this matrix to be A A + L,, where A, is the term resulting from discretization of (I) with Neumann boundary conditions, and L3 is additional the term resulting from the Robin-type -tpp condition, as in:
Dirichlet: u (.)),,,
Neumann: ñ Vu =gx(x),
Robin: au + π Vu 9R 9R(x),
(3)
where ñř is the outward normal to the edge on the boundary and denotes inhomogeneous data.
The matrix L2 has the same dimensions as A,, so that A is well-defined. However, it has significantly different sparsity pattern, with nonzero entries only in rows/columns corresponding to nodes on the boundary of subdomain D. In practice we identify a cycle or path in the graph corresponding to A, with the property that every node in the cycle is on the boundary of D but not the boundary of D (the discretized domain), and then restrict the nonzeros in L; to the entries corresponding to the edges in this cycle/path (including self-edges, corresponding to entries on the diagonal of T.:).
Using this notation, the ORAS preconditioner can be written as
MORAS = Σ R
A
R
(4)
Our aim is to learn the values in the matrices Tup, aiming to outperform the classical choices of these values predicted by Fourier analysis, in the cases where those values are known, and to learn

---

suitable values for cases, such as finite-element discretizations on unstructured grids, where no known optimized Schwarz parameters exist. We optimize the values for the case of stationary (Richardson) iteration, but evaluate the performance of the resulting methods both as stationary iterations and as preconditioners for FGMRES.
Graph Neural Networks (GNNs): Applying learning techniques to graph structured data necessi- tates stepping beyond multilayer perceptron (MLP) and conventional convolutional neural networks (CNN) to a type of network that leverages the underlying graph nature of the problem, namely graph convolutional neural networks (GCNNs). GCNNs are typically divided into two categories: spectral and spatial [19]. Spectral GCNNs, first introduced by Bruna et al. [20), consider graph Laplacian cigenbasis and define convolution as diagonal operators. As such, spectral GCNN methods suffer from time complexity problems due to the necessity for the eigendecomposition of the Laplacian matrix. Nevertheless, in follow-up work 21 [2]] remedies have been proposed to mitigate this. Unlike spectral methods, spatial GCNNs consider local propagation of information in graph domains as a convolution graph. One popular framework is the message passing neural network (MPNN) (23). which is based on sharing information among neighbor nodes in each round of convolution pass. This has been generalized [24] by introducing a configurable graph network block structure consisting of node and edge convolution modules and global attribute. In an effort alleviate computational complexity of GCNNs [25], topology adaptive graph convolution networks (TAGCN) can be con- structed by defining learnable filters. This is not only computationally simpler, but also allows for adapting to the topology of the graphs when scanning them for convolution.
3 Method
3.1 Optimization problem and the loss function
Throughout this paper, we use |- to denote the e norm of matrix or vector, | for the Frobenius norm of A, and p(A) as the spectral radius of A. The optimization problem that we seek to solve is to find optimal values for the entries in the matrices L, constrained by given sparsity patterns, to minimize p(T'), where T =1 MORAS is the error-propagation operator for the stationary iteration corresponding MORAS defined in (4). The spectral radius p(T) corresponds to the asymptotic convergence factor of the stationary iteration, giving a bound on the asymptotic convergence of the method. Formally defined in terms of the extremal eigenvalue of T (since T is not symmetric), direct minimization of p(T) is difficult since backpropagation of an eigendecomposition is numeri- cally unstable (26). To overcome this, Luz et al. [12] propose to relax the minimization of p(T) (for a similar AMG preconditioner) to minimizing the Frobenius norm of T, [|]] In our experiments, however, we find that this is insufficient, leading preconditioners that do not scale. One reason is that while the Frobenius norm is an upper bound on the spectral radius, it does not appear to be a suitably "tight" bound for use in this context (see Section 4.2 and Figure O). Instead, we use relaxation inspired by Gelfand's 's formula, that ((T)
lim
and the common bound that
K K
*
|T x01| |
P(T) S | |
sup 2y0
=
sup (1T x=I) *
(5)
-=||>||-1
for some finite K E N. This results in the optimization problem
min 1,111122 ............ x: 11:1ax
sup
|T x||.
(6)
|=1
sparsity of L.
3.2 Numerical evaluation of the loss function
We denote the action of evaluating the GNN by f(()) (where f represents the network parameters), and consider a discretized problem with DoF set D, of size n. The set D can be decomposed into subdomains either by using fixed geometric choices of the subdomain (e.g., for Umite-difference discretizations), using the METIS graph partitioner (27), or a k-means-basse based clustering algorithm (best known as Lloyd's algorithm which has O(n) time complexity) (28 29]. For unstructured problems, we use a k-means-har based algorithm, decomposing to subdomains D. for = 1,2, ......... S
with overlap ò; see Supplementary Materials for details. The GNN then takes D and its decomposition

---

as ipputs, as well as sparsity constraints on the matrices L, for i = 1,2,..... S, and outputs values for these matrices:
1 I ...... L + f()
Using the learned subdomain interface matrices, we then obtain the modified MLORAS (Machine
Learning Optimized Restricted Additive Schwarz) operator, M ORAS. simply using Aة A, (,(0)
in We denote the associated eπror propagation operator by T = =I- - (8) MRSS
While Gelfand's formula and the associated upper bound in (5) are valid in any norm, is natural
to consider them with respect to the e norm this setting. However, this raises the same issue as encountered in 122 that generally requires eigendecomposition compute the norm. To avoid
this, we use a stochastic sampling of (T(0)) [|]
generated by the sample set X E RR*m for some
m E N. given as
X [1 T2,..... ml I V 2j ~ R" uniformly, ||aا|| = 1.
(8)
Here, we randomly select m points uniformly on unit sphere in R", which can be done using the method in [30]] We then define
K (TT))) K 01 |.
(T(0) K
y-(1)
x2
taking each column of X as the initial guess to the solution of the homogeneous problem Az and taking K steps of the stationary algorithm to generate (T(()) æij. Since we normalize each column
K of to have |ر 1, the value of (T(") zj serves as a lower bound for (IT1))* Thus.
taking the maximum of the values in Y provides a practical loss function that we use below, defining
( (a)
(10)
We note similarities between this loss function and that used in 31.. but that we are able to use the maximum of (0) (giving better approximation to the norm) our context instead of averaging.
Ultimately, the cost of our algorithm depends strongly on the chosen values of m and K. For sufficiently large values of m, we now show that the maximum value in y(0) is an arbitrarily good
approximation (I)() K 11 (in the statistical sense).
Theorem 1. For any nonzero matrix T, E and 8 < there exist M, K.K. N such that for any > M, if one chooses m points, Tj. uniformly at random from {or E IR"", ||||| = 1), then Y = [| x1
|| xell....., T xm satisfies
(ITTT) - max ma) >1
(11)
Proof. According to Gelfand's theorem, there exists I.
E
N such that Vl I,
p(T) sup a=||-||=1
|x 5 Take any K L and let ē
Since R is finite-
2 T |
dimensional, there exists an x" E", ||**|| such that sup
x"1 || Denote the | x = |
-=||>/||1 volume of the surface of the n-dimensional sphere of unit radius around the origin IR" by Cior. and the volume of the region this sphere within radius EK of x by Cz.K.
Given 8 1, let M satisfy M
log(5)
Then,
log
Ce.K
Cias
P P(|=" zill 2K for all i) (1
Ce.K Cior
M
< 6.
(12)
Thus, with probability of least 8, we expect at least one point from selection of m M points uniformly distributed the sphere of radius unit radius to be in C.,K. Let that point be Zr,
5

---

giving
| x²1 - -| 2,
< (IT* 2²1 -11T xr
(13)
< ( [| (x" - -x+)]l) *
(14)
c < | * || ---- I = |T
(15)
2
using Lemma (3 (see Supplementary Materials) and the reverse triangle inequality. Finally, by the triangle inequality, with probability of least 1 8, we have:
((T)) - max k s ((T)
17 x k + sup -=||-||=1
x max(Y) k S.. (16)
sup
-=||=||=1
Remark. According to [32, since the optimal interface values are cowringcommputable. if the depth of the GNN is at least the diameter of the graph, and a TAGConv layer followed by a feature encoder Turing-complett the optimal interface values can be learned. In our setting, for problem on structured grid of size N x N with two identical rectangular subdomains, this implies that the GNN will be able to learn the optimal interface values given, if and only the GNN has depth at least 2N, has deep enough feature encoders, and the width of the layers is unbounded.
Remark. Theorem 1 guarantees convergence of the loss function to the spectral radius, in the limits of many samples and many stationary iterations. To the best of our knowledge, such a guarantee is not known for the previous loss functions used in the area [72 Moreover, there are substantial improvements the numerical results using the new loss function comparison to that of (12), as shown in F
Theorem Assuming bounded subdomain size, the time complexity to evaluate the optimal Schwarz parameters using our method is (n)), where n is the number of nodes in the grid.
Proof. Given bounded Lloyd subdomain size and fixed number of Lloyd aggregation cycles, sub- domain generation has (nn) time complexity (28) (see the Supplementary Material). To evaluate each TAGConv layer, computes y Σi= Gevz + bl,, where L is the number of node features, G. E R is the graph filter, b is learnable bias, and E IR are the node features. Moreover, the graph filter is polynomial in the adjacency matrix M of the graph: G. = Σド= 9t.j It. where is constant and gt.j are the coefficients of filter polynomial. Since the graph has bounded node degrees, it implies that M is sparse and the action of M3 has O(n) cost, and therefore, the full TAGCony evaluation also has ((n) cost. Moreover, the cost of edge feature and the feature networks are ((n). resulting in overall (n)) cost.
3.3 Training Details
We use graph neural network based on four TAGConv layers and ResNet node feature encoders consisting of cight blocks after each TAGConv layer, see Supplementary Materials for more details on the structure of the GNN, and Section 4.3 for an ablation study. The training set in our study consists of 1000 unstructured grids with piecewise linear finite elements and with grids ranging from 90-850 nodes (and an average of 310 nodes). The grids are generated by choosing either a regular grid (randomly selected 60% of the time) or randomly generated convex polygon; pygmsh [33 is used to generate the mesh on the polygon interior. A sample of the training grids is shown Figure We train the GNN for four epochs with mini batch size of 25 using the ADAM optimizer [34] with a fixed learning rate of 10-4 For the numerical evaluation of the loss function (10) we use K iterations and 500 samples. The code_ was implemented using Py PyTorch Geometric (5)
PyAMG [36), and NetworkX (37] All training was performed on an 8-core i9 Macbook Pro using CPU only.
The first is "Brute Force" network, which is a single-layer neural network without an activation function, trained only on a structured
Moreover, we train two special networks for use in Section 4.1
All code and data for this paper is https: http:///hhhaaaaccccc////////////
(MIT licensed).

---

Figure 2: Example grids from the training set.
Jacobi
AS
RAS
ORAS 000
ORAS 002
MLORAS (Ours)
Brute Force
10-1
10-5
Rosidual
10-9
2 Stationary Iteration
FGMRES Iteration
Figure 3: Results for the Helmholtz problem on 10 x 10 structured grid (left) with two identical subdomains, with convergence plots for the methods used as preconditioners for FGMRES (center) and as stationary algorithms (right).
10 10 grid with two identical subdomains using the ADAM optimizer. The purpose of training this to obtain the optimal interface values as benchmark for comparison, including against our method. The second is the same network used for method, MLORAS, but overtrained only the problem to understand the learning capabilities of the method and choice GNN architecture.
Section 4.1
4 Results
4.1 Two-subdomain structured grids
We first consider rectangular structured grids with two subdomains. Although restrictive, these problems allow us to directly compare existing OSM parameters, which are only available for structured grids and exactly two subdomains. We follow 7 and consider Helmholtz problem (T) on the unit square with 7 1 and homogeneous Dirichlet boundary conditions. We discretize the problem on an x N rectangular grid with L/(N + 1) grid spacing, using the standard five-point finite difference stencil.
Figure 3 shows the results of different methods on a 10 10 structured grid with identical subdomains. We see that the overtrained MLORAS network learns interface parameters that out- perform the Restrictive Additive Schwarz (RAS) 6 method: more significantly, MLORAS also outperforms the existing Optimized-RAa RAS (ORAS) methods that were analytically derived for this specific problem (zeroth-order 000 and second-order 002 from 7 Moreover, these results, the MLORAS network almost reaches the performance of the "Brute Force" network (obtained by directly optimizing all interface values for this single structured grid), indicating that using a GNN to encode the optimal interface values does not significantly restrict the search space.
7

---

finite Brute Force
MLORAS (Ouzs)
ORAS 002
2
250
250
500
Figure 4: Interface values for the x 10 structured grid with two identical subdomains. From left
right: brute force optimization of the interface values, MLORAS, and second-order ORAS (002) [
1018 nodes
200 nodes
Resicual porm
10-2
Jacobi
10
RAS
ML.ORAS (Ours)
10-10
30
100
200
250
300
350
100
150
250
300
350
FGMRES Iteration
Ineration
Figure 5: Example convergence smaller (left) and larger (right) unstructured grids.
To understand the performance of the methods in more detail, Figure 4 plots the interface values output by the brute force optimization, the MLORAS network, and the 002 ORAS algorithm. This shows that the MLORRS network is choosing interface values very close those directly optimized by the brute force method, unlike those selected by the O02 method.
4.2 Unstructured grids
To evaluate the performance of the MLORAS network unstructured grids, consider 16 unstructured triangular grids in 2D with sizes ranging from about 90 to 40 000 nodes. These grids are defined convex subsets of (0, 1) (0, 1); solve the Helmholtz problem (I) these domains with = and homogeneous Dirichlet boundary conditions, discretized with piecewise- elements.
Example convergence plots are shown in Figure where our method (MLORAS) compared to RAS (Restricted Additive Schwarz (6]), AS (Additive Schwarz), and Jacobi methods as preconditioner for FGMRES. Here we see that the MLORAS network able to learn optimized interface parameters scale to problems that much
for unstructured grids that outperform RAS, and that MLORAS larger than those in the training set, which are all below n = 1000 nodes. Importantly, MLORAS retains advantage over RAS even as the grid size increases.
Figure 6 shows the performance of three different methods across unstructured test grids. The three methods are: (1) "F-norm optimized"" the same network as MLORAS. but trained to directly optimize the Frobenius norm of the error propagation matrix, || F. (2) the RAS 6 method, and (3) the ML method. We do not show ORAS results here, as it cannot be applied to unstructured
grids. Figure reveals two important facts. First, MILORAS consistently outperforms RAS over the entire testing set, both an FGMRES preconditioner and stationary algorithm, and this

---

F-noru oplimized
RAS
MLORAS (Ous)
PGMItSS steps
Statlonary erroz
0.75
0.50
Probenins
0.25
24
0.00
10a
10a
10
10a
10a
10
10a
10a
10
Grid size
Grid size
Grid size
Figure 6: Convergence all unstructured grids in the testing set. Left: the number of preconditioned FGMRES steps required to solve the problem to within relative eπror of 10-12 12 Center: error reduction by the stationary iteration after 10 iterations. Right: Frobenius norm for each method on each test problem.
graph cony layers
Residual block length
33.0
Avg PGMRES steps
0.15
Avg stationary eITOY
32.5
0.14
32.0
0.13
31.5
0.12
31.0
0.11
0
8
2
6
8
10
12
14
16
Residual block length
Number of graph convolution layers
Figure 7: Ablation study results. The left panel varies the residual block length while keeping the number of TAGConv layers fixed at 4. The right panel varies the number of TAGConv layers while keeping the residual block length fixed 8. The solid blue lines (left axis) show the average number of FGMRES steps needed to reduce the relative error below 10 10-12, while the dashed green lines (right axis) show the average stationary algorithm error reduction after 10 iterations. The red circles mark the results for the network with the best performance (residual block length of 8 and 4 TAGCony layers), which was used for all other studies in this paper. Error bars show one standard error of the mean.
advantage maintained even for large grids. Second, we see that the Frobenius norm is indeed worse choice for the loss function than our new loss in (6). We see this from the fact that MLORAS has worse Frobenius In addition, when we explicitly optimize the Frobenius norm (the "F-norm optimized" method), we see that we do obtain the lowest Frobenius norm, but this translates to the worst convergence rate.
than either of the other
methods, but it has the best convergence rate.
4.3 Ablation study
To understand the impact of the network structure on performance, we conduct an ablation study by varying the ResNet block length and the number of TAGConv layers. In each case, we train a network on five different training sets, each consisting of 1000 grids generated as described in Section The trained networks are tested on 50 unstructured grids, each with 2400 to 2600 nodes. For each architecture, the mean performance is computed, together with error bars from the fivefold repetition, with results shown in Figure 7 We see that a higher residual block length is always better, but that optimal number of TAGConv layers 4, and using more layers actually decreases performance. This phenomenon has been observed in other contexts (e.g., (38) and (3) and can be attributed to over-smoothing GNNs. The use of residual blocks in our architecture is, thus, important to allow us to increase network depth without needing to increase the number of GNN layers.

---

Conclusion
We propose an unsupervised learning method based on Graph Convolutional Neural Networks (GC- NNs) for extending optimized restricted additive Schwarz (ORAS) methods to multiple subdomains and unstructured grid cases. Our method is trained with novel loss function, stochastically minimiz- ing the spectral radius of the error propagation operator obtained using the learned interface values. The time complexity evaluating the loss function, as well as obtaining the interface values using our neural network are all linear in problem size. Moreover, the proposed method ís able to outperform ORAS, both as a stationary algorithm and preconditioner for FGMRES on structured grids with two subdomains, considered in the conventional ORAS literature. On more general cases, such as unstructured grids with arbitrary subdomains of bounded size, our method outperforms RAS consistently, both stationary algorithm and preconditioner for FGMRES. The main limitations of the current work are that the method was studied for two PDE cases, namely the Helmholtz problem in the main paper and the non-uniform Poisson problem in Appendix D We also defer the study of nonlinear or time-dependent PDEs to future work.
Acknowledgement
The work of S.P.M. was partially supported by an NSERC Discovery Grant. The authors thank the referees for their insightful comments. The authors have no competing interests to declare.
10

---

References
[1] Andrea Toselli and Olof Widlund. Domain decompasition wethods-algoritmms and theory, volume 34 of Springer Series in Computational Mathematics. Springer-Veri Berlin, 2005.
[2] Alfio Quarteroni and Alberto Valli. Domain decomposition methods for partial differential equations. Numerical Mathematics and Scientific Computation. The Clarendon Press, Oxford University Press, New York, 1999. Oxford Science Publications.
[3] Victorita Dolean, Pierre Jolivet, and Frédéric Nataf. An introduction to domain decomposition methods. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2015.
[4] M.J. Gander, L. Halpern, and F. Nataf. Optimized Schwarz methods. In Proceedings of the 12th International Conference on Domain Decomposition, pages 15-27. ddm.org. 2000.
[5] Martin J. Gander and Felix Kwok. Optimal interface conditions for an arbitrary decomposition into subdomains. In Domain decomposition methods in science and engineering XIX, volume 78 of Lect. Notes Comput. Sci. Eng., pages 101-108. Springer, Heidelberg, 2011.
[6] Xiao-Chuan Cai and Marcus Sarkis. A restricted additive Schwarz preconditioner for general sparse linear systems. Siam Journal on Scientific Computing. 21(2):792-797, 1999.
[7] Amik St-Cyr, Martin J Gander, and Stephen J Thomas. Optimized multiplicative, additive, and restricted additive Schwarz precoodiiincnnnniiiiinngg SIAM Journal on Scientific Computing, (9(6)24002 2425, 2007.
[8] Frédéric Magoulès, Daniel Szyld. and Cédric Venet. Asynchronous optimized Schwarz methods with and without overlap. Numer: Math., 137(1):199-22)) 227, 2017.
[9] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics- informed neural networks:
A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
[10] Zongyi Li, Nikola Kovachki, Kamyar Arizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differen- tial equations. arXiv preprint ar arXiv:2010.0885 2020.
[11] Daniel Greenfeld. Meirav Galun, Ronen Basri, Irad Yavnch, and Ron Kimmel. Learning to optimize multigrid PDE solvers. In International Conference on Machine Learning, pages 2423. PMLR, 2019.
[12] Ilay Luz, Meirav Galun, Haggai Maron, Ronen Basri, and Irad Yavneh. Learning algebraic multigrid using graph neural networks. In International Conference on Machine Learning, pages 6489-6499, PMLR, 2020.
[13] Ali Taghibakhshi, Scott MacLachlan, Luke Olson, and Matthew West. Optimization based algebraic multigrid coarsening using reinforcement learning. Advances in Neural Information Processing Systems, 34, 2021.
[14] Alexander Heinlein, Axel Klawonn, Martin Lanser, and Janine Weber. Combining machine learning and domain decomposition methods for the solution of partial differential equations- review. GAMM-Miiteillngenn
2000:202200001, 2021.
[15) Alexander Heinlein, Axel Klawonn, Martin Lanser, and Janine Weber. Machine learning in adaptive domain decomposition methods methods--preiecting the geometric location of constraints. SIAM Journal on Scientific Computing, 41 (166)
2019.
[16] Ke Li, Kejun Tang, Tianfan Wu, and Qifeng Liao. D3M: A deep domain decomposition method for partial differential equations. IEEE Access, 8:5283-5294, 2019.
[17] Wuyang Li, Xueshuang Xiang, and Yingxiang Xu. Deep domain decomposition method: Elliptic problems. In Mathematical and Scientific Machine Learning, pages 269-286. PMLR, 2020.
11

---

[18] Valentin Mercier, Serge Gratton, and Pierre Boudier. A coarse space acceleration of deep-DDM, arXiv preprint arXiv:2112.03732, 2021.
[19] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and Yu Philip. A comprehensive survey graph neural networks. IEEE Transactions Neural Networks and Learning Systems, 2020.
[20] Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun. Spectral networks and locally connected networks on graphs. arXiv preprint arXiv: 2013.
[21] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst. Convolutional neural networks on graphs with fast localized spectral filtering. arXiv preprint arXi:: 09375, 2016.
[22] Thomas N Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. arXiv preprint arXiv: ar.iv:6609.0...00 wetp://09900000 2016.
[23] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl. Neural message passing for quantum chemistry. In International Conference on Machine Learning, pages 1263-1272, PMLR, 2017.
[24] Peter W Battaglia, Jessica B Hamrick, Victor Bapst, Alvaro Sanchez-Gonzalez, Vinicius Zambaldi. Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, al. Relational inductive biases, deep learning, and graph networks. arXiv preprint arXiv: 1806.0166. 2018.
[25] Jian Du, Shanghang Zhang, Guanhang Wu, José MF Moura, and Soummya Kar. Topology adaptive graph convolutional networks. arXiv preprint arXiv: 1710.10700 2017.
[26] Wei Wang, Zheng Dang, Yinlin Hu, Pascal Fua, and Mathieu Salzmann. Backpropagation friendly eigendecomposition. Advances in Neural Information Processing Systems, 32, 2019.
[27] George Karypis and Vipin Kumar. fast and high quality multilevel scheme for partitioning irregular graphs. SIAM Journal scientific Computing, 20(1):359-33) -929,999911111 1998.
[28] William N Bell. Algebraic multigrid for discrete differential forms. PhD thesis, University of Illinois at Urbana- Urbana-Champaa Champaign, 2008
[29] Stuart Lloyd. Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2):129-137) 1982.
[30] George EP Box. A note on the generation of random normal deviates. Ann. Math. Statist., 29:610-611, 1958.
[31] Alexandr Katrutsa, Talgat Daulbaev, and Ivan Oseledets. Black-box learning of multigrid parameters. Journal of Computational and Applied Mathematics, 368: 12524, 2020.
[32] Andreas Loukas. What graph neural networks cannot learn: depth width. arXiv preprint 1rXiv:190..03199. 2019.
[33] Nico Schlömer. pygmsh: A Python frontend for Gmsh, 2021. (GPL-3.0 License).
[34] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6688, 2014.
[35] Matthias Fey and Jan E. Lenssen. Fast graph representation learning with PyTorch Geometric. In ICLR Workshop on Representation Learning on Graphs and Manifolds, 2019. (code is MIT licensed).
[36] Nathan Bell, Luke N. Olson, and Jacob Schroder. PyAMG: Algebraic multigrid solvers in python. Journal of Open Source Software, 7172)4114,, 2022. (code is MIT licensed).
[37] Aric Hagberg, Pieter Swart, and Daniel S Chult. Exploring network structure, dynamics, and function using networkx. Technical report, Los Alamos National Lab. (LANII), Los Alamos, NM (United States), 2008. (code is BSD licensed).
12

---

[38] Lingxiao Zhao and Leman Akoglu. Pairnorm: Tackling oversmoothing in GNNs. arXiv preprint arXiv:1909.122 12223, 2019.
[39] Deli Chen, Yankai Lin, Wei Li, Peng Li, Jie Zhou, and Xu Sun. Measuring 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-343880 3488, 2020.
[40] Thomas H Cormen, Charles Leiserson, Ronald Rivest, and Clifford Stein. Introduction to algorithmss (Second edition). MIT Press, 2001.
13

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the pa- per's contributions and scope? [Yes] See the abstract and the last paragraph of the introduction.
(b) Did you describe the limitations of your work? [Yes] See Section 5
(c) Did you discuss any potential negative societal impacts of your work? [N/A] This paper only concerns learning interface values for domain decomposition methods and thus has no direct negative societal impacts.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See the statements of Theorems [ and 2
(b) Did you include complete proofs of all theoretical results? [Yes] See the proofs of Thhorrems[aan 2
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material as URL)? [Yes] URL for the repository including all the code and data is provided as a footnote in the first paragraph of Section [X
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] See Section 3.3
(c) Did you report error bars (c.g.. with respect the random seed after running experi- ments multiple times)? [Yes] See Figure 7
(d) Did you include the total amount of compute and the type of resources used (e.g.. type of GPUs, internal cluster, cloud provider)? [Yes] See Section 3.3
4. If you using existing assets (e.g., code, data, models) curating/releusing
assets... (a) If your work uses existing assets, did you cite the creators? [Yes] We cited references
155
36
57
the first paragraph of Section 33
(b) Did you mention the license of the assets? [Yes] We mentioned the licenses for
35
36
37
in the references.
(c) Did you include any assets either in the supplemental material URL? [Yes] Provided in Checklist Question 3(a).
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A] No such data was used.
(e) Did you discuss whether the data you are using/curating contains personally identifiable information offensive content? [N/A] No such data was used.
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A] No participants.
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A] No participants.
(c) Did you include the estimated hourly wage paid participants and the total amount spent on participant compensation? [N/A] No participants.
14


Reviews:
[{"review_id": "maQVmc61Wxn", "reviewer": "Reviewer_7Eu3", "summary": "Learning Interface Conditions in Domain Decomposition Solvers\n\nIn this work, the authors use machine learning tools to enhance the domain decomposition solvers. It uses a graph convolutional network to learn the preconditioner. Specifically, the GNN takes the domain D and its decomposition as inputs, and outputs the boundary matrices L.\n\nNumerically, the paper mainly considers a Helmholtz equation with constant coefficient and inhomogeneous boundary conditions. The experiments show the proposed method has a better convergence rate compared to the state-of-the-art solvers. ", "questions": "There could be more definitions or an example of the matrices L. What is its shape and what do the entries represent? \n", "limitations": "The authors discuss that the study is limited to the Helmholtz equation. \nBesides, in the reviewer's opinion, the decomposition method as an iterative method is limited to elliptic PDEs. It will be interesting to see how the decomposition work for time-dependent equation and discontinuous problem. Do we need to change the decomposition over time?", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: \n- The purposed method can be used to solve problems on larger domains (>10000 nodes). \n- The method is theoretical guarantees inherited from the Domain Decomposition solver.\n- And it has a better convergence rate than the numerical solvers.\n- Linear complexity\n\nWeakness: \n\n- In the reviewer’s opinion, the testing problem is somehow too simple. The paper mentioned non-linear equations, but only a linear Helmholtz equation is considered. It will be also interesting to study the time-dependent equations. Can the iterative solver be used to solve hyperbolic PDEs?\n- The new loss function appears somewhat increment to the reviewer.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "Ut3emhkB2_J", "reviewer": "Reviewer_4BC2", "summary": "This paper presents a learning-based approach to find the optimal boundary values for domain-decomposed PDE solvers. The major contribution is the application of TAGConv to solving PDEs, and a new numerical method to compute the loss defined on the spectral radius of the difference matrix.", "questions": "None. The paper is well-written and comprehended. ", "limitations": "The author only tested their method on Helmholtz equation, which is the major limitation. But the authors also have mentioned this limitation and further validation remains for future works.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The method looks solid, and is well-validated. The results look promising, and the application of TAGConv to domain-decomposition is novel. \n\nHowever, the paper doesn't have much contribution to machine learning itself, and thus I'm not sure it fits the scope of NeurIPS. I would recommend the author resubmit this work to a journal for numerical PDE solvers.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "OF2_BZVPG8W", "reviewer": "Reviewer_H1jk", "summary": "Domain decomposition solvers are effective for large-scale PDE problems and optimized Schwarz domain decomposition methods are commonly used for structured-grid problems. In this work, the authors proposed to learn the parameters in Schwarz preconditioner that generalize it to unstructured grid and enable faster convergence. The main difference compared to the prior work, which works with a different preconditioner, is the proposed new loss function inspired by Gelfand's formula; and a practical approximation of the loss function via stochastic sampling.  Experiments on the 2D Helmholtz equation have been done to demonstrate the effectiveness of the approach.   ", "questions": "1. how are boundary nodes handled for the input graph to the GNN?\n2. it's a bit unclear to be for the input/output of the GNN and how is that done. In line 139, it is noted that the input $i$ the $D$ (total domain) and its decomposition, and the sparsity constraint $L_i$, and predict the matrices. How exactly the \"decomposition\", and \"sparsity constraints\" are encoded in the input graph? According to the supplement, the output is the \"learned interface edge features\" and that will give us the predicted matrices, is this understanding correct? ", "limitations": "There is only one experiment setting with the Helmholtz equation. It would be more appealing if we have another experiment with a different PDE and boundary condition. \n", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "originality - the proposed loss function is different from the prior works and the stochastic sampling technique used to approximate the loss function is new to be applied to this problem, to the best of my knowledge. \n\nquality/clarity - the figures are clear and the article is well-written. More detailed captions could be added to facilitate reading. It would be nice to add a list of contributions by the end of the introduction to make the contributions clear. \n\nsignificance - the proposed approach demonstrates its effectiveness in the 2D Helmholtz equation problem. The approach scale linearly with the problem made is potentially be used for larger-scale problems. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ["I don’t know"]}, "extra_text": {}}]

Author rebuttal:

Decision: Accept