{'fraction_non_alphanumeric': 0.03041399810008274, 'fraction_numerical': 0.01147611313700855, 'mean_word_length': 4.767674089784376, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coarse-Grid Selection Using Simulated Annealing

Multilevel techniques are efficient approaches for solving the large linear systems that arise from discretized partial differential equations and other problems. While geometric multigrid requires detailed knowledge about the underlying problem and its discretization, algebraic multigrid aims to be less intrusive, requiring less knowledge about the origin of the linear system. A key step in algebraic multigrid is the choice of the coarse/fine partitioning, aiming to balance the convergence of the iteration with its cost. In work by MacLachlan and Saad, a constrained combinatorial optimization problem is used to define the"best"coarse grid within the setting of a two-level reduction-based algebraic multigrid method and is shown to be NP-complete. Here, we develop a new coarsening algorithm based on simulated annealing to approximate solutions to this problem, which yields improved results over the greedy algorithm developed previously. We present numerical results for test problems on both structured and unstructured meshes, demonstrating the ability to exploit knowledge about the underlying grid structure if it is available.


INTRODUCTION
Multigrid and other multilevel methods are well-established algorithms for the solution of large linear systems of equations that arise in many areas of computational science and engineering. Multigrid methods arise from the observation that basic iterative methods, such as the (weighted) Jacobi and Gauss-Seidel iterations, effectively eliminate only part of the error in an approximation to the solution of the problem, and that the complementary space of errors can be corrected from a coarse representation of the problem. While geometric multigrid has been shown to be highly effective for problems where a hierarchy of discrete models can be built directly, many problems benefit from the use of algebraic multigrid (AMG) techniques, where graph-based algorithms and other heuristics are used to define the multigrid hierarchy directly from the matrix of the finest-grid problem.
First introduced in the early 1980's 2, 3 , AMG is now widely used and available in standard libraries, such as hypre 4, 5 , PETSc 6, 7 , Trilinos 8,9 , and PyAMG 10 . While there are many variants of AMG (as discussed below), the common features of AMG algorithms are the use of graph-based algorithms to construct a coarse grid from the given fine-grid discretization matrix (possibly with some additional information, such as geometric locations of the degrees of freedom for elasticity problems 11 ) followed by construction of an algebraic interpolation operator from the coarse grid to the fine grid. Both of these components have received significant attention in the literature, with an abundance of schemes for creating the coarse grid and for determining the entries in interpolation. In this paper, we focus on the coarse-grid correction process, adopting a combinatorial optimization viewpoint on the selection of coarse-grid points in the classical AMG setting.

arXiv:2105.13280v3 [math.NA] 7 Apr 2022
One of the primary differences between so-called classical (Ruge-Stüben) AMG and (smoothed) aggregation approaches is in how the coarse grid itself is defined 12 . In aggregation-based approaches, coarse grids are defined by clustering fine-grid points into aggregates (or subdomains), while, in classical AMG, a subset of the fine-grid degrees of freedom (points) is selected in order to form the coarse grid. In the original method 3,13 , this was accomplished using a greedy algorithm to construct a maximal independent set over the fine-grid degrees of freedom as a "first pass" at forming the coarse-grid set, which is then augmented by a "second pass" algorithm that ensures suitable interpolation can be defined from the final coarse grid. Many alternative strategies to forming the coarse grid have been proposed over the years, particularly for the parallel case, to overcome the inherent sequentiality of the original algorithm 5,[14][15][16][17] . Other approaches, for example based on redefining the notion of strength of connection [18][19][20][21] or use of compatible relaxation principles 22,23 , have also been proposed.
Much of this work has been motivated by the failure of the classical AMG heuristics to yield robust solvers for difficult classes of problems, such as the Helmholtz equation, convection-dominated flows, or coupled systems of PDEs. While much research has tried to generalize these heuristics to give sensible algorithms for broader classes of problems, another approach is to consider methods that abandon these heuristics in favour of algorithms that are directly tied to rigorous multigrid convergence theory. In recent years, two main classes of AMG algorithms have arisen in this direction. One class of methods arises in the aggregation setting, where well-chosen pairwise aggregation algorithms coupled with suitable relaxation can yield convergence guarantees for symmetric, diagonally dominant M-matrices 24,25 . The second class of methods are those based on reductionbased AMG principles, where there is a direct connection between properties of the fine-grid submatrix and the guaranteed convergence rate of the scheme 1, 26-28 . Reduction-based multigrid methods 29 are a generalization of cyclic reduction 30 , in which conditions on the coarse/fine partitioning are relaxed so that while the exact Schur complement may not be sparse, it can be accurately approximated by a sparse matrix. This notion is made rigorous by MacLachlan et al. 26 , who propose a theoretical framework to analyse convergence based on diagonal dominance of the fine-grid submatrix (in contrast to the fine-grid submatrix being diagonal in the classical setting of cyclic reduction), and extended to Chebyshev-style relaxation schemes by Gossler and Nabben 28 . While the theoretical guarantees offered by such methods are attractive, these approaches suffer from two key limitations. First, the classes of problems for which convergence rates are guaranteed are limited (requiring diagonal dominance and/or M-matrix structure), Second, the guaranteed convergence rates of stationary iterations are generally poor in comparison to methods based on classical heuristics. An active area of research (disjoint from the goals of this paper) is addressing the first approach 31 , while the second can be ameliorated by the use of Krylov subspace methods to accelerate stationary convergence. We note that there are other families of theoretical analysis of AMG algorithms 32 , including theory tailored to heuristic approaches, such as compatible relaxation 23 . Much of this theory, however, is of a confirmational and not a predictive nature, i.e., the convergence bounds rely on properties of the output of heuristics for choosing coarse grids and interpolation operators, and not on the inputs to these processes, such as the system matrix and algorithmic parameters. Thus, while it provides guidance in the development of heuristic methods, it does not provide the concrete convergence bounds offered by pairwise-aggregation or AMGr methods.
While both methods 24,25 provide practical algorithms to generate coarse grids that satisfy their convergence guarantees, this is a notable omission in the work of MacLachlan et al. 26 and much of the work on AMGr. MacLachlan and Saad 1,27 identify that the choice of optimal coarse grids for AMGr can be quantified as the solution of an NP-complete integer linear programming problem. They use this formulation to guide construction of a family of greedy coarsening algorithms that aim to maximize the size of the fine-grid set subject to maintaining a fixed level of diagonal dominance in the fine-grid submatrix (and, consequently, the resulting convergence bound on the AMG method). In this paper, we follow the theory and framework for AMGr 1, 26 , but aim to overcome some limitations of the underlying greedy coarsening algorithm. In particular, we develop a coarsening algorithm based on simulated annealing to partition the coarse and fine points. Numerical results show that this algorithm achieves smaller coarse grids (than the greedy approach) that satisfy the same diagonal-dominance criterion. Hence, these grids lead to moderately more efficient (two-level) algorithms, by reducing the size of the coarse-grid problem without changing the convergence bound guaranteed by the theory. We emphasize that this work is presented more as a proof-of-concept than as a practical coarsening algorithm, due to the high cost of simulated annealing. In related work, we have shown that the "online" cost of the approach proposed here can be traded for a significant "offline" cost using a reinforcement learning approach 33 . However, a necessary next step in the work proposed here is in investigating more cost-effective alternatives to these methods.
This paper is arranged as follows. In Section 2, AMG coarsening and the existing greedy coarsening algorithm for AMGr are discussed. The new coarsening algorithm based on simulated annealing is outlined in Section 3. Numerical experiments are given in Section 4, demonstrating performance of this approach on both isotropic and anisotropic problems, discretized on both structured and unstructured meshes. Concluding remarks and potential future research are discussed in Section 5.

AMG COARSENING
Geometric multigrid is known to be highly effective for many problems discretized on structured meshes. However, it is naturally more difficult to make effective coarse grids (and effective coarse-grid operators) for problems discretized on unstructured meshes. AMG was developed specifically to address this, automating the formation of coarse-grid matrices without any direct mesh information. The first coarsening approach in AMG, often called classical or Ruge-Stüben (RS) coarsening, was developed by Brandt, McCormick, and Ruge 2, 3 and later by Ruge and Stüben 13 . We summarize this approach below, primarily to allow comparison with the reduction-based AMG method 26,29 and the greedy coarsening approach proposed by MacLachlan and Saad 1 that is the starting point for the research reported herein.
AMG algorithms make use of a two-stage approach, where a setup phase precedes the cycling in the solve phase. Common steps in AMG setup phases include recursively choosing a coarse grid and defining intergrid transfer operators. The details of these processes depend on the specifics of the AMG algorithm under consideration, with both point-based and aggregation-based approaches to determining the coarse grid, and a variety of interpolation schemes possible for both of these.

Classical coarsening
As in geometric multigrid, the coarse grid in AMG is selected so that errors not reduced by relaxation can be accurately approximated on coarse grids. In an effective scheme, these errors are interpolated accurately from coarse grids that have substantially fewer degrees of freedom than the next finer grid, thus significantly reducing the cost of solving the coarse-grid residual problem.
In classical AMG, for an × matrix , the index set Ω = {1, … , } is split into sets and , with Ω = ∪ and ∩ = ∅. Each degree of freedom (DoF) or point ∈ is denoted a -point and a point ∈ is denoted an -point. This splitting is constructed by considering the so-called strong connections in the graph of matrix between the fine-level variables. Given a threshold value, 0 < ≤ 1, the variable is said to strongly where is the coefficient of in the th equation. The set of points that strongly influence , denoted by , is defined as the set of points on which point strongly depends. The set of points that are strongly influenced by is denoted by . Two heuristics are followed in RS coarsening to select a coarse grid: H1 : For every -point, , every point ∈ should either be a coarse-grid point or should strongly depend on at least one point in that also strongly influences .
H2 : The set of -points should be a maximal subset of all points, where no -point strongly depends on another -point. In practice, strong enforcement of both of H1 and H2 is not always possible; the classical interpolation formula relies on H1 being strongly enforced, while H2 is used only to encourage the selection of small, sparse coarse grids.

Greedy coarsening and underlying optimization
While the classical coarsening algorithm has proven effective for many problems when coupled with suitable construction of the interpolation operator 13 , the process provides few guarantees in practice. Indeed, there is an inherent disconnect between the selection of any single coarse-grid point and the impact on the quality of the resulting interpolation operator. This has motivated coupled approaches to selecting coarse grids and defining interpolation, including reduction-based AMG, or AMGr.
Reduction-based multigrid was proposed by Ries et al. 29 , building on earlier work aiming to improve multigrid convergence for the standard finite-difference (FD) Poisson problem. The fundamental idea of reduction-based multigrid lies in defining the multigrid components to approximate those of cyclic-reduction algorithms 30 . In particular, one way to interpret cyclic reduction is as a two-level multigrid method with idealized relaxation, interpolation, and restriction operators. Suppose that the coarse/fine partitioning is already determined, and consider the reordering of the linear system = according to the partition, writing An exact algorithm for the solution of = in this partitioned form is given by This can be turned into an iterative method for solving = in the usual way, by replacing the right-hand side vector, , by the evolving residual, and computing updates to a current approximation, ( ) , giving 1.
In this form, this remains an exact algorithm: given any initial guess, (0) , we have the solution = (1) . A truly iterative method results from approximating the three instances of −1 in the above algorithm, namelŷ leading to the following iteration: 1. ( In multigrid terminology, the first step can be seen as a so-called -relaxation step, where we approximate −1 using some simple approximation, such as with weighted Jacobi or Gauss-Seidel. The second step is a coarse-grid correction step, where the residual after -relaxation is restricted by injection, and the coarse-grid correction, , is computed by solving a linear system with an approximation,̂ , of the true Schur complement, − −1 . The final two steps correspond to an interpolation of the coarse-grid correction, writing the interpolation operator = , for some approximation of the ideal interpolation operator, ≈ −1 . If the approximation in the first step is sufficiently poor that a significant residual is expected to remain at the -points after relaxation, it is also reasonable to augment the restriction in the second step, using either the transpose of interpolation or some better approximation to an ideal restriction operator than just injection. As written above, the algebraic form of reduction-based multigrid retains the key disadvantage of classical AMG -there is little connection between the algorithmic choices (of the coarse/fine partitioning, and the approximations of , of the Schur complement, and of the idealized interpolation and restriction operators) with the resulting convergence of the algorithm. To address these limitations, MacLachlan et al. 26 proposed a reduction-based AMG algorithm (AMGr) that attempts to connect convergence with properties of . In particular, they presented a two-level convergence theory with a convergence rate quantified by the approximation of ≈ , with the expectation that application of −1 to both a vector and to is computationally feasible. In what follows, we use the standard notation that matrices ⪰ when ⪰ for all vectors . Then the multigrid cycle with pre-relaxation sweeps, coarse-level correction, and post-relaxation sweeps has error propagation operator MG 2 = ⋅ ⋅ which satisfies If a partitioning and approximation, , of are found that satisfy the assumptions given above, then this theorem establishes existence of an interpolation operator, , giving multigrid convergence with a direct tie to the approximation parameter, . In this theory, tight spectral equivalence between and is required to ensure good performance of the solver. This leads to the goal of constructing the fine-grid set so that there is a guarantee of tight spectral equivalence between and .
To meet that goal, MacLachlan and Saad 1 proposed to partition the rows and columns of the matrix in order to ensure the diagonal dominance of , so that could be chosen as a diagonal matrix. In particular, a diagonal dominance factor is introduced for each row , defined as With this, is said to be -diagonally dominant if ≥ for all ∈ , where > 1∕2 measures the diagonal dominance of . If is -diagonally dominant, then it is shown that the diagonal matrix, , with ( ) = (2 − 1 ) for all ∈ leads to 0 ⪯  ⪯ 2−2 2 −1 , giving a -dependent convergence bound if the other assumptions of Theorem 1 are satisfied. Furthermore, if is symmetric, positive-definite, and diagonally dominant, then this prescription for guarantees that all conditions of Theorem 1 are satisfied, so long as is -diagonally dominant. In addition to establishing this connection between the diagonal dominance parameter and the convergence parameter, , MacLachlan and Saad 1 posed the partitioning algorithm as an optimization problem for given > 1∕2, asking for the largest F-set such that ≥ for every ∈ . Such a -dominant guarantees good equivalence between a diagonal matrix, , and , and the largest such -set would make the coarse-grid problem smallest. This leads to an optimization problem of the form The greedy algorithm 1 for (3) selects rows and columns of the matrix for in a greedy manner to ensure the diagonal dominance of , as described in Algorithm 1. Here, the set contains all the DoFs that are unpartitioned (not yet assigned to be coarse or fine points). Initially, all degrees of freedom are assigned to , while the and sets are empty. The rows that directly satisfy the diagonal dominance criterion are initially added to the -set (and removed from ). If there are no diagonally dominant DoFs in the set, then the point having the least diagonal dominance is made a -point. This selection may make some other points in diagonally dominant, whereupon these DoFs are added to the -set and removed from , and the process is repeated until the set becomes empty. while ≠ ∅ do 10: = argmin ∈ {̂ } 11: if̂ ≥ then 15: end if 17: end for 18: end while 19: Return , 20: end function While the convergence bound for AMGr is attractive, it has several shortcomings. First and foremost is the strict assumption on diagonal dominance needed for the convergence guarantee to be valid -we explore cases where this is not the case below in Sections 4.3 and 4.4, and see that these problems are, indeed, more difficult to handle in this framework. Second, we emphasize that while the convergence bound obtained in Theorem 1 depends only on (and, thus, is independent of mesh size if the coarse grid is generated by Algorithm 1 or the techniques introduced below with constant ). These bounds are generally worse than the observed convergence for standard AMG approaches. Using = 0.56, as we do for all examples below, gives = 22∕3, and the convergence bound for = 1 in (2) is 0.977. While this can be improved by using either more relaxation per cycle or by accelerating convergence with a Krylov method, we consider only stationary cycles here, accepting poorer convergence factors in exchange for a direct relationship to the convergence theory summarized above. We note that this is also similar to the stationary convergence bound for the method in Napov and Notay 24 , which is 0.93 for problems that satisfy the assumptions therein.

SIMULATED ANNEALING
The optimization problem in (3) is a combinatorial optimization problem. Because such problems arise in many areas of computational science and engineering, significant research effort has been devoted to developing algorithms for their solution, both of general-purpose type (e.g., branch and bound techniques) and for specific problems (e.g., the travelling salesman problem) 34 . For many problems, the size of the solution space makes exhaustive brute-force algorithms infeasible; for some such problems, branch and bound techniques may be successful in paring down the solution space to a more manageable size. In many cases, however, there are no feasible exact algorithms, and stochastic search algorithms, such as simulated annealing (SA), genetic algorithms, or tabu search methods, can be employed 35 . In this work, we apply SA algorithms to approximate the global minimum of the optimization problem in (3).
Simulated annealing (SA) is a probabilistic method used to find global optima of cost functions that might have a large number of local optima. The SA algorithm randomly generates a state at each iteration and the cost function is computed for that state. The value of the cost function for a state determines whether the state is an improvement. If the current state improves the value of the objective function, it is accepted to exploit the improved result. The current state might also be accepted, with some probability less than one, even if it is worse than the previous state, however the probability of accepting a bad state decreases exponentially with the "badness" of the state. The purpose of (sometimes) accepting inferior states, known as "exploration" is to avoid being trapped in local optima, and the inferior intermediate states are considered in order to give a pathway to a globally better solution. The total number of iterations of SA depends on an initial "temperature" and the rate of decrease of that temperature. The temperature also affects the probability of accepting a bad state, with the exploration phase of the algorithm becoming less probable as the temperature decreases, to ensure convergence to a global optimum.

Idea and generic algorithm
Consider the set of fine points to be the current state. To find that maximizes a fitness function ( ), such as ( ) = | | in (3), simulated annealing proceeds as shown in Algorithm 2. The temperature starts at an initial value and decays by a factor ∈ (0, 1) at each iteration (the "cooling schedule"). A key choice when using simulated annealing is a method for choosing a neighbor statẽ that is close to the current state . In the case of optimizing the choice of a subset ⊂ Ω, a straightforward choice is to choosẽ by randomly adding or removing an element from . Finally, the function ( ,̃ , ) is the probability of accepting a move from a current state with fitness to the new state with fitness̃ . The standard acceptance function is This probability function always accepts transitions that raise the fitness, and sometimes accepts transitions that decrease the fitness. Occasionally accepting fitness-decreasing transitions is essential to escape local maxima in the energy landscape, with the chance of such transitions being controlled by the current temperature. By analogy with the physical annealing of metals, at high temperatures we will accept almost all fitness-decreasing transitions, allowing exploration of the fitness landscape, but as the temperature cools we will accept fewer of these transitions, until we eventually become trapped near a fitness maximum.
Algorithm 2 Generic simulated annealing (SA) algorithm 1: Initialize to a random state and to an initial temperature 2: for steps iterations do 3: Randomly pick a neighbor statẽ of 4: if rand(0, 1) < ( ( ), (̃ ), ) then Much work has been devoted to the design of optimal transition probability functions 36 and cooling schedules 37 in simulated annealing algorithms. While it can be shown that simulated annealing converges to the global maximum with probability one as the cooling time approaches infinity 38 , in practice the performance depends significantly on the selection of the neighbors̃ of the current state . A common heuristic for choosing neighbors is to select states with similar fitness, which is more efficient because we are less likely to reject such transitions, although this is in tension with the desire to escape from steep local maxima.
There are also different algorithmic possibilities for the handling of constraints on the state . Given a constraint subset of allowable states, one common method of enforcing the constraint is to pick only neighbor states̃ that are in the constraint subset. This method is appropriate if the constraint subset is connected and constraint-satisfying neighbors can always be found, and it is easy to see that this algorithm inherits all theoretical properties of the unconstrained version. An alternative method, which we will use in Section 3.2, is to allow constraint-violating neighbors to be selected but keep a record of the highestfitness constraint-satisfying state visited. This method is advantageous when constraint-violating paths in state space allow rapid transitions between (possibly disconnected) areas in the constraint subset. This algorithm also preserves the global guarantees of convergence under the condition that the global maximum is constraint satisfying.

Adaptation to coarse/fine partitioning
While it is possible to apply SA directly to the optimization problem in (3), preliminary experiments show that this is inefficient, due to the global coupling of the degrees of freedom in the optimization problem. To overcome this, we consider a domain decomposition approach to solving the optimization problem, where we first divide the discrete set of degrees of freedom into (non-overlapping) subdomains (the construction of which is considered in detail in the following subsection) and apply SA to the subdomain problems. These subdomain problems are not independent, therefore information is exchanged about the tentative partitioning on adjacent subdomains; this is accomplished in either an additive (Jacobi-like) or multiplicative (Gauss-Seidellike) manner. For sufficiently small subdomains, however, SA can efficiently find near-optimal solutions to localized versions of the optimization in (3), and we focus on this process below.
SA is run on each subdomain to partition its DoFs into (local) -and -sets. As the suitability of this partitioning in a global sense naturally depends on decisions being made on adjacent subdomains, we must be careful to consider what happens to DoFs in the global mesh adjacent to those in the subdomain. If these subdomains already have their own tentative ∕ partitions computed (as is expected to be the case on all but the first sweep through the domain), then this partitioning is considered fixed, and used to guide decisions on the current subdomain. If no tentative partitioning has yet been computed on a neighbouring subdomain, then the points in the neighbouring subdomain that are connected to some DoF of the current subdomain are considered to be -points, while all other DoFs in the neighbouring subdomain are considered to be -points. These assumptions are used only for the purposes of computing the partitioning on the current subdomain, to apply hard constraints on the partitioning.
Specifically, if Ω = {1, 2, … , } is the global set of degrees of freedom, partitioned into disjoint subdomains as Ω = ∪ =1 Ω , then the optimization problem to be solved on subdomain is given as where information from other subdomains enters in the constraint, both as additional points for which the constraint must be satisfied and in the right-hand side of the constraint where we sum over ∈ (and not ∈ ). Here, is the current set of points in Ω that are in the -set, and we define = Ω ⧵ to be the complementary -set. Further, we usē and to denote the sets of tentative -points in Ω and Ω, respectively, where the standard "closure" notation, Ω , denotes the set of points , such that either ∈ Ω or ≠ 0 for some ∈ Ω . In the localized optimization problem above, the set̄ contains both the points in and any point ∈ Ω ⧵ Ω such that either is in the -set in a neighbouring subdomain, Ω , that has a tentative partition or ∈ Ω for a neighbouring subdomain, Ω that does not yet have a tentative partition.
Within an annealing step on a given subdomain, Ω , we take the current tentative set as the initial guess for the partitioning, with the exception of the first step, where we take = ∅ for all ∈ Ω (to ensure we start from a configuration that satisfies the constraint). As a benchmark for the annealing process, we initialize ( ) as the largest size of a constraint-satisfying -set seen so far on Ω (taken to be zero on the first iteration), and to be the number of constraint-satisfying -points in the current set . At each annealing step on Ω , we swap points in and out of , either increasing, decreasing, or maintaining its size, with equal probability. The basic algorithm for these swaps is given in Algorithm 3, where we take the sets and as input, along with values and giving the numbers of points to swap from to and vice-versa. We note that there are two possible ways to do this swap, either selecting the elements from and independently and then moving the elements, or first moving the selected elements from and, then, selecting the elements from the updated set to move. We follow the first way as preliminary experiments suggested that it gives slightly better results. The main annealing step, as given in Algorithm 4, then takes the form of selecting whether to increase, decrease, or maintain the size of , checking if the selected action is possible, and performing it if it is. To increase the size of , we first check that has sufficient entries to move. If so, we increase the size of by selecting + entries of to move to and entries of to move to , for pre-determined values of and (typically = 1, = 0, although these values could also be drawn from a suitable distribution). If the decision is made to swap points, we check that both and have sufficient points to swap, then swap points from each set into the other (typically = 1). Finally, if the decision is made to decrease the size of , we check that it has points to remove, then move + points from to and points from to (ensuring + > 0; typically = 1, = 0). Finally, we measure the fitness of the resulting tentative -set and decide whether or not to accept it before decrementing the temperature by the relative factor and the number of further annealing steps to take by one.
The fitness score of a given (tentative) partition over Ω is directly calculated as the number of points in the set that satisfy the diagonal dominance criterion, as outlined in Algorithm 5. In the acceptance step of the algorithm, given as Algorithm 6, we compute the fitness score for the tentativē and compare it to that of the current (last accepted) set̄ . If the fitness score increases or remains same, then we automatically accept the step and update , , and . In this case, we additionally check if all points in the current̄ -set are constraint satisfying and, if so, if this set increases the value of ( ) . If so, we update the value of ( ) (and update the global -set). If̃ < , then the step is accepted with a probability that decreases with temperature and −̃ , but the additional check need not be done.

Localization and Gauss-Seidel variants
In Section 4, we consider both structured-grid and unstructured-grid problems; consequently, we consider both geometric and algebraic decompositions of Ω into {Ω } =1 . An important advantage of algebraic partitioning, however, is that it can be used to generate deeper multigrid hierarchies, since geometric partitioning can only be used to generate a single coarse level (which has unstructured DoF locations and, thus, does not naturally lead to further geometric partitioning). Here, we outline both strategies.   For structured grids, we can consider geometric decomposition of the fine grid into subdomains. For problems with (eliminated) Dirichlet boundary conditions, we typically satisfy the diagonal dominance criterion already at all points adjacent to a Dirichlet boundary, so these points are taken as -points from the beginning and not included in the decomposition into subdomains. A natural strategy is to subdivide the remaining points into square or rectangular subdomains of equal size (modulo boundary/corner cases). We consider this in Section 4 for both finite-difference (FD) and bilinear finite-element (FE) discretizations on uniform meshes. When the number of points (in one dimension) to be decomposed is not evenly divided by the given subdomain size (in one dimension), the right-most and bottom-most subdomains on the mesh are of smaller size, given by the remainder in that division. On structured grids, we can consider either a lexicographical Gauss-Seidel iteration over the subdomains or a four-colored Gauss-Seidel iteration. In preliminary experiments, the four-colored Gauss-Seidel iteration outperformed the lexicographical strategy by a small margin, so we use this.
On both structured and unstructured grids, we also consider algebraic decomposition using Lloyd aggregation to define the subdomains. Lloyd aggregation, proposed by Bell 39 , is a natural application of Lloyd's algorithm 40 to subdivide the DoFs of a matrix into well-shaped subdomains. Given a desired number of subdomains (or average size per subdomain), the unit-distance graph of the matrix is constructed and one "center" point for each subdomain is randomly selected. Each (tentative) subdomain is then selected as the set of points that are closer to the subdomain center point than to any other subdomain center, using a modified form of the Bellman-Ford algorithm. Then, for each subdomain, the center point is reselected (again using a modified form of the Bellman-Ford algorithm), as the current centroid of the subdomain (a DoF having maximum distance from the subdomain boundary, breaking ties arbitrarily). This process is repeated (reforming subdomains around the new centers, then reassigning centers) until the assignment to subdomains has converged. In the algebraic case, a multicoloured Gauss-Seidel iteration strategy is slightly more complicated, since we would need to compute a colouring of the subdomains; hence, we simply use lexicographical Gauss-Seidel to iterate between subdomains, noting that the advantage of the four-color iteration in the structured case is quite small.

Benchmark results
A natural comparison is with that of the greedy strategy 1 , which we provide in the following numerical results. For the structuredgrid discretizations, we have also explored "hand-tuned" optimization strategies to try and maximize the size of the global F-set while satisfying the constraint.
For the five-point finite-difference stencil of the Laplacian on a uniform 2D mesh, for any value 1 2 < < 4 7 , the diagonal dominance criterion will be satisfied if every -point has at least one -neighbour. An optimal strategy 33 for this case arises by dividing the fine mesh into "X-pentominos", sets of five grid points with one center point and its four cardinal neighbours, plus edge/corner cases where only a subset of an X-pentomino is needed. Then, the -set can be selected as the center points of each X-pentomino (or subset of one), with the remaining points assigned to the -set. In an ideal case (e.g., with periodic boundary conditions, or minimal edge cases), this leads to an -set with size equal to 4∕5 of the size of Ω. Such a coarsening is shown in the left panel of Figure 1. Note that some points adjacent to Dirichlet BCs are naturally chosen as coarse points in this strategy, despite their inherent diagonal dominance, because they are center points of an X-pentomino that includes just one point in the region away from the boundary. These could be equally well treated by making the interior point in these X-pentominos a -point and leaving the center point on the boundary as an -point, but there is no advantage to doing so. For the nine-point bilinear finite-element discretization of the Laplacian on a uniform 2D mesh, for any value of less than 4∕7, the diagonal dominance criterion will be satisfied if every -point has at least two -neighbours. Here, we partition the points (again except those adjacent to a Dirichlet boundary) into square subdomains of size 3 × 3, and select two points consistently in each subdomain as -points. In an ideal case, where we can fill the domain with 3 × 3 "bricks", this leads to an -set with 7∕9 of the size of Ω. Such a coarsening is shown in the right panel of Figure 1. When the domain cannot be filled perfectly with 3 × 3 bricks, the remaining square/rectangular regions still require one or two points to be selected, reducing the optimal size of the resulting set.

RESULTS
In the results below, we set the initial (global) temperature at = 1, and compute the reduction rate, , so that ts = 0.1, where ts is the total number of annealing steps to be attempted. Thus, all experiments end when the global temperature reaches 0.1. In addition, for each problem, we determine ts by fixing a total number of steps per DoF in the system, and report results based on the alloted number of steps per DoF. This work is further subdivided into Gauss-Seidel sweeps by fixing values of the number of annealing steps per DoF per sweep, with the number of sweeps determined as the ratio of total steps per DoF to steps per DoF per sweep. We first consider structured-grid experiments for both the finite-difference (five-point) and bilinear finite element (nine-point) discretizations of the Laplacian, using a 32 × 32 mesh to explore how much work is needed to determine near-optimal partitions using both geometric and algebraic divisions into subdomains. Using these experiments to identify "best practices", we then explore how the coarsening algorithm performs as we change the problem size, move from structured to unstructured finiteelement meshes, and isotropic to anisotropic operators.

Structured-grid discretizations with geometrically structured subdomains
We start by considering the five-point finite-difference stencil on a fixed (uniform) 32 × 32 mesh, and consider the effects of changing both the size of the Gauss-Seidel subdomains and the distribution of work in the algorithm. As a measure of quality of the results, we consider the maximum of the ratio | |∕|Ω| over all constraint-satisfying -sets generated in a single run of the annealing algorithm. As a comparison, for this problem, the best "by-hand" optimization of the size of the -set yields | |∕|Ω| = 0.8047, as depicted by the blue lines in Figure 2, while the greedy algorithm 1 yields | |∕|Ω| = 0.561 (not depicted because it is far from the data shown here). We vary three algorithmic parameters in Figure 2, the total number of SA steps per DoF, with values ranging from 5000 to 2 000 000, the number of SA steps per DoF in a single sweep of Gauss-Seidel on each subdomain, and the size of the subdomains used in the Gauss-Seidel sweeps. Of these, the number of SA steps per DoF per sweep of Gauss-Seidel is correlated least with the quality of results. For the smallest subdomain sizes (in particular, the 2 × 2 subdomains at the top left of Figure 2), we see a sharp improvement in performance when increasing this parameter from 1 to about 10, but plateauing results beyond that. For larger subdomains, there is again some advantage in most cases to choosing a value slightly larger than 1, but little advantage in choosing large values (and, indeed, we see sharp degradation in performance with larger values of this parameter on the 6 × 6 mesh at bottom right). Thus, in further results, we generally fix this value to be relatively small, either 1 or 5.
For small subdomain sizes -e.g., the 2 × 2 and 3 × 3 subdomain cases shown in the top row of Figure 2 -we observe little benefit of increasing the total amount of work in the algorithm, with comparable results for all numbers of total annealing steps per DoF. In contrast, for the 6×6 subdomains (bottom right of Figure 2), we see strong stratification of the results with increasing total number of SA steps. The 4 × 4 case also shows improving behaviour with increasing number of total SA steps, albeit less-so than with larger subdomains. Overall, we see this consistently in our numerical results: for small subdomain sizes, the quality of coarsening achieved quickly stagnates, with each subdomain having too little freedom to make adjustments into better global configurations while still satisfying local constraints. For larger subdomains, we see both better overall configurations and improvements with larger work budgets; however, we also see increasing numbers of SA steps per DoF needed to achieve the best performance.
A key question is how to balance the parameters in the SA algorithm to achieve reasonable performance at an acceptable cost. To examine this, we fix the total number of annealing steps per DoF and consider the best partitioning achieved by varying the subdomain size in the left panel of Figure 3. Here, we run for a number of different values of the number of SA steps per DoF per Gauss-Seidel sweep, and take the best partitioning observed over these grids for each subdomain size (noting that the optimal choice varies somewhat with subdomain size, typically being larger for small subdomain sizes and smaller for large subdomain sizes). When using 200 000 SA steps per DoF, there is a clear maximum in the graph for 4 × 4 subdomains, although the relative difference in quality in not substantial; however, when using 2 000 000 SA steps per DoF, we see continued improvement in the quality of the partitioning up to the 6 × 6 subdomain case. In the right panel of Figure 3, we look more closely at the convergence of the results with increasing numbers of SA steps per DoF for the case of 6 × 6 subdomains, again taking the best results obtained for different values of the number of SA steps per DoF per Gauss-Seidel sweep. Here, we see a clear improvement in the results up to (10 5 ) SA steps per DoF, and continued improvement up to 2 000 000 SA steps per DoF. We recall that we fix the SA temperature reduction rate, , so that → 1 as the number of SA steps increases, yielding the same total reduction in for each experiment. Thus, the results in Figure 3 are consistent with the expected behaviour of SA, that we can achieve results arbitrarily close to the global maximizer of our functional but only if we take many steps and slowly "cooling" the SA iteration. For a more practical algorithm, we emphasize the behaviour at lower numbers of SA steps / DoF, emphasizing that we achieve results within 5% of the best-known solution already with only 3000 steps / DoF, and within 2% at around 50 000 steps / DoF.
A natural question that arises from the right-hand panel of Figure 3 is whether the best meshes obtained occur early or late in the annealing process. That is, while we clearly see benefit from slow "cooling" of the annealing process (many SA steps per DoF), it is important to identify when the optimal results are obtained during the annealing process. Figure 4 shows the annealing history for a sample run, on the 32 × 32 uniform grid five-point finite-difference stencil, with 2 000 000 SA steps per DoF (for a total of almost 2 × 10 9 annealing steps). At right, we see the temperature decay, following an exponential curve from = 1.0 at the first step to = 0.1 at the final step. At left, we plot the changes in the ratio | |∕|Ω| compared to the by-hand optimal mesh for this grid (also shown). For comparison, the ratio from the greedy algorithm is also shown. We see that after an initial rapid improvement in the value of the ratio, there is a secondary period where performance improves notably but steadily, up to between 5 × 10 8 and 10 9 total annealing steps. This demonstrates that many annealing steps are, indeed, needed to reach the best partitioning seen here, although a good partitioning is still found after many fewer steps.
We next address the nature of the grids generated by annealing, and whether they resemble grids that could be selected geometrically for this problem. Figure 5 shows two of the grids generated, along with the subdomains used in their generation. These represent the "best" grids found by the annealing procedure, with ratios of | |∕|Ω| of around 0.795. While these grids yield competitive ratios, there is no clear global geometric pattern, nor obvious relationship to the best by-hand optimized grid shown in Figure 1. Furthermore, there are no clear improvements of these grids that could be readily made, such as single coarse points that could be omitted without leading to constraint violations. This suggests that the energy landscape for this problem is likely dominated by locally optimal configurations that are separated by states with constraint violations and/or sharp changes in energy.  Grid partitioning for 32 × 32 uniform grid with five-point FD stencil using two million SA steps per DoF. At left, the partitioning is generated using primarily 4 × 4 subdomains, with 25 SA steps per DoF per cycle. At right, the partitioning is generated using 6 × 6 subdomains, with 5 SA steps per DoF per cycle. The grid at left has 814 -points, while that at right has 816. Figure 6 shows how the performance of the annealing algorithm scales with problem size, for both the geometric choice of subdomains for the finite-difference discretization discussed so far, and for the finite-element discretization. In addition, an algebraic selection (discussed below) is added for comparison. All methods (including the by-hand optimization) perform relatively well for small meshes, but degrade as the mesh size increases. The amount of work needed to achieve these results with the annealing algorithm also increases with grid size. For the 8 × 8 mesh, using a single (global) subdomain, equally-good partitions to the by-hand optimization can be found with just 2000 annealing steps per DoF (and fewer when more subdomains are used). For the 16 × 16 mesh, more work and larger subdomains are needed to achieve such performance. With small (2 × 2 or 3 × 3) subdomains, even using 200 000 annealing steps per DoF did not achieve performance equal to the by-hand partitioning on the 16 × 16 mesh. For larger subdomains, the algorithm did equal the results of the by-hand partitioning, but with increasing work as subdomain size increased: for 4 × 4 subdomains, 10 000 annealing steps per DoF are needed, while 50 000 annealing steps per DoF are needed for 5 × 5 subdomains. For 6 × 6 subdomains, even using 200 000 annealing steps per DoF, we could not recover results matching the by-hand partitioning, although we speculate that this would have occurred with even more work invested. For larger domain sizes, the results shown in Figure 6 represent the best results found for a given grid over all runs with varying subdomain sizes, total SA steps per DoF, and SA steps per DoF per GS sweep. While these best results are, in general, achieved with the largest allocations of SA steps per DoF tried in our experiments, the marginal benefit of considering such large amounts of work to generate the coarsenings are quite low. Here, in all cases where we invested the computational effort to explore the question, we found less than a 1% improvement in | |∕|Ω| when increasing beyond (10 5 ) SA steps per DoF (and, in many cases, the improvement was only by one or two points).
For the nine-point finite-element stencil, we see similar results to those reported above and, consequently, do not include figures detailing the individual experiments in as much detail. Most notably, the best partitioning that we achieve by hand is a slightly worse in this case (as detailed in Section 3.4), and the partitioning using the greedy algorithm is better than in the FD case, yielding | |∕|Ω| = 0.752. At left of Figure 7, we show an analogous figure to Figure 2 for the case of 5 × 5 subdomains. Here, we observe the same stratification. However, we are able to recover the same quality of partitioning as the best by-hand partitioning using 5 × 5 subdomains and 2 000 000 SA steps per DoF. Also note that we see the same mild dependence on the number of SA steps per DoF per sweep as we did in the FD case. The right panel of Figure 7 shows how the quality of partitioning changes with subdomain size and total work budget. Similarly, there is an improvement in performance with additional work for larger subdomain sizes, but that improvement stagnates with additional work for smaller subdomain sizes. Sample grids generated for the finite-element case, including colouring to indicate subdomain choice, are shown in Figure 8.
Finally, we verify that the result two-level AMGr algorithms are at least as effective as predicted in theory, by measuring asymptotic convergence factors ( ), grid complexities ( grid , equal to the ratio of sum of the number of DoFs on each level of the hierarchy to that on the finest level), and operator complexities ( op , equal to the ratio of the sum of the number of nonzero entries in the system matrix on each level of the hierarchy to that on the finest level). We approximate by running the cycle with   Grid partitioning for 32 × 32 uniform grid with nine-point FE stencil using two million SA steps per DoF. At left, the partitioning is generated using 5 × 5 subdomains, with one SA step per DoF per cycle. At right, the partitioning is generated using 6 × 6 subdomains, also with one SA step per DoF per cycle. The grid at left has 814 -points, while that at right has 811. a random initial guess, (0) , and zero right-hand side, and then compute ‖ ‖ , where ( ) is the approximation (to the true solution, which is the zero vector) after cycles. Because the convergence factors of AMGr are somewhat larger than those for classical AMG, we take = 800 to ensure that we sample the asymptotic behaviour suitably. Table 1 reports this data for two-level cycles for both the FD and FE discretizations, for both the geometric subdomain choice considered here and the algebraic subdomain choice discussed next. Considering the geometric subdomain choice, we see grid-independent convergence factors of about 0.9 for the FD case and 0.7 for the FE case. While these are notably worse than are observed for typical multigrid methods for these problems, they are consistent with the existing results for AMGr and, in particular, conform with the convergence rate bound from Theorem 1 of 0.977 for = 0.56 and = 1. Notably, neither the convergence factor nor the measured complexities degrade substantially with problem size.

Structured-grid discretizations with algebraically chosen subdomains
We next consider partitioning the structured-grid problems using an algebraic choice of the subdomains based on Lloyd aggregation. Figure 9 shows the change in the maximum number of -points (scaled by the total number of DoFs) with the change of the subdomain size for both the FD and FE discretizations (left and right, respectively). Similarly to the case of geometric partitioning, we see relatively poor performance for small subdomain sizes, regardless of the work allocated to the SA process. For larger subdomain sizes, we see improving results with number of SA steps per DoF, as seen above. Also as seen above, the optimal subdomain size varies with total work allocation, increasing as we increase the amount of work per DoF, but even with 2 000 000 SA steps per DoF, we do not see the best performance with largest subdomains for the FE discretization. Considering variation in problem size, Figure 6 shows that, as in the geometric subdomain case, we see some decrease in performance as problem size grows, but that this decrease seems to (mostly) plateau at larger problem sizes. In comparison with the geometric subdomain choice, we see some small degradation in performance with algebrically chosen subdomains, particularly in the FD case, but it is small in comparison with the optimality gap between the by-hand optimized solutions and those generated by SA.
Performance of the resulting two-grid cycles is tabulated in Table 1, where we see very comparable convergence factors, as well as grid and operator complexities, as in the geometric subdomain case. Table 2 tabulates convergence factors, grid complexities, and operator complexities for three-level V-and W-cycles. For the FD case, we see some degradation in convergence when using V-cycles, while W-cycles give convergence factors similar to those in the two-level case. For the FE problem, there is less degradation for V-cycles, but still W-cycles are required to get convergence factors comparable to the two-level case. As we have increased the depth of the multigrid hierarchy, the grid and operator complexities in Table 2 are naturally larger than those in Table 1, but the growth in these complexity measures is consistent with optimally scaling multigrid methods.

Anisotropic problems on structured grids
Next For 0 < < ∕2, the axis of the small diffusion coefficient in the problem rotates clockwise from being in the positive -direction for small to the positive -direction for ≈ ∕2. Anisotropic problems cause difficulty for the greedy coarsening algorithm 1 , where large grid and operator complexities and poor algorithmic performance were overcome by augmenting the coarse grids with the second pass of the Ruge-Stüben coarsening algorithm and using classical AMG interpolation in place of the AMGr interpolation operator. We emphasize that both the FD and FE discretizations of this problem result in non-diagonally dominant system matrices and, consequently, the convergence guarantee from Theorem 1 does not apply in this setting. We include these problems, nonetheless, both to stress-test our approach and highlight the need for further research in this direction.
As an example of a problem in this class, we fix = 10 −6 and = ∕3. Figure 10 shows the coarse-grid points selected for the bilinear finite-element discretization of the problem on a uniform 32 × 32 mesh. At left, we see that the partitioning generated by the SA algorithm correctly detects that this is an anisotropic problem with strong coupling primarily in thedirection and weak coupling primarily in the -direction, producing grids that are consistent with semi-coarsening, with some

FIGURE 10
Grid partitioning for an anisotropic diffusion problem with = 10 −6 and = ∕3, discretized on a 32 × 32 mesh using the bilinear FE stencil. The initial partitioning, at left, is generated using algebraic subdomains averaging 20 points per subdomain, using 2 000 000 SA steps per DoF and 1 SA step per DoF per GS sweep, and has 664 -points and 360 -points. At right is the grid augmented using the second pass of the classical AMG algorithm, resulting in 343 -points and 681 -points.
regions of coarsening along diagonals. Unfortunately, however, this grid leads to poor convergence using AMGr interpolation. As in the original experiments by MacLachlan and Saad 1 , each fine-grid point has multiple coarse-grid neighbours (leading to an increase in operator complexity, due to the nonzero structure of −1 ), but since is no longer diagonally dominant in the anisotropic case, the assumption that − − is positive semi-definite is violated, and the resulting performance is poor. To overcome this failure, we augment the coarse-grid set using the second-pass algorithm from Ruge-Stüben AMG, using the classical strength of connection parameter of 0.30 to determine strong connections in the graph. At right of Figure 10, we show the resulting graph, which now satisfies the requirement that every pair of strongly connected -points has a common -neighbour (which is clearly not satisfied in the partitioning at left). Unfortunately, this comes at a heavy cost, as the grid at right now has many more -points than we would like. Below, we verify that these grids lead to effective multigrid hierarchies, when coupled with classical AMG interpolation, but note the increased grid and operator complexities in these hierarchies. A key question for future work is whether we can determine algebraic interpolation operators for grids such as those at left that lead to effective AMGr performance. Tables 3 and 4 present convergence factors and grid and operator complexities for twoand three-level cycles, respectively. For geometric partitioning, we use 200 000 SA steps per DoF with one or two SA steps per DoF per cycle and 5 × 5 or 6 × 6 subdomains (depending on what worked best for a given problem). Similar choices are made using algebraic partitioning, although we see slight improvements using 2 000 000 SA steps per DoF for some problems. For the three-level tests, we uniformly use 200 000 SA steps per DoF, with 1 SA step per DoF per GS sweep, and algebraic subdomains with average size of 36 points. In Table 3, we again see that there is relatively little difference in results generated using geometric and algebraic choices of the Gauss-Seidel subdomains. Notably, both choices lead to effective cycles for both the finite-difference and finite-element discretizations, albeit at the cost of increased grid and operator complexities. Table 4 again shows degradation in performance moving from two-grid to three-grid V-cycles, but that three-grid W-cycles mostly recover comparable convergence to the two-level case, with some notable degradation for the finite-element operator.

Discretizations on unstructured grids
We next consider two diffusion problems with homogeneous Dirichlet boundary conditions on unstructured triangulations of square domains, discretized using linear finite elements. First, we consider an isotropic diffusion operator (−∇ ⋅ ∇ = ) on [−1, 1] 2 , generated by taking an initially unstructured mesh, refining it a set number of times, and then smoothing the resulting mesh. We consider three levels of refinement for this example, resulting in meshes with 1433, 5617, and 22 241 DoFs. The SAbased partitioning scheme appears to perform similarly well in this setting, with a splitting found for the smallest mesh shown in Figure 11. While we no longer have hand-generated estimates of the optimal coarsening factors, we note that the SA-based partitioning scheme outperforms the greedy algorithm, generating -sets with 1106, 4241, and 16 604 points, for these three  grids, respectively, in comparison to -sets of 1024, 3916, and 15 079 points. Two-and three-level convergence factors, as well as grid and operator complexities for these systems are given in Table 5. Here, we observe that these results are quite similar to those seen for the structured-grid discretizations above, with some degradation in convergence seen for the three-level V-cycle results, but not in the W-cycle results. Here, the three-level results are computed using 200 000 SA steps per DoF, with 1 SA step per DoF per GS sweep, on subdomains with an average of 36 points per subdomain. Slight modifications to these parameters yield small improvements in two-level results; as a result, we use these choices in the table. We next consider the anisotropic diffusion operator (again with Dirichlet boundary conditions) considered above, with = 0.01 and = ∕3, on an unstructured triangulation of the unit square taken from Brannick and Falgout 23 , matching the problem labelled 2D-M2-RLap in that paper. As in Brannick and Falgout 23 , we consider three levels of refinement of the mesh, with 798, 3109, and 12 273 DoFs, respectively. Figure 12 shows two different partitions for this mesh. At left, we give the partitioning generated using the SA-based partitioning algorithm proposed here, and at right we give the splitting after a second pass of the Ruge-Stüben coarsening algorithm where strength of connection is computed using the classical strength parameter of 0.55. As with the structured-grid case above, we found the second pass is needed to achieve good convergence factors for the anisotropic problem, but that it does so at the expense of coarsening at a much slower rate. Convergence factors, grid  complexities, and operator complexities for these problems are reported in Table 6. As above, the complexities are higher for the anisotropic operators than the isotropic (due to the use of the second pass). Here, we observe degradation in convergence with grid refinement, although again with less degradation for W-cycles than V-cycles. The three-level results are computed with the same parameters as for the isotropic problem above, with more SA steps per DoF used for the two-level results, as this yielded slight improvements. While the performance reported in Table 6 is far from optimal, we note that the convergence factors for the larger two meshes are better than those reported for compatible relaxation 23 , while the operator complexities are comparable to those reported therein for BoomerAMG 5 . We again emphasize that these anisotropic diffusion problems do not satisfy the convergence bounds stated above and, consequently, are included to stress-test the framework presented here.

Convection-diffusion problems on structured grids
For our final tests, we consider the upwind finite-difference discretization of two singularly perturbed convection-diffusion equations. While these problems lead to non-symmetric discretization matrices (and, as such, Theorem 1 no longer applies), they represent a plausible set of test problems for reduction-based methods due to their "nearly triangular" M-matrix structure 31 .
We discretize these problems on uniform × meshes with the standard 5-point finite-difference for the diffusion term, and a first-order upwind finite-difference discretization for the convection terms. In all experiments, we use 200 000 SA steps per DoF to generate the partitioning, and generate the AMGr interpolation operator by computing  Table 7 reports measured asymptotic convergence factors, along with grid and operator complexities for these problems, showing insensitivity to both problem size and the singular perturbation parameter, . Figure 13 shows two sample partitionings, one for each convection direction, illustrated at = 10 −5 for = 32.

CONCLUSIONS AND FUTURE WORK
While existing heuristics for coarse-grid selection within AMG offer reliable performance for a wide class of problems, there remains much interest in both improving our understanding of AMG convergence and developing new approaches that offer guarantees of convergence for even wider classes of problems. We believe a promising, yet underexplored, area of research is in

FIGURE 13
Partitioning for convection-diffusion problems on uniform grids, for the grid-aligned case, at left, and the nongrid-aligned case, at right. In both cases, we depict the partitioning for = 10 −5 and = 32, generated using 200 000 SA steps per DoF.
the extension of reduction-based AMG approaches, that have been shown to be effective for both isotropic diffusion problems 1,26 and more interesting classes of problems, such as some hyperbolic PDEs 31 . In this paper, we propose a new coarsening algorithm for AMGr, based on applying simulated annealing to the optimization problem first posed by MacLachlan and Saad 1 . For isotropic problems on both structured and unstructured meshes, we find that the SA-based partitioning approach outperforms the original greedy algorithm, sometimes dramatically, while sharing some of the existing limitations of the AMGr framework, particularly for anisotropic problems. Nonetheless, we see this as an important proof-of-concept, showing that randomized search and other derivative-free optimization algorithms can be successfully applied to the combinatorial optimization problems in AMG coarsening.
The main drawback of this approach is the high computational cost of simulated annealing. For example, computing a single coarse grid for a 32 × 32 mesh with 200 000 SA steps per DoF required around 2.5 hours on a modern workstation, primarily for evaluating the fitness functional in the optimization, with expected scaling in problem size and number of SA steps per DoF.
Two key questions that this raises are whether the fitness functional can be approximated in a more efficient manner (e.g., using a value neural network) and whether other optimization techniques that rely on fewer samplings of the fitness functional can be applied. Both of these are topics of our current research.
This research also exposes a known weakness of the AMGr methodology (and, to our knowledge, one of all AMG-like algorithms with guarantees on convergence rates) when applied to problems that are not diagonally dominant (such as finiteelement discretizations of anisotropic diffusion equations). Another current research direction is identifying whether the diagonal choice of can be generalized to yield better convergence (while retaining a guaranteed convergence rate), and whether such changes can be accommodated within the optimization problems solved here.