Convergence analysis of a two-grid method for nonsymmetric positive definite problems

Multigrid is one of the most efficient numerical methods for solving large-scale systems of linear equations that arise from discretized partial differential equations. It has been shown to be a powerful solver, with linear or near-linear computational complexity, for a large class of linear systems; see, e.g., [9, 3, 21, 22]. The fundamental module of multigrid methods is a two-grid scheme, which consists of two complementary error-reduction processes: smoothing and coarse-grid correction. The smoothing process is typically a simple iterative method, such as the (weighted) Jacobi and Gauss–Seidel methods. One limitation of these classical methods is that they are only efficient at eliminating high-frequency (i.e., oscillatory) error in general. The remaining low-frequency (i.e., smooth) error will be further reduced by the coarse-grid correction, which is essentially a projection-type method. These two processes will be applied iteratively in multigrid algorithms until a desired residual tolerance is achieved.

Much of multigrid analysis within the literature relies on the problem matrix, $A$, being symmetric positive definite (SPD). For such problems, the convergence theory of multigrid methods has been well developed; see, e.g., [24, 4, 5, 22, 25, 27, 28]. When $A$ is SPD, the convergence analysis of a symmetric two-grid method (i.e., pre- and postsmoothing steps are performed in a symmetric way) amounts to analyzing the eigenvalues of its iteration matrix (also called the error propagation matrix), because the $A$-norm of the iteration matrix, called the convergence factor, coincides with its spectral radius. In the case of Galerkin coarse solver, an identity has been established to characterize the convergence factor of two-grid methods [24, 5, 29], which has been widely used to analyze two-grid methods; see, e.g., [5, 25, 1, 26]. Nevertheless, it is often too costly to solve the Galerkin coarse-grid system exactly, especially when its size is relatively large. Recently, Xu and Zhang [27, 28] extended the two-grid convergence theory to the case of non-Galerkin coarse solvers. They proposed two theoretical frameworks for analyzing the convergence of inexact two-grid methods, which have been successfully applied to the convergence analysis of multigrid methods.

Unlike the SPD case, the convergence theory of multigrid methods for nonsymmetric problems is still not mature. Most existing algorithms are based largely on heuristic or incomplete theory. In general, a nonsymmetric matrix fails to induce a norm and the error propagation matrix of the coarse-grid correction is an oblique projection. Thus, two fundamental questions need to be answered when studying two-grid methods for nonsymmetric problems. The first one is how to choose a suitable convergence measure. Ideally, the selected measure can facilitate us to derive an identity for characterizing the convergence of two-grid methods, as in the SPD case. The simplest choice is the spectral radius of iteration matrix [17, 23, 15, 6]. However, it describes only the asymptotic convergence, which may be a completely misleading indicator [8, 14, 18], especially when the iterations are not carried out enough times. Some other measures (e.g., the $\sqrt{A^{T}A}$- and $\sqrt{AA^{T}}$-norms) and the associated convergence estimates can be found in [2, 12, 14, 13, 16, 18] and the references therein.

The second question is how to build the intergrid transfer operators $R\in\mathbb{R}^{n_{\rm c}\times n}$ (restriction) and $P\in\mathbb{R}^{n\times n_{\rm c}}$ (prolongation) such that the coarse-grid correction will not increase error [16], where $n_{\rm c}$ is the number of coarse variables and $n$ is the number of fine ones ($A\in\mathbb{R}^{n\times n}$). The error propagation matrix of the correction process takes the form $I-P(RAP)^{-1}RA$, provided that $RAP\in\mathbb{R}^{n_{\rm c}\times n_{\rm c}}$ is nonsingular. Note that $I-P(RAP)^{-1}RA$ is a projection (i.e., an idempotent operator). For any operator norm $\|\cdot\|$, it holds that

If $\|I-P(RAP)^{-1}RA\|>1$, then the correction process may increase error, which is contrary to the core idea of multigrid methods. Indeed, the coarse-grid correction is key to the fast convergence of multigrid algorithms. Therefore, for a fixed $R$, one needs to design a prolongation $P$ such that $\|I-P(RAP)^{-1}RA\|=1$ in order to ensure that the correction process will not increase error.

In this paper, we are concerned with the convergence theory of a two-grid method for nonsymmetric positive definite problems (e.g., linear systems arising from discretized convection-diffusion equations). In the two-grid method, an SPD matrix, $M\in\mathbb{R}^{n\times n}$, is used as a smoother and the induced $M$-norm is used for measuring two-grid convergence. In addition, $M^{-1}A^{T}R^{T}$ is used as a prolongation operator in order to ensure that the coarse-grid correction will not increase error. Such a prolongation can be viewed as a multiplicative perturbation of the classical choice $R^{T}$ (since $M^{-1}A^{T}\approx I$). Unlike the existing estimates, we here establish a succinct identity for the convergence factor of the two-grid method. It can be conveniently used to derive a class of optimal restriction operators and analyze how the convergence factor is influenced by the row space of $R$. Compared with the exact case, the convergence theory of inexact two-grid methods is of more practical significance, while, for nonsymmetric problems, such theory is scarce in the literature. Motivated by this observation, we present some convergence estimates for an inexact variant of the two-grid method (both linear and nonlinear coarse solvers are considered), which lay a theoretical foundation for developing new algorithms.

The rest of this paper is organized as follows. In Section 2, we present a two-grid method for solving nonsymmetric positive definite linear systems. In Section 3, we establish an identity for the convergence factor of the two-grid method, followed by two applications of the identity. In Section 4, we propose and analyze an inexact variant of the two-grid method. In Section 5, we give some concluding remarks.

For convenience, we first list some notation used in the subsequent discussions.

• –

  $I_{n}$ denotes the $n\times n$ identity matrix (or $I$ when its size is clear from context).

• –

  $\lambda_{\min}(\cdot)$, $\lambda_{\min}^{+}(\cdot)$, and $\lambda_{\max}(\cdot)$ stand for the smallest eigenvalue, the smallest positive eigenvalue, and the largest eigenvalue of a matrix, respectively.

• –

  $\lambda_{i}(\cdot)$ denotes the $i$th smallest eigenvalue of a matrix.

• –

  $\lambda(\cdot)$ denotes the spectrum of a matrix.

• –

  $\|\cdot\|_{2}$ denotes the spectral norm of a matrix.

• –

  $\|\cdot\|_{M}$ denotes the norm induced by an SPD matrix $M\in\mathbb{R}^{n\times n}$: if $\mathbf{v}\in\mathbb{R}^{n}$, then $\|\mathbf{v}\|_{M}=\sqrt{\mathbf{v}^{T}M\mathbf{v}}$; if $X\in\mathbb{R}^{n\times n}$, then $\|X\|_{M}=\max_{\mathbf{v}\in\mathbb{R}^{n}\backslash\{0\}}\frac{\|X\mathbf{v}\|_{M}}{\|\mathbf{v}\|_{M}}$.

• –

  $\mathbb{E}[\cdot]$ denotes the expectation of a random variable.

Consider solving the linear system

where $A\in\mathbb{R}^{n\times n}$ is nonsymmetric but positive definite (namely, $\mathbf{v}^{T}A\mathbf{v}>0$ for all $\mathbf{v}\in\mathbb{R}^{n}\backslash\{0\}$, which implies that $A$ is nonsingular), $\mathbf{u}\in\mathbb{R}^{n}$, and $\mathbf{f}\in\mathbb{R}^{n}$. A natural splitting of $A$ is given by

where $A_{\rm sym}=\frac{1}{2}(A+A^{T})$ and $A_{\rm skew}=\frac{1}{2}(A-A^{T})$ are the symmetric and skew-symmetric parts of $A$, respectively. It is easy to check that $A$ is positive definite if and only if the symmetric matrix $A_{\rm sym}$ is positive definite.

Several basic assumptions involved in the analysis of two-grid methods are summarized as follows.

• •

  Let $M\in\mathbb{R}^{n\times n}$ be an SPD smoother such that $\|I-M^{-1}A\|_{M}\leq 1$.

• •

  Let $R\in\mathbb{R}^{n_{\rm c}\times n}$ be a restriction matrix of rank $n_{\rm c}$, where $n_{\rm c}\,(<n)$ is the number of coarse variables.

• •

  Let $P\in\mathbb{R}^{n\times n_{\rm c}}$ be a prolongation (or interpolation) matrix of rank $n_{\rm c}$.

• •

  Assume that the coarse-grid matrix $A_{\rm c}:=RAP\in\mathbb{R}^{n_{\rm c}\times n_{\rm c}}$ is nonsingular.

With the above assumptions, an exact two-grid method for solving (2.1) can be described by Algorithm 1, in which $A_{\rm c}^{-1}$ is used as a coarse solver.

In fact, Algorithm 1 involves the following two error-reduction processes:

where

Then

Recall that $\|I-M^{-1}A\|_{M}\leq 1$. The initial error $\big{\|}\mathbf{u}-\mathbf{u}^{(0)}\big{\|}_{M}$ will not be amplified by (2.3a).

It can be seen from (2.4) that $\varPi_{A}$ is a projection along (or parallel to) $\operatorname*{null}(RA)$ onto $\operatorname*{range}(P)$, which leads to

with equality if and only if $\varPi_{A}=M^{-1}\varPi_{A}^{T}M$. To ensure that the error $\big{\|}\mathbf{u}-\mathbf{u}^{(1)}\big{\|}_{M}$ will not be amplified by (2.3b), it suffices to choose a prolongation $P$ such that

Observe that $\varPi_{A}$ is a projection along $\operatorname*{null}(RA)$ onto $\operatorname*{range}(P)$, while the right-hand side of (2.5) is a projection along $\operatorname*{null}(P^{T}M)$ onto $\operatorname*{range}(M^{-1}A^{T}R^{T})$. Hence, (2.5) entails that

which can be satisfied by

for any $n_{\rm c}\times n_{\rm c}$ nonsingular matrix $W$. Taking $W=I_{n_{\rm c}}$ yields the prolongation

For any coarse vector $\mathbf{v}_{\rm c}\in\mathbb{R}^{n_{\rm c}}$, $P_{\star}\mathbf{v}_{\rm c}$ can be easily computed, because the action of $M^{-1}$ is always available. Furthermore, it is possible that $P_{\star}$ has a sparse structure. For instance, if $A$ is sparse, $M=\omega\operatorname*{diag}(A)$, and $R=\big{(}0\,\ I_{n_{\rm c}}\big{)}$, then $P_{\star}$ is sparse.

Combining (2.3a) and (2.3b), we obtain

where

which is called the iteration matrix (or error propagation matrix) of Algorithm 1. Then, it holds that

with

The quantity (2.8) is referred to as the $M$-convergence factor of Algorithm 1, which will be analyzed in the next section.

In this section, we analyze the convergence of Algorithm 1 with prolongation $P_{\star}$. A succinct identity for characterizing the convergence factor $\|E_{\rm TG}\|_{M}$ is established, followed by discussions on how the identity can be used to derive a class of optimal restriction operators and analyze the influence of $\operatorname*{range}(R^{T})$ on $\|E_{\rm TG}\|_{M}$.

We start with a technical lemma, which provides a sufficient and necessary condition for $\|E_{\rm TG}\|_{M}<1$ (see the proof of Theorem 3.3 for details).

Using (2.8) and Lemma 3.1, we can derive the following identity.

The identity (3.3) provides a straightforward approach to analyzing the convergence properties of Algorithm 1. Of particular interest is a restriction operator that minimizes the convergence factor $\|E_{\rm TG}\|_{M}$, provided that $M$ is preselected.

The following inequality, known as the Poincaré separation theorem (see, e.g., [11, Corollary 4.3.37]), will be used to analyze the optimal restriction.

Based on (3.3) and Lemma 3.4, we can establish the following optimal restriction theory.

To analyze the influence of $\operatorname*{range}(R^{T})$ on $\|E_{\rm TG}\|_{M}$, we need the following lemma; see, e.g., [11, Corollary 4.3.5].

The following theorem shows that $\|E_{\rm TG}\|_{M}$ decreases when expanding the row space of $R$.