The Role of Energy Minimization in Algebraic Multigrid Interpolation

TL;DR

This paper reviews algebraic multigrid convergence theory, proposes a weighted energy minimization approach for interpolation, and explores how energy minimization influences AMG performance through theoretical and numerical analysis.

Contribution

It introduces a general weighted energy minimization functional for AMG interpolation and discusses its practical implications and limitations.

Findings

Weighted energy minimization offers new insights into AMG interpolation.

Numerical results highlight the role of energy minimization and constraint vectors.

The proposed preconditioner aids in solving Sylvester- or Lyapunov-type equations.

Abstract

Algebraic multigrid (AMG) methods are powerful solvers with linear or near-linear computational complexity for certain classes of linear systems, Ax=b. Broadening the scope of problems that AMG can effectively solve requires the development of improved interpolation operators. Such development is often based on AMG convergence theory. However, convergence theory in AMG tends to have a disconnect with AMG in practice due to the practical constraints of (i) maintaining matrix sparsity in transfer and coarse-grid operators, and (ii) retaining linear complexity in the setup and solve phase. This paper presents a review of fundamental results in AMG convergence theory, followed by a discussion on how these results can be used to motivate interpolation operators in practice. A general weighted energy minimization functional is then proposed to form interpolation operators, and a novel `diagonal' preconditioner for Sylvester- or Lyapunov-type equations developed simultaneously. Although results based on the weighted energy minimization typically underperform compared to a fully constrained energy minimization, numerical results provide new insight into the role of energy minimization and constraint vectors in AMG interpolation.

Equations

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+M^{-1}(\mathbf{b}-A{\mathbf{x}}_k),$$

$$A=\left[\begin{array}{cc}{A}_{ff}& A_{fc}\\ A_{cf}& A_{cc}\end{array}\right].$$

$$E_{TG}=(I-M^{-T}A)(I-{\pi }_A)(I-M^{-1}A).$$

$$c_1{\mathbf{v}}^{T}X\mathbf{v}\le {\mathbf{v}}^T\tilde{M}\mathbf{v}\le c_2{\mathbf{v}}^{T}X\mathbf{v}.$$

$$\Vert (I-{\pi }_X)\mathbf{v}{\Vert }_X^2$$

$$\Vert (I-{\pi }_{\tilde{M}})\mathbf{v}{\Vert }_{\tilde{M}}^2$$

$$⟨A,B{⟩}_F$$

$$\mathrm{tr}(ABC)$$

$$\Vert E_{TG}{\Vert }_A=\Vert (I-{\pi }_A)(I-M^{-1}A){\Vert }_A^2=\Vert (I-M^{-T}A)(I-{\pi }_A){\Vert }_A^2.$$

$$\Vert (I-M^{-T}A)\mathbf{v}{\Vert }_A^2\le \Vert \mathbf{v}{\Vert }_A^2-\delta \Vert (I-{\pi }_A)\mathbf{v}{\Vert }_A^2,$$

$$\Vert (I-M^{-T}A)(I-{\pi }_A)\mathbf{v}{\Vert }_A^2$$

$$\le \Vert (I-{\pi }_A)\mathbf{v}{\Vert }_A^2-\alpha /\beta \Vert (I-{\pi }_A{)}^2\mathbf{v}{\Vert }_A^2$$

$$\le (1-\alpha /\beta )\Vert \mathbf{v}{\Vert }_A^2.$$

$$\Vert \mathbf{v}-P{\mathbf{v}}_c{\Vert }_A^2\le \beta \Vert A\mathbf{v}{\Vert }_X^2.$$

$$\Vert (I-{\pi }_A)\mathbf{v}{\Vert }_A^2=⟨A(I-{\pi }_A)\mathbf{v},(I-{\pi }_A)\mathbf{v}-P{\mathbf{v}}_c⟩\le \Vert A(I-{\pi }_A)\mathbf{v}{\Vert }_X\Vert (I-{\pi }_A)\mathbf{v}-P{\mathbf{v}}_c{\Vert }_{X^{-1}}.$$

$$\Vert \mathbf{v}-P{\mathbf{v}}_c{\Vert }_{X^{-1}}^2\le \beta \Vert \mathbf{v}{\Vert }_A^2.$$

$$\frac{\Vert \mathbf{v}-P{\mathbf{v}}_c{\Vert }_{\tilde{M}}^2}{\Vert \mathbf{v}{\Vert }_A^2}\le K.$$

$$K_{TG}=\underset{\mathbf{v}\ne 0}{\max }\frac{\Vert (I-{\pi }_{\tilde{M}})\mathbf{v}{\Vert }_{\tilde{M}}^2}{\Vert \mathbf{v}{\Vert }_A^2},$$

$$c_1\underset{\mathbf{v}\ne 0}{\max }\frac{\Vert (I-{\pi }_X)\mathbf{v}{\Vert }_X^2}{\Vert \mathbf{v}{\Vert }_A^2}$$

$$\le K_{TG}\le \underset{\mathbf{v}\ne 0}{\max }\frac{\Vert (I-{\pi }_X)\mathbf{v}{\Vert }_{\tilde{M}}^2}{\Vert \mathbf{v}{\Vert }_A^2}\le c_2\underset{\mathbf{v}\ne 0}{\max }\frac{\Vert (I-{\pi }_X)\mathbf{v}{\Vert }_X^2}{\Vert \mathbf{v}{\Vert }_A^2}.$$

$${\lambda }_{\text{min}}(\tilde{M})\underset{\mathbf{v}\ne 0}{\max }\frac{\Vert (I-Q_P)\mathbf{v}{\Vert }^2}{\Vert \mathbf{v}{\Vert }_A^2}\le K_{TG}\le {\lambda }_{\text{max}}(\tilde{M})\underset{\mathbf{v}\ne 0}{\max }\frac{\Vert (I-Q_P)\mathbf{v}{\Vert }^2}{\Vert \mathbf{v}{\Vert }_A^2},$$

$$\Vert (I-Q_P)\mathbf{v}{\Vert }^2\le \frac{K_{TG}}{\Vert A\Vert }\Vert \mathbf{v}{\Vert }_A^2.$$

$${\kappa }_s\le {\lambda }_{\text{min}}(X_s^{-1}{A}_s)\le {\lambda }_{\text{max}}(X_s^{-1}{A}_s)\le c_2.$$

$$\Vert A(I-{\pi }_A)\mathbf{v}{\Vert }_X\le \Vert I-{\pi }_A^T{\Vert }_X\Vert A\mathbf{v}{\Vert }_X,$$

$$\Vert E_{MG}{\Vert }_A=1-\frac{1}{K_{MG}},$$

$${\Vert \left(I-{\left(M^{(k)}\right)}^{-T}{A}^{(k)}\right){\mathbf{v}}^{(k)}\Vert }_{A^{(k)}}^2\le {\Vert {\mathbf{v}}^{(k)}\Vert }_{A^{(k)}}^2-\frac{\alpha }{\Vert A^{(k)}\Vert }{\Vert A^{(k)}{\mathbf{v}}^{(k)}\Vert }^2,$$

$${\Vert {\mathbf{v}}^{(k)}-P^{(k)}{\mathbf{v}}^{(k+1)}\Vert }_{A^{(k)}}^2\le \frac{\beta }{\Vert A^{(k)}\Vert }{\Vert A^{(k)}{\mathbf{v}}^{(k)}\Vert }^2,$$

$${\Vert {\mathbf{v}}^{(k)}-P^{(k)}{\mathbf{v}}^{(k+1)}\Vert }^2\le \frac{\beta }{\Vert A^{(k)}\Vert }{\Vert {\mathbf{v}}^{(k)}\Vert }_{A^{(k)}}^2.$$

$${\Vert (I-{\pi }_k){\mathbf{v}}^{(k)}\Vert }^2\le \frac{\beta }{\Vert A^{(k)}{\Vert }^2}{\Vert A^{(k)}{v}^{(k)}\Vert }^2,$$

$${\Vert (I-Q_k){\mathbf{v}}^{(k)}\Vert }^2$$

$${\Vert {\mathbf{v}}^{(k)}-P^{(k)}{\hat{\mathbf{v}}}^{(k+1)}\Vert }_{A^{(k)}}^2\le \Vert A^{(k)}\Vert {\Vert {\mathbf{v}}^{(k)}-P^{(k)}{\hat{\mathbf{v}}}^{(k+1)}\Vert }^2\le \frac{\beta }{\Vert A^{(k)}\Vert }{\Vert A^{(k)}{v}^{(k)}\Vert }^2.$$

$$\underset{P}{\min }\underset{\mathbf{v}}{\max }\frac{\Vert (I-PR)\mathbf{v}{\Vert }^2}{\Vert \mathbf{v}{\Vert }_A^2}=\left[\begin{array}{c}-A_{ff}^{-1}{A}_{fc}\\ I\end{array}\right]:=P_{ideal},$$

$$[A_{ff}^{-1}{]}_{ij}\le C{q}^{|i-j{|}_G-1},$$

Full text

{AMS}

AMS subject classifications

1 Introduction

Algebraic multigrid (AMG) was designed as a solver for large, sparse linear systems, typically M-matrices, resulting from the discretization of elliptic PDEs. For many such problems, AMG has been shown to achieve fast convergence, and scale in parallel to hundreds of thousands of processors [1]. Such convergence and scaling properties are desirable for solvers, and substantial work has been devoted to broadening the applicability of AMG. Continued research has made for a rich theoretical basis for AMG [15, 16, 31, 32, 19], as well as many numerical implementations and variations that are either robust for a larger class of linear systems [9, 10, 21], or effective at solving a specific problem such as linear elasticity [2] or Hemholtz [25]. Nevertheless, the one-size-fits-all AMG solver remains elusive, in part because many of the theoretical results are difficult to use in a practical setting.

A novel feature of AMG in contrast to many other linear solvers is that the setup and solve complexity in terms of floating point operations (FLOPs) are both typically linear or log-linear in the total number of degrees of freedom (DOFs). This is fundamental to good scaling of time to solution with increasing problem size, but also limits the options in algorithm design, particularly when trying to directly use theoretical results on convergence. Two common aspects seen in AMG convergence theory are the use of orthogonal projections onto subspaces and requiring a given approximation property to hold for all vectors. In both cases, namely constructing an orthogonal projection or enforcing a constraint for $n$ basis vectors, the complexity of explicitly enforcing such requirements is at least quadratic in $n$ and, thus, not feasible in keeping with the desired linear complexity of AMG. Furthermore, AMG can only have linear complexity when all operators are sparse, including the coarse-grid operators constructed for a multilevel algorithm. In the abstract setting of convergence theory, such sparsity constraints are not accounted for, adding an additional barrier to the direct use of convergence theory in practical methods.

Here, we review the tension between theory and practice in AMG and propose a new variant of AMG that aims to directly address these complications. An overview of AMG convergence theory is given in Section 2. Fundamental results on two-grid and multigrid convergence theory are presented in a simple and consistent manner, to clarify what is required of interpolation operators for effective AMG convergence, and the so-called “optimal” and “ideal” interpolation operators are introduced. Section 3 proposes a discussion of AMG interpolation operators used in practice, and how they relate to theoretical results, along with an examination of how different theoretical results can be approximated in linear complexity. This leads to the introduction of a general weighted functional to be minimized in forming interpolation operators in Section 4, which is shown to have a unique solution for a fixed interpolation sparsity pattern. A conjugate gradient method is developed to approximate the solution, with a novel preconditioner that is applicable to general equations with a Sylvester- or Lyapunov-like form (Section 5). Numerical results demonstrate that a constrained energy minimization [8, 20, 26, 33] consistently outperforms a weighted energy minimization. Although this may seem intuitive for AMG researchers, a number of other interesting results also come up that lead to open questions on interpolation in AMG:

• •

Enforcing one constraint vector to be (almost) exactly in the range of interpolation is fundamental to good AMG convergence. However, adding additional constraint vectors that are not effectively reduced by the current AMG hierarchy does not necessarily improve convergence, which is at odds with motivation of traditional adaptive AMG methods.

• •

Energy-minimization applied to columns of $P$ (while maintaining constraints) is also fundamental to a convergent AMG method for some more difficult model problems. However, although further iterations of energy minimization continue to reduce the associated residual, convergence of the resulting AMG solver does not improve after a small number of iterations.

• •

Using a diagonal preconditioner for energy-minimization iterations applied to columns of $P$ can offer significant improvement in convergence of the resulting AMG solver.

2 Theoretical framework

Multilevel solvers come in various forms, including geometric multigrid (GMG), AMG, finite element algebraic multigrid (AMGe), algebraic multilevel iterations (AMLI), and the method of subspace corrections. In this work, we focus on AMG as a general method to solve a linear system $A\mathbf{x}=\mathbf{b}$ using only “algebraic” information, in contrast to GMG and AMGe, which require additional information, on the underlying grid or finite element stiffness matrices, respectively. Subspace corrections are presented in a more general yet framework than AMG, but analysis of subspace correction can also be applied to AMG [31].

The basis for AMG as an iterative method to solve $A\mathbf{x}=\mathbf{b}$ is in reducing error through two processes: “relaxation” and “coarse-grid correction.” If designed properly, these processes are complementary in the sense that they are effective on different error modes and, together, effectively reduce all types of error. Relaxation refers to a general iterative method of the form

[TABLE]

and is often chosen to be a simple method such as Jacobi or Gauss-Seidel. This process is typically efficient at removing “high-frequency error,” or error associated with large eigenvalues of $A$. Convergence of a relaxation scheme in the “energy norm” or “$A$-norm,” $\|\mathbf{v}\|_{A}^{2}=\langle A\mathbf{x},\mathbf{x}\rangle$, is equivalent to bounding the error-propagation matrix in the $A$-norm, namely $\|I-M^{-1}A\|_{A}<1$. Furthermore, $\|I-M^{-1}A\|_{A}<1$ if and only if $M+M^{T}-A$ is symmetric positive definite (SPD) [32, Theorem 2.3.1]. In this case, we say that $M$ is an $A$-convergent relaxation operator.

Multigrid originated in the geometric setting, where high-frequency error actually has a high physical frequency. Since standard relaxation schemes such as Jacobi and Gauss-Seidel are able to capture this error well, the natural way to capture the converse, low-frequency error, is to recursively coarsen the underlying grid so that low-frequency modes on the initial fine grid appear high-frequency on coarser grids. Relaxation on coarser grids will then reduce this error, and the results can be interpolated back to the fine grid. The algebraic concept is much the same but low-frequency refers to algebraically smooth modes, corresponding to large eigenvalues of $I-M^{-1}A$ or, typically, small eigenvalues of $A$, and vice-versa for high-frequency. Because there is no explicit grid in the algebraic setting, algebraic coarsening is based on choosing a coarse subspace which can capture algebraically smooth error from the fine grid.

For an $n\times n$ SPD matrix $A$, consider an $\ell^{2}$-orthogonal decomposition of $\mathbb{R}^{n}$, where any $\mathbf{x}$ can be decomposed as $\mathbf{x}=R^{T}\mathbf{y}_{c}+S\mathbf{y}_{f}$, where $RS=0$. Here, $S$ corresponds to the space on which relaxation is effective, and $R$ defines the coarse space on which a coarse-grid correction is constructed. Defining the interpolation operator, $P$, we assume that $PR$ is a projection, which requires $RP=I$. A Galerkin coarse-grid operator is formed, $A_{c}=P^{T}AP$, and an exact coarse-grid correction (in the $A$-norm) given by the $A$-orthogonal projection onto $\operatorname{Range}(P)$, $\pi_{A}=P(P^{T}AP)^{-1}P^{T}A$. In this work, a CF-style splitting will be used (or in the case of aggregation-based coarsening, a root-node approach, where one node in each aggregate is declared a C-point and the rest F-points [21]). A CF-splitting has the useful property that coarse-grid nodes are a subset of nodes on the current grid, allowing for $A$ to be written in the block form

[TABLE]

In this case, splitting operators take the form $R=\begin{bmatrix}0&I\end{bmatrix}$, $S=\begin{bmatrix}I&0\end{bmatrix}^{T}$, and $P=\begin{bmatrix}W^{T}&I\end{bmatrix}^{T}$, where $RAR^{T}=A_{cc}$ and $S^{T}AS=A_{ff}$. Together, the two-grid error-propagation matrix operator for AMG, with an $A$-symmetric relaxation scheme based on $M^{-1}$ and $M^{-T}$, is given by

[TABLE]

Convergence of $E_{TG}$ is generally considered in the $A$-norm, where each iteration reduces error in the $A$-norm by at least a factor of $\|E_{TG}\|_{A}$. Noting that $E_{TG}$ is symmetric in the $A$-norm, it follows that eigenvectors of $E_{TG}$ are $A$-orthogonal and $\|E_{TG}\|_{A}=\rho(E_{TG})$. Thus, optimizing the AMG convergence rate can be viewed equivalently as minimizing $\|E_{TG}\|_{A}$ or $\rho(E_{TG})$. Two-grid convergence can also be considered in terms of the spectral equivalence between a preconditioner, $B_{TG}$, and $A$, where $E_{TG}=I-B_{TG}^{-1}A$ [32, Proposition 5.1.2]; however, here we bound convergence in terms of $\|E_{TG}\|_{A}$ for consistency. A multilevel method is implemented and analyzed as a two-grid method with an inexact coarse-grid solve, where the coarse-grid “solve” recursively calls a two-grid method on the coarse-grid problem. Multigrid convergence is also considered in the $A$-norm, where we want to bound $\|E_{MG}\|_{A}\leq K<1$.

Let $M$ and $M^{T}$ be $A$-convergent relaxation operators and define the symmetrized relaxation operator as $\widetilde{M}=M^{T}(M+M^{T}-A)^{-1}M$, so that $I-\widetilde{M}^{-1}A=(I-M^{-1}A)(I-M^{-T}A)$. This symmetrizes the action of $M$ and is used primarily as a theoretical tool (as $\widetilde{M}$ is rarely easily computable). Common bounds on two-grid convergence factors come from considering various orthogonal projections onto the range of the interpolation operator, $P$. Define $\pi_{X}:=P(P^{T}XP)^{-1}P^{T}X$ as the unique $X$-orthogonal projection onto $\operatorname{Range}(P)$ for some nonsingular operator $X$, e.g. $A$ or $\widetilde{M}$, and $Q_{P}:=P(P^{T}P)^{-1}P$ as the $l^{2}$-orthogonal projection onto $\operatorname{Range}(P)$. For computable bounds, assume that $X$ is spectrally equivalent to $\widetilde{M}$, denoted $X\simeq\widetilde{M}$; that is, there exists $0<c_{1}\leq c_{2}$ such that

[TABLE]

It follows from (1) and the definition of orthogonal projections that

[TABLE]

for all $\mathbf{v}$.

Finally, note the following identities with respect to the Frobenius inner product and trace that are used regularly in this work:

[TABLE]

and let $A\circ B$ denote the Hadamard product, defined as the element-wise multiplication of two matrices $A,B\in\mathbb{R}^{m\times n}$.

2.1 Two-grid convergence

Substantial work has been devoted to understanding convergence theory of AMG in the two-grid setting [18, 24, 32, 31, 15, 16, 19, 35]. By Lemma 4.1 of [22], we can analyze $\|E_{TG}\|_{A}$ either directly, or by considering the variants with only pre- or post-relaxation, as

[TABLE]

One of the simplest two-grid convergence bounds is given by Lemma 2.3 of [23]:

Theorem 2.1**.**

If there is a $\delta>0$ such that

[TABLE]

for all $\mathbf{v}$, then $\|E_{TG}\|_{A}\leq 1-\delta$.

Following [27, 19], sufficient conditions for two-grid convergence are given in the following theorem

Theorem 2.2**.**

Let symmetric and positive-definite matrix $X$ be given, and assume that there exist $\alpha,\beta>0$ such that $\|(I-M^{-T}A)\mathbf{v}\|_{A}^{2}\leq\|\mathbf{v}\|_{A}^{2}-\alpha\|A\mathbf{v}\|_{X}^{2}$ and $\|(I-\pi_{A})\mathbf{v}\|_{A}^{2}\leq\beta\|A\mathbf{v}\|_{X}^{2}$ for all $\mathbf{v}$. Then, $\|E_{TG}\|_{A}^{2}\leq 1-\alpha/\beta$.

Proof 2.3**.**

We prove the bound on $\|E_{TG}\|_{A}$ by proving the corresponding bound on $\|(I-M^{-T}A)(I-\pi_{A})\|_{A}$. For any $\mathbf{v}$,

[TABLE]

The first assumption on Theorem 2.2 is commonly referred to as the smoothing property, since it assumes that relaxation effectively reduces the error in an approximation when the residual associated with that error is large (when measured in the $X$-norm). The second assumption on Theorem 2.2 is referred to as the strong approximation property, since it assumes that coarse-grid correction is effective at reducing errors when the associated residuals are small; this is equivalent to assuming that such errors are well-approximated within $\operatorname{Range}(P)$. This assumption is termed the strong approximation property as it can clearly be replaced by a weaker one, that $\|(I-\pi_{A})\mathbf{v}\|_{A}^{2}\leq\beta\|A(I-\pi_{A})\mathbf{v}\|_{X}^{2}$ for all vectors, $\mathbf{v}$, stating that coarse-grid correction is effective at reducing errors for which the residual after coarse-grid correction is small. This latter assumption is commonly referred to as the weak approximation property. The difference between the weak and strong approximation properties comes up in the multilevel setting, and is discussed in Section 2.2.

In practice, the weak and strong approximation properties are typically considered in slightly altered forms. For the strong approximation property, an equivalent statement is that for any $\mathbf{v}$, there exists a $\mathbf{v}_{c}$ such that

[TABLE]

While a similar equivalence could be derived for the weak approximation property, a more typical bound arises by noting that, for any $\mathbf{v}$ and $\mathbf{v}_{c}$,

[TABLE]

Thus, a sufficient condition for the weak approximation property to hold is that for any $\mathbf{v}$ there exists a $\mathbf{v}_{c}$ such that

[TABLE]

Note that this is trivially true for $\mathbf{v}\in\operatorname{Range}(P)$, but the weak approximation property is implied by this condition for $\mathbf{v}\in\operatorname{Range}(P)^{\perp}$. It is in this form that much of the recent two-grid AMG theory has been developed.

Of particular interest is the result obtained taking $X=\widetilde{M}^{-1}$, so that the smoothing property in Theorem 2.2 trivially holds with $\alpha=1$. In this case, the two-grid convergence bound in Theorem 2.2 is determined entirely by the constant in the weak approximation property. Indeed, in this setting a sharp bound on convergence is possible [16, 32].

Theorem 2.4** (Weak approximation property).**

Let $A$ be SPD, $\widetilde{M}=M^{T}(M+M^{T}-A)^{-1}M$ for some relaxation scheme $M$, and $P$ the interpolation operator for a two-grid method. Suppose $\exists$ $K$ such that for any $\mathbf{v}\neq 0$, there exists a $\mathbf{v}_{c}$ such that

[TABLE]

Then the two-grid method converges uniformly, and $\|E_{TG}\|_{A}\leq 1-\frac{1}{K}$. Furthermore, the best (minimal) constant $K$ over all $P$ is given by

[TABLE]

in which case $\|E_{TG}\|_{A}=1-\frac{1}{K_{TG}}$.

Equation (5) gives a sharp bound on two-grid convergence, but can be generalized to any matrix $X$ and corresponding $X$-orthogonal projection onto $\operatorname{Range}(P),\pi_{X}$. Let $X$ be spectrally equivalent to $\widetilde{M}$ as in equation (1). Then, from equations (2) and (3),

[TABLE]

In the case of $X=I$, (6) simplifies to considering interpolation error in the $l^{2}$-norm [35]

[TABLE]

motivating the often-used simpler form of the WAP,

[TABLE]

The necessarily complementary role of relaxation and coarse-grid correction in AMG is accounted for in the WAP by requiring interpolation accuracy with respect to $\widetilde{M}$, that is the coarse-grid correction must account for low-eigenvalue modes of $\widetilde{M}$, which are not effectively reduced through relaxation with $\widetilde{M}$. A bound on two-grid convergence can also be formulated as two independent constraints based on coarse-grid selection and bounding the energy in the range of $P$ as follows [15, 32].

Lemma 2.5** (Two-grid energy-stability).**

Let $X$ be spectrally equivalent to $\widetilde{M}$ as in (1), and define $A_{s}=S^{T}AS,X_{s}=S^{T}XS$, where

[TABLE]

If $PR$ is bounded in energy, $\|PR\|_{A}^{2}\leq C$ for some $C$, a WAP in the $X$-norm, with projection $PR$ is satisfied. These are sufficient conditions for uniform two-grid convergence, with $K_{TG}\leq\frac{c_{2}}{\kappa_{s}}\|PR\|_{A}^{2}$.

Lemma 2.5 can be seen as an energy-stability constraint coupled with a compatibility measure of the fine and coarse grids. Equation (7) measures how well relaxation, $\widetilde{M}$, or the spectrally equivalent $X$, captures information about the fine-grid operator, $A_{s}$. This is based on the idea of compatible relaxation [4, 7], which ensures that the relaxation scheme is able to effectively reduce error on the fine grid. Then, assuming a compatible choice of grids, interpolation must be stable in energy, where $\|PR\mathbf{v}\|_{A}^{2}\leq C\|\mathbf{v}\|_{A}^{2}$. Note, Lemma 2.5 is equivalent to Theorem 2.4, however the differing explicit conditions make for different approaches to constructing multigrid hierarchies. One important difference is that in this case, the lemma is formulated in terms of a general operator norm as opposed to a constraint for all $\mathbf{v}$. This distinction is discussed in more detail in Section 3.

2.2 Multilevel convergence

Now that two-grid convergence theory has been introduced, let us continue by considering when two-level convergence can be extended to the multilevel setting. For two-level convergence, either the weak or strong approximation property provide sufficient conditions for convergence; however, the same is not true when considering the multilevel case. In the multilevel setting, if the smoothing and strong approximation properties hold on all levels of the multigrid hierarchy, with constants that are uniformly bounded (independently of the level in the hierarchy), then multilevel convergence of the multigrid V-cycle can be proven [27]. Even if the weak approximation property holds uniformly, though, multilevel convergence still cannot be guaranteed. Since $A(I-\pi_{A})=(I-\pi_{A}^{T})A$, we can easily derive the bound

[TABLE]

showing that if $\|I-\pi_{A}^{T}\|_{X}$ is not uniformly bounded across the levels in the hierarchy, then a uniform weak approximation property does not imply a uniform strong approximation property.

The standard example of this is the use of a piecewise constant interpolation operator, $P$, for any standard discretization of the Poisson problem. If $\mathbf{v}$ is a smooth vector, then $\|A\mathbf{v}\|_{X}$ will be small for many reasonable choices of $X$, such as $X=D^{-1}$, where $D$ is the diagonal of the system matrix, $A$. After coarse-grid correction, $(I-\pi_{A})\mathbf{v}$ will have jumps induced by the piecewise-constant interpolation, so $\|(I-\pi_{A})\mathbf{v}\|_{A}$ will be large, reflecting the high-frequency character of $(I-\pi_{A})\mathbf{v}$. With this, the strong approximation property can only be achieved with a large constant, $\beta$. In contrast, since $\|A(I-\pi_{A})\mathbf{v}\|_{X}$ will also be large, the weak approximation property can be fulfilled with a moderate constant, $\beta$. As is well-known, piecewise constant interpolation is sufficient for good two-level convergence, but not multilevel, consistent with the theoretical results.

Let $E_{MG}$ be the error-propagation matrix for a $V(1,1)$-cycle with a full multigrid hierarchy. The resulting convergence factor is bounded by the $A$-norm of $E_{MG}$, which can take the form

[TABLE]

for some $K_{MG}\geq 1$. Multigrid with an arbitrary number of levels can be thought of as a recursive use of two-grid methods with inexact coarse-grid solves, which is typically how convergence theory is formulated in the multilevel setting. The standard multilevel convergence result and some equivalent or sufficient conditions are stated below. In all cases, we assume that the multigrid hierarchy is specified by matrices $A^{(k)}$ and interpolation operators $P^{(k)}$, with the convention that $P^{(k)}$ is the interpolation operator from level $k+1$ to level $k$, and $A^{(k+1)}=\left(P^{(k)}\right)^{T}A^{(k)}P^{(k)}$. Furthermore, we assume that on each level, a relaxation scheme is specified that satisfies a consistent smoothing property,

[TABLE]

for all vectors, $\mathbf{v}^{(k)}$, with $\alpha$ independent of $k$.

Theorem 2.6** (Strong approximation property).**

If, for every $\mathbf{v}^{(k)}$, there exists a $\mathbf{v}^{(k+1)}$ such that

[TABLE]

for some $\beta$ independent of $k$. Then the multilevel $V(1,1)$-cycle converges uniformly, and $K_{MG}\leq\beta/\alpha$.

Note that the strong approximation property (SAP), is similar to the multilevel generalization of the WAP of equation (4) when $X=\frac{1}{\|A\|}I$,

[TABLE]

In this form of the WAP, interpolation of an eigenvector, $\mathbf{v}^{(k)}$, must be accurate in the $l^{2}$-norm to the order of its corresponding eigenvalue, with constant $\frac{\beta}{\|A^{(k)}\|}$. A stronger statement is required by the SAP, namely that interpolation of an eigenvector in the $A^{(k)}$-norm must be accurate to the order of its corresponding eigenvalue. A detailed look at the SAP can be found in Theorem 5.6.1 and Chapter 6 of[32].

Two sufficient conditions for the strong approximation property are stated below. Both rely on the smoothing property stated above holding uniformly across all levels. To simplify notation, we define $\pi_{k}=P^{(k)}\left(A^{(k+1)}\right)^{-1}\left(P^{(k)}\right)^{T}A^{(k)}$ as the $A^{(k)}$-orthogonal projection on level $k$ of the hierarchy, and $Q_{k}=P^{(k)}\left(\left(P^{(k)}\right)^{T}P^{(k)}\right)^{-1}\left(P^{(k)}\right)^{T}$ as the $\ell^{2}$-orthogonal projection on level $k$ of the hierarchy.

Corollary 2.7** ($l^{2}$-boundedness of $\pi_{k}$).**

If, for every $\mathbf{v}^{(k)}$,

[TABLE]

for some $\beta$ independent of $k$, then the multilevel $V(1,1)$-cycle converges uniformly, and $K_{MG}\leq\beta/\alpha$.

Corollary 2.8** (WAP$(A^{2})$).**

If, for every $\mathbf{v}^{(k)}$,

[TABLE]

for some $\beta$ independent of $k$, then the SAP holds with constant $\beta$ and $K_{MG}\leq\beta/\alpha$.

Proof 2.9**.**

For $\hat{\mathbf{v}}^{(k+1)}=\left(\left(P^{(k)}\right)^{T}P^{(k)}\right)^{-1}\left(P^{(k)}\right)^{T}\mathbf{v}^{(k)}$,

[TABLE]

2.3 “Optimal” and “ideal” interpolation

Returning to the two-level case, alongside bounds on two-grid convergence factors, specific interpolation operators have been derived as the best interpolation operator in certain contexts. Let $P=\begin{bmatrix}W\\ I\end{bmatrix}$ and $R=\begin{bmatrix}0&I\end{bmatrix}$, and consider relaxing the numerator of (5) from the $\widetilde{M}$-norm to the following problem:

[TABLE]

where $P_{ideal}$ is so-called “ideal interpolation” [15]. Due to the inverse of $A_{ff}$, $P_{ideal}$ is often a difficult operator to compute directly and may be dense, neither of which are compatible with the goals of AMG (see Section 3). However, denoting the graph distance between nodes $i$ and $j$ in $A_{ff}$ by $|i-j|_{G}$, then the following decay property is well-known:

[TABLE]

for some constant $C$ and $q<1$, where $q\approx\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$, for condition number, $\kappa$, of $A_{ff}$ [8, 13]. Thus, for sparse, well-conditioned $A_{ff}$ (as is expected with a proper choice of coarse grid), coefficients of $A_{ff}^{-1}$ decay exponentially fast away from the diagonal. Under this assumption, $A_{ff}^{-1}$ can be approximated well with a sparse matrix, at least in a Frobenius sense. Approaches to approximating $P_{ideal}$ can be found in [21, 8], each of which contain some variation of the following result:

Lemma 2.10** (Ideal interpolation).**

Let $P_{ideal}$ be as in (11) and let $P$ take the form $P=\begin{bmatrix}W\\ I\end{bmatrix}$, restricted to a fixed nonzero sparsity pattern. Then minimizing the difference between columns of $P$ and $P_{ideal}$ in the $A$-norm is equivalent to minimizing each column of $P$ in the $A$-norm. Furthermore, the solution of this minimization is unique.

However, considering that (11) does not provide a sharp bound on convergence, ideal interpolation typically does not provide optimal (two-grid) convergence factors over all $P$. The optimal $P$ with respect to two-grid convergence is given in the following lemma [6, Lemma 1].

Lemma 2.11** (Optimal interpolation).**

Let $0<\lambda_{1}\leq...\leq\lambda_{n}$ and $\mathbf{v}_{1},...,\mathbf{v}_{n}$ denote the eigenvalues and eigenvectors, respectively, of the generalized eigenvalue problem

[TABLE]

Then the minimal convergence rate of the two-grid method $\|E_{TG}(P)\|_{A}$ over all $P$ with dim$(P)=n_{c}$ is given by

[TABLE]

with corresponding optimal interpolation matrix given by

[TABLE]

Although results in [6] suggest that, at times, a sparse approximation to $P_{\textnormal{opt}}$ may be feasible, it is certainly more difficult to develop a cheap, sparse approximation to $P_{\textnormal{opt}}$ compared with $P_{\textnormal{ideal}}$. That being said, Lemma 2.11 does corroborate the general AMG approach of including eigenvectors of $A$ (actually of $\widetilde{M}^{-1}A$) associated with small eigenvalues in $\operatorname{Range}(P)$. In fact, it follows from Lemma 2.11 that if the first $n_{c}+1$ eigenvalues of $A\mathbf{v}=\lambda\widetilde{M}\mathbf{v}$ are all approximately zero, AMG cannot achieve strong convergence factors. This highlights the importance of the distribution of eigenvalues on the performance of AMG.

3 Interpolation in practice

AMG is a popular solver largely because of its linear complexity in the setup and solve phase. Let $A$ be SPD and consider forming a multigrid hierarchy to solve $A\mathbf{x}=\mathbf{b}$. Interpolation operators in AMG methods are often (implicitly) constructed with the goal of controlling or minimizing some functional with a theoretical relation to convergence, such as (5), (8), (9), or (11). However, there are two important factors that must be considered in practice and are generally absent from theory – (i) the process used to form interpolation operators must remain linear in complexity in keeping with the desired $O(n)$ total cost for AMG methods, and (ii) interpolation operators must remain sparse in order to construct a sparse coarse-grid matrix that can be used in a recursive process.

These constraints can prove difficult to achieve when designing AMG methods, and make approximating some of the bounds in Section 2 more tractable than others. In particular, the operators used in the convergence theory, such as $\pi_{A}$, $P_{ideal}$, and $P_{opt}$, generally cannot be easily computed in any practical setting. Furthermore, these convergence results are typically required to hold for all $\mathbf{v}$ or, equivalently, for some basis for the space such as the eigenvectors of $A$. Constructing interpolation or coarsening based on a full basis of vectors is generally not tractable in linear complexity and, thus, two forms of approximation are often used, (i) work with a candidate set of $k$ vectors, where $k\ll n$, or (ii) work in an operator norm, which is a supremum over all vectors.

The first approach is to directly satisfy conditions of a theorem but only for a set of candidate vectors of dimension $k\ll n$. This is a standard approach for satisfying the WAP or SAP, and classical AMG [27] can be viewed as doing this for only the constant vector, while smoothed aggregation (SA) [30] may use a larger basis, such as the rigid-body modes for elasticity problems. Most adaptive multigrid methods [10, 9, 5] can also be viewed in this way, where the candidate vectors arise from the adaptive process. As in (5), the (two-grid) convergence rate is bounded by the maximum of $K_{TG}$ over all $\mathbf{v}\neq\mathbf{0}$. Noting the denominator of $\|\mathbf{v}\|_{A}^{2}$, the maximum will generally occur for $\mathbf{v}$ associated with small eigenvalues of $A$, where $\|\mathbf{v}\|_{A}^{2}$ is very small compared to $\|\mathbf{v}\|_{\widetilde{M}}^{2}$. For differential operators, it is common to have a zero or near-zero row sum, making the constant a good representation of low-energy modes. Developing and using additional candidate vectors is the basis for adaptive approaches, which are designed for difficult linear systems beyond the scope of classical SA or AMG [5, 9, 10]. In such solvers, an adaptive process is used to develop a set of target vectors representative of low-energy modes of $A$. Interpolation is then constrained to interpolate these modes either exactly or nearly so, and the process is repeated on coarse grids.

An alternative approach is to formulate the minimization over all vectors $\mathbf{v}\neq\mathbf{0}$. Consider the energy-stability constraint in Lemma 2.5, based on controlling $\|PR\|_{A}^{2}$. As an induced $A$-norm, $\|PR\|_{A}$ is defined via a supremum over all $\mathbf{v}\neq\mathbf{0}$. However, $\|PR\|_{A}$ can be bounded using the Frobenius norm, which gives an indirect approach to bounding $\|PR\|_{A}$. Let $R=\begin{bmatrix}0&I\end{bmatrix}$, then

[TABLE]

where $S_{A}^{-1}:=RA^{-1}R^{T}=(A_{cc}-A_{cf}A_{ff}^{-1}A_{fc})^{-1}$, is the inverse of the Schur complement of $A$ in $A_{cc}$. Although an interesting equivalence, the Schur complement is difficult to form in practice. A more tractable approach is obtained by using an intermediate bound,

[TABLE]

Given the form of $R=(\mathbf{0},I)$, this simplifies to

[TABLE]

Minimizing $\operatorname{tr}(P^{T}AP)$ was proposed in this form in [8], and is equivalent to minimizing columns of $P$ in the $A$-norm. This approach has been used in smoothed aggregation (SA) [30], root-node AMG [21, 28], and the general energy-minimization framework proposed in [26]. Recall from Lemma 2.10 that minimizing energy in columns of $P$ is also equivalent to minimizing the difference between columns of $P$ and $P_{ideal}$ in the $A$-norm.

It is worth considering the leading constant in (12), $\|A^{-1}\|=\frac{1}{\lambda_{\text{min}}(A)}$, as this is likely large and could lead to a poor bound on $\|PR\|_{A}^{2}$. Note that the energy constraint in Lemma 2.5 can also be formulated as

[TABLE]

for all vectors $\mathbf{v}$ with $\mathbf{v}_{c}=R\mathbf{v}$ [16, Theorem 5.2]. The factor of $\|A^{-1}\|$ in (12) accounts for the possibility that $P$ is chosen so that $\mathbf{v}_{c}^{T}A_{c}\mathbf{v}_{c}$ is an $O(1)$ quantity when $\mathbf{v}$ corresponds to the smallest eigenvalue of $A$ (corresponding to a bad choice of $P$). In practice, we make the heuristic assumption that the choice of $P$ will not be so bad and, thus, minimizing $\operatorname{tr}(P^{T}AP)=\sum_{i}\lambda_{i}(A_{c})$ is an effective way to control $\|PR\|_{A}$.

4 Trace-minimization

As discussed in Section 3 and can be seen in other AMG methods, AMG interpolation operators are often constructed based on some combination of ensuring that a given set of candidate vectors is interpolated exactly, while ensuring energy stability of the coarse-grid operator. In this direction, we now propose to form $P$ through minimizing a general weighted functional combining these two approaches,

[TABLE]

for $\tau\in[0,1)$ and candidate vector $\mathbf{v}$. If multiple candidate vectors, $\{\mathbf{v}_{i}\}$, are available a priori, for example, the rigid body modes in elasticity, then we minimize over the maximum $\mathbf{w}\in$ Span$\{\mathbf{v}_{i}\}$:

[TABLE]

This is a complementary approach, focusing on achieving accurate interpolation of the low-energy modes in the candidate set as well as energy stability on the coarse grid. It is also complementary in the sense that the first term is defined over a candidate set of vectors, $\{\mathbf{v}_{i}\}$, while the second term is defined over $P$, and should improve interpolation regardless of the provided candidate vectors.

Let $P$ take the form $P=\begin{bmatrix}W\\ I\end{bmatrix}$, and consider minimizing (13). Define a set of $n_{B}$ candidate vectors as $A$-orthonormalized columns of a matrix $B=\begin{bmatrix}B_{f}\\ B_{c}\end{bmatrix}$, and let $X\simeq\widetilde{M}$ as in (1). Then, consider minimizing $K_{TG}$ from (5), restricted to unit linear combinations of $\mathbf{v}\in\mathcal{B}=\{\mathbf{w}\in\operatorname{Range}(B)\mid\|\mathbf{w}\|_{A}=1\}$:

[TABLE]

This approximates the WAP in the $X$-norm using an $l^{2}$-projection onto $\operatorname{Range}(P)$ (as opposed to the optimal $\pi_{\widetilde{M}}$-orthogonal projection). Recall the second term in (13) corresponds to minimizing the columns of $P$ in the $A$-norm. Expanding $\operatorname{tr}(P^{T}AP)$ gives

[TABLE]

Substituting equations (14) and (15) into (13) gives a functional of $W$ to minimize in forming $P$. Dropping terms independent of $W$ and pulling out a factor of two for a more familiar form, define

[TABLE]

Observe that (16) is a quadratic functional in $W$. Define a bounded linear operator, $\mathcal{L}$, and right-hand-side, $\mathcal{B}$, as

[TABLE]

in which case $\mathcal{F}(W)=\tfrac{1}{2}\langle\mathcal{L}W,W\rangle_{F}-\langle W,\mathcal{B}\rangle_{F}$. Note that if $A_{ff}$ and $X_{ff}$ are symmetric and positive definite, then $\mathcal{L}$ is self-adjoint and positive definite in the Frobenius norm:

[TABLE]

Using the symmetry of $\mathcal{L}$, the first and second Frechét derivative of $\mathcal{F}$ are given by:

[TABLE]

Since $\mathcal{L}$ is self-adjoint and positive definite, $\mathcal{F}^{\prime\prime}(W)[V]\geq 0$ $\forall$ $V$. Thus, the minimum of $\mathcal{F}$ in $W$ is achieved at $W$ such that $\mathcal{F}^{\prime}(W)[V]=0$, and $\mathcal{F}^{\prime}(W)=0$ $\forall$ $V$ if and only if $\mathcal{L}W=\mathcal{B}$. This has a unique solution, $W=\mathcal{L}^{-1}\mathcal{B}$. However, it is likely that $\mathcal{L}^{-1}\mathcal{B}$ is dense and not practical, motivating a constrained sparsity pattern for $W$.

4.1 Constrained sparsity pattern

In practice, the sparsity pattern of $W$ must be fixed a priori in order to control the operator complexity of $W$ and $A_{c}$. Define a vector space

[TABLE]

for a set of indices $\mathcal{N}$ denoting a fixed sparsity pattern for $W$. A Hilbert space $\mathcal{H}$ can be defined over $\mathcal{X}$ with the Frobenius inner product, $\langle A,B\rangle_{F}=\sum_{ij}A_{ij}B_{ij}$. It is easily verified that $\mathcal{X}$ is complete over the norm induced by $\langle\cdot,\cdot,\rangle_{F}$ due to the completeness of $\mathbb{R}$. Now define the bounded linear functional $\hat{\mathcal{L}}:\mathcal{H}\to\mathcal{H}$ as

[TABLE]

and a corresponding bilinear form

[TABLE]

A quadratic form as in (16) restricted over $\mathcal{N}$ can then be defined as

[TABLE]

where $\hat{\mathcal{B}}\in\mathcal{H}$ is $\mathcal{B}$ restricted to $\mathcal{N}$. Note that in $\mathcal{H}$, $\langle W,\mathcal{B}\rangle_{F}=\langle W,\hat{\mathcal{B}}\rangle_{F}$. A similar derivation as shown for $\mathcal{L}$ confirms that $\hat{\mathcal{L}}$ is self-adjoint and $a(W,V)$ symmetric. Then, observe that for $W\in\mathcal{H},W\neq 0$, $\hat{\mathcal{L}}$ and $a(W,V)$ are positive:

[TABLE]

The following standard lemma of functional analysis can then be invoked to find a minimizer of $\hat{\mathcal{F}}(W)$ (18).

Lemma 4.1**.**

Let $a(x,y)$ be a bounded, symmetric, positive-definite bilinear form on a Hilbert space $\mathcal{H}$, and $\mathcal{G}(x)$ be a bounded linear functional on $\mathcal{H}$. Then the following are equivalent

[TABLE]

Furthermore, there exists a unique solution $x\in\mathcal{H}$ satisfying (19), (20).

Based on Lemma 4.1, we seek the unique solution to

[TABLE]

which can be iterated towards using the preconditioned conjugate gradient method introduced in Section 5.

Remark 4.2**.**

A conceptual limiting case of the proposed weighted energy minimization is to interpolate candidate vectors exactly and minimize energy based on that constraint. However, this does not directly fit into the framework of (13). Constrained energy-minimization has been proposed in various forms [8, 20, 26, 33], and was used as a basis for root-node AMG in [21]. Defining the affine space $\mathcal{A}=\{W\text{ : }W\in\mathcal{H}\text{ and }WB_{c}=B_{f}\}$, the constrained minimization problem is given by

[TABLE]

Since the linear operator now consists of normal matrix multliplication, $A_{ff}W$ as opposed to the left and right multiplication in (17), the existence and uniqueness of a solution to (22) can be shown in a linear algebra setting (see [26]) along with a CG implementation based on projecting into $\mathcal{A}$.

5 Preconditioned conjugate gradient

Because $\hat{\mathcal{L}}$ is self-adjoint and positive in $\mathcal{H}$, conjugate gradient (CG) in the Hilbert space setting is a competitive approach to solving (21) in an iterative fashion. It is generally advisable to precondition CG iterations for optimal convergence. Here, we construct a diagonal preconditioner for (21) to make iterations more robust when $\hat{\mathcal{L}}$ is poorly conditioned at a marginal increase of computational cost.

Unlike with matrices, however, it is not clear what the “diagonal” of $\hat{\mathcal{L}}$ is. Let $W\in\mathbb{R}^{N_{f}\times N_{c}}$, and define the operator $\overline{(W)}$ as the columns of $W$ stacked in a column-vector. Note that $\overline{(W^{T})}$ then gives the rows of $W$ stacked as a column-vector. Let $Y$ be the permutation matrix such that $\overline{(W)}=Y\overline{(W^{T})}$ and $YY^{T}=Y^{T}Y=I$, which can be thought of as a mapping of $W$ from row-major format to column-major format. First note the following lemma with regards to Kronecker products and the action of $Y$.

Lemma 5.1**.**

Let $Y$ be a permutation matrix mapping $W\in\mathbb{R}^{N_{f}\times N_{c}}$ from row-major format to column-major format, that is, $\overline{(W)}=Y\overline{(W^{T})}$. Then, for any $P\in\mathbb{R}^{N_{f}\times N_{f}}$ and $Q\in\mathbb{R}^{N_{c}\times N_{c}}$,

[TABLE]

Proof 5.2**.**

First consider the structure of $Y$. Note the following relations between $W,\overline{(W)}$, and $\overline{(W^{T})}$, i.e. $W$ stored as a standard dense matrix, a column-major matrix, and a row-major matrix, respectively,

[TABLE]

Defining $Y$ such that $Y\overline{(W^{T})}=\overline{(W)}$, it follows that

[TABLE]

and the action of $YAY^{T}$ is then given as

[TABLE]

Now consider the element-wise Kronecker products of $P$ and $Q$:

[TABLE]

for $i,k\in[0,N_{f}],j,l\in[0,N_{c}]$. Combining (23), (24), and (25) gives

[TABLE]

It follows that $Y(P\otimes Q)Y^{T}=Q\otimes P$.

Remark 5.3**.**

Lemma 5.1 is a known result that we arrived at inadvertently, where $Y$ is known as the “Perfect Shuffle” matrix [11]. Its relation to row-major and column-major storage of matrices is, to our knowledge, a new contribution to the literature.

Now consider finding the diagonal of $\mathcal{L}$ by looking at $\overline{\mathcal{L}}$ as an operator on $\overline{(W)}$. To do so, represent the action of $A_{ff}W$ through $(I_{N_{c}}\otimes A_{ff})\overline{(W)}$, where $(I_{N_{c}}\otimes A_{ff})$ gives a block diagonal matrix of $N_{c}$ $A_{ff}$’s, each to be multiplied by one column of $W$. Recalling the identity $(A\otimes B)(C\otimes D)=(AB\otimes CD)$ and Lemma 5.1,

[TABLE]

This derivation can be naturally extended to $\mathcal{H}$, where $W\in\mathcal{H}$ has a specified sparsity pattern, $\mathcal{N}$, by setting the $k$th row and column of $\overline{\mathcal{L}}$ equal to zero for all $k$ such that $\overline{(W)}_{k}:=W_{ij}$, and $(i,j)\not\in\mathcal{N}$. Because $\overline{\mathcal{L}}$ is a block operator with block size $N_{f}\times N_{f}$, it follows that there is a distinct “diagonal” in $\mathcal{L}$ corresponding to each $j$th column of $W$,

[TABLE]

A diagonal preconditioning for $\hat{\mathcal{L}}$ is then given by taking the Hadamard product with $\mathcal{D}\in\mathcal{H}$, where the $j$th column of $\mathcal{D}$ is given by the element-wise inverse of (26):

[TABLE]

In the case of $A_{ff}$ having a constant or near-constant diagonal, and letting $X_{ff}$ be the diagonal of $A_{ff}$ (a common practical choice), $\mathcal{D}$ is constant or near-constant. In practice, preconditioning with $\mathcal{D}$ is important for problems in which diagonal elements of $A$ or target vectors $B$ consist of a wide range of values.

5.1 Sylvester and Lyapunov equations

In fact, (27) can be used to define a preconditioner for general systems of the Sylvestor- or Lyapunov-type:

[TABLE]

for solution matrix $W$, where $A,B,C$ and $D$ need not be symmetric (of course an appropriate Krylov solver must be chosen based on properties of the functional). A diagonal preconditioner for (28) is given by taking the Hadamard product with

[TABLE]

Systems of the form in (28) arise often in the context of optimal control theory. Letting $B=C=I$, (28) is a Sylvester equation; letting $B=A^{T}$, $C=-I$, and $D=I$, (28) is a discrete Lyapunov equation; and letting $B=C=I$ and $D=A^{T}$, (28) is a continuous Lyapunov equation. There have been many efforts at developing Krylov methods and preconditioners for such systems; for example, see [12, 14, 17, 29, 34]. Here we develop a simple preconditioner for problems of the form (29), that is easy to construct and apply.

6 Numerical results

In this section, we present numerical results for a variety of problems, comparing a weighted energy minimization and constrained energy minimization, and analyzing the choice of constraint vector. The method proposed here is implemented in the PyAMG library [3]; AMG methods such as strength-of-connection, coarsening, etc., follow that of [21], and the reader is referred there for details. In figures, RN refers to a constrained energy minimization using root-node AMG [21] and TM${}_{10^{k}}$ refers to weighted energy minimization proposed here with weight $\tau=10^{k}$. The test problems considered are:

Anisotropic diffusion: 2-dimensional rotated anisotropic diffusion, discretized with linear finite elements, on an unstructured triangular mesh:

[TABLE]

where

[TABLE]

Due to the unstructured mesh, all angles $\theta\in(0,\nicefrac{{\pi}}{{2}})$ are effectively equivalent from a solver perspective; thus, moving forward we (arbitrarily) let $\theta=\nicefrac{{3\pi}}{{16}}$. Mesh spacing is taken to be $h\approx 1/1000$, resulting in approximately 1.25M DOFs. 2. 2.

Diffusion with an oscillatory coefficient: 2-dimensional diffusion problem (as in equations (30-31)), discretized with linear finite elements on a structured, regular triangular mesh with $2N^{2}$ elements, with a piecewise linear coefficient that oscillates at every other grid point, regardless of mesh size:

[TABLE]

where

[TABLE]

Taking $h=1/N$, we recognize $Nx=x/h$ as the (integer) index of a mesh point in the $x$-direction, with a similar interpretation of $Ny$. For points $(x,y)$ not on the mesh, $f(x,y)$ is interpolated linearly (as a function in the finite-element space). This results in a number of coefficient oscillations that grows proportionally with the number of mesh points, resulting in a checkerboard-like pattern with alternating large and small coefficients. Diffusion problems with large and frequent coefficient changes, such as this one, traditionally make difficult test problems for multigrid methods.

6.1 Determining $P$

Here we look at AMG convergence as a function of the number of iterations of CG used to determine $P$ and how constraint vectors are enforced, either exactly, in a constrained energy minimization, or weighted by $\tau\in(0,1)$, in a weighed energy minimization. For all results, a V-cycle is applied with two iterations of Jacobi pre- and post-relaxation as a preconditioner for CG. Unless otherwise specified, the constraint vector is chosen as the constant vector with several Jacobi smoothing iterations applied; weighted energy-minimization uses the diagonal preconditioning of Section 5; and constrained energy-minimization uses the diagonal preconditioning of [26].

6.1.1 Anisotropic Diffusion

Figure 1 shows the work-per-digit-of-accuracy (WPD) as a function of the number of iterations of CG used to determine $P$, for variations in energy minimization applied to anisotropic Poisson (Problem # 1), with anisotropy $\epsilon\in\{1,0.001,0\}$. WPD is defined as

[TABLE]

where $C$ is the cycle-complexity of the multigrid solver and $\rho$ is the average convergence rate of the solver over all iterations. This metric measures how many work-units, defined as the floating point operations to perform a single matrix-vector multiply, are required to reduce the residual by one order of magnitude. This metric is particularly useful for cross-comparisons of solvers with differing sparsity structures. For more detail, see, for instance, [21].

Interpolation is fixed to use a degree-four sparsity pattern, that is, the sparsity pattern for each column of $P$ reaches out to neighbors within graph distance four from the corresponding C-point (see [21, 28]). This wider sparsity pattern often leads to better convergence rates for difficult problems [28], but also requires more iterations of energy-minimization. Essentially, wider sparsity patterns create more interpolation coefficients in $P$, which are then determined through energy-minimization.

Several immediate results follow from Figure 1. First, there is a limit at which additional iterations to determine $P$ no longer improve convergence. For the isotropic case $(\epsilon=1)$, the best convergence rates are obtained by simply enforcing the constraint with a single constrained smoothing pass; additional energy-minimization steps do not improve convergence. As the level of anisotropy increases ($\epsilon\to 0$), the number of iterations of CG required to achieve the best performance increases. However, convergence of the AMG solver based on a given constraint vector and coarsening scheme remains bounded below, regardless of further energy minimization of $P$. Second, it is clear that enforcing the constraint exactly or near-exactly is fundamental to good convergence, even for the simplest isotropic problem. Although theory tells us that interpolating low-energy modes is necessary for good convergence, the fact that this cannot be achieved through weighted energy minimization is slightly non-intuitive. Energy-minimization reduces the columns of $P$ in the $A$-norm, which should thus build $P$ to include low-energy modes in its range. Heuristically, it seems that after a handful of CG iterations, the range of $P$ would contain sufficient low-energy modes for good convergence. However, it is clear in Figure 1 that even in the isotropic case, using a large $\tau=0.1$ to focus on energy minimization over constraints leads to very poor performance.

Together, these points underline the role of energy minimization in AMG convergence as an acceleration technique. For some difficult problems, energy minimization is critical to achieving scalable convergence. Strongly anisotropic diffusion is one such example that typically proves difficult for standard AMG methods, but can be solved effectively with constrained energy minimization [21]. Nevertheless, regardless of energy minimization, strong convergence cannot be obtained without enforcing or nearly-enforcing an appropriate constraint vector (Figure 1).

6.1.2 Diffusion with an oscillatory coefficient

Figure 2 shows the WPD as a function of the number of iterations of CG used to determine $P$ for variations in energy minimization applied to the oscillating coefficient problem (Problem # 2), with coefficient oscillations of $K=10^{6}$ and $K=10^{3}$. For $K=10^{6}$, Figure 2(a) shows results for diagonal preconditioning of energy-minimization and Figure 2(b) shows the case of no preconditioning. Figure 2(c) shows the case of diagonal preconditioning with $K=10^{3}$. Comparing Figures 2(a) and 2(b), we see that using preconditioning in weighted energy minimization reduces the number of iterations necessary to achieve good convergence. Moreover, preconditioning appears to actually improve the best achievable AMG convergence factor in practice. For constrained energy-minimization, energy minimization iterations without preconditioning increases the WPD by $3-5\times$ within a reasonable number of iterations on $P$ (of course, asymptotically the preconditioned and non-preconditioned results are equivalent but, in practice, only $O(1)$ iterations are done.) This raises an interesting question as to if better preconditioners for energy minimization can actually improve the AMG solver’s performance in a way that additional iterations with a diagonal preconditioner cannot in practical time; however, this is a topic for future study.

Focusing on the more practical solvers in Figure 2(a), we also see that the results mirror those in Figure 1. Overall, the constrained energy-minimization case performs best, with weighted energy-minimization able to approach the constrained case only for the right $\tau$ values and enough energy-minimization iterations on $P$. Again, there is a limit beyond which additional energy-minimization iterations no longer improve AMG convergence. For constrained energy-minimization, relatively few iterations are needed. Lastly, enforcing the constraint exactly or near-exactly is fundamental to good convergence. Using energy-minimization with larger $\tau$ values leads to poor performance.

The effects of the oscillating coefficient $K$ can be seen by comparing Figures 2(a) and 2(c). Interestingly, the larger $K$ value leads to a need for smaller $\tau$ values for the weighted case, (compare the curves for $\tau=10^{-7}$). Overall, apart from changing the size of beneficial $\tau$ values, the size of the coefficient oscillation does not noticeably affect either the weighted or constrained energy-minimization.

A final note of interest is that larger interpolation sparsity patterns do not help here. Thus, a moderate sparsity pattern of degree three is chosen for these results.

Remark 6.1**.**

We did not find tracking the CG residual norm during energy-minimization to be useful and, hence, omit plots of this information. The key difficulty is that it is not clear how to connect the residual norm to the eventual multigrid convergence rate. In other words, it is not clear how to use the residual norm to halt the energy-minimization process. For instance, taking the cases of constrained energy-minimization from Figures 1 and 2, it is clear that at most five iterations of energy-minimization are needed. However, the residual norm continues to decrease monotonically by multiple orders of magnitude from iteration five to iteration 19. Yet, this extra residual reduction does not speed up convergence of the resulting multigrid solver. In practice, the number of iterations needed typically equals the degree of the sparsity pattern of $P$ plus some small number, usually two or three. This number of iterations is required to first fill the allowed sparsity pattern, and then to provide two or three iterations of additional improvement.

6.2 Constraint vectors and adaptivity

In Section 6.1, we learned two things: (i) for good convergence, it is important that $P$ exactly or almost exactly interpolates an appropriate constraint vector, and (ii) coupled with a good constraint, energy minimization can improve convergence, but only by a fixed amount. This leads to the natural idea of adding an additional constraint vector when further energy minimization of $P$ no longer improves convergence. Such an approach is the basis of adaptive multigrid methods, where a set of constraint vectors are developed that are then included or approximately included in the range of $P$ [5, 9, 10]. There are multiple ways to generate constraint vectors; here we take the simple approach of generating a random vector $\mathbf{x}_{0}$ and applying some form of improvement iterations (either relaxation or V-cycles) to reduce $\|\mathbf{x}_{0}\|_{A}$. Table 3 shows results for constrained energy minimization AMG applied to the anisotropic Poisson problem, with varying numbers of improvement iterations and varying numbers of constraint vectors.

Several interesting things follow from the results in Table 3. First, the difference in convergence factor between two-grid and multigrid is very small. This indicates that we are solving our coarse-grid problem well using V-cycles, and that convergence is limited by how “good” the coarse-grid problem is, and not how accurately we are solving it. Moreover, naively adding constraint vectors that were not accounted for in the range of $P$ does not improve convergence and, in fact, degrades convergence in all cases, while increasing the setup complexity. Although more involved processes have been developed for adaptive multigrid methods, these simple tests give insight that improving convergence is not as simple adding new constraint vectors.

7 Conclusions

This paper explores the role of energy minimization in AMG interpolation from a theoretical and practical perspective. The eventual goal is to develop improved interpolation techniques that are more robust than current state-of-the-art, without the significant overhead setup cost of fully adaptive methods. A minimization framework is developed based on a weighted combination of interpolating known low-energy modes with a global energy minimization over $P$. On one hand, accurately interpolating the constraint vectors proves to be of fundamental importance to good convergence, as observed where constrained energy minimization consistently performs best, and weighted energy minimization performs best with the relative weight of interpolating constraints $\gg 0.99$. However, convergence generally does not improve when additional constraint vectors are added beyond the first. This either means, for these test problems, (i) accurately interpolating one constraint vector leads to convergence factors close to the optimal rate for the given coarse grid [6], or (ii) there are other factors fundamental to convergence of AMG that are not being addressed in this framework. Results here do not suggest the newly proposed algorithm is superior to existing methods, but do provide insight on the convergence of AMG in the practical setting, as well as the relation to AMG convergence theory.

Acknowledgments

The authors gratefully acknowledge the contributions of Ludmil Zikatanov to the work presented here.