2.1 Model elliptic PDE operators

We consider the boundary value problems

where $\unicode[STIX]{x1D6FC}:\unicode[STIX]{x1D6FA}\mapsto \mathbb{R}^{d\times d}$ is an SPD matrix function satisfying

for some positive constants $\unicode[STIX]{x1D6FC}_{0}$ and $\unicode[STIX]{x1D6FC}_{1}$ . Here $d=1,2,3$ , and $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{d}$ is a bounded domain with boundary $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ .

A variational formulation for (2.1) is as follows: find $u\in V$ such that

Here

and $V$ is a Sobolev space that can be chosen to address different boundary conditions accompanying equation (2.1). One case is that of mixed boundary conditions

where $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{D}\cup \unicode[STIX]{x1D6E4}_{N}$ . The pure Dirichlet problem is when $\unicode[STIX]{x1D6E4}_{D}=\unicode[STIX]{x1D6E4}$ , while the pure Neumann problem is when $\unicode[STIX]{x1D6E4}_{N}=\unicode[STIX]{x1D6E4}$ . We thus have $V$ as

When we consider a pure Dirichlet problem, $\unicode[STIX]{x1D6E4}_{D}=\unicode[STIX]{x1D6E4}$ , we denote the space by $V=H_{0}^{1}(\unicode[STIX]{x1D6FA})$ . In addition, for pure Neumann boundary conditions, the following condition is added to ensure the existence of the solution to (2.3):

One of the most commonly used model problems is when

which corresponds to the Poisson equation

This simple problem provides a good representative model for isotropic problems.

There are two other cases that are of special interest. The first case is when $\unicode[STIX]{x1D6FC}$ is a scalar and has discontinuous jumps such as

The second case is when $\unicode[STIX]{x1D6FC}$ is a diagonal matrix such as (for $d=2$ )

which corresponds to the operator

In both cases we assume that $\unicode[STIX]{x1D716}$ is sufficiently small to investigate the robustness of algorithms with respect to discontinuous jumps and anisotropy.

Figure 2.1.

(a) Regular (uniform) triangulations for the unit square, and (b) unstructured mesh approximating the unit disk.

2.2 Examples of finite difference and finite element discretizations

As an illustrative example, we consider a finite difference discretization of the Poisson equation (2.8) with pure Dirichlet boundary conditions on the unit square $\unicode[STIX]{x1D6FA}=(0,1)\times (0,1)$ . We consider a uniform triangulation of $\unicode[STIX]{x1D6FA}$ (see Figure 2.1(a)), and we set

We use the standard central difference approximation to the Laplacian operator

The finite difference scheme is then given by

where

and $u_{i,j}\approx u(x_{i},y_{j})$ . The approximations $u_{i,j}$ are found by solving a linear system. We order the points $(x_{i},y_{j})$ lexicographically, and for $k=1,\ldots ,n^{2}$ we have

We can then write (2.12) as

where

A slightly different scheme is obtained using more of the neighbouring points of $(x_{i},y_{j})$ . We can build an approximation using 8 points $(x_{i\pm 1},y_{j\pm 1})$ together with the ‘centre’ point $(x_{i},y_{j})$ . As a result we have the following 9-point finite difference scheme:

Again, if we order $(x_{i},y_{j})$ lexicographically, then (2.17) is the linear system (2.15) corresponding to the 9-point finite difference discretization of the Laplace equation, where

We now give an example of finite element discretization. Given a triangulation ${\mathcal{T}}_{h}$ for $\unicode[STIX]{x1D6FA}$ , such as that given in Figure 2.1, let $V_{h}\subset V$ be a finite element space consisting of piecewise linear (or higher-order) polynomials with respect to the triangulation ${\mathcal{T}}_{h}$ . The finite element approximation of the variational problem (2.3) is: find $u_{h}\in V_{h}$ such that

Assume $\{\unicode[STIX]{x1D719}_{i}\}_{i=1}^{N}$ is the nodal basis of $V_{h}$ , namely $\unicode[STIX]{x1D719}_{i}(x_{j})=\unicode[STIX]{x1D6FF}_{ij}$ for any nodes $x_{j}$ . We write $u_{h}(x)=\sum _{j=1}^{N}\unicode[STIX]{x1D707}_{j}\unicode[STIX]{x1D719}_{j}(x)$ . Equation (2.19) is then equivalent to

which is a linear system of equations:

The matrix $A$ is known as the stiffness matrix of the nodal basis $\{\unicode[STIX]{x1D719}_{i}\}_{i=1}^{N}$ .

For $d=2$ and the special uniform triangulation shown in Figure 2.1, this stiffness matrix for the Laplacian operator turns out to be exactly the one given by (2.12). This special case is an example of the close relationship between finite difference and finite element methods.

We note that the finite element method is based on the variational formulation (2.3), whereas the finite difference method is not. In the development of the AMG method, however, the variational method is also used to derive coarse level equations for finite difference methods.

For any $T\in {\mathcal{T}}_{h}$ , we define

We say that the mesh ${\mathcal{T}}_{h}$ is shape-regular if there exists a uniformly bounded constant $\unicode[STIX]{x1D70E}\geq 1$ such that

We call $\unicode[STIX]{x1D70E}$ the shape-regularity constant.

Let $h=\max _{T\in {\mathcal{T}}_{h}}\overline{h}_{T}$ , with $\overline{h}_{T}$ defined as in (2.21). We say that the mesh ${\mathcal{T}}_{h}$ is quasi-uniform if there exists a uniformly bounded constant $C>0$ such that

2.3 Spectral properties

We now discuss the spectral properties of the partial differential operator ${\mathcal{L}}$ given in (2.1).

We recall the well-known Courant–Fischer min-max principle (Courant and Hilbert Reference Courant and Hilbert1924) for eigenvalues of symmetric matrices.

Next, we recall Theorem 1 in Fan (Reference Fan1949), which is known as the Ky–Fan trace-minimization principle.

Finally, following Xu (Reference Xu1992), we use the notation $a\lesssim b$ to represent the existence of a generic positive constant $C$ , which is independent of salient parameters, such as problem size or anisotropic ratio, and such that $a\leq Cb$ . Furthermore, we write $a{\displaystyle\mathop{=}\limits_{{\sim}}}b$ if and only if $a\lesssim b$ and $b\lesssim a$ .

We have the well-known Weyl estimate for the asymptotic behaviour of the Laplacian operator (Weyl Reference Weyl1911, Reference Weyl1912, Reed and Simon Reference Reed and Simon1978).

Next, we extend Weyl’s law to discretized PDE operators. The following theorem gives a discrete version of Weyl’s law for finite element discretization. Further details of this result can be found in Xu, Zhang and Zikatanov (Reference Xu, Zhang and Zikatanov2016b ).

2.4 Properties of finite element matrices

The main algebraic property for the stiffness matrices given by (2.20) is that it is sparse with $O(N)$ non-zeros, symmetric positive definite (for both Dirichlet and mixed boundary conditions) and semi-definite for pure Neumann boundary conditions. Its eigenvalues satisfy a discrete Weyl’s law.

For simplicity, we will only consider pure Dirichlet boundary conditions in the rest of this section.

We next discuss some more refined spectral properties of finite element stiffness matrices from uniform grids for the unit square domain $\unicode[STIX]{x1D6FA}=(0,1)\times (0,1)$ for $d=2$ . We begin with the Poisson equation. It is easy to derive a closed-form solution of the eigenpairs of $A$ given by (2.15), and we have

and

Now consider the important case of anisotropic problem (2.11). We order the vertices of the triangulation lexicographically, and, as before, denote them by $\{(ih,jh)\}_{i,j=0}^{n}$ . The stiffness matrix is then

Obviously,

and it is easily verified that

which leads to the expression

for the eigenvalues, and

for the corresponding eigenvectors.

3.1 Dual and inner product

For convenience of exposition, we will assume that $V$ is equipped with an inner product $(\cdot ,\cdot )$ . By the Riesz representation theorem, for any $f\in V^{\prime }$ , there is a unique $u\in V$ such that

It is via the Riesz representation theorem that we identify the dual space $V^{\prime }$ with $V$ . In the rest of this paper, for convenience, we will always assume that $V^{\prime }=V$ for any finite-dimensional vector space $V$ . As a result, we have

Thanks to the identification $V^{\prime }=V$ via (3.5), the identities in (3.6) are clear.

We would like to point out that, in an abstract discussion of all iterative methods for problem (3.1), it suffices to use the abstract dual pairing $\langle \cdot ,\cdot \rangle$ without having to introduce an inner product $(\cdot ,\cdot )$ on $V$ . However, we find that using an inner product is convenient for exposition, as we shall see later. We further point out that we will not use an inner product to identify $V^{\prime }=V$ for any infinite-dimensional vector space in this paper.

If $\{\unicode[STIX]{x1D719}_{i}\}_{i=1}^{N}$ is a basis of $V$ , we will always choose a basis $\{\unicode[STIX]{x1D713}_{i}\}_{i=1}^{N}$ of $V^{\prime }$ that is dual to the basis of $V$ , namely

Such a dual basis will be used only for theoretical considerations and not for the actual implementation of any algorithms.

We will consider two kinds of linear vector spaces. The first kind is $V=\mathbb{R}^{n}$ and the inner product is just the dot product

A canonical basis of $\mathbb{R}^{N}$ is formed by the column vectors of the identity matrix, $\{e_{i}\}_{i=1}^{N}$ . It is easy to see that the dual basis of $(\mathbb{R}^{N})^{\prime }=\mathbb{R}^{N}$ is just the original basis $\{e_{i}\}_{i=1}^{N}$ itself.

The second kind is a finite-dimensional functional subspace of $L^{2}(\unicode[STIX]{x1D6FA})$ for a given domain $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{d}$ ( $1\leq d\leq 3$ ), equipped with the $L^{2}$ -inner product

One commonly used linear vector space is a finite element space $V_{h}$ , and the nodal basis functions $\{\unicode[STIX]{x1D719}_{i}\}_{i=1}^{N}$ are often used as a basis. In this case, the dual basis functions $\{\unicode[STIX]{x1D713}_{i}\}_{i=1}^{N}$ of $V^{\prime }$ are no longer the original nodal basis functions but rather a set of functions (which are usually globally supported) that satisfy (3.7). This set of dual basis functions is usually needed to derive the matrix representation of operators between various spaces and their duals, but they are not needed for the actual implementation of relevant algorithms.

For a linear operator

its adjoint,

is defined as follows:

Since $V$ plays the role of both $V$ and its dual $V^{\prime }$ , (3.8) and (3.9) can have four different meanings:

1. (i) if $L:V\mapsto V$ , then $L^{\prime }:V^{\prime }\mapsto V^{\prime }$ ,

2. (ii) if $L:V\mapsto V^{\prime }$ , then $L^{\prime }:V\mapsto V^{\prime }$ ,

3. (iii) if $L:V^{\prime }\mapsto V$ , then $L^{\prime }:V^{\prime }\mapsto V$ ,

4. (iv) if $L:V^{\prime }\mapsto V^{\prime }$ , then $L^{\prime }:V\mapsto V$ .

Thanks to the identification we made between $V^{\prime }$ and $V$ via (3.5), definition (3.10) is applicable to all four cases.

If $V=\mathbb{R}^{n}$ and $(u,v)=(u,v)_{\ell ^{2}}$ , then $L^{\prime }=L^{T}$ , namely the matrix transpose. We say that an operator $A:V\mapsto V^{\prime }$ is symmetric positive definite (SPD) if

When $A$ is SPD, it defines another inner product $(\cdot ,\cdot )_{A}$ on $V$ ,

and a corresponding norm,

We use the superscript ‘*’ for the adjoint operator with respect to $(\cdot ,\cdot )_{A}$ , that is,

It is easy to see that

and $(BA)^{\ast }=BA$ if and only if $B^{\prime }=B$ .

3.2 Matrix representation

Let $V_{c}\subset V$ be a subspace and consider the inclusion operator $\imath _{c}:V_{c}\mapsto V$ . Assume that $\{\unicode[STIX]{x1D719}_{i}^{c}\}_{i=1}^{n_{c}}$ and $\{\unicode[STIX]{x1D719}_{i}\}_{i=1}^{n}$ are basis functions of $V_{c}$ and $V$ respectively, and the matrix representation of $\imath _{c}$ is the matrix

Equation (3.12) is shorthand for the expansion of the basis in $V_{c}$ via the basis in $V$ :

What is the matrix representation of $\imath _{c}^{\prime }:V^{\prime }\mapsto V_{c}^{\prime }$ ? Although we have $V^{\prime }=V$ and $V_{c}^{\prime }=V_{c}$ , we need to use dual bases $\{\unicode[STIX]{x1D713}_{i}^{c}\}_{i=1}^{n_{c}}\subset V_{c}^{\prime }$ and $\{\unicode[STIX]{x1D713}_{i}\}_{i=1}^{n}\subset V^{\prime }$ respectively. With respect to these dual bases, the matrix representation of $\imath _{c}^{\prime }$ is simply $P^{T}$ (the transpose of $P$ ), since it is easy to verify that

Now consider a linear operator

There are two different ways to get a matrix representation of $A$ , because $V$ plays two roles here: $V$ is $V$ itself, and $V=V^{\prime }$ . For the first case, we use the same basis $\{\unicode[STIX]{x1D719}_{i}\}$ for $V$ as the domain of $A$ and $V$ as the range of $A$ . In this case, the matrix representation of $A$ is the matrix

In the second case, we use the basis $\{\unicode[STIX]{x1D719}_{i}\}$ for $V$ as the domain space of $A$ , but use the dual basis $\{\unicode[STIX]{x1D713}_{i}\}$ for $V^{\prime }=V$ as the range space of $A$ . In this case, the matrix representation of $A$ is the matrix

It is easy to see that

and

The matrix $\tilde{A}$ in (3.17) is often called the stiffness matrix of $A$ and the matrix $M$ in (3.18) is called the mass matrix.

In the early multigrid literature, a discrete inner product equivalent to the $L^{2}$ -inner product was often introduced for finite element spaces so that the corresponding mass matrix becomes diagonal. However, if we view the underlying finite element operator (3.14) in a slightly different way,

we will then see easily that the introduction of the discrete $L^{2}$ -inner product is not necessary.

If $V=\mathbb{R}^{n}$ and we choose the canonical basis $\{e_{i}\}$ for $V$ , we would not encounter the mass matrix problem as in the functional space case, since in this case $\{e_{i}\}$ is also the dual basis of $V^{\prime }$ . This is certainly convenient, but the convenience tends to mask some subtle but important difference between various vectors and matrices in a given problem and the objects (functions) that they represent.

Given a matrix $A\in \mathbb{R}^{n\times n}$ , we can view it as either

or

As it turns out, when $A$ is obtained from the discretization of partial differential equations, (3.21) is more informative than (3.20). Hence we write a matrix equation

It is sometimes helpful to regard $x$ and $b$ as living in two ‘different’ spaces,

3.3 Eigenvalues and eigenvectors

Let us briefly discuss eigenvalues and eigenvectors for the symmetric operator $T:V\mapsto V$ . If $(\unicode[STIX]{x1D706},\unicode[STIX]{x1D719})$ is an eigenpair of $T$ ,

then it is easy to see that $(\unicode[STIX]{x1D706},\tilde{\unicode[STIX]{x1D719}})$ is an eigenpair of the matrix representation $\tilde{T}$ of $T$ :

Here $\tilde{\unicode[STIX]{x1D719}}\in \mathbb{R}^{n}$ is the vector representation of $\unicode[STIX]{x1D719}$ :

For an operator $A$ defined in (3.14), we need to be cautious when talking about eigenvalues of $A$ . Although we identify $V^{\prime }=V$ via (3.5), $V^{\prime }$ and $V$ play two different roles, so $A$ is essentially a mapping between two ‘different’ spaces $V$ and $V^{\prime }$ , and the spectrum of $A$ should be defined carefully. However, if we consider a symmetric operator

then $RA:V\mapsto V$ is an operator that is symmetric with respect to the $A$ -inner product. In this case, if $(\unicode[STIX]{x1D706},\unicode[STIX]{x1D719})$ is an eigenpair of $RA$ , then $(\unicode[STIX]{x1D706},\tilde{\unicode[STIX]{x1D719}})$ is an eigenpair of $\tilde{R}\tilde{A}$ (i.e. equal to the matrix representation of $RA$ ).

We consider the trivial identification operator

namely $\jmath v=v$ for all $v\in V^{\prime }=V$ . It is easy to see that the matrix representation of $\jmath$ is the inverse of the mass matrix $M=((\unicode[STIX]{x1D719}_{j},\unicode[STIX]{x1D719}_{i}))$ , namely

Using this identification operator, the operator $\jmath A:V\mapsto V$ is a symmetric operator from $V$ to $V$ . We can then talk about its spectrum. For example, if $(\unicode[STIX]{x1D706},\unicode[STIX]{x1D719})$ is an eigenpair of $\jmath A$ , then $(\unicode[STIX]{x1D706},\tilde{\unicode[STIX]{x1D719}})$ satisfies

This is often the generalized eigenvalue problem appearing in finite element analysis.

Although, for all $v\in V$ , $\jmath Av=Av$ because of the identification introduced above, $\jmath A$ and $A$ are, strictly speaking, two different operators: they have two different ranges and their matrix representations are different.

The discussions above, although simple, may seem a little confusing at first glance, but an unambiguous understanding and clarification of these concepts and the underlying subtleties will be helpful for the presentation of algebraic multigrid methods in the rest of this paper. For more detailed discussions on relevant topics, we refer to Xu (Reference Xu2016).

3.4 Bibliographical notes

For a general reading on the basic linear algebra materials used here, we refer to Halmos (Reference Halmos1974) and Xu (Reference Xu1992, Reference Xu2016). In particular, for a more detailed discussion related to dual spaces and matrix representations, we refer to Xu (Reference Xu1992, Reference Xu2016).

4.1 Basic iterative methods

A general linear iterative method for solving (4.1) can be written as follows. Given $u^{0}\in V$ , we define

where $B:V^{\prime }\mapsto V$ is a linear operator which can be thought of as an approximate inverse of $A$ .

Sometimes it is more desirable for the iterator $B$ to be symmetric. If $B$ is not symmetric, there is a natural way to symmetrize it. Consider the following iteration:

The symmetrized iteration (4.3) can be written as

where

which satisfies

Obviously, $\unicode[STIX]{x1D70C}(I-\bar{B}A)<0$ $\Longleftrightarrow \bar{B}>0$ $\Longleftrightarrow$ $G\equiv (B^{\prime })^{-1}+B^{-1}-A>0$ .

4.2 Jacobi and Gauss–Seidel methods

For $A=(a_{ij})\in \mathbb{R}^{n\times n}$ , we write

where $D$ is the diagonal of $A$ , and $L$ and $U$ are the strict lower and upper triangular parts of $A$ respectively.

Given $\unicode[STIX]{x1D714}>0$ , the (modified) Jacobi method can be written as (4.2) with

as shown in Algorithm 1. The (modified) Gauss–Seidel method can be written as (4.2) with

as shown in Algorithm 2.

The next result, which follows easily from Theorem 4.1, is well known.

In practice, it is often easy to choose $\unicode[STIX]{x1D714}$ to satisfy (4.9) so that the modified Jacobi method is guaranteed to converge. In the rest of this paper, we may always assume that such a choice of $\unicode[STIX]{x1D714}$ is made. For the Gauss–Seidel method, we will always choose $\unicode[STIX]{x1D714}=1$ (optimal successive over-relaxation is not generally used in multigrid methods). The Jacobi and Gauss–Seidel methods together with their convergence theory can be extended to block matrices in a straightforward fashion.

4.3 The subspace correction method

We consider a sequence of spaces $V_{1},\ldots ,V_{J}$ . These spaces, which will be known as auxiliary spaces, are not necessarily subspaces of $V$ , but each of them is related to the original space $V$ by a linear operator

Our very basic assumption is that the following decomposition holds:

This means that for any $v\in V$ there exists $v_{i}\in V_{i}$ (which may not be unique) such that

Furthermore, we assume that each $V_{i}$ is equipped with an energy inner product $a_{i}(\cdot ,\cdot )$ . We define

by

Let $\unicode[STIX]{x1D6F1}_{i}^{\prime }:V^{\prime }\mapsto V_{i}^{\prime }$ be the adjoint of $\unicode[STIX]{x1D6F1}_{i}$ :

Let $P_{i}=\unicode[STIX]{x1D6F1}_{i}^{\ast }:V\mapsto V_{i}$ be the adjoint of $\unicode[STIX]{x1D6F1}_{i}$ with respect to the A-inner products:

The following identity holds:

If $u$ is the solution of (4.1), by (4.14), we have

where

This equation may be regarded as the restriction of (4.1) to $V_{i}$ . We assume that each such $A_{i}$ has an approximate inverse or preconditioner:

The parallel subspace correction (PSC for short) method is (4.2) with $B=B_{\operatorname{psc}}$ given by

The successive subspace correction (SSC for short) method is defined in Algorithm 3.

The algorithm is equivalent to (4.2) with $B=B_{\operatorname{ssc}}$ given by

where

The Jacobi and Gauss–Seidel method can be interpreted as PSC and SSC based on the decomposition

with exact subspace solves such that

By Theorem 4.5,

By Lemma 4.6,

We note that

where

In the above presentation, most results are for SPD problems. We would like to point out that all of these results can be extended to a more general class of problems, namely symmetric positive semi-definite problems. When $A$ is a matrix, we further assume that all the diagonals of $A$ are non-zero and hence positive. When the subspace correction method is used for a more general symmetric positive semi-definite operator $A$ , we further assume that each $A_{i}$ is SPD. We shall use the acronym SPSD to denote matrices or operators that satisfy the aforementioned properties.

4.4 Bibliographical notes

The general notion of subspace corrections by means of space decompositions was described by Xu (Reference Xu1992), based on Bramble, Pasciak, Wang and Xu (Reference Bramble, Pasciak, Wang and Xu1991b , Reference Bramble, Pasciak, Wang and Xu1991c ). It is an abstract point of view encompassing the theory and practice of a large class of iterative algorithms such as multigrid and domain decomposition methods. In the past two decades a great deal of effort has been put into the investigation of theoretical and practical issues related to these methods. General results, applicable in many cases, in the theory of additive and multiplicative methods in Hilbert space can be found in Xu and Zikatanov (Reference Xu and Zikatanov2002). For a literature review and basic results we refer the reader to the monographs and survey articles by Hackbusch (Reference Hackbusch1985), Bramble (Reference Bramble1993), Vassilevski (Reference Vassilevski2008), Xu (Reference Xu1989, Reference Xu and Ainsworth1997), Xu and Zou (Reference Xu and Zou1998), Yserentant (Reference Yserentant1993), Toselli and Widlund (Reference Toselli and Widlund2005), Griebel and Oswald (Reference Griebel and Oswald1995) and Smith, Bjørstad and Gropp (Reference Smith, Bjørstad and Gropp1996). For detailed studies of classical iterative methods, we refer to the monographs by Young (Reference Young1971), Hackbusch (Reference Hackbusch1994), Varga (Reference Varga2000) and Saad (Reference Saad2003).

Note that in this section we have considered SPSD matrices, and according to Lee, Wu, Xu and Zikatanov (Reference Lee, Wu, Xu and Zikatanov2007, Reference Lee, Wu, Xu and Zikatanov2008) and Ayuso de Dios, Brezzi, Marini, Xu and Zikatanov (Reference Ayuso de Dios, Brezzi, Marini, Xu and Zikatanov2014), all the results in this section are valid for SPSD problems with semi-norms. Relations between the auxiliary space method and subspace correction methods is drawn in Chen (Reference Chen2011). In the classical multigrid literature (Brandt et al. Reference Brandt, Mccormick and Ruge1982, Brandt, McCormick and Ruge Reference Brandt, Mccormick and Ruge1985, Ruge and Stüben Reference Ruge, Stüben and McCormick1987, Trottenberg, Oosterlee and Schüller Reference Trottenberg, Oosterlee and Schüller2001), the notions of algebraically smooth (low) frequencies and algebraic high frequencies play an important role. They are also instrumental in the design of new AMG methods. As indicated by the convergence estimates, for a given smoother, the desirable coarse spaces should capture or well approximate the lower end of the spectrum of the relaxed matrix $\bar{R}A$ or $D^{-1}A$ . This is usually referred to as near-null space (Treister and Yavneh Reference Treister and Yavneh2015, Lai and Olson Reference Lai and Olson2011, Xu Reference Xu2009, Brezina et al. Reference Brezina, Falgout, Maclachlan, Manteuffel, Mccormick and Ruge2006).

5.1 A two-level method

A two-level method typically consists of the following components:

1. (i) a smoother, $R:V^{\prime }\mapsto V$ ,

2. (ii) a coarse space $V_{c}$ , which may or may not be a subspace of $V$ , linked to $V$ via a prolongation operator,

   $$\begin{eqnarray}P:V_{c}\mapsto V,\end{eqnarray}$$

3. (iii) a coarse space solver, $B_{c}:V_{c}^{\prime }\mapsto V_{c}$ .

In the discussion below we need the inner product

and the accompanying norm $\Vert \cdot \Vert _{\bar{R}^{-1}}$ . Here we recall that the definition of $\bar{R}$ is analogous to that in (4.5).

We always assume that $\bar{R}$ is SPD, and hence the smoother $R$ is always convergent. Furthermore,

The restriction of (4.1) is then

where

The coarse space solver $B_{c}$ is often chosen to be the exact solver, namely $B_{c}=A_{c}^{-1}$ , for analysis, but in a multilevel setting, $B_{c}$ is recursively defined and it is an approximation to $A_{c}^{-1}$ . We distinguish two different cases in choosing $B_{c}$ :

When $A$ is semi-definite, we use $N$ to denote the kernel of $A$ and we always assume that $N\subset V_{c}$ . Let

where the orthogonality is understood with respect to the $(\cdot ,\cdot )_{\bar{R}^{-1}}$ -inner product. Let $Q_{1}:V\mapsto W$ be the orthogonal projection with respect to the $(\cdot ,\cdot )_{\bar{R}^{-1}}$ -inner product

Here $A_{c}$ is semi-definite on $V_{c}$ but invertible on $W_{c}$ . We denote the restriction of $A_{c}$ on $W_{c}$ by $\hat{A}_{c}$ , and define the pseudo-inverse of $A_{c}$ by

With a slight abuse of notation, we will still use $A_{c}^{-1}$ to denote the pseudo-inverse of $A_{c}$ , namely

We will use similar notation for the pseudo-inverse of other relevant singular operators and matrices in the rest of the paper.

We choose to define an AMG algorithm in terms of an operator $B:V^{\prime }\mapsto V$ , which can be considered as an approximate inverse or a preconditioner of $A$ . A typical two-level MG method is shown in Algorithm 4.

In the rest of this section, we take $B_{c}=A_{c}^{-1}$ .

There are usually two different (and mathematically equivalent) ways to choose $V_{c}$ , $V$ , $P$ and $R$ . The first, known as the operator version, is such that

In this case $P=\imath _{c}$ , where

is the natural inclusion of $V_{c}$ into $V$ . In the application to finite element discretization for second-order elliptic boundary value problems, $V_{c}$ and $V$ are just the finite element subspaces of $H^{1}(\unicode[STIX]{x1D6FA})$ . This type of notation is convenient for analysis. However, this is not an algorithm that can be directly used for implementation.

The second, known as the matrix version, is such that

and

is the prolongation matrix.

These two different notations are related via the use of basis functions of $\{\unicode[STIX]{x1D719}_{i}^{c}\}_{i=1}^{n_{c}}\subset V_{c}$ and $\{\unicode[STIX]{x1D719}_{i}\}_{i=1}^{n}\subset V$ . As noted earlier, the prolongation matrix $P$ given in (3.12)–(3.13) is simply the matrix representation of $\imath _{c}$ given in (5.9), and we have

The following observation is clear.

5.2 An optimal two-level AMG theory

The aim of the AMG method is to balance the interplay between smoother $R$ and the coarse space $V_{c}$ . Most existing AMG methods first fix a smoother, which is often given by the Jacobi or Gauss–Seidel methods (or their combinations and variations), and then optimize the choice of coarse spaces. This is the approach we will mainly discuss in this paper. However, we also comment on a different approach by first fixing the coarse space and then trying to optimize the choice of the smoother. It is also possible to try to make an optimal choice of smoother and coarse space simultaneously, but we will not address this approach here.

Let $Q_{c}:V\mapsto V_{c}$ be the orthogonal projection with respect to the $(\cdot ,\cdot )_{\bar{R}^{-1}}$ -inner product

By the definition of $W$ and $W_{c}$ , we have that $(\cdot ,\cdot )_{A}$ is an inner product on $W$ , $\Vert \cdot \Vert _{A}$ is a norm on $W$ , and the projection $\unicode[STIX]{x1D6F1}_{c}:V\mapsto W_{c}$ is well defined:

The two-level convergence rate is obtained in the following theorem.

Theorem 5.3 can be stated as follows using the matrix representation introduced in Section 3.

In the case of two subspaces, we have the following theorem giving the precise convergence rate of the corresponding SSC method.

Using the matrix representation introduced in Section 3, we state the matrix version of Theorem 5.6 below. For simplicity, with abuse of notation, we still use $A$ to denote the matrix representation of operator $A$ .

The following theorem is important in motivating most AMG algorithms.

Since the eigenvalues of $\bar{R}A$ are expensive to compute, the practical value of Theorem 5.6 is limited. However, it provides useful guidance in the design of practical AMG methods.

For finite element discretizations we can use Weyl’s law combined with Theorem 5.6 to prove an estimate for the convergence rate of a two-grid method with an optimal coarse space.

Theorem 5.3 provides an explicit estimate for the convergence of a two-level method in terms of $K(V_{c})$ . For a given method, a smaller bound on $K(V_{c})$ means a faster convergence rate. In particular, the two-level AMG method is uniformly convergent if $K(V_{c})$ is uniformly bounded with respect to mesh parameters.

5.3 Quasi-optimal theories

We now look at the necessary and sufficient condition for uniform convergence of a two-level method as proved in Section 5.2 (see Theorem 5.3):

Here $K(V_{c})$ is the smallest constant for which (5.28) holds for all $v\in V$ . The space $V_{c}$ which minimizes $K(V_{c})$ is $V_{c}^{\text{opt}}$ . A similar argument was used in the proof of Theorem 5.6 in Section 5.2. Here we generalize the result to semi-definite  $A$ .

For a given smoother $R$ , one basic strategy in the design of AMG is to find a coarse space such that $K(V_{c})$ is as small as possible in practice. There are many cases, however, in which the operator $\bar{R}^{-1}$ in the definition of $K(V_{c})$ is difficult to work with.

One commonly used approach is to replace $\bar{R}^{-1}$ with a simpler but spectrally equivalent SPD operator. More specifically, we assume that $D:V\mapsto V^{\prime }$ is an SPD operator such that

That is,

where

Examples of such equivalent norms for Schwarz smoothers are given in (4.20) and (4.27). As a rule, the norm defined by $\bar{R}$ corresponding to the symmetric Gauss–Seidel method, that is, $R$ defined by the pointwise Gauss–Seidel method, can be replaced with the norm defined by the diagonal of $A$ (i.e., by the Jacobi method, which, while not always convergent as a relaxation, provides an equivalent norm).

In terms of this operator $D$ , we introduce the quantity

where $Q_{D}:V\mapsto V_{c}$ is the $(u,v)_{D}$ -orthogonal projection.

By (5.16), (5.31) and (5.30), we have

The proof of Theorem 5.11 is straightforward and indicates that, if $c_{D}$ and $c^{D}$ are ‘uniform’ constants, the convergence rate of the two-level method is ‘uniformly’ dictated by the quantity $K(V_{c},D)$ .

We say that $V_{c}$ is quasi-optimal if the following inequality holds:

with a constant $\unicode[STIX]{x1D6FE}>0$ independent of the size of the problem.

The construction of an approximation to the optimal coarse space $V_{c}^{\text{opt}}$ that is used in most AMG algorithms relies on two operators $A_{M}$ and $D_{M}$ which satisfy

with constants $c_{1}$ and $c_{2}$ independent of the problem size. Here, on the right-hand side, we have a seminorm $\Vert \cdot \Vert _{A_{M}}$ , because sometimes $A_{M}$ is only semi-definite. We point out that here $A_{M}$ and $D_{M}$ are analogues of ${\displaystyle\mathop{A}\limits_{{\sim}}}_{W}$ and ${\displaystyle\mathop{D}\limits_{{\sim}}}$ defined in (6.6) and (6.7), respectively, in the general framework of Section 6; the assumptions in (5.36) are analogous to Assumption 6.4 which we made in the general AMG framework in Section 6.

5.4 Algebraically high and low frequencies

In geometric MG, algebraically smooth error is also smooth in the usual geometric sense. However, in AMG settings, smooth error can be geometrically non-smooth. In order to make this distinction, we use the term algebraically smooth error when we refer to the error in the AMG setting that is not damped (eliminated) by the smoother $R$ . In general, good interpretation of the algebraically smooth error leads to an efficient and robust AMG algorithm. Careful characterization of the algebraically smooth error is needed, since in this case we can try to construct a coarser level which captures these error components well.

Here is a more formal definition of an algebraically smooth error.

The set of algebraically smooth vectors will be denoted by

We point out that this is a set of vectors (a ball, or rather an ellipsoid) and not a linear vector space in general. It is then clear that the elements of this set need to be approximated well by elements from the coarse space.

The rationale of Definition 5.13 can be seen from the following simple result (it is always true that $(\bar{R}Av,v)\leq \Vert v\Vert _{A}^{2}$ ).

We can easily show that this definition is equivalent to saying that the algebraically smooth error components are the components for which the smoother converges slowly. Indeed, inequality (5.42) is clearly equivalent to $((I-\overline{R}A)v,v)_{A}\geq (1-\unicode[STIX]{x1D700})(v,v)_{A}$ , namely

The property (5.43) is often referred to as the smoothing property.

Thanks to (7.9), we have the following result.

On the other hand, note that by the definition of $\Vert \cdot \Vert _{\bar{R}^{-1}}$ , we always have $\Vert v\Vert _{\bar{R}^{-1}}\geq C\Vert v\Vert _{A}$ , with constant $C$ independent of the parameters of interest. Drawing from analogy with the geometric multigrid method, we introduce the notion of algebraic high frequency as follows.

The set of algebraically high-frequency vectors will be denoted by

The concept of algebraic high and low frequencies will be used in the two-level AMG theory in Section 5.3, and also in the design of classical AMG in Section 12.1.

We now introduce the notion of near-null space.

5.5 Smoothing properties of Jacobi and Gauss–Seidel methods

The essence of multigrid methods is that simple iterative methods such as Jacobi and Gauss–Seidel methods have a special property known as the smoothing property. As an illustration, we apply the Gauss–Seidel method to

with $A$ given by (2.16) for isotropic problems and (2.34) for anisotropic problems, respectively. We first choose $\unicode[STIX]{x1D707}$ randomly (as shown in Figure 5.1) and then compute $A\unicode[STIX]{x1D707}$ for both (2.16) and (2.34) to compute right-hand sides $b=A\unicode[STIX]{x1D707}$ respectively. We then apply the Gauss–Seidel method to both equations with initial guess $\unicode[STIX]{x1D707}^{0}=0$ (see Figure 5.2).

Figure 5.1.

Initial error for both (2.16) and (2.34).

Figure 5.2.

(a) Error after applying five Gauss–Seidel iterations to (2.16). (b) Error after applying five Gauss–Seidel iterations to (2.34). (c) Error after applying one block Gauss–Seidel iteration to (2.34).

Note that, for $A$ given by (2.34) when $\unicode[STIX]{x1D716}\ll 1$ , we have

The corresponding eigenfunctions, which can be viewed as ‘algebraic low frequencies’, can be highly oscillatory in the $x$ -direction.

As an illustration of the difference between algebraic high/low frequencies and geometric high/low frequencies, we consider the linear system given by (2.34) for anisotropic problems. Clearly, $A$ can be written as

We define the vector $\unicode[STIX]{x1D707}\in R^{N}$ as

Then it is easy to compute that

We have

and

Letting $D$ be the diagonal of $A$ , we then have

This shows that

which implies that $\unicode[STIX]{x1D707}$ is an algebraic low frequency if $\unicode[STIX]{x1D716}$ is sufficiently small.

On the other hand, if we denote the nodal basis functions corresponding to the uniform finite element mesh by $\{\unicode[STIX]{x1D719}_{ij}:1\leq i,j\leq n\}$ , that is, $\unicode[STIX]{x1D719}_{ij}$ is a piecewise linear function such that

then we define

If we consider the finite element function corresponding to $\unicode[STIX]{x1D707}$ , namely, the function defined by

then, in the geometric point of view, this function is highly oscillatory in the $x$ -direction, which is a geometric high frequency (see Figure 5.3).

Figure 5.3.

Graph of function defined in (5.49). This function is highly oscillatory in the $x$ -direction.

5.6 Convergence theory in view of algebraic high and low frequencies

We next present a convergence theory based on algebraic high and low frequencies. We first prove the following lemma.

5.7 Bibliographical notes

One of the first results on two-level convergence of AMG methods can be found in papers by Brandt et al. (Reference Brandt, Mccormick and Ruge1982) and Ruge and Stüben (Reference Ruge, Stüben and McCormick1987). There has been much research on extending the MG theory to algebraic settings (Maitre and Musy Reference Maitre and Musy1983, Bank and Douglas Reference Bank and Douglas1985, Mandel Reference Mandel1988) and the algebraic variational approach to the two-level MG theory (McCormick Reference Mccormick1985, Reference Mccormick1984, McCormick and Ruge Reference Mccormick and Ruge1982).

For two-grid convergence, sharper results, including two-sided bounds, are given by Zikatanov (Reference Zikatanov2008) and also considered by Falgout and Vassilevski (Reference Falgout and Vassilevski2004) and Falgout, Vassilevski and Zikatanov (Reference Falgout, Vassilevski and Zikatanov2005). These two-level results are more or less a direct consequence of the abstract theory provided by Bramble et al. (Reference Bramble, Pasciak, Wang and Xu1991c ), Xu (Reference Xu1992) and Xu and Zikatanov (Reference Xu and Zikatanov2002). A survey of these and other related results can be found in a recent article by MacLachlan and Olson (Reference Maclachlan and Olson2014). The two approaches included in this section were recently developed by Xu, Zhang and Zikatanov (Reference Xu, Zhang and Zikatanov2016a , Reference Xu, Zhang and Zikatanov2016c ).

Theorem 5.3 can be found in Zikatanov (Reference Zikatanov2008), and can be viewed as a consequence of the XZ identity (Xu and Zikatanov Reference Xu and Zikatanov2002) in the special case of two subspaces from the general framework of the subspace correction method. The original proof of this theorem in Zikatanov (Reference Zikatanov2008) was based on the XZ identity. The proof here is new and more direct.

Multilevel results are difficult to establish in general algebraic settings, and most of them are either based on unrealistic assumptions or they use geometric grids to prove convergence. We refer to Vaněk, Mandel and Brezina (Reference Vaněk, Mandel and Brezina1996) and Brezina and Vassilevski (Reference Brezina, Vassilevski, Lirkov, Margenov and Wasniewski2011) for results in this direction. Rigorous multilevel results for finite element equations can be derived using the auxiliary space framework, which was developed by Xu (Reference Xu1996) for quasi-uniform meshes. More recently, Chen, Nochetto and Xu (Reference Chen, Nochetto and Xu2012) showed that multilevel convergence results for adaptively refined grids are optimal. A multilevel convergence result on shape-regular grids using AMG based on quad-tree coarsening (in two dimensions) and oct-tree coarsening (in three dimensions) is shown in Grasedyck, Wang and Xu (Reference Grasedyck, Wang and Xu2015).

Finally, we point out that the notation used in parts of this section originates in Bank and Dupont (Reference Bank and Dupont1980), Bramble and Pasciak (Reference Bramble and Pasciak1987) and Bramble, Pasciak and Xu (Reference Bramble, Pasciak and Xu1990), and is convenient for the analysis, especially when finite element equations are considered.

6.1 Bibliographical notes

Matsokin and Nepomnyashchikh (Reference Matsokin and Nepomnyashchikh1985) first proved the fictitious space lemma. The work of Xu (Reference Xu1996) on the auxiliary space method is related. An additive version of Lemma 6.1 can be found in Xu and Zikatanov (Reference Xu and Zikatanov2002), and the most general case (including multiplicative preconditioners) is found in Xu (Reference Xu2016).

In later sections we show how the general theory used here can be applied to various AMG algorithms, for example classical AMG (Brandt et al. Reference Brandt, Mccormick and Ruge1982, Ruge and Stüben Reference Ruge, Stüben and McCormick1987), smoothed aggregation AMG (Míka and Vaněk Reference Míka and Vaněk1992a , Reference Míka and Vaněk1992b ), spectral AMGe (Chartier et al. Reference Chartier, Falgout, Henson, Jones, Manteuffel, Mccormick, Ruge and Vassilevski2003, Efendiev, Galvis and Vassilevski Reference Efendiev, Galvis, Vassilevski, Huang, Kornhuber, Widlund and Xu2011), and other algorithms.

Many of the works in the domain decomposition (DD) literature also use techniques for defining coarse spaces, which, to a large extent, have similar aims to the coarse space constructions for AMG outlined in this section. We refer to Toselli and Widlund (Reference Toselli and Widlund2005), Widlund (Reference Widlund1994, Reference Widlund2009), Dohrmann, Klawonn and Widlund (Reference Dohrmann, Klawonn and Widlund2008), Spillane et al. (Reference Spillane, Dolean, Hauret, Nataf, Pechstein and Scheichl2014) and the references therein for more details on using local eigenspaces for constructing coarse spaces in DD methods.

7.1 Sparse matrix and its adjacency graph

An undirected graph (or simply a graph) ${\mathcal{G}}$ is a pair $({\mathcal{V}},{\mathcal{E}})$ , where ${\mathcal{V}}$ is a finite set of points called vertices, and ${\mathcal{E}}$ is a finite set of pairs of vertices, known as edges. We often write ${\mathcal{V}}=\{1,\ldots ,n\}$ for some fixed $n$ . We will not consider directed graphs in this paper because the graphs corresponding to the symmetric sparse matrices are undirected.

An edge $e\in {\mathcal{E}}$ is an unordered pair $(j,\,k)$ , where $j,\,k\in {\mathcal{V}}$ . The vertices $j$ and $k$ are said to be adjacent if $(j,k)\in {\mathcal{E}}$ . A path from a vertex $j$ to a vertex $k$ is a sequence $(j_{0},\,j_{1},\,j_{2},\ldots ,\,j_{l})$ of vertices where $j_{0}=j$ , $j_{l}=k$ , and $(j_{i},\,j_{i+1})\in {\mathcal{E}}$ for all $i=0,\,1,\ldots ,\,l-1$ . A vertex $j$ is connected to a vertex $k$ if there is a path from $j$ to $k$ . ${\mathcal{G}}=({\mathcal{V}},\,{\mathcal{E}})$ is connected if every pair of vertices is connected by a path; otherwise it is said to be disconnected. A graph ${\mathcal{G}}_{0}=({\mathcal{V}}_{0},{\mathcal{E}}_{0})$ is called a subgraph of ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$ if ${\mathcal{V}}_{0}\subset {\mathcal{V}}$ and ${\mathcal{E}}_{0}\subset {\mathcal{E}}$ .

The neighbourhood $N(i)$ is the set of vertices adjacent to the vertex $i$ . The degree or valency of a vertex is the number of edges that connect to it. These are defined as

A path connecting two vertices $i$ and $j$ is a sequence of edges $(k_{0},k_{1}),(k_{1},k_{2})$ , $\ldots ,(k_{m-1},k_{m})$ in ${\mathcal{E}}$ such that $k_{0}=i$ and $k_{m}=j$ . The length of the path is the number of edges in it. The distance between two vertices $i$ and $j$ is the length of the shortest path connecting $i$ and $j$ , and we denote it by $\operatorname{dist}(i,j)$ . If $i$ , $j$ are not connected, then $\operatorname{dist}(i,j)=\infty$ . The diameter of a graph is the largest distance between two vertices, that is, $\operatorname{diam}({\mathcal{G}})=\max _{(i,j)\in {\mathcal{E}}}\operatorname{dist}(i,j)$ . An independent set is a set of vertices in which no two are adjacent. A maximal independent set is an independent set such that adding any other vertex to the set forces the set to contain an edge.

Given a symmetric matrix $A\in \mathbb{R}^{n\times n}$ , the adjacency graph of $A$ is an undirected graph, denoted by ${\mathcal{G}}(A)$ , ${\mathcal{G}}=({\mathcal{V}},\,{\mathcal{E}})$ with ${\mathcal{V}}=\{1,\,2,\ldots ,\,n\}$ . The edges ${\mathcal{E}}$ are defined as

A matrix $A$ is called irreducible if its adjacency graph ${\mathcal{G}}(A)=({\mathcal{V}},{\mathcal{E}})$ is connected. Otherwise, $A$ is called reducible.

An example of a symmetric matrix is shown in Figure 7.1(a) and the corresponding graph is shown in Figure 7.1(b). Pictorial representation of a graph is often not available, and a graph can be drawn in different ways with different coordinates of the vertices. As a general rule, sparse matrices do not provide any geometric information for the underlying graph and only the combinatorial or topological properties of ${\mathcal{G}}(A)$ .

Figure 7.1.

(a) A sparse symmetric matrix, and (b) its associated graph.

Given

we define

We say that $X$ has sparsity pattern given by $S$ if and only if $X\in \mathbb{R}_{S}^{n\times m}$ .

Often, the sparsity pattern of a matrix is determined in advance and the set $S$ is determined by a given matrix. For $Y\in \mathbb{R}^{n\times m}$ , we let

We now consider the graphs associated with finite element or finite difference stiffness matrices. In the case of finite elements, this is the space of finite element functions $V_{h}$ , and for finite difference discretizations this is the space of mesh-functions $V_{h}$ , which can be identified with $\mathbb{R}^{N}$ .

We assume that we have a finite-dimensional space $V_{h}^{N}$ which we use to discretize the Neumann problem. We also have the finite element space for the Dirichlet (or mixed boundary condition) problem, and we assume that the following inclusions hold: $V_{h}=V_{h}^{D}\subset V_{h}^{N}$ . Equivalently, we have a subspace of the variables which vanish on $V_{h}$ : for example, the values of the finite element or finite difference solution at the nodes on the Dirichlet boundary vanish.

Let $A^{{\mathcal{N}}}$ be the matrix corresponding to the finite element or finite difference discretization of the model second-order elliptic equation with Neumann boundary conditions. Clearly we have the identity

Here, the sum is over all edges ${\mathcal{E}}$ of the graph ${\mathcal{G}}(A^{{\mathcal{N}}})=({\mathcal{V}},{\mathcal{E}})$ , $\unicode[STIX]{x1D6FF}_{e}v=v_{i}-v_{j}$ , if ${\mathcal{E}}\ni e=(i,j)$ , $i<j$ . Also, $\unicode[STIX]{x1D714}_{e}=-(A^{{\mathcal{N}}})_{ij}$ are the off-diagonal entries of $A^{{\mathcal{N}}}$ . Note that since we consider the Neumann problem, the bilinear form defined by $A^{{\mathcal{N}}}$ vanishes for $u$ (resp.  $v$ ) such that $u_{i}=1$ (resp.  $v_{i}=1$ ) for all $i$ . For both 5-point and 9-point stencils we have that $\unicode[STIX]{x1D714}_{e}=1$ for all $e$ .

We consider the stiffness matrix $A^{{\mathcal{N}}}$ , corresponding to the model problem (2.1) with Neumann boundary conditions on a bounded domain $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{d}$ , that is, we have the boundary condition

where $\boldsymbol{n}$ is the unit normal vector to $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ pointing outward. It is easy to derive the stiffness matrices corresponding to the Dirichlet or mixed boundary condition problem: we just restrict the bilinear form defined by $A^{{\mathcal{N}}}$ to the subspace given by

7.2 $M$ -matrix relatives of finite element stiffness matrices

A symmetric matrix $A\in \mathbb{R}^{n\times n}$ is called an $M$ -matrix if it satisfies the following three properties:

As the first step in creating a space hierarchy, most AMG algorithms for $Au=f$ with positive semi-definite $A$ use a simple filtering of the entries of  $A$ , and construct an $M$ -matrix which is then used to define crucial AMG components.

A few remarks are in order.

1. (i) We have used the term $M$ -matrix to denote semi-definite matrices, and we are aware that this is not the precise definition. It is, however, much more convenient to refer to $M$ -matrices, and we decided to relax the definition here in the hope that the inaccuracy would pay off by appealing to the reader.

2. (ii) We point out that the restricted $M$ -matrix relatives are instrumental in the definition of coarse spaces and also in the convergence rate estimates. This is clearly seen later in Section 5.2 where we present the unified two-level theory for AMG.

3. (iii) Often, the case is that the one-sided inequality in (7.9) is in fact a spectral equivalence.

By definition, we have the following simple but important result.

This result means that we only need to work on the $M$ -matrix relative of $A$ in order to get the estimate (7.10).

In this section we show how to construct the $M$ -matrix relative to the matrix resulting from a finite element discretization of the model problem (2.1) with linear elements. We first consider an isotropic problem with Neumann boundary condition (7.4) and isotropic $\unicode[STIX]{x1D6FC}=a(x)I$ . Construction of $M$ -matrix relatives in the case of anisotropic tensor $\unicode[STIX]{x1D6FC}(x)$ in (2.1) is postponed to Section 14.2.

In the rest of this section we make the following assumptions on the coefficient and the geometry of $\unicode[STIX]{x1D6FA}$ .

• ∙ The domain $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{d}$ is partitioned into simplices $\unicode[STIX]{x1D6FA}=\cup _{T\in {\mathcal{T}}_{h}}T$ .

• ∙ The coefficient $a(x)$ is a scalar-valued function, and its discontinuities are aligned with the partition ${\mathcal{T}}_{h}$ .

• ∙ We consider the Neumann problem, and hence the bilinear form (2.3) is

  (7.12) $$\begin{eqnarray}\int _{\unicode[STIX]{x1D6FA}}a(x)\unicode[STIX]{x1D6FB}v\cdot \unicode[STIX]{x1D6FB}u=\mathop{\sum }_{(i,j)\in {\mathcal{E}}}(-a_{ij})\unicode[STIX]{x1D6FF}_{e}u\unicode[STIX]{x1D6FF}_{e}v=\mathop{\sum }_{e\in {\mathcal{E}}}\unicode[STIX]{x1D714}_{e}\unicode[STIX]{x1D6FF}_{e}u\unicode[STIX]{x1D6FF}_{e}v.\end{eqnarray}$$

• ∙ It is well known that the off-diagonal entries of the stiffness matrix $A$ are given by

  $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D714}_{e} & = & \displaystyle -(\unicode[STIX]{x1D719}_{j},\unicode[STIX]{x1D719}_{i})_{A}=\mathop{\sum }_{T\supset e}\unicode[STIX]{x1D714}_{e,T},\nonumber\\ \displaystyle \displaystyle \unicode[STIX]{x1D714}_{e,T} & = & \displaystyle \displaystyle \frac{1}{d(d-1)}\overline{a}_{T}|\unicode[STIX]{x1D705}_{e,T}|\cot \unicode[STIX]{x1D6FC}_{e,T},\quad \overline{a}_{T}=\displaystyle \frac{1}{|T|}\int _{T}a(x)\text{d}x.\nonumber\end{eqnarray}$$
  Here, $e=(i,j)$ is a fixed edge with end-points $x_{i}$ and $x_{j}$ , $T\supset e$ is the set of all elements containing $e$ , $|\unicode[STIX]{x1D705}_{e,T}|$ is the volume of the $(d-2)$ -dimensional simplex opposite to $e$ in $T$ , and $\unicode[STIX]{x1D6FC}_{e,T}$ is the dihedral angle between the two faces in $T$ not containing $e$ .

• ∙ Let ${\mathcal{E}}$ denote the set of edges in the graph defined by the triangulation and let ${\mathcal{E}}^{-}$ be the set of edges where $a_{ij}\geq 0$ , $i\neq j$ . The set complementary to ${\mathcal{E}}^{-}$ is ${\mathcal{E}}^{+}={\mathcal{E}}\setminus {\mathcal{E}}^{-}$ . Then, with $\unicode[STIX]{x1D714}_{e}=-a_{ij}$ and $\unicode[STIX]{x1D6FF}_{e}u=(u_{i}-u_{j})$ , $e=(i,j)$ , we have

  (7.13) $$\begin{eqnarray}\int _{\unicode[STIX]{x1D6FA}}a(x)\unicode[STIX]{x1D6FB}v\cdot \unicode[STIX]{x1D6FB}u=\mathop{\sum }_{e\in {\mathcal{E}}^{+}}\unicode[STIX]{x1D714}_{e}\unicode[STIX]{x1D6FF}_{e}u\unicode[STIX]{x1D6FF}_{e}v-\mathop{\sum }_{e\in {\mathcal{E}}^{-}}|\unicode[STIX]{x1D714}_{e}|\unicode[STIX]{x1D6FF}_{e}u\unicode[STIX]{x1D6FF}_{e}v.\end{eqnarray}$$

• ∙ We also assume that the partitioning is such that the constant function is the only function in the null space of the bilinear form (7.12). This is, of course, the case when $\unicode[STIX]{x1D6FA}$ is connected (which is true, as $\unicode[STIX]{x1D6FA}$ is a domain).

The non-zero off-diagonal entries of $A$ may have either positive or negative sign, and usually ${\mathcal{E}}^{-}\neq \emptyset$ . The next theorem shows that the stiffness matrix $A$ defined via the bilinear form (7.12) is spectrally equivalent to the matrix $A_{+}$ defined as

Thus, we can ignore any positive off-diagonal entries in $A$ , or equivalently we may drop all $\unicode[STIX]{x1D714}_{e}$ for $e\in {\mathcal{E}}^{-}$ . Indeed, $A_{+}$ is obtained from $A$ by adding to the diagonal all positive off-diagonal elements and setting the corresponding off-diagonal elements to zero. This is a stronger result that we will need later, because it gives spectral equivalence with the $M$ -matrix relative of  $A_{+}$ .

A simple corollary, which we will use later to prove estimates for the convergence rate, is as follows.

Corollary 7.5 together with Lemma 7.3 provides a theoretical foundation for using $M$ -matrix relatives to design AMG for finite element matrices.

7.3 Bibliographical notes

We have introduced some standard notions from graph theory. For the reader interested in more detailed descriptions, we refer to the classical textbooks by Diestel (Reference Diestel2010) and Gibbons (Reference Gibbons1985) as a general introduction to graph theory, and to Saad (Reference Saad2003) and Varga (Reference Varga2000) for considerations linking graphs, sparse matrices and iterative methods.

Our results on $M$ -matrix relatives are related to the some of the works on preconditioning by $Z$ -matrices and $L$ -matrices (Kraus and Schicho Reference Kraus and Schicho2006, Kraus Reference Kraus2008). They are implicitly used in most of the AMG literature (Ruge and Stüben Reference Ruge, Stüben and McCormick1987) where the classical connection strength definition gives an $M$ -matrix. We point out that the $M$ -matrix property and the existence of an $M$ -matrix relative is often not sufficient to achieve even a two-level uniform convergence of AMG. A typical example is a matrix which has been rescaled and the constant is no longer in the kernel of the discrete operator. In this case, the standard AMG application may fail, and the near kernel needs to be recovered by different means, such as the adaptive AMG processes considered in Section 15 and the references given therein.

8.1 Basic idea and strength function

If a subset of vertices on which an algebraically smooth vector, say $v$ , changes very slowly, we only need to keep one degree of freedom to represent $v$ in this subset. In other words, we can either aggregate this subset together or keep one vertex and delete the rest of the vertices in this subset. We say that the vertices in this subset are strongly connected to each other. Connection strength is a concept introduced to identify strongly connected pairs of vertices. Roughly speaking, we say that $i$ and $j$ are strongly connected if $v_{i}\approx v_{j}$ .

We imagine coarsening the graph in two stages: the first step is to remove some edges and the second step is to remove some vertices. The second step is the goal. Let us examine how the second step is carried out: we either (i) aggregate some neighbouring vertices together or (ii) pick a maximal independent set, denoted by ${\mathcal{C}}$ , in the filtered graph and then remove all the remaining vertices. Using the argument above, (i) each aggregate should only consist of strongly connected vertices, or (ii) each of the deleted vertices should be strongly connected to some point in ${\mathcal{C}}$ . To guarantee either of these two situations, we then have to remove all the weakly connected edges in the first step.

Let us further use some heuristic arguments to motivate how the connection strength should be defined. Let $v$ be an algebraically smooth (5.40), namely

Let $u=v/\Vert v\Vert _{\bar{R}^{-1}}$ . We then have

Thanks to Lemma 5.16, we can assume that $A$ is an $M$ -matrix, namely $-a_{ij}=|a_{ij}|$ . From (8.1), we have the following observations.

1. (i) A larger $|a_{ij}|$ means a smaller $(u_{i}-u_{j})^{2}$ .

2. (ii) An algebraically smooth error varies more slowly in the direction where $|a_{ij}|$ is larger.

This observation leads to the following definition of connection strength. Given a threshold $\unicode[STIX]{x1D703}>0$ , we say that the vertex $j$ of the adjacency graph is $\unicode[STIX]{x1D703}$ -strongly connected to vertex $i$ if

Note that by the definition above we may have $j$ strongly connected to $i$ while $i$ is not strongly connected to $j$ . As a result, the adjacency graph corresponding to the matrix of strong connections may not be symmetric.

However, our theoretical framework is given in terms of the symmetrized operator $\bar{R}$ , regardless of whether the original smoother $R$ is symmetric. This is due to the fact that we used the energy norm, namely the $A$ -norm, to measure the convergence rate, and the resulting convergence rate is given in terms of $\bar{R}$ . This is the best convergence theory we have, and we will use this theory to study AMG algorithms. As a result, we will only consider strength functions that are symmetric. A strength function

associated with an SPSD matrix is symmetric if $s_{c}(i,j)=s_{c}(j,i)$ .

Given a threshold $\unicode[STIX]{x1D703}>0$ , we say that $i$ and $j$ are $\unicode[STIX]{x1D703}$ -strongly connected if

We then define the strength matrix:

Note that $S$ is a Boolean matrix with entries equal to 0 or 1 depending on the connection strength.

Given the non-overlapping decomposition

we extend the definition of strength function to

and we assume that $s_{c}$ is symmetric.

Given a threshold $\unicode[STIX]{x1D703}>0$ , we say that ${\mathcal{A}}_{i}$ and ${\mathcal{A}}_{j}$ are $\unicode[STIX]{x1D703}$ -strongly connected to each other if

An example of such a strength function is given in (8.13).

In the AMG literature, a number of heuristics have been proposed for identifying strong connections, particularly when considering discretizations of anisotropic equations. In general, connection strength is a notion that is difficult to address theoretically or to relate to the convergence rate of an algorithm. We refer to the classical papers and monographs mentioned in Section 9.6 for further discussions on related issues. Current trends in AMG development aim to re-evaluate the role of the classical definition of connection strength.

We would finally like to comment that connection strength is used to define the sparsity of $P$ , and it is crucial, for example, in the proof of the convergence of the two-level method for discretizations of elliptic equations with jump coefficients. Choosing the ‘right’ sparsity of $P$ is crucial, as a denser $P$ would lead to a better approximation from a coarser space, and a sparser $P$ would lead to a less expensive algorithm.

In the rest of the section, we discuss different definitions of connection strength:

1. (i) classical AMG,

2. (ii) lean AMG,

3. (iii) local-optimization based AMG.

8.2 Classical AMG

With the above motivation, we define the strength function as follows:

Definition (8.7) is a symmetrized version of the strength function used in the classical AMG literature (see (8.2)).

The following definition is also commonly used in classical AMG algorithms:

Again, this is a symmetrized version of strength functions used in the AMG literature.

Finally, we may also have the following two definitions, which are based on Cauchy–Schwarz for SPSD matrices:

and

Note that definition (8.7) (which is mostly associated with classical AMG) ignores all the non-negative entries of the stiffness matrix $A=(a_{ij})$ .

8.3 Local-optimization strength function

Suppose that each pair of the index $\{i,j\}\subset \{1,\ldots ,n\}$ is associated with a space $V_{ij}$ , which is not necessarily a subspace of $V$ . Assume now that we have two operators: $A_{ij}:V_{ij}\mapsto V_{ij}^{\prime }$ , which is symmetric positive semi-definite, and $D_{ij}:V_{ij}\mapsto V_{ij}^{\prime }$ , which is symmetric positive definite.

Given a number $k_{ij}<\dim V_{ij}$ , we define a coarse space $V_{ij}^{c}\subset V_{ij}$ as

where $\unicode[STIX]{x1D701}_{ij}^{(k)}$ is the eigenvector corresponding to the $k$ th smallest eigenvalue of $D_{ij}^{-1}A_{ij}$ .

We let $Q_{ij}:V_{ij}\mapsto V_{ij}^{c}$ denote the orthogonal projection with respect to $(\cdot ,\cdot )_{D_{ij}}$ . Motivated by (6.13), we define the strength function $s_{c}$ as follows:

A special case of the above definition is introduced in aggregation-based algebraic multigrid (AGMG). Suppose now that we have a set of aggregates $\{{\mathcal{A}}_{1},\ldots ,{\mathcal{A}}_{J}\}$ . We fix a pair $\{i,j\}\subset \{1,\ldots ,J\}$ , and define

We then let $V_{ij}$ define the restriction of $V$ to $G$ , that is,

where

We use $A_{ij}$ and $D_{ij}$ , respectively, to denote the restriction of $A$ and $D$ to  $G$ . We then choose $k_{ij}=1$ and define the local coarse space $V_{ij}^{c}$ by

and the orthogonal projection $Q_{ij}:V_{ij}\mapsto V_{ij}^{C}$ with respect to $(\cdot ,\cdot )_{D_{ij}}$ , namely,

The strength function based on aggregation is defined as

Another example is choosing $V_{ij}=\mathbb{R}^{2}$ and

We choose $k_{ij}=1$ and define the coarse space $V_{ij}^{c}\subset V_{ij}$ as

By direct computation, the strength function defined in (8.11) is given by

We note that $s_{1}\geq s_{2}$ and hence

We point out that the strength function given by (8.14) is obtained by using the theory in Section 6, while the other strength functions such as (8.7) are obtained by heuristic considerations.

8.4 Lean AMG

Instead of using the absolute value of matrix entries as the criterion to determine if two points are strongly coupled, lean AMG uses affinity to measure the connection strength, which is based on the following heuristic observation: Given a vector $v$ , if $(i,j)$ is a strong connected pair of vertices, after several relaxation on $v$ , namely

the values of $v_{i}$ and $v_{j}$ should be close.

In lean AMG, we generate $K$ test vectors. Each test vector is the result of applying $\unicode[STIX]{x1D708}$ Gauss–Seidel relaxation sweeps to $Ax=0$ , starting from randomly generated vectors $x^{(1)},\ldots ,x^{(K)}\in \mathbb{R}^{n}$ . Further, we denote

Here $X_{i}^{T}$ is the $i$ th row of $X_{n\times K}$ . The strength function for lean AMG is then defined as

8.5 Bibliographical notes

Classical algorithms for determining connection strength can be found in Brandt et al. (Reference Brandt, Mccormick and Ruge1982), Stüben (Reference Stüben1983), Brandt et al. (Reference Brandt, Mccormick and Ruge1985), Ruge and Stüben (Reference Ruge, Stüben and McCormick1987) and Briggs, Henson and McCormick (Reference Briggs, Henson and Mccormick2000). The original measure of connection strength given by Ruge (Reference Ruge1983, Reference Ruge1985), McCormick and Ruge (Reference McCormick and Ruge1989) and Brandt et al. (Reference Brandt, Mccormick and Ruge1985) is non-symmetric, but for theoretical considerations, which only depend on the symmetrized smoothers, it suffices to use the slightly more restrictive but symmetric versions of connection strength.

Some extensions of these classical algorithms for defining strong connections are based on different measures for connectivity and distance such as measure of importance and algebraic distance. Details can be found in Ruge and Stüben (Reference Ruge, Stüben and McCormick1987) and Trottenberg et al. (Reference Trottenberg, Oosterlee and Schüller2001, Appendix A).

Connection strength functions have had little theoretical backing in the past. The results developed in this section, such as local Poincaré inequalities and especially the strength functions (8.15), show that such heuristics are reasonable, and their choice can be motivated by theoretical results.

For aggregation AMG, typically a symmetric connection strength function is used, as defined by Vaněk et al. (Reference Vaněk, Mandel and Brezina1996). Some recent aggregation algorithms also define connection strength based on sharp theoretical results, and use the local two-level convergence rate (8.11) as a measure (Notay Reference Notay2010, Napov and Notay Reference Napov and Notay2012).

For aggregations based on matching (aggregates of size $2$ ), the ‘heavy edge’ matching algorithms in Karypis and Kumar (Reference Karypis and Kumar1998) correspond to a connection strength function selecting aggregates depending on edge weights in the adjacency graph. Some recent works (e.g. Livne and Brandt Reference Livne and Brandt2012) use connection strength functions based on the size of the entries in the Gram matrix formed by a set of smoothed test vectors.