Algebraic multilevel preconditioners for the graph Laplacian based on matching in graphs
Abstract.
This paper presents estimates of the convergence rate and complexity of an algebraic multilevel preconditioner based on piecewise constant coarse vector spaces applied to the graph Laplacian. A bound is derived on the energy norm of the projection operator onto any piecewise constant vector space, which results in an estimate of the twolevel convergence rate where the coarse level graph is obtained by matching. The twolevel convergence of the method is then used to establish the convergence of an Algebraic Multilevel Iteration that uses the twolevel scheme recursively. On structured grids, the method is proven to have convergence rate and complexity for each cycle, where denotes the number of unknowns in the given problem. Numerical results of the algorithm applied to various graph Laplacians are reported. It is also shown that all the theoretical estimates derived for matching can be generalized to the case of aggregates containing more than two vertices.
1991 Mathematics Subject Classification:
65N30, 65N151. Introduction
Algebraic Multigrid (AMG) attempts to mimic the main components of Geometric Multigrid in an algebraic fashion, that is, by using information from the coefficient matrix only to construct the multilevel solver. The basic algorithm uses a setup phase to construct a nested sequence of coarse spaces that are then used in the solve phase to compute the solution. The two main approaches to the AMG setup algorithm are classical AMG [3, 4] and (smoothed) aggregation AMG [13, 17, 9, 18, 11], which are distinguished by the type of coarse variables used in the construction of AMG interpolation.
In the classical AMG algorithm, the coarse variables are chosen using a coloring algorithm which is designed to find a suitable maximal independent subset of the fine variables. Then, given the coarse degrees of freedom, a row of interpolation is constructed for each fine point from its neighboring coarse points. In contrast, the aggregationbased AMG setup algorithm partitions the fine variables into disjoint subdomains, called aggregates. Then, a column (or several columns as in [18]) of interpolation is associated to each aggregate, which has nonzero entries only for the unknowns belonging to this aggregate. The focus of this paper is on the development and, in particular, the analysis of the latter aggregationtype methods.
The idea of aggregating unknowns to coarsen a system of discretized partial differential equations dates back to work by Leont’ev in 1959 [12]. Simon and Ando developed a related technique for aggregating dynamic systems in 1961 [13] and a twogrid aggregationbased scheme was considered in the context of solving Markov chain systems by Takahashi in 1975 [15]. Aggregationbased methods have been studied extensively since and numerous algorithms and theoretical results have followed [9, 18, 11]. Vanék introduced an extension of these methods known as smoothed aggregation multigrid in which smoothing steps are applied to the columns of the aggregationbased interpolation operator to accelerate twolevel convergence and a modification of this twolevel algorithm with overcorrection is presented in [17]. A multilevel smoothed aggregation algorithm and its convergence analysis are found in [16] and, in [19], an improved convergence theory of the method is presented. The latter theory is then extended to allow for aggressive coarsening, provided an appropriate polynomial smoother is used [8]. A further generalization known as adaptive smoothed aggregation is developed in [7]. Variants of the above approaches continue to be developed for use in scientific computing and have been developed for higher order partial differential equations [18], convection diffusion problems [11], Markov chains [14, 2], and the Dirac equation in quantum chromodynamics [5].
In this paper, an aggregation based Algebraic Multigrid method for the graph Laplacian is presented. The approach constructs the sequence of coarse graphs recursively using a pairwise aggregation, or matching, form of interpolation. However, it is demonstrated here, that the convergence rate of a twolevel method based on such a construction is uniformly bounded for the graph Laplacian on general graphs and, thus, can be used within an Algebraic Multilevel Iteration (AMLI) [1, 19] as a preconditioner to the Conjugate Gradient iteration to obtain a nearly optimal solver. A noteworthy feature of the approach is its simplicity, which makes it possible to analyze the convergence and complexity of the method with few assumptions and without any geometric information.
The remainder of the paper is organized as follows. In Section 2, we introduce the graph Laplacian problem and discuss some of its applications. In Section 3, we introduce a graph matching algorithm and demonstrate that the energy norm of the projection onto the coarse space is a key quantity in deriving convergence and complexity estimates of the method. Additionally, we introduce an approach computing an approximation of the energy norm of this projection operator. In Section 4, we present an analysis on the twolevel method for the graph Laplacian operator. In Section 5, we consider the convergence and complexity of the resulting AMLI method, and in Section 6 we provide numerical results and address some practical issues of the method.
2. Problem formulation and notation
Consider an unweighted connected graph , where denotes the set of vertices and denotes the set of edges of . The variational problem considered here is as follows: Given an satisfying , where is a constant vector, find a , where denotes the cardinality of the set of vertices , such that
(2.1) 
where
(2.2) 
Define the discrete gradient operator such that
where and are standard Euclidean bases. The operator is named graph Laplacian since
The operator is symmetric and positive semidefinite and its kernel is the space spanned by the constant vector. These properties can also be verified by the matrix form of defined in the following way
where is the degree of the th vertex.
Efficient multilevel graph Laplacian solvers are important in numerous application areas, including finite element and finite difference discretizations of elliptic partial differential equations (PDE), data mining, clustering in images, and as a preconditioner for weighted graph Laplacians. Moreover, the theory developed here for multilevel aggregation solvers applied to graph Laplacians should provide insights on how to design a solver for more general weighted graph Laplacians, which cover also anisotropic diffusion problems. Two generalizations of the graph Laplacian systems are as follows.

Weighted graph Laplacians: Assume that the graph is weighted and the th edge is assigned a weight , then the corresponding bilinear form of is
Define as a diagonal matrix whose th diagonal entry is equal to , then the matrix can be decomposed as
Finite element and finite difference discretizations of elliptic PDEs with Neumann boundary conditions results in such weighted graph Laplacians.

Positive definite matrices: Assume that is defined in terms of the bilinear form
By introducing a Lagrange multiplier , the system can be rewritten as an augmented linear system
These graph Laplacians with lower order terms are similar to discretized PDE with Dirichlet boundary conditions. A solution for this augmented linear system directly results to the solution of .
The present paper focuses on designing a multilevel preconditioner that is constructed by applying recursively a space decomposition based on graph matching. The aim is to analyze the matching AMLI solver for the graph Laplacian in detail as a first step in gaining an indepth understanding of a multilevel solver elliptic PDEs. The extension of the proposed algorithm to general graph problems is also a subject of current research.
3. Space decomposition based on matching
In this section, an outline of the basic idea of matching is provided, a commutative diagram which can be used to estimate the energy norm of the projection onto the piecewise constant coarse vector space resulting from a matching in a graph is given, and some auxiliary results that are needed later on in the convergence analysis are discussed.
3.1. Subspaces by graph partitioning and graph matching
A graph partitioning of is a set of subgraphs such that
In this paper, all subgraphs are assumed to be non empty and connected. The simplest non trivial example of such a graph partitioning is a matching, i.e, a collection (subset ) of edges in such that no two edges in are incident.
For a given graph partitioning, subspaces of are defined as
Note that each vertex in corresponds to a connected subgraph of and every vertex of belongs to exactly one such component. The vectors from are constants on these connected subgraphs. Of importance is the orthogonal projection on , which is denoted by , and defined as follows:
(3.1) 
Given a graph partitioning, the coarse graph is defined by assuming that all vertices in a subgraph form an equivalence class, and that and are the quotient set of and under this equivalence relation. That is, any vertex in corresponds to a subgraph in the partitioning of , and the edge exists in if and only if the th and th subgraphs are connected in the graph . Figure 1 is an example of matching of a graph and the resulting coarse graph.
As mentioned above, the reason to focus on matching is that it simplifies the computation of several key quantities used in the upcoming estimates derived for a perfect matching and it is possible to show that a matching which is not perfect can be analyzed in a similar way.
3.2. Commutative diagram
Let be the discrete gradient of a graph Laplacian , as defined in (2.1), and be defined as in (3.1). Assume that there exists an operator such that the following commutative diagram holds true:
The proof of this assumption is provided later on. From the commutative relation it follows that
(3.2) 
Thus, an estimate on the seminorm of amounts to an estimate of the norm of . In the next subsection, an explicit form of is constructed and an estimate of its norm is derived.
Remark 3.1.
A more general approach for weighted graph Laplacians is to assume that the weight matrix , therefore the bound on the norm becomes
where can have some negative weights, which results in a matrix this is complex valued. A detailed analysis in such a setting and the application of this idea to anisotropic diffusion problems are discussed in [6].
3.3. Construction of in case of piecewise constant spaces
Here, we proceed with an explicit construction and norm estimate of the operator .
For any graph partitioning in which the subgraphs are connected, a given edge belongs to the set of “internal edges”, whose vertices belong to the same subgraph, or to the set of “external edges”, whose vertices belong to two distinct subgraphs. For example, let and denote the subgraphs 1 and 2 in Fig. 2, then is an internal edge and is an external edge.
Since the vector has the same value on the two endpoints of the edge , we have that . Accordingly, all entries in , the th row of , are set to zero:
For the external edge , it follows that satisfies
(3.3) 
for every . The following Lemma is useful in computing explicitly the entries of .
Lemma 3.2.
Let be a positive semidefinite operator and let be a basis in . Assume that the null space of is one dimensional, namely there exist a nonzero vector such that , and for every integer we have . We then have:

For any , the operator with is invertible.

The following identity holds for all :
Proof.
To establish (i) it suffices to show that implies . Assuming that for some it follows that:
Note that both terms on the right side of the above identity are nonnegative and, hence, their sum can be zero if and only if both terms are zero. Since is positive semidefinite by assumption with one dimensional null space, from we conclude that for some . For the second term we have that , and since for all , it follows that and hence . This proves (i).
Remark 3.3.
A special case is given by taking and , which denote the standard Euclidean bases. Then, it follows that
(3.4) 
in which denotes the average value of .
Next, denote by the restriction of to a subgraph , and set . Then, can be expressed as the th component of for in (3.5).
(3.5) 
where are the local indices of the vertex set , the operator and the term are the averaging operator and the constant vector restricted on the subgraph , and maps the global edge indices to the local edge indices.
Applying this formula for and , gives the row of the operator on the edge that connects and as follows
(3.6) 
Here is given by
which then makes the summation in (3.6) valid. The vector is defined in a similar way.
Assume that the global indices of the vertices in and are ordered consecutively as decreasing integers starting at and increasing integers starting at . Then, the th row of can be expressed as
(3.7) 
where the number is on the th position in this row of . Note that from (3.5) it follows that the property (3.3) holds for this construction of , since by definition of
where , and and , both in local indices, are the incident vertices of .
4. A twolevel method
In this section, the orthogonal projection given in (3.1) based on a matching is proven to be stable assuming the maximum degree of the graph is bounded. Then, a twolevel preconditioner is derived and the condition number of the system preconditioned by this twolevel method is proven to be uniformly bounded (under the same assumption).
4.1. Twolevel stability
The construction of for a matching proceeds as follows. First, note that all rows of that correspond to an edge are identically zero. On the other hand, if the edge , then it is an external edge and, thus, by (3.7), the th row of is
Hence,
(4.1) 
The alternative way of describing the entries in is by showing that,
(4.2) 
Formula (4.1) implies that, the th row of can be a zero row if , or a row with 3 nonzero entries if , which results to
Formula (4.2) implies that, the th column of can have exactly 1 nonzero entry if , or nonzeros entries whose values are if . Here is the number of edges satisfying and for any given , thus is bounded by , where is the maximum degree of the graph, since an edges can have at most neighboring edges. This leads to
On a graph whose maximal degree is larger or equal to 2, the estimates on the infinity norm and norm of result to the following estimate on :
(4.3) 
Remark 4.1.
Applying Gerschgorin’s theorem directly to the matrix leads to a sharper estimate: .
Formula (4.3) implies directly the following lemma.
Lemma 4.2.
On any graph whose maximum degree is 2 (e.g. such graph is a path), the operator defined in (4.1) satisfies and the following estimate holds
Numerical tests show that this is a sharp estimate on the seminorm and that leads to fast convergent and reliable AMG methods.
4.2. A twolevel preconditioner
Here, using an estimate of the stability of the matching projection (i.e, the norm , where is defined via the matching) twolevel convergence is established. Assume that for a graph Laplacian a perfect matching is given and consider the matrix whose th column is given by
(4.4) 
where and is the th edge in . Further, define to be the projection from to , i.e.,
Similar to the definition of , define as the matrix whose columns are given by
(4.5) 
where and is the th edge in . Then, the matrix is orthogonal and the columns of and form a hierarchical bases, which can be used to relate the twolevel method to a block factorization as follows.
Given , , and , define
A direct calculation then shows that
where
(4.6) 
is the Schur complement and
(4.7) 
Next, define as the unweighted coarse graph and denote by the graph Laplacian of . In contrast to most of the existing AMG methods, here , except in special cases, e.g., for 1 dimensional problems. Let, be a positive constant such that
(4.8) 
Then, the fact that all weights in the graph corresponding to are larger than or equal to one implies , and
Consider the twolevel preconditioner which uses the coarse graph Laplacian by
Let be a preconditioner for , and be a preconditioner for the graph Laplacian . Then, a twolevel preconditioner is defined by
(4.9) 
where
As observed in [10] and [19], this gives a block matrix representation of the twolevel method
where the pseudoinverse operator denoted by is used since the graph Laplacian is semidefinite. The convergence of the twolevel method can now be estimated by comparing and the preconditioner .
The remainder of this section is dedicated to establishing a spectral equivalence between and for the twolevel matching algorithm. The proof uses the following Lemma.
Lemma 4.3.
For any the Schur complement as given in (4.6) satisfies
(4.10) 
Proof.
Note that
because here, is the orthogonal projection of onto the space spanned by the columns of and, thus, minimizes the distance (in norm) between and this space. Hence,
Let be a vector satisfying , then the following lemma now holds.
Lemma 4.4.
Let , where is defined as in (4.8). Then for any , such that , we have
(4.11) 
Proof.
By Lemma 4.3 we have
Furthermore,
Note that the only difference between the preconditioners and is that the former matrix uses , whereas the latter uses to define the 22 block. The spectral equivalence constant between the operators and is obtained as follows:
which implies
Hence, for any and
which is equivalent to (4.11) since is nonsingular. ∎
Since the twolevel method requires exact solvers for and the graph Laplacian , the convergence rate of a method that uses which is defined by replacing these exact solves with approximate ones is of interest. Combining Lemma 4.4 and the twolevel convergence estimate (Theorem 4.2 in [10]), yields the following result.
Theorem 4.5.
If the preconditioners and are spectrally equivalent to and such that
then
(4.12) 
Note that this estimate reduces to (4.11) when and .
4.3. Convergence estimate for matching
We here show the sharpness of the estimation given by Theorem 4.5 when the graph Laplacian corresponds to a structured grid, and the coarse space is given by aligned matching.
Define an dimensional hypercubic grid as the graph Laplacian such that the following conditions are satisfied.

A vertex corresponds to an vector , and , . Here is an Euclidean basis and are given positive integers that represent the numbers of vertices along all dimensions.

An edge is in the edge set if and only and .
Then the energy norm can be estimated for aligned matching on a hypercubic grid .
Lemma 4.6.
Let be an dimensional hypercibic grid and is a fixed dimension. Assume that is an even number. The matching along the th dimension is defined as
Let be the projection onto the piecewise constant space resulting from the matching . Then satisfies .
Proof.
Define the set be the collection of all edges along the th dimension, as
Also define and the graph Laplacians and , derived from and respectively. ∎
The graphs in the set are paths, whose maximum degree is 2, and is a also matching on these paths. Therefore by Lemma 4.2 it is true that
(4.13) 
On the other hand, the matching is aligned on the set , meaning that any two matched pairs are connected through 0 or 2 edges in , thus the edges in set can then be subdivided into many sets of edges of the same type, one of which is shown in Fig. 3.
Notice that in this figure, the edge and are in , while and are in . Using the definition of , the energy norm of is estimated on the the subset of indicated by Fig. 3, by
Thus implies that
(4.14) 
Remark 4.7.
A similar estimate follows for aligned partitionings consisting of line segments of size . Namely, in this case it can be shown that holds. This estimate in turn agrees with properties of Chebyshev polynomials, suggesting the use of an AMLI method equipped with certain Chebyshev polynomials. Comparing this result with the result from Theorem 4.6 suggests that using a more shape regular partitioning rather than one consisting of lines is more appropriate since this gives smaller values of the seminorm .
A bound on the constant follows by using that is well conditioned and that its condition number depends on the degree of the graph, but not on the size of the graph.
Lemma 4.8.
Let be the perfect matching on a graph maximum whose degree is , and let be defined as in (4.5), then we have
Proof.
The norm of the vector is computed by definition:
We also have
From the Lemma it follows that for any there exists a smoother such that the bound on the constant in Theorem 4.5 is
This result in turn implies that an efficient solver for can be constructed by applying a few Conjugate Gradient iterations with an overall cost that is linear with respect to the size of .
The constant in (4.8) can be estimated by checking the weights of the graph for the graph Laplacian . Taking any two distinct subgraphs (edges) in the matching, say the th and th such that , it follows that the corresponding entry is equal to the number of exterior edges that connect to these subgraphs. For an aligned matching aligned a fixed dimension of a hypercubic grid, these weights are bounded by . For any general graph , the weights in are bounded by , since there are at most distinct edges that connect to any other 2 distinct edges. Then, letting to denote the unweighted graph Laplacian on the graph defined by , and noting that all offdiagonal entries of are equal to , it follows that
Remark 4.9.
These estimates can be generalized to other subgraph partitionings in a similar way. As an example, consider again a graph for a hypercubic grid of any dimension. Then, for line aggregates of size (aligned with the grid) the following estimate holds
Such estimates give insight into the design of a nearly optimal multilevel method. Moreover, the bounds are sharp enough, namely, the corresponding multilevel method can be proven to have convergence rate and complexity.
5. Algebraic multilevel iteration (AMLI) based on matching
In this section, a multilevel method that uses recursively the twolevel matching methods from Section 4.2 in combination with a polynomial stabilization, also known as Algebraic Multilevel Iteration (AMLI) cycle is analyzed. Here, the focus is on proving a nearly optimal convergence rate, that is, one which is nearly independent of the number of unknowns , and at the same time has low computational complexity.
5.1. Multilevel hierarchy
Assume that is an graph Laplacian matrix where . For define the matching and the prolongation operator according to (4.4), then compute the graph Laplacian of the coarse graph (Recall that,