# Symbol-based Convergence Analysis in Block Multigrid Methods with
Applications for Stokes Problems111This work was partially supported by INdAM-
GNCS

Matthias Bolten<EMAIL_ADDRESS>School of Mathematics and Natural
Sciences, University of Wuppertal, Wuppertal, Germany  Marco Donatelli
<EMAIL_ADDRESS>Department of Science and High Technology.
University of Insubria, Como, Italy.  Paola Ferrari
<EMAIL_ADDRESS>Isabella Furci<EMAIL_ADDRESS>
###### Abstract

The main focus of this paper is the study of efficient multigrid methods for
large linear systems with a particular saddle-point structure. Indeed, when
the system matrix is symmetric, but indefinite, the variational convergence
theory that is usually used to prove multigrid convergence cannot be directly
applied.

However, different algebraic approaches analyze properly preconditioned
saddle-point problems, proving convergence of the Two-Grid method. In
particular, this is efficient when the blocks of the coefficient matrix
possess a Toeplitz or circulant structure. Indeed, it is possible to derive
sufficient conditions for convergence and provide optimal parameters for the
preconditioning of the saddle-point problem in terms of the associated
generating symbols. In this paper, we propose a symbol-based convergence
analysis for problems that have a hidden block Toeplitz structure. Then, they
can be investigated focusing on the properties of the associated generating
function $\mathbf{f}$, which consequently is a matrix-valued function with
dimension depending on the block size of the problem. As numerical tests we
focus on the matrix sequence stemming from the finite element approximation of
the Stokes problem. We show the efficiency of the methods studying the hidden
9-by-9 block multilevel structure of the obtained matrix sequence. Moreover we
propose an efficient algebraic multigrid method with convergence rate
independent of the matrix size. Finally, we present several numerical tests
comparing the results with state-of-the-art strategies.

###### keywords:

Multigrid methods, block matrix sequences, spectral symbol, Stokes equation

###### MSC:

[2010] 65N55, 15B05, 35P20, 34M40

††journal: Applied Numerical Mathematics

[raggedentrytext]toclinesection,subsection,subsubsection,paragraph,subparagraph

##  1 Introduction

Many works have focused on efficient iterative methods for the solution of
large linear systems with coefficient matrix in saddle-point form [3, 19, 5].
Indeed, many preconditioned Krylov subspace methods have been developed based
on the study of the spectral properties of the coefficient matrix [3, 4, 6,
23]. When we deal with multilevel and multilevel block systems the performance
of some known preconditioners can deteriorate [20]. Multigrid methods (MGM)
are in general a good choice in terms of computational cost and optimal
convergence rate. Indeed, optimal convergence rates are obtained by
constructing and employing a proper sequence of linear systems of decreasing
dimensions. Moreover, when the coefficient matrix has a particular multilevel
(block) Toeplitz structure, the convergence analysis for such methods can be
obtained in compact form exploiting the concept of symbols [15, 16].

The knowledge and the use of the symbol has been successfully employed for
designing very efficient preconditioning and multigrid techniques in various
settings and several applications [11, 12, 14, 27]. However, even in the case
where the coefficient matrix is somehow related to a Toeplitz/Toeplitz-like
structure, some saddle-point problems are intrinsically indefinite. This
limits the possible convergence results that can be exploited.

In many works [3, 4, 11] a common choice for problems with saddle-point
structure is that of constructing preconditioning for two-by-two block systems
involving the Schur complement. In many cases stemming from the approximation
of Partial Differential Equations (PDEs), even though the whole problem is
indefinite, the latter procedure is based on the solution of a linear system
with a positive definite coefficient matrix.

In this paper we consider the approach suggested by [19] which provides a
point-wise smoother and a coarse grid correction to the linear system after
preconditioning from the left and the right using lower and upper triangular
matrices, respectively.

This strategy is particularly effective in the case where the blocks of the
matrix possess a scalar unilevel Toeplitz and circulant structure. Indeed, in
[7] sufficient conditions for the Two-Grid Method (TGM) convergence and “level
independency” for the W-cycle have been obtained in terms of the generating
functions of the blocks. Studying the latter functions, also the optimal
parameters for the preconditioning of the saddle-point problem and for the
point smoother were found. Clearly, in the scalar setting, the involved
generating functions were univariate or multivariate scalar trigonometric
polynomials. Consequently, the exploited MGM theory is that for structured
unilevel or multilevel scalar matrices [1, 2]. When discretizing a PDE with
higher order approximation techniques, the coefficient matrix - or its blocks
- possess a (multilevel) block structure.

In this paper we focus on saddle-point problems where the blocks of the
coefficient matrix have a (hidden) multilevel block structure. Consequently,
we combine the method given in [19] with the recent results on TGM convergence
analysis for block structured matrices [8, 12]. Moreover, we show that the
structure can be preserved on the coarse level, allowing for a recursive
application of the method, after proper symmetrization. The global algorithm
is then a 3-step procedure whose main body is formed by a symbol-based block
multigrid method.

The outline of the paper is as follows. In Subsections 1.1– 1.4 we provide the
notation and the necessary results on multigrid methods for both saddle-point
and symmetric definite positive (SPD) matrices. Regarding SPD matrices we
focus especially on the case where the multigrid is applied to Toepliz
matrices with both scalar or matrix-valued generating functions. At the end of
the Section 1 we report the main results of the TGM when applied to saddle-
point problems. In Section 2 we provide TGM convergence results for saddle-
point systems with blocks that are multilevel block structured matrices. In
Section 3 we introduce the Stokes problem in 2D and we focus on a Finite
Element approximation in Subsection 3.1. The discrete problem involves solving
a saddle-point linear system with block multilevel Toeplitz blocks. We
theoretically show that such a problem is within the setting and fulfils the
assumptions of the TGM convergence results described in Section 2. In Section
4, we introduce an iterative 3-step procedure that can be exploited for
solving recursively linear systems with structure described as in Subsection
3.1. Precisely, we focus on the structure and on the spectral features of the
involved matrix sequences in order to prove the efficiency of the multigrid
strategy with the proposed projection operators given in Subsection 4.2.
Finally, we conclude with Section 5 presenting selected numerical tests that
confirm the efficiency of the proposed multigrid method in comparison with a
state-of-the-art multigrid strategy for saddle-point structured problems.

###  1.1 Notation and preliminary results

This subsection is devoted to fix the notation and present the key concepts on
multigrid methods and Toeplitz matrices.

We report the main definitions and properties of Toeplitz and circulant
matrices generated by a function $f$. The domain of $f$ can be either one-
dimensional $[-\pi,\pi]$ or $d$-dimensional $[-\pi,\pi]^{d}$ and this gives
either unilevel or $d$-level structured matrices, respectively. The co-domain
of $f$ can be either the complex field or the linear space of $d\times d$
complex matrices, generating scalar or block structured matrices,
respectively. We will denote the generating function by bold $\mathbf{f}$ when
it is a matrix-valued function. When we want to highlight that the generating
function is multivariate, we specify its argument
$\mathbf{f}(\boldsymbol{\theta})$ in bold, where
$\boldsymbol{\theta}=(\theta_{1},\theta_{2},\dots,\theta_{d})\in[-\pi,\pi]^{d}$.
The structure and the size of the Toeplitz matrix is clearly related to its
generating function.

In the whole paper we adopt the following notation on norms and relations
between Hermitian Positive Definite (HPD) matrices. If
$X\in\mathbb{C}^{N\times N}$ is a HPD matrix, then
$\|\cdot\|_{X}=\|X^{1/2}\cdot\|_{2}$ denotes the Euclidean norm weighted by
$X$ on $\mathbb{C}^{N}$. If $X$ and $Y$ are Hermitian matrices, then the
notation $X\leq Y$ means that $Y-X$ is a nonnegative definite matrix. Finally,
given a matrix-valued function $\mathbf{f}$ defined on the cube
$[-\pi,\pi]^{m}$ we denote by $\|\mathbf{f}\|_{\infty}={\rm
ess\,sup}_{\boldsymbol{\theta}\in[-\pi,\pi]^{m}}\|\mathbf{f}(\boldsymbol{\theta})\|_{2}$.

In what follows, we report the formal definition of Toeplitz and circulant
matrices.

######  Definition 1.

Given a Lebesgue integrable function
$\mathbf{f}:[-\pi,\pi]^{d}\rightarrow\mathbb{C}^{s\times s}$, its Fourier
coefficients are given by

$\hat{\mathbf{f}}_{\mathbf{k}}=\frac{1}{(2\pi)^{d}}\int_{[-\pi,\pi]^{d}}\mathbf{f}(\boldsymbol{\theta}){\rm
e}^{-\iota\left\langle{\mathbf{k}},\boldsymbol{\theta}\right\rangle}\mathrm{d}\boldsymbol{\theta}\in\mathbb{C}^{s\times
s},\qquad\mathbf{k}=(k_{1},\ldots,k_{d})\in\mathbb{Z}^{d},$

where $\boldsymbol{\theta}=(\theta_{1},\ldots,\theta_{d})$,
$\left\langle\mathbf{k},\boldsymbol{\theta}\right\rangle=\sum_{i=1}^{d}k_{i}\theta_{i}$,
and the integrals of matrices are computed elementwise.

The $\mathbf{n}$th Toeplitz matrix, where $\mathbf{n}=(n_{1},\ldots,n_{d})$,
associated with $\mathbf{f}$ is the matrix of dimension $s\cdot
n_{1}\cdot\ldots\cdot n_{d}$ given by

$T_{\mathbf{n}}(\mathbf{f})=\sum_{\mathbf{1}-\mathbf{n}\leq\mathbf{k}\leq\mathbf{n}-\mathbf{1}}J_{n_{1}}^{k_{1}}\otimes\cdots\otimes
J_{n_{d}}^{k_{d}}\otimes\hat{\mathbf{f}}_{\mathbf{k}},$

where $\mathbf{1}$ is the vector of all ones, $J^{k_{\xi}}_{n_{\xi}}$ is the
$n_{\xi}\times n_{\xi}$ matrix whose $(i,h)$th entry equals 1 if
$(i-h)=k_{\xi}$ and $0$ otherwise, and where $\mathbf{s}\leq\mathbf{t}$ means
that $s_{j}\leq t_{j}$ for any $j=1,\ldots,d$.

The $\mathbf{n}$th circulant matrix associated with $\mathbf{f}$ is the matrix
of dimension $s\cdot n_{1}\cdot\ldots\cdot n_{d}$ given by

$\mathcal{C}_{\mathbf{n}}(\mathbf{f})=\sum_{\mathbf{1}-\mathbf{n}\leq\mathbf{k}\leq\mathbf{n}-\mathbf{1}}Z_{n_{1}}^{k_{1}}\otimes\cdots\otimes
Z_{n_{d}}^{k_{d}}\otimes\hat{\mathbf{f}}_{\mathbf{k}},$

where $Z^{k_{\xi}}_{n_{\xi}}$ is the $n_{\xi}\times n_{\xi}$ matrix whose
$(i,h)$th entry equals 1 if $(i-h)\mod n_{\xi}=k_{\xi}$ and $0$ otherwise.

######  Theorem 1 ([10]).

Let $\textbf{f}\in L^{1}([-\pi,\pi]^{d},s)$ be a matrix-valued function with
$d\geq 1,s\geq 2$. Then, the following (block-Schur) decomposition of
$\mathcal{C}_{\bf n}(\textbf{f})$ is valid:

$\mathcal{C}_{\bf n}(\textbf{f})=(F_{\bf n}\otimes I_{s})D_{\bf
n}(\textbf{f})(F_{\bf n}\otimes I_{s})^{H},$ (1)

where

$D_{\bf n}(\textbf{f})=\underset{{{\bf 0}\leq{\bf k}\leq{\bf n}-{\bf 1}}}{{\rm
diag}}\left(S_{\bf n}(\textbf{f})\left(\theta_{\bf k}^{({\bf
n})}\right)\right),$ (2)

$F_{\bf n}$ is the d-level Fourier matrix, $F_{\bf
n}=F_{n_{1}}\otimes\dots\otimes F_{n_{d}}$, $I_{s}$ is the identity matrix of
size $s$ and

$\theta_{\bf k}^{({\bf n})}=2\pi\frac{{\bf k}}{{\bf n}},\quad{{\bf k}={\bf
0},\dots,{\bf n}-{\bf 1}}.$ (3)

Moreover, $S_{\bf n}(\textbf{f})(\cdot)$ is the $\bf n$-th Fourier sum of f
given by

$S_{\bf
n}(f)(\boldsymbol{\theta})=\sum_{j_{1}=1-n_{1}}^{n_{1}-1}\cdots\sum_{j_{d}=1-n_{d}}^{n_{d}-1}\hat{f}_{\bf
j}{\rm e}^{{\iota}\left\langle{\bf j},\boldsymbol{\theta}\right\rangle},\ \ \
\ \left\langle{\bf j},\theta\right\rangle=\sum_{t=1}^{d}j_{t}\theta_{t}.$ (4)

######  Remark 1.

The latter results states that the eigenvalues of $\mathcal{C}_{\bf
n}(\textbf{f})$ are given by the evaluations of $\lambda_{t}(S_{\bf
n}(\textbf{f})(\cdot))$, $t=1,\ldots,s$, at the grid points (3).

Moreover, if f is a trigonometric polynomial of fixed degree (with respect to
$\bf n$), then it is worth noticing that $S_{\bf
n}(\textbf{f})(\cdot)=\textbf{f}(\cdot)$ for $\bf n$ large enough. Therefore,
in such a setting, the eigenvalues of $\mathcal{C}_{\bf n}(\textbf{f})$ are
the evaluations of $\lambda_{t}(\textbf{f}(\cdot))$, $t=1,\ldots,s$, at the
very same grid points. The latter implies that if $\textbf{f}\geq 0$
trigonometric polynomial and $\lambda_{\min}(\textbf{f})$ does not vanish at
any of the grid points (3), then $\mathcal{C}_{\bf n}(\textbf{f})$ is HPD.

######  Remark 2.

In Definition 1, we considered square matrix-valued functions. In the general
case where $\mathbf{f}:[-\pi,\pi]^{d}\rightarrow\mathbb{C}^{s_{1}\times
s_{2}}$, when we write $\mathcal{C}_{\mathbf{n}}(\mathbf{f})$, we mean that

$\mathcal{C}_{\mathbf{n}}(\mathbf{f})=\boldsymbol{\Pi}^{T}_{s_{1}n^{d}\times
s_{1}n^{d}}\begin{bmatrix}\mathcal{C}_{\mathbf{n}}({f}^{(1,1)})&\mathcal{C}_{\mathbf{n}}({f}^{(1,2)})&\dots&\mathcal{C}_{\mathbf{n}}({f}^{(1,s_{2})})\\\
\mathcal{C}_{\mathbf{n}}({f}^{(2,1)})&\mathcal{C}_{\mathbf{n}}({f}^{(2,2)})&\dots&\mathcal{C}_{\mathbf{n}}({f}^{(2,s_{2})})\\\
\vdots&&\ddots&\vdots\\\
\mathcal{C}_{\mathbf{n}}({f}^{(s_{1},1)})&\mathcal{C}_{\mathbf{n}}({f}^{(s_{1},2)})&\dots&\mathcal{C}_{\mathbf{n}}({f}^{(s_{1},s_{2})})\end{bmatrix}\boldsymbol{\Pi}_{s_{2}n^{d}\times
s_{2}n^{d}},$ (5)

where ${f}^{(j,k)}$ is the $(j,k)$th component of the function $\mathbf{f}$
and, for $s\in\\{s_{1},s_{2}\\}$, the permutation matrices
$\boldsymbol{\Pi}_{sn^{d}\times sn^{d}}$ are defined by

$\boldsymbol{\Pi}_{sn^{d}\times sn^{d}}=[\Pi_{1}|\Pi_{2}|\dots|\Pi_{n^{d}}],$
(6)

with

$[\Pi_{k}]_{i,j}=\begin{cases}1&\mbox{if }i=(j-1)n^{d}+k,\\\
0&\mbox{otherwise},\end{cases},\qquad i=1,\dots,sn^{d},\quad j=1,\dots,s.$

This approach is consistent with Definition 1, however, we avoided to include
non-square matrices in the definition of circulants because in the next
sections we will exploit the circulant algebra structure, which is natural in
the square case.

######  Remark 3.

In the non-square case it is still possible to exploit symbols to perform
products between circulant matrices as follows. Suppose
$\mathbf{f}:[-\pi,\pi]^{d}\rightarrow\mathbb{C}^{s_{1}\times s_{2}}$,
$\mathbf{g}:[-\pi,\pi]^{d}\rightarrow\mathbb{C}^{s_{2}\times s_{3}}$ and
denote the components of the product $\mathbf{f}\mathbf{g}$ with
${h}^{(i,j)}$, then

$\displaystyle\mathcal{C}_{\mathbf{n}}(\mathbf{f})\mathcal{C}_{\mathbf{n}}(\mathbf{g})=$
$\displaystyle\boldsymbol{\Pi}^{T}_{s_{1}n^{d}\times
s_{1}n^{d}}\begin{bmatrix}\mathcal{C}_{\mathbf{n}}({f}^{(i,j)})\end{bmatrix}_{\begin{subarray}{c}i=1,\dots,s_{1}\\\
j=1,\dots,s_{2}\end{subarray}}\boldsymbol{\Pi}_{s_{2}n^{d}\times
s_{2}n^{d}}\cdot$ $\displaystyle\cdot\boldsymbol{\Pi}^{T}_{s_{2}n^{d}\times
s_{2}n^{d}}\begin{bmatrix}\mathcal{C}_{\mathbf{n}}({g}^{(i,j)})\end{bmatrix}_{\begin{subarray}{c}i=1,\dots,s_{1}\\\
j=1,\dots,s_{2}\end{subarray}}\boldsymbol{\Pi}_{s_{3}n^{d}\times s_{3}n^{d}}$
$\displaystyle=$ $\displaystyle\boldsymbol{\Pi}^{T}_{s_{1}n^{d}\times
s_{1}n^{d}}\begin{bmatrix}\mathcal{C}_{\mathbf{n}}({f}^{(i,j)})\end{bmatrix}_{\begin{subarray}{c}i=1,\dots,s_{1}\\\
j=1,\dots,s_{2}\end{subarray}}\begin{bmatrix}\mathcal{C}_{\mathbf{n}}({g}^{(i,j)})\end{bmatrix}_{\begin{subarray}{c}i=1,\dots,s_{1}\\\
j=1,\dots,s_{2}\end{subarray}}\boldsymbol{\Pi}_{s_{3}n^{d}\times s_{3}n^{d}}$
$\displaystyle=$ $\displaystyle\boldsymbol{\Pi}^{T}_{s_{1}n^{d}\times
s_{1}n^{d}}\begin{bmatrix}\mathcal{C}_{\mathbf{n}}({h}^{(i,j)})\end{bmatrix}_{\begin{subarray}{c}i=1,\dots,s_{1}\\\
j=1,\dots,s_{2}\end{subarray}}\boldsymbol{\Pi}_{s_{3}n^{d}\times s_{3}n^{d}},$
$\displaystyle=$
$\displaystyle\mathcal{C}_{\mathbf{n}}(\mathbf{f}\mathbf{g}).$

The generating function plays an important role for investigating some
analytic properties of the related Toeplitz matrix. Under particular
hypotheses, $f$ can provide important spectral information on $T_{n}(f)$. In
this case we say that $f$ is the spectral symbol of the matrix sequence
$\\{T_{n}(f)\\}_{n}$. In the following, we report an important localization
result on the eigenvalues of $T_{\bf n}(\mathbf{f})$ generated by a Hermitian
$d\times d$ matrix-valued function.

######  Theorem 2 ([24, 25]).

Let $\mathbf{f}\in L^{1}([-\pi,\pi]^{m})$ be a Hermitian $d\times d$ matrix-
valued function with $m\geq 1,d\geq 2$. Let $m_{1}$ be the essential infimum
of the minimal eigenvalue of $\mathbf{f}$, $M_{1}$ be the essential supremum
of the minimal eigenvalue of $\mathbf{f}$, $m_{d}$ be the essential infimum of
the maximal eigenvalue of $\mathbf{f}$, and $M_{d}$ be the essential supremum
of the maximal eigenvalue of $\mathbf{f}$.

  1. 1.

If $m_{1}<M_{1}$ then all the eigenvalues of $T_{\bf n}(\mathbf{f})$ belong to
the open set $(m_{1},M_{d}]$ for every ${\bf n}\in\mathbb{N}^{m}$. If
$m_{d}<M_{d}$ then all the eigenvalues of $T_{\bf n}(\mathbf{f})$ belong to
the open set $[m_{1},M_{d})$ for every ${\bf n}\in\mathbb{N}^{m}$.

  2. 2.

If $m_{1}=0$ and $\boldsymbol{\theta}_{0}$ is the unique zero of
$\lambda_{\min}(\mathbf{f})$ such that there exist positive constants $c,C,t$
for which

$c\|\boldsymbol{\theta}-\boldsymbol{\theta}_{0}\|^{t}\leq\lambda_{\min}(f(\boldsymbol{\theta}))\leq
C\|\boldsymbol{\theta}-\boldsymbol{\theta}_{0}\|^{t},$

then the minimal eigenvalue of $T_{\bf n}(\mathbf{f})$ goes to zero as
$(n_{1}n_{2}\dots n_{m})^{-t/m}$.

###  1.2 Multigrid methods

In this subsection we briefly report the relevant results concerning the
convergence theory of algebraic multigrid methods [22] especially when the
coefficient matrix is a block-Toeplitz matrix generated by a matrix-valued
trigonometric polynomial [8].

In the general case, we are interested in solving a linear system
$A_{N_{0}}x_{N_{0}}=b_{N_{0}}$ where $A_{N_{0}}$ is a Hermitian positive
definite matrix. Assume $N_{0}>N_{1}>\dots>N_{\ell}>\dots N_{\ell_{\min}}$ and
define two sequences of full-rank rectangular matrices
$P_{N_{\ell},N_{\ell+1}}\in\mathbb{C}^{N_{\ell}\times N_{\ell+1}}$ and
$R_{N_{\ell+1},N_{\ell}}\in\mathbb{C}^{N_{\ell+1}\times N_{\ell}}$. Let us
consider the stationary iterative methods $\mathcal{V}_{N_{\ell},\rm{pre}}$,
with iteration matrix ${V}_{N_{\ell},\rm{pre}}$, and
$\mathcal{V}_{N_{\ell},\rm{post}}$, with iteration matrix
${V}_{N_{\ell},\rm{post}}$.

An iteration of an algebraic Multigrid Method (MGM) is given by Algorithm 1.

Algorithm 1
MGM$(A_{N_{\ell}},\mathcal{V}_{N_{\ell},\rm{pre}}^{\nu_{\rm{pre}}},\mathcal{V}_{N_{\ell},\rm{post}}^{\nu_{\rm{post}}},R_{N_{\ell+1},N_{\ell}},P_{N_{\ell},N_{\ell+1}},b_{N_{\ell}},x_{N_{\ell}}^{(j)},\ell)$

0\. if $\ell=\ell_{\min}$

1\.
$\,\quad{\rm\mathbf{solve}}\,A_{N_{\ell}}x_{N_{\ell}}^{(j+1)}=b_{N_{\ell}}$

2\. else

3\.
$\,\quad\tilde{x}_{N_{\ell}}=\mathcal{V}_{N_{\ell},\rm{pre}}^{\nu_{\rm{pre}}}(A_{N_{\ell}},{b}_{N_{\ell}},x_{N_{\ell}}^{(j)})$

4\. $\,\quad r_{N_{\ell}}=b_{N_{\ell}}-A_{N_{\ell}}\tilde{x}_{N_{\ell}}$

5\. $\,\quad r_{N_{\ell+1}}=R_{N_{\ell+1},N_{\ell}}r_{N_{\ell}}$

6\. $\,\quad
A_{N_{\ell+1}}=R_{N_{\ell+1},N_{\ell}}A_{N_{\ell}}P_{N_{\ell},N_{\ell+1}}$

7\.  for $i=1:\gamma$

8\.
MGM$(A_{N_{\ell+1}},\mathcal{V}_{N_{\ell+1},\rm{pre}}^{\nu_{\rm{pre}}},\mathcal{V}_{N_{\ell+1},\rm{post}}^{\nu_{\rm{post}}},R_{N_{\ell+2},N_{\ell+1}},P_{N_{\ell+1},N_{\ell+2}},b_{N_{\ell+1}},x_{N_{\ell+1}}^{(j)},\ell+1)$

9\. end

10\.
$\quad\hat{x}_{N_{\ell}}=\tilde{x}_{N_{\ell}}+P_{N_{\ell},N_{\ell+1}}y_{N_{\ell+1}}$

11\. $\quad
x_{N_{\ell}}^{(j+1)}=\mathcal{V}_{N_{\ell},\rm{post}}^{\nu_{\rm{post}}}(A_{N_{\ell}},{b}_{N_{\ell}},\hat{x}_{N_{\ell}})$

The steps in lines $3$ and $11$ consist, respectively, in applying
$\nu_{\rm{pre}}$ times a pre-smoother and $\nu_{\rm{post}}$ times a post-
smoother of the given iterative methods. The steps in lines $4$ to $10$ define
the “coarse grid correction”. When $\gamma$ and $\ell_{\min}$ are equal to 1
the Algorithm 1 represents one iteration of the Two-Grid method (TGM). In this
case a linear system with coefficient matrix $A_{N_{1}}$ has to be solved.
Consequently, if the coefficient matrix has a large size, $\ell_{\min}$ should
be taken such that $N_{\ell_{\min}}$ becomes small enough for solving cheaply
step 1. In the latter situation, the choice $\gamma=1$ defines the V-cycle
method and the choice $\gamma=2$ defines the W-cycle method.

The next Theorem provides the conditions which are sufficient for the
convergence of the TGM in the particular case where the matrices
$A_{N_{\ell+1}}$ at the coarse levels are computed according to the Galerkin
approach ($R_{N_{\ell+1},N_{\ell}}=P_{N_{\ell},N_{\ell+1}}^{H}$). Moreover,
the theorem guarantees the linear convergence of the method, which means that
the number of iterations needed by the algorithm to reach a given accuracy
$\epsilon$ is bounded from above by a constant independent of the size of the
problem. In this case we will say that the method is convergent and optimal.

######  Theorem 3.

([22]) Let $A_{N_{\ell}}$ and $V_{N_{\ell},{\rm post}}$ be defined as in the
TGM algorithm and assume that no pre-smoothing is applied
$(V_{N_{\ell},\rm{pre}}^{\nu_{\rm{pre}}}=I_{N_{\ell}})$. Let us define

${\rm
TGM}(A_{N_{\ell}},I_{N_{\ell}},V_{N_{\ell},\rm{post}}^{\nu_{\rm{post}}},P_{N_{\ell},N_{\ell+1}})\\\
=V_{N_{\ell},\rm{post}}^{\nu_{\rm{post}}}\left[I_{N_{\ell}}-P_{N_{\ell},N_{\ell+1}}\left(P_{N_{\ell},N_{\ell+1}}^{H}A_{N_{\ell}}P_{N_{\ell},N_{\ell+1}}\right)^{-1}P_{N_{\ell},N_{\ell+1}}^{H}A_{N_{\ell}}\right],$

the iterative matrix of the TGM. Assume that

  * (a)

the smoother fulfills the smoothing property:

$\exists\sigma_{\rm{post}}>0\,:\;\|V_{N_{\ell},\rm{post}}x_{N_{\ell}}\|_{A_{N_{\ell}}}^{2}\leq\|x_{N_{\ell}}\|_{A_{N_{\ell}}}^{2}-\sigma_{\rm{post}}\|x_{N_{\ell}}\|_{A_{N_{\ell}}^{2}}^{2},\quad\forall
x_{N_{\ell}}\in\mathbb{C}^{N_{\ell}}.$

  * (b)

The grid transfer operator fulfills the approximation property:

$\exists\gamma>0\,:\;\min_{y\in\mathbb{C}^{N_{\ell+1}}}\|x_{N_{\ell}}-P_{N_{\ell},N_{\ell+1}}y\|_{2}^{2}\leq\gamma\|x_{N_{\ell}}\|_{A_{N_{\ell}}}^{2},\quad\forall
x_{N_{\ell}}\in\mathbb{C}^{N_{\ell}}.$

Then $\gamma\geq\sigma_{\rm{post}}$ and

$\displaystyle\|{\rm
TGM}(A_{N_{\ell}},I_{N_{\ell}},V_{N_{\ell},\rm{post}}^{\nu_{\rm{post}}},P_{N_{\ell},N_{\ell+1}})\|_{A_{N_{\ell}}}\leq\sqrt{1-\sigma_{\rm{post}}/\gamma}.$

###  1.3 Multigrid methods for block circulant and block Toeplitz systems

If the coefficient matrix is of the form
$\mathcal{C}_{\mathbf{n}}(\mathbf{f})$ or $T_{\mathbf{n}}(\mathbf{f})$, with
$\mathbf{f}$ being a matrix-valued trigonometric polynomial, $\mathbf{f}\geq
0$, sufficient conditions for the convergence and optimality of the TGM can be
expressed in relation to $\mathbf{f}$.

In particular, in [8] the authors derive sufficient conditions that should be
fulfilled by the grid transfer operators in order to validate the
approximation property (b) of Theorem 3. In the following we report a brief
description of the result [8] in the general multilevel block Toeplitz
setting.

Suppose $\mathbf{f}$ is a $s\times s$ matrix-valued trigonometric polynomial,
$\mathbf{f}\geq 0$, such that there exist unique
$\boldsymbol{\theta}_{0}\in[0,2\pi)^{d}$ and $\bar{\jmath}\in\\{1,\dots,s\\}$
such that

$\left\\{\begin{array}[]{ll}\lambda_{j}(\mathbf{f}(\boldsymbol{\theta}))=0&\mbox{for
}\boldsymbol{\theta}=\boldsymbol{\theta}_{0}\mbox{ and }j=\bar{\jmath},\\\
\lambda_{j}(\mathbf{f}(\boldsymbol{\theta}))>0&{\rm
otherwise}.\end{array}\right.$ (7)

The latter means that only one eigenvalue function of $\mathbf{f}$ has exactly
one zero in $\boldsymbol{\theta}_{0}$ and $\mathbf{f}$ is positive definite in
$[0,2\pi)^{{d}}\backslash\\{\boldsymbol{\theta}_{0}\\}$. Define
$q_{\bar{\jmath}}(\boldsymbol{\theta}_{0})$ as the eigenvector function of
$\mathbf{f}(\boldsymbol{\theta}_{0})$ associated with
$\lambda_{\bar{\jmath}}(\mathbf{f}(\boldsymbol{\theta}_{0}))=0$. Moreover,
define
$\Omega(\boldsymbol{\theta})=\left\\{\boldsymbol{\theta}+\pi\boldsymbol{\eta},\,\boldsymbol{\eta}\in\\{0,1\\}^{{d}}\right\\}$.
Let us define a matrix-valued trigonometric polynomial $\mathbf{p}(\cdot)$
such that

  * 1.

$\sum_{\boldsymbol{\xi}\in\Omega(\boldsymbol{\theta})}\mathbf{p}(\boldsymbol{\xi})^{H}\mathbf{p}(\boldsymbol{\xi})>0,\quad\forall\,\boldsymbol{\theta}\in[0,2\pi)^{{d}},$
(8)

which implies that the trigonometric function

$\mathbf{s}(\boldsymbol{\theta})=\mathbf{p}(\boldsymbol{\theta})\left(\sum_{\boldsymbol{\xi}\in\Omega(\boldsymbol{\theta})}\mathbf{p}(\boldsymbol{\xi})^{H}\mathbf{p}(\boldsymbol{\xi})\right)^{-1}\mathbf{p}(\boldsymbol{\theta})^{H}$

is well-defined for all $\boldsymbol{\theta}\in[0,2\pi)^{{d}}$.

  * 2.

$\mathbf{s}(\boldsymbol{\theta}_{0})q_{\bar{\jmath}}(\boldsymbol{\theta}_{0})=q_{\bar{\jmath}}(\boldsymbol{\theta}_{0}).$
(9)

  * 3.

$\limsup_{\boldsymbol{\theta}\rightarrow\boldsymbol{\theta}_{0}}\lambda_{\bar{\jmath}}(\mathbf{f}(\boldsymbol{\theta}))^{-1}(1-\lambda_{\bar{\jmath}}(\mathbf{s}(\boldsymbol{\theta})))=c,$
(10)

where $c\in\mathbb{R}$ is a constant.

In order to have optimal convergence of the TGM it is sufficient compute a
grid transfer operator as

$P_{N_{\ell},N_{\ell+1}}=T_{\mathbf{n}_{\ell}}(\mathbf{p})\left(K_{N_{\ell},N_{\ell+1}}\otimes
I_{s}\right),$ (11)

where $K_{N_{\ell},N_{\ell+1}}$ is a $N_{\ell}\times N_{\ell+1}$ multilevel
cutting matrix obtained as the Kronecker product of proper unilevel cutting
matrices and $T_{\mathbf{n}_{\ell}}(\mathbf{p})$ is a multilevel block-
Toeplitz matrix generated by $\mathbf{p}$.

###  1.4 TGM for saddle-point systems

In the previous section we focused on the case where the coefficient matrix is
a block structured and non-negative definite matrix. For matrices in saddle-
point form, Rüge-Stuben theory and Theorem 3 cannot be directly exploited. In
[19] a general TGM procedure has been presented for a system of the form

$\mathcal{A}=\begin{bmatrix}A&B^{T}\\\ B&-C\end{bmatrix}$ (12)

where $A$ is an $n\times n$ HPD matrix, $C$ is an $\hat{n}\times\hat{n}$ non-
negative definite matrix and $B$ is such that its has rank is $\hat{n}\leq n$
or that $C$ is positive definite on the null space of $B^{T}$. Indeed, we have
the following result.

######  Theorem 4 ([19]).

Let $\mathcal{A}$ be as in (12). Let $\alpha$ be a positive number such that
$\alpha<2(\lambda_{\max}(D_{A}^{-1}A))^{-1}$ and define

$\hat{\mathcal{A}}=\mathcal{LAU},\qquad\mathcal{L}=\begin{bmatrix}I_{n}\\\
\alpha
BD_{A}^{-1}&-I_{\hat{n}}\end{bmatrix},\qquad\mathcal{U}=\begin{bmatrix}I_{n}&-\alpha
D_{A}^{-1}B^{T}\\\ &I_{\hat{n}}\end{bmatrix},$ (13)

and

$\hat{C}=C+B(2\alpha D_{A}^{-1}-\alpha^{2}D_{A}^{-1}AD_{A}^{-1})B^{T}.$ (14)

Let $P_{A}$ and $P_{\hat{C}}$ be, respectively, $n\times k$ and
$\hat{n}\times\ell$ matrices of rank $k<n$ and $\ell<\hat{n}$, define the
prolongation

$\mathcal{P}=\begin{bmatrix}P_{A}&\\\ &P_{\hat{C}}\end{bmatrix}$ (15)

for the global system involving $\mathcal{\hat{A}}$ and suppose that the pairs
$(A,P_{A})$ and $(\hat{C},P_{\hat{C}})$ fulfill the approximation property as
in Theorem 3. Then, the spectral radius of the TGM iteration matrix using one
iteration of the damped Jacobi method with relaxation parameter $\omega$ as
post smoother satisfies

$\rho({\rm
TGM}(\mathcal{\hat{A}},\mathcal{P},\omega))\leq\max\left(1-\frac{\omega}{{\kappa}({A},P_{A})},\,1-\frac{\omega}{{\kappa}({\hat{C}},P_{\hat{C}})},\,\omega\hat{\gamma}_{A}-1,\,\omega\hat{\gamma}_{\hat{C}}-1,\,\sqrt{1-\frac{\omega(2-\omega\tilde{\gamma})}{\tilde{\kappa}}}\right),$
(16)

where

$\displaystyle\hat{\gamma}_{A}=\left(\alpha\left(2-\alpha\lambda_{\max}(D_{A}^{-1}A)\right)\right)^{-1},$
$\displaystyle\qquad\hat{\gamma}_{\hat{C}}=\lambda_{\max}\left(D_{\hat{C}}^{-1}(C+BA^{-1}B^{T})\right),$
$\displaystyle\tilde{\gamma}=\frac{2\hat{\gamma}_{A}\hat{\gamma}_{\hat{C}}}{\hat{\gamma}_{A}+\hat{\gamma}_{\hat{C}}},$
$\displaystyle\qquad\tilde{\kappa}=\frac{2{\kappa}({A},P_{A}){\kappa}({\hat{C}},P_{\hat{C}})}{{\kappa}({A},P_{A})+{\kappa}({\hat{C}},P_{\hat{C}})}.$

##  2 TGM for saddle-point systems in the multilevel block setting

The aim of the current section is to provide the theoretical background for
the development and analysis of a TGM for the Stokes problem discretized with
the technique of Subsection 3.1. The starting point is Theorem 4, which we
want to combine with the theory for multilevel block structured matrices of
Subsection 1.3. The derivations are an extension of the results in [7], where
the authors consider a general unilevel circulant setting and a multilevel
circulant setting which arises from the discretization of Stokes-like
equations. The generalization that we consider in what follows consists in
dealing with matrix-valued symbols instead of scalar-valued ones and this
forces us to give the conditions on the symbols in a different form, owing to
the non-commutativity of the matrix-matrix product.

The following results concern the properties of a matrix of the form

$\mathcal{A}=\begin{bmatrix}\mathcal{C}_{\mathbf{n}}(\mathbf{f}_{A_{x}})&O&\mathcal{C}_{\mathbf{n}}^{T}(\mathbf{f}_{B_{x}})\\\
O&\mathcal{C}_{\mathbf{n}}(\mathbf{f}_{A_{y}})&\mathcal{C}_{\mathbf{n}}^{T}(\mathbf{f}_{B_{y}})\\\
\mathcal{C}_{\mathbf{n}}(\mathbf{f}_{B_{x}})&\mathcal{C}_{\mathbf{n}}(\mathbf{f}_{B_{y}})&-\mathcal{C}_{\mathbf{n}}(\mathbf{f}_{C})\end{bmatrix}\in\mathbb{R}^{(2s_{a}+s_{c})n^{2}\times(2s_{a}+s_{c})n^{2}}$

with multilevel block circulant blocks generated by multivariate matrix-valued
functions. For the sake of readability, we focus on the case $d=2$,
$\mathbf{n}=(n,n)$ and, moreover, we fix $s_{a}\geq s_{c}$. An extension to
larger dimensional problems is straightforward, in this case additional blocks
will be part of $\mathcal{A}$, e.g., for $d=3$ we will have an additional
block $A_{z}$ on the block diagonal and additionally $B_{z}$ and $B_{z}^{T}$
in the last block row and column, respectively.

######  Lemma 1.

Let $\tilde{A}=\begin{bmatrix}A_{x}&O\\\
O&A_{y}\end{bmatrix}\in\mathbb{R}^{2sn^{2}\times 2sn^{2}}$,
$A_{x}=\mathcal{C}_{\mathbf{n}}(\mathbf{f}_{A_{x}})$,
$A_{y}=\mathcal{C}_{\mathbf{n}}(\mathbf{f}_{A_{y}})$ with
$\mathbf{f}_{A_{x}}$, $\mathbf{f}_{A_{y}}$ $s\times s$ HPD matrix-valued bi-
variate trigonometric polynomials.

If $\alpha$ is a positive number such that

$\alpha<2\left(\max_{z\in\\{x,y\\}}\left\\{\max_{j=1,\dots,s}\left\\{\left(\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{z}})\right)^{-1}\right\\}\left\|\mathbf{f}_{A_{z}}\right\|_{\infty}\right\\}\right)^{-1}$

then the matrix $2\alpha
D^{-1}_{\tilde{A}}-\alpha^{2}D^{-1}_{\tilde{A}}\tilde{A}D^{-1}_{\tilde{A}}$ is
HPD, where $D_{\tilde{A}}=\operatorname*{diag}(\tilde{A})$.

###### Proof.

Since $\tilde{A}$ is HPD, then $D_{\tilde{A}}$ is HPD. Hence, writing

$2\alpha
D^{-1}_{\tilde{A}}-\alpha^{2}D^{-1}_{\tilde{A}}\tilde{A}D^{-1}_{\tilde{A}}=D^{-\frac{1}{2}}_{\tilde{A}}\left(2\alpha
I_{s}-\alpha^{2}D^{-\frac{1}{2}}_{\tilde{A}}\tilde{A}D^{-\frac{1}{2}}_{\tilde{A}}\right)D^{-\frac{1}{2}}_{\tilde{A}}$

we see that we need to prove that the Hermitian matrix $2\alpha
I_{s}-\alpha^{2}D^{-\frac{1}{2}}_{\tilde{A}}\tilde{A}D^{-\frac{1}{2}}_{\tilde{A}}$
is HPD. The eigenvalues of $2\alpha
I_{s}-\alpha^{2}D^{-\frac{1}{2}}_{\tilde{A}}\tilde{A}D^{-\frac{1}{2}}_{\tilde{A}}$
are of the form
$2\alpha-\alpha^{2}\lambda_{j}\left(D^{-\frac{1}{2}}_{\tilde{A}}\tilde{A}D^{-\frac{1}{2}}_{\tilde{A}}\right)$,
where
$\lambda_{j}\left(D^{-\frac{1}{2}}_{\tilde{A}}\tilde{A}D^{-\frac{1}{2}}_{\tilde{A}}\right)$
are the positive eigenvalues of
$D^{-\frac{1}{2}}_{\tilde{A}}\tilde{A}D^{-\frac{1}{2}}_{\tilde{A}}$. Recalling
that $\alpha>0$ by assumption, we have

$\lambda_{\min}\left(\alpha
I_{s}-\alpha^{2}D^{-\frac{1}{2}}_{\tilde{A}}\tilde{A}D^{-\frac{1}{2}}_{\tilde{A}}\right)=2\alpha-\alpha^{2}\lambda_{\max}\left(D^{-\frac{1}{2}}_{\tilde{A}}\tilde{A}D^{-\frac{1}{2}}_{\tilde{A}}\right),$

which implies that $\alpha
I_{s}-\alpha^{2}D^{-\frac{1}{2}}_{\tilde{A}}\tilde{A}D^{-\frac{1}{2}}_{\tilde{A}}$
is HPD if

$0<\alpha<\frac{2}{\lambda_{\max}\left(D^{-\frac{1}{2}}_{\tilde{A}}\tilde{A}D^{-\frac{1}{2}}_{\tilde{A}}\right)}.$
(17)

By construction, $D^{-1}_{\tilde{A}}=\begin{bmatrix}D^{-1}_{A_{x}}&O\\\
O&D^{-1}_{A_{y}}\end{bmatrix}$, with

$D^{-1}_{A_{x}}=I_{n^{2}}\otimes\operatorname*{diag}_{1\leq j\leq
s}\left(\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{x}})\right)^{-1},\quad
D^{-1}_{A_{y}}=I_{n^{2}}\otimes\operatorname*{diag}_{1\leq j\leq
s}\left(\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{y}})\right)^{-1},$

where
$\hat{a}_{0}^{(1,1)}(\mathbf{f}_{A_{x}}),\dots,\hat{a}_{0}^{(s,s)}(\mathbf{f}_{A_{x}})$
and
$\hat{a}_{0}^{(1,1)}(\mathbf{f}_{A_{y}}),\dots,\hat{a}_{0}^{(s,s)}(\mathbf{f}_{A_{y}})$
are the positive diagonal entries of the 0th block diagonal Fourier
coefficient of $\mathbf{f}_{A_{x}}$ and $\mathbf{f}_{A_{y}}$, respectively.
Consequently, we write

$\displaystyle\lambda_{\max}$
$\displaystyle\left(D^{-\frac{1}{2}}_{\tilde{A}}\tilde{A}D^{-\frac{1}{2}}_{\tilde{A}}\right)=$
$\displaystyle=\max_{z\in\\{x,y\\}}\lambda_{\max}\left(D^{-\frac{1}{2}}_{A_{z}}\tilde{A}D^{-\frac{1}{2}}_{A_{z}}\right)$
$\displaystyle=\max_{z\in\\{x,y\\}}\lambda_{\max}\left(\mathcal{C}_{\mathbf{n}}\left(\operatorname*{diag}_{1\leq
j\leq
s}\left(\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{z}})\right)^{-\frac{1}{2}}\mathbf{f}_{A_{z}}\operatorname*{diag}_{1\leq
j\leq
s}\left(\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{z}})\right)^{-\frac{1}{2}}\right)\right)$
$\displaystyle\leq\max_{z\in\\{x,y\\}}\left\|\operatorname*{diag}_{1\leq j\leq
s}\left(\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{z}})\right)^{-\frac{1}{2}}\mathbf{f}_{A_{z}}\operatorname*{diag}_{1\leq
j\leq
s}\left(\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{z}})\right)^{-\frac{1}{2}}\right\|_{\infty}$
$\displaystyle\leq\max_{z\in\\{x,y\\}}\left\\{\max_{j=1,\dots,s}\left\\{\left(\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{z}})\right)^{-1}\right\\}\left\|\mathbf{f}_{A_{z}}\right\|_{\infty}\right\\}.$

Hence, if

$0<\alpha<\frac{2}{\max_{z\in\\{x,y\\}}\left\\{\max_{j=1,\dots,s}\left\\{\left(\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{z}})\right)^{-1}\right\\}\left\|\mathbf{f}_{A_{z}}\right\|_{\infty}\right\\}}$

the parameter $\alpha$ fulfils condition (17) and the thesis of the Lemma
follows. ∎

######  Remark 4.

In order to follow the procedure analysed by Theorem 4, we need to study the
matrix $\hat{C}$ defined in (14). Theorem 4 requires that either $B$ is full-
rank or $C$ is HPD on the null space of $B^{T}$. In the next theorem, we will
require that $\hat{C}$ is HPD, but this implies that the latter condition on
$B$ and $C$ is fulfilled. Indeed, if $\hat{C}$ is HPD and there exists
$v\in\mathbb{R}^{s_{c}n^{2}}$ such that $v^{T}B=0$, then

$\displaystyle 0<v^{T}\hat{C}v<v^{T}Cv+v^{T}B\left(2\alpha
D_{A}^{-1}-\alpha^{2}D_{A}^{-1}AD_{A}^{-1}\right)B^{T}v<v^{T}{C}v.$

Conversely, note that if $B$ is full-rank or $C$ is HPD on the null space of
$B^{T}$, then $\hat{C}$ is HPD as long as $\alpha$ is chosen according to
Lemma 1.

We are now ready to state the main theorem concerning the convergence of a TGM
of the form of Theorem 4 in the multilevel block circulant setting.

######  Theorem 5.

Let us define the matrix

$\mathcal{A}=\begin{bmatrix}\tilde{A}&B^{T}\\\
B&-C\end{bmatrix}\in\mathbb{R}^{(2s_{a}+s_{c})n^{2}\times(2s_{a}+s_{c})n^{2}}$
(18)

having the blocks with the following structure:

  * 1.

$\tilde{A}=\begin{bmatrix}A_{x}&O\\\ O&A_{y}\end{bmatrix}$ where
$A_{x}=\mathcal{C}_{\mathbf{n}}(\mathbf{f}_{A_{x}})$,
$A_{y}=\mathcal{C}_{\mathbf{n}}(\mathbf{f}_{A_{y}})$ with
$\mathbf{f}_{A_{x}}$, $\mathbf{f}_{A_{y}}$ being $s_{a}\times s_{a}$ matrix-
valued bivariate trigonometric polynomials.

  * 2.

$C=\mathcal{C}_{\mathbf{n}}(\mathbf{f}_{C})$ with $\mathbf{f}_{C}$
$s_{c}\times s_{c}$ non-negative bivariate trigonometric polynomial.

  * 3.

$B=[B_{x},B_{y}]$ where $B_{x}=\mathcal{C}_{\mathbf{n}}(\mathbf{f}_{B_{x}})$,
$B_{y}=\mathcal{C}_{\mathbf{n}}(\mathbf{f}_{B_{y}})$ with $\mathbf{f}_{B_{x}}$
and $\mathbf{f}_{B_{y}}$ being $s_{a}\times s_{c}$ matrix-valued bivariate
trigonometric polynomials such that

$\mathbf{f}_{B_{x}}\left(\mathbf{f}_{A_{x}}\right)^{-1}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\left(\mathbf{f}_{A_{y}}\right)^{-1}\mathbf{f}_{B_{y}}^{H}\in
L^{\infty}\left([-\pi,\pi]^{2},s_{c}\right).$ (19)

Let $\alpha$ be a positive number such that

$\alpha<2\left(\max_{z\in\\{x,y\\}}\left\\{\max_{j=1,\dots,s}\left\\{\left(\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{z}})\right)^{-1}\right\\}\left\|\mathbf{f}_{A_{z}}\right\|_{\infty}\right\\}\right)^{-1}$
(20)

and define $\mathcal{\hat{A}}=\mathcal{LAU}$ where

$\mathcal{L}=\begin{bmatrix}I_{2s_{a}n^{2}}\\\ \alpha
BD_{\tilde{A}}^{-1}&-I_{s_{c}n^{2}}\end{bmatrix},\quad\mathcal{U}=\begin{bmatrix}I_{2s_{a}n^{2}}&-\alpha
D_{\tilde{A}}^{-1}B^{T}\\\ &I_{s_{c}n^{2}}\end{bmatrix}.$ (21)

Let $P_{\tilde{A}}$ and $P_{\hat{C}}$ be the full rank matrices defining the
prolongation $\mathcal{P}$ for the global system $\mathcal{\hat{A}}$. That is,

$\mathcal{P}=\begin{bmatrix}P_{\tilde{A}}&\\\ &P_{\hat{C}}\end{bmatrix},$

where $\hat{C}=C+B\left(2\alpha
D_{\tilde{A}^{-1}}-\alpha^{2}D_{\tilde{A}^{-1}}\tilde{A}D_{\tilde{A}^{-1}}\right)B^{T}=C_{\mathbf{n}}(\mathbf{f}_{\hat{C}})$,

$P_{\tilde{A}}=\begin{bmatrix}\mathcal{C}_{\mathbf{n}}(\mathbf{p}_{A_{x}})(K_{\mathbf{n}}^{T}\otimes
I_{s_{a}})&\\\
&\mathcal{C}_{\mathbf{n}}(\mathbf{p}_{A_{y}})(K_{\mathbf{n}}^{T}\otimes
I_{s_{a}})\end{bmatrix},\quad
P_{\hat{C}}=\mathcal{C}_{\mathbf{n}}(\mathbf{p}_{\hat{C}})(K_{\mathbf{n}}^{T}\otimes
I_{s_{c}})$

with $\mathbf{p}_{A_{x}}$, $\mathbf{p}_{A_{y}}$, and $\mathbf{p}_{\hat{C}}$
bi-variate matrix-valued trigonometric polynomials. Consider a TGM associated
with a basic scheme using only a single step of damped Jacobi post-smoothing
with relaxation parameter $\omega$. If

  1. 1.

the trigonometric polynomials $\mathbf{f}_{A_{x}}$, $\mathbf{f}_{A_{y}}$, and
$\mathbf{f}_{\hat{C}}$ are in the setting of (7) and they are full rank at all
grid points defined in (3);

  2. 2.

the pairs $(\mathbf{p}_{A_{x}},\mathbf{f}_{A_{x}})$,
$(\mathbf{p}_{A_{y}},\mathbf{f}_{A_{y}})$, and
$(\mathbf{p}_{\hat{C}},\mathbf{f}_{\hat{C}})$ fulfil relations (8)–(10);

  3. 3.

the relaxation parameter $\omega$ is a positive number such that

$\omega<2\min\left\\{2\alpha-\alpha^{2}\max_{j}\left\\{\frac{\|\mathbf{f}_{A_{x}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{x}})}},\frac{\|\mathbf{f}_{A_{y}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{y}})}}\right\\},\right.\\\
\left.\max_{j}\hat{a}_{0}^{(j,j)}(\mathbf{f}_{\hat{C}})\left\|\mathbf{f}_{C}+\mathbf{f}_{B_{x}}\mathbf{f}_{A_{x}}^{-1}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\mathbf{f}_{A_{y}}^{-1}\mathbf{f}_{B_{y}}^{H}\right\|_{\infty}^{-1}\right\\},$
(22)

then $\rho({\rm TGM}(\mathcal{\hat{A}},\mathcal{P},\omega))<1.$

###### Proof.

As a consequence of Remark 1, the assumptions on $\mathbf{f}_{A_{x}}$,
$\mathbf{f}_{A_{y}}$ and $\mathbf{f}_{\hat{C}}$ imply that $\tilde{A}$ is HPD.
Moreover, we see from Remark 4 that either $B$ is full-rank or $C$ is HPD on
the null space of $B^{T}$.

The assumptions on $\mathbf{p}_{A_{x}}$, $\mathbf{p}_{A_{y}}$,
$\mathbf{p}_{\hat{C}}$ and $\mathbf{f}_{A_{x}}$, $\mathbf{f}_{A_{y}}$,
$\mathbf{f}_{\hat{C}}$ ensure the validation of the approximation property by
[8, Section 6]. That is, the pairs $(\tilde{A},P_{\tilde{A}})$ and
$({\hat{C}},P_{\hat{C}})$ fulfil the approximation property in the 2D block
setting. Therefore, we can apply Theorem 4 and by (16) we have

$\rho({\rm
TGM}(\mathcal{\hat{A}},\mathcal{P},\omega))\leq\max\left(1-\frac{\omega}{{\kappa}(\tilde{A},P_{\tilde{A}})},1-\frac{\omega}{{\kappa}({\hat{C}},P_{\hat{C}})},\omega\hat{\gamma}_{\tilde{A}}-1,\omega\hat{\gamma}_{\hat{C}}-1,\sqrt{1-\frac{\omega(2-\omega\tilde{\gamma})}{\tilde{\kappa}}}\right).$

In conclusion, in order to prove that $\rho({\rm
TGM}(\mathcal{\hat{A}},\mathcal{P},\omega))<1$, we need to prove that each of
the quantities in the maximum is bounded from above by a constant strictly
smaller that 1.

The quantities $1-\frac{\omega}{{\kappa}({A},P_{A})}$ and
$1-\frac{\omega}{{\kappa}({\hat{C}},P_{\hat{C}})}$ are strictly smaller than 1
because the approximation property constants ${{\kappa}({A},P_{A})}$ and
${{\kappa}({\hat{C}},P_{\hat{C}})}$ are finite and positive.

Concerning the two terms $\omega\hat{\gamma}_{\tilde{A}}-1$ and
$\omega\hat{\gamma}_{\hat{C}}-1$ in the maximum, we estimate

$\begin{split}\hat{\gamma}_{\tilde{A}}&=\left(\alpha\left(2-\alpha\lambda_{\max}(D_{\tilde{A}}^{-1}\tilde{A})\right)\right)^{-1}=\\\
&\left(2\alpha-\alpha^{2}\lambda_{\max}\left(\begin{bmatrix}D^{-1}_{A_{x}}A_{x}&O\\\
O&D^{-1}_{A_{y}}A_{y}\end{bmatrix}\right)\right)^{-1}\leq\\\
&\left(2\alpha-\alpha^{2}\max_{j}\left\\{\frac{\|\mathbf{f}_{A_{x}}\|_{\infty}}{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{x}})},\frac{\|\mathbf{f}_{A_{y}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{y}})}}\right\\}\right)^{-1},\end{split}$
(23)
$\begin{split}\hat{\gamma}_{\hat{C}}=&\lambda_{\max}\left(D_{\hat{C}}^{-1}(C+B\tilde{A}^{-1}B^{T})\right)\leq\max_{j}\frac{1}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{\hat{C}})}}\lambda_{\max}\left(\mathcal{C}_{n}\left(\mathbf{f}_{C}+\mathbf{f}_{B_{x}}\mathbf{f}_{A_{x}}^{-1}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\mathbf{f}_{A_{y}}^{-1}\mathbf{f}_{B_{y}}^{H}\right)\right)\\\
\leq&\max_{j}\frac{1}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{\hat{C}})}}\left\|\mathbf{f}_{C}+\mathbf{f}_{B_{x}}\mathbf{f}_{A_{x}}^{-1}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\mathbf{f}_{A_{y}}^{-1}\mathbf{f}_{B_{y}}^{H}\right\|_{\infty},\end{split}$
(24)

where in the latter inequality we are using assumption (19) and the sub-
additivity of norms. Hence with the choice of $\omega$ in assumption 3., we
have

$\displaystyle\omega\hat{\gamma}_{\tilde{A}}-1<$ $\displaystyle
2\left(2\alpha-\alpha^{2}\max_{j}\left\\{\frac{\|\mathbf{f}_{A_{x}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{x}})}},\frac{\|\mathbf{f}_{A_{y}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{y}})}}\right\\}\right)\cdot$
$\displaystyle\cdot\left(2\alpha-\alpha^{2}\max_{j}\left\\{\frac{\|\mathbf{f}_{A_{x}}\|_{\infty}}{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{x}})},\frac{\|\mathbf{f}_{A_{y}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{y}})}}\right\\}\right)^{-1}-1$
$\displaystyle=$ $\displaystyle 1,$

and

$\displaystyle\omega\hat{\gamma}_{\hat{C}}-1<$ $\displaystyle
2\max_{j}\hat{a}_{0}^{(j,j)}(\mathbf{f}_{\hat{C}})\left\|\mathbf{f}_{C}+\mathbf{f}_{B_{x}}\mathbf{f}_{A_{x}}^{-1}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\mathbf{f}_{A_{y}}^{-1}\mathbf{f}_{B_{y}}^{H}\right\|_{\infty}^{-1}\cdot$
$\displaystyle\cdot\max_{j}\frac{1}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{\hat{C}})}}\left\|\mathbf{f}_{C}+\mathbf{f}_{B_{x}}\mathbf{f}_{A_{x}}^{-1}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\mathbf{f}_{A_{y}}^{-1}\mathbf{f}_{B_{y}}^{H}\right\|_{\infty}-1$
$\displaystyle=$ $\displaystyle 1.$

Finally, in order to prove that
$\sqrt{1-\frac{\omega(2-\omega\tilde{\gamma})}{\tilde{\kappa}}}<1$, we prove
that $2-\omega\tilde{\gamma}>0$. From hypothesis 3., we have

$\begin{split}\omega<&2\alpha-\alpha^{2}\max_{j}\left\\{\frac{\|\mathbf{f}_{A_{x}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{x}})}},\frac{\|\mathbf{f}_{A_{y}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{y}})}}\right\\}+\\\
&+\max_{j}\hat{a}_{0}^{(j,j)}(\mathbf{f}_{\hat{C}})\left\|\mathbf{f}_{C}+\mathbf{f}_{B_{x}}\mathbf{f}_{A_{x}}^{-1}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\mathbf{f}_{A_{y}}^{-1}\mathbf{f}_{B_{y}}^{H}\right\|_{\infty}^{-1}.\end{split}$
(25)

Moreover, using the estimations (23) and (24), we have

$\begin{split}\tilde{\gamma}&=\frac{2\hat{\gamma}_{A}\hat{\gamma}_{\hat{C}}}{\hat{\gamma}_{A}+\hat{\gamma}_{\hat{C}}}\\\
&\leq
2\left(\frac{\left(2\alpha-\alpha^{2}\max_{j}\left\\{\frac{\|\mathbf{f}_{A_{x}}\|_{\infty}}{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{x}})},\frac{\|\mathbf{f}_{A_{y}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{y}})}}\right\\}\right)^{-1}\max_{j}\frac{\left\|\mathbf{f}_{C}+\mathbf{f}_{B_{x}}\mathbf{f}_{A_{x}}^{-1}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\mathbf{f}_{A_{y}}^{-1}\mathbf{f}_{B_{y}}^{H}\right\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{\hat{C}})}}}{\left(2\alpha-\alpha^{2}\max_{j}\left\\{\frac{\|\mathbf{f}_{A_{x}}\|_{\infty}}{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{x}})},\frac{\|\mathbf{f}_{A_{y}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{y}})}}\right\\}\right)^{-1}+\max_{j}\frac{\left\|\mathbf{f}_{C}+\mathbf{f}_{B_{x}}\mathbf{f}_{A_{x}}^{-1}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\mathbf{f}_{A_{y}}^{-1}\mathbf{f}_{B_{y}}^{H}\right\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{\hat{C}})}}}\right),\end{split}$
(26)

where we used the fact that the function $(x,y)\mapsto 2\frac{xy}{x+y}$ is
increasing in $(0,+\infty)\times(0,+\infty)$ with respect to the component-
wise partial ordering in $\mathbb{R}^{2}$. Then, combining (25) and (26), we
obtain $\omega\tilde{\gamma}<2$. ∎

######  Remark 5.

In the previous theorem, some of the assumptions concern the pair
$\left({\hat{C}},P_{\hat{C}}\right)$. The matrix $\hat{C}$ is computed from
$A$, $B$ and $C$, so the required assumptions can be stated involving
explicitly the latter three matrices. For the sake of brevity, we discuss this
point here in a particular setting, which is relevant from the application
point of view. Indeed, we focus on the case where the matrix $\mathcal{A}$ in
(18) is such that $C$ is the null matrix and $s_{C}=1$, so that
${f}_{\hat{C}}$ is the scalar–valued function

${f}_{\hat{C}}=\mathbf{f}_{B_{x}}\mathbf{g}_{x}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\mathbf{g}_{y}\mathbf{f}_{B_{y}}^{H},$

where we denoted by $\mathbf{g}_{z}$ the generating function of $2\alpha
D^{-1}_{{A}_{z}}-\alpha^{2}D^{-1}_{{A_{z}}}{A_{z}}D^{-1}_{{A}_{z}}$, for
$z\in\\{x,y\\}$. If we choose $\alpha$ according to (20), then the functions
$\mathbf{g}_{z}$ are HPD at all points. Assumption 1. of Theorem 5 then
corresponds to requiring that there exists a $\boldsymbol{\theta}_{0}$ such
that $\mathbf{f}_{B_{x}}(\boldsymbol{\theta})$ and
$\mathbf{f}_{B_{y}}(\boldsymbol{\theta})$ are both the null vectors if and
only if $\boldsymbol{\theta}=\boldsymbol{\theta}_{0}$. Regarding the choice of
${p}_{\hat{C}}$, it is straightforward to see that in this scalar setting
relations (8)–(10) are fulfilled by the pair $({p}_{\hat{C}},{f}_{\hat{C}})$
if we have (8) and

$\limsup_{\boldsymbol{\theta}\rightarrow\boldsymbol{\theta}_{0}}\frac{\sum_{\boldsymbol{\xi}\in\Omega(\boldsymbol{\theta})\setminus\\{\boldsymbol{\theta}\\}}\left|{p}(\boldsymbol{\xi})\right|^{2}}{\sum_{z\in\\{x,y\\}}\mathbf{f}_{B_{z}}(\boldsymbol{\theta})\mathbf{g}_{z}(\boldsymbol{\theta})\mathbf{f}_{B_{z}}^{H}(\boldsymbol{\theta})}<\infty.$

Since the enumerator is a scalar, this second condition is the finiteness of
the following limit superior:

$\limsup_{\boldsymbol{\theta}\rightarrow\boldsymbol{\theta}_{0}}\frac{1}{\sum_{z\in\\{x,y\\}}\left(\mathbf{f}_{B_{z}}(\boldsymbol{\theta})/\sum_{\boldsymbol{\xi}\in\Omega(\boldsymbol{\theta})\setminus\\{\boldsymbol{\theta}\\}}\left|{p}(\boldsymbol{\xi})\right|\right)\mathbf{g}_{z}(\boldsymbol{\theta})\left(\mathbf{f}_{B_{z}}^{H}(\boldsymbol{\theta})/\sum_{\boldsymbol{\xi}\in\Omega(\boldsymbol{\theta})\setminus\\{\boldsymbol{\theta}\\}}\left|{p}(\boldsymbol{\xi})\right|\right)}.$

Recalling again that the matrices $\mathbf{g}_{z}(\boldsymbol{\theta})$ are
HPD for all $\boldsymbol{\theta}$, the latter condition corresponds to
requiring that

$\limsup_{\boldsymbol{\theta}\rightarrow\boldsymbol{\theta}_{0}}\frac{\mathbf{f}_{B_{x}}(\boldsymbol{\theta})}{\sum_{\boldsymbol{\xi}\in\Omega(\boldsymbol{\theta})\setminus\\{\boldsymbol{\theta}\\}}\left|{p}(\boldsymbol{\xi})\right|}\neq\mathbf{0}\qquad\mbox{or}\qquad\limsup_{\boldsymbol{\theta}\rightarrow\boldsymbol{\theta}_{0}}\frac{\mathbf{f}_{B_{y}}(\boldsymbol{\theta})}{\sum_{\boldsymbol{\xi}\in\Omega(\boldsymbol{\theta})\setminus\\{\boldsymbol{\theta}\\}}\left|{p}(\boldsymbol{\xi})\right|}\neq\mathbf{0}.$
(27)

######  Remark 6.

We expect that the convergence result in Theorem 5 holds also for multilevel
block Toeplitz matrices generated by trigonometric polynomials. Indeed, they
are a low rank correction of multilevel block circulant matrices with the same
generating function. The convergence could slightly deteriorate, but the
combination of the proposed multigrid method with Krylov methods kills the
outliers and guarantees the efficiency proven in the circulant case.

##  3 Saddle-point matrices stemming from the Stokes equation

In this section we are interested in applying the multigrid method described
before to large linear systems stemming from the finite element approximation
of the Stokes equation. The Stokes equation has the form

$\left\\{\begin{array}[]{rll}\Delta\mathbf{u}+\mathop{\mathrm{grad}}p&=-\Phi&\text{in}\,\Omega,\\\
\mathop{\mathrm{div}}\mathbf{u}&=0&\text{in}\,\Omega,\\\
\mathbf{u}&=\mathbf{g}_{D}&\text{on}\,\Gamma_{D},\\\
\frac{\partial\mathbf{u}}{\partial\nu}&=\mathbf{g}_{N}&\text{on}\,\Gamma_{N},\end{array}\right.$
(28)

where $\mathbf{u}$ is the velocity and $p$ is the pressure, by $\Gamma_{D}$ we
denote the Dirichlet boundary part of $\Gamma:=\partial\Omega$ and by
$\Gamma_{N}$ the Neumann part. For the construction and analysis of the system
matrices we consider the domain $\Omega=(0,1)^{2}$, homogenous Dirichlet
boundary conditions are attained for $x_{1}=0$ and $x_{2}=1$, homogenous
Neumann boundary conditions for $x_{1}=1$ and $x_{2}=0$.

###  3.1 A TGM for a Finite Element discretization of the Stokes problem

The weak form of (28) is given by

$\begin{array}[]{rll}a(\mathbf{u},\mathbf{v})-b(\mathbf{v},q)&=(\mathbf{\Phi},\mathbf{v}),&\text{for
all $\mathbf{v}$},\\\ b(\mathbf{u},q)&=0,&\text{for all
$\mathbf{q}$},\end{array}$

where

$a(\mathbf{u},\mathbf{v})=\int\limits_{\Omega}\nabla\mathbf{u}\cdot\nabla\mathbf{v}\
dx\quad\text{and}\quad
b(\mathbf{u},q)=\int\limits_{\Omega}\mathop{\mathrm{div}}\mathbf{u}\ q\ dx.$

Higher order elements are needed for the discretization of the velocities
$\mathbf{u}$ than for the discretization of the pressure $p$, e.g., Q2-Q1
Taylor-Hood-elements. Alternatively, linear elements can be used for the
displacement but using a subdivision of the element, corresponding to a
2-times finer—or once refined—mesh for the displacements. We choose the latter
approach that is often called Q1-iso-Q2/Q1 [9]. The properties of the system
matrix are similar to those of the system matrix produced by Q2-Q1 elements.
The advantage of using Q1-iso-Q2/Q1 is that the discrete operator representing
$a$ is more sparse than that for Q2-Q1 and thus cheaper to apply. The
resulting algebraic system of equations given by

${\mathcal{A}}_{N}\mathbf{x}=\left[\begin{array}[]{cc|c}{A_{x}}&O&B_{x}\\\
O&A_{y}&B_{y}\\\\[-1.00006pt] \hline\cr
B_{x}^{H}&B_{y}^{H}&O\end{array}\right]\begin{bmatrix}\mathbf{u}\\\
\mathbf{p}\end{bmatrix}=\begin{bmatrix}\mathbf{b_{u}}\\\
\mathbf{b_{p}}\end{bmatrix}=\mathbf{b},$ (29)

where $A_{x},A_{y}\in\mathbb{R}^{4n^{2}\times 4n^{2}}$ and
$B_{x},B_{y}\in\mathbb{R}^{4n^{2}\times n^{2}}$, so $N=9n^{2}$. In particular,
the matrices $A_{x}$ and $A_{y}$ are non-negative and they are a permutation
of $4\times 4$ block bi-level Toeplitz matrices with partial dimension $n$.
Precisely, let $\mathbf{e}_{i}$ be the $i$th column of the identity matrix of
size $4n^{2}$, we can define a proper $4n^{2}\times 4n^{2}$ permutation
matrix, $\Pi_{4}=[P_{1}|P_{2}|\dots|P_{4}]$, $P_{l}\in\mathbb{R}^{4n^{2}\times
n^{2}},\,l=1,\dots,4$, such that the $k$th column of $P_{l},\,l=1,\dots,4$, is
$e_{l+4(k-1)}$. The matrix $\Pi_{4}$ transforms the matrix ${A_{x}}$ into the
following a $4\times 4$ block bi-level Toeplitz matrix

$T_{\mathbf{n}}(\mathbf{f}_{A_{x}}(\theta_{1},\theta_{2}))=\Pi_{4}{A_{x}}\Pi_{4}^{T},$

where

$\mathbf{f}_{A_{x}}(\theta_{1},\theta_{2})=-\frac{1}{3}\begin{bmatrix}-8&1+{\rm
e}^{-\iota\theta_{1}}&1+{\rm
e}^{-\iota\theta_{2}}&\varphi(-\theta_{1},-\theta_{2})\\\ 1+{\rm
e}^{\iota\theta_{1}}&-8&\varphi(\theta_{1},-\theta_{2})&1+{\rm
e}^{-\iota\theta_{2}}\\\ 1+{\rm
e}^{\iota\theta_{2}}&\varphi(-\theta_{1},\theta_{2})&-8&1+{\rm
e}^{-\iota\theta_{1}}\\\ \varphi(\theta_{1},\theta_{2})&1+{\rm
e}^{\iota\theta_{2}}&1+{\rm e}^{\iota\theta_{1}}&-8\end{bmatrix},$ (30)

with $\varphi(\theta_{1},\theta_{2})=1+{\rm e}^{\iota\theta_{1}}+{\rm
e}^{\iota\theta_{2}}+{\rm e}^{\iota(\theta_{1}+\theta_{2})}$. An analogous
transformation can be applied to ${A_{y}}$ and the generating function is
$\mathbf{f}_{A_{y}}(\theta_{1},\theta_{2})=\mathbf{f}_{A_{x}}(\theta_{2},\theta_{1})$.

We consider an analogous $N\times N$ permutation matrix,
$\Pi_{9}=[P_{1}|P_{2}|\dots|P_{9}]$, $P_{l}\in\mathbb{R}^{N\times
n^{2}},\,l=1,\dots,9$, such that the $k$th column of $P_{l},\,l=1,\dots,9$, is
$e_{l+9(k-1)}$, with $e_{i}$ being the $i$th column of the identity matrix of
size $9n^{2}$. We can transform the global matrix $\mathcal{A}_{N}$ in a
$9\times 9$ block bi-level system. This property will be used later when we
will exploit the information of the permuted matrix
$G_{N}=\Pi_{9}\mathcal{A}_{N}\Pi_{9}^{T}=T_{\mathbf{n}}(\mathbf{f}),$ where
$T_{\mathbf{n}}(\mathbf{f})$ is the bi-level $9\times 9$ block Toeplitz
$T_{\bf{n}}(\mathbf{f})=\left[\hat{\mathbf{f}}_{\mathbf{i}-\mathbf{j}}\right]_{\mathbf{i},\mathbf{j}=\mathbf{1}}^{\bf
n}\in\mathbb{C}^{N\times N}$ generated by
$\mathbf{f}:[-\pi,\pi]^{2}\rightarrow\mathbb{C}^{9\times 9}$. In particular,

$\mathbf{f}(\theta_{1},\theta_{2})=\begin{bmatrix}\mathbf{f}_{A_{x}}(\theta_{1},\theta_{2})&O&\mathbf{f}_{B_{x}}(\theta_{1},\theta_{2})\\\
O&\mathbf{f}_{A_{y}}(\theta_{1},\theta_{2})&\mathbf{f}_{B_{y}}(\theta_{1},\theta_{2})\\\
\mathbf{f}_{B_{x}}^{H}(\theta_{1},\theta_{2})&\mathbf{f}_{B_{y}}^{H}(\theta_{1},\theta
2)&0\end{bmatrix},$ (31)

where

$\mathbf{f}_{B_{x}}(\theta_{1},\theta_{2})=\begin{bmatrix}\frac{1}{48}(-\iota\sin\theta_{1})(1+2\cos\theta_{2})\\\
\frac{1}{24}(1-{\rm e}^{\iota\theta_{1}})(5+\cos\theta_{2})\\\
\frac{1}{8}(\iota\sin\theta_{1})(1+{\rm e}^{\iota\theta_{2}})\\\
\frac{1}{8}(1-{\rm e}^{\iota\theta_{1}})(1+{\rm
e}^{\iota\theta_{2}})\end{bmatrix}\\!\\!,\,\mathbf{f}_{B_{y}}(\theta_{1},\theta_{2})=\begin{bmatrix}\frac{1}{48}(1+2\cos\theta_{1})(-\iota\sin\theta_{2})\\\
\frac{1}{8}(1+{\rm e}^{\iota\theta_{1}})(-\iota\sin\theta_{2})\\\
\frac{1}{24}(5+\cos\theta_{1})(1-{\rm e}^{\iota\theta_{2}})\\\
\frac{1}{8}(1+{\rm e}^{\iota\theta_{1}})(1-{\rm
e}^{\iota\theta_{2}})\end{bmatrix}\\!\\!.$ (32)

To solve efficiently the system (29) with TGM satisfying the assumptions of
Theorem 5, we investigate the structure of the coefficient matrix
${\mathcal{A}}_{N}$. It is a block matrix of the form

${\mathcal{A}}_{N}=\left[\begin{array}[]{ccc}\tilde{A}&B^{T}\\\
B&-C\end{array}\right],$

with $C=O$, $B=[B_{x},B_{y}]$ and $\tilde{A}=\begin{bmatrix}{A_{x}}&O\\\
O&{A_{y}}\end{bmatrix}$.

###  3.2 Spectral analysis of $\tilde{A}$ and choice of $P_{\tilde{A}}$

Lemma 1 and Theorem 5 highlight the importance of the block $\tilde{A}$ for
the choice of the parameter $\alpha$ and the construction of an efficient
operator $P_{\tilde{A}}$. In practice we can just consider the matrix $A_{x}$
and the related generating function $\mathbf{f}_{A_{x}}$, since $A_{y}$ and
$\mathbf{f}_{A_{y}}$ enjoy the same spectral properties. For this reason in
the following we will avoid the specification $x$ or $y$. We already know that
$A$ is similar to a bi-level $4\times 4$ block Toeplitz
$T_{\bf{n}}(\mathbf{f}_{A})$. The following proposition provides us with
important spectral information on $\mathbf{f}_{A}$, and consequently
$T_{\bf{n}}(\mathbf{f}_{A})$.

######  Proposition 1.

Let $\mathbf{f}_{A}$ be the $4\times 4$ matrix-valued function defined in (30)
and $T_{\bf{n}}(\mathbf{f}_{A})=\Pi_{4}A\Pi_{4}^{T}$ the associated Toeplitz
matrix. Then the four eigenvalue functions of $\mathbf{f}_{A}$ are

$\begin{split}&\lambda_{1}(\mathbf{f}_{A})=3-\frac{1}{3}{\rm
e}^{-\iota\left(\frac{\theta_{1}+\theta_{2}}{2}\right)}\left({\rm
e}^{\iota\theta_{1}}+{\rm e}^{\iota\frac{\theta_{1}}{2}}+1\right)\left({\rm
e}^{\iota\theta_{2}}+{\rm e}^{\iota\frac{\theta_{2}}{2}}+1\right);\\\
&\lambda_{2}(\mathbf{f}_{A})=3-\frac{1}{3}{\rm
e}^{-\iota\left(\frac{\theta_{1}+\theta_{2}}{2}\right)}\left({\rm
e}^{\iota\theta_{1}}-{\rm e}^{\iota\frac{\theta_{1}}{2}}+1\right)\left({\rm
e}^{\iota\theta_{2}}-{\rm e}^{\iota\frac{\theta_{2}}{2}}+1\right);\\\
&\lambda_{3}(\mathbf{f}_{A})=3+\frac{1}{3}{\rm
e}^{-\iota\left(\frac{\theta_{1}+\theta_{2}}{2}\right)}\left({\rm
e}^{\iota\theta_{2}}-{\rm e}^{\iota\frac{\theta_{2}}{2}}+1\right)\left({\rm
e}^{\iota\theta_{1}}+{\rm e}^{\iota\frac{\theta_{1}}{2}}+1\right);\\\
&\lambda_{4}(\mathbf{f}_{A})=3+\frac{1}{3}{\rm
e}^{-\iota\left(\frac{\theta_{1}+\theta_{2}}{2}\right)}\left({\rm
e}^{\iota\theta_{1}}-{\rm e}^{\iota\frac{\theta_{1}}{2}}+1\right)\left({\rm
e}^{\iota\theta_{2}}+{\rm e}^{\iota\frac{\theta_{2}}{2}}+1\right).\\\
\end{split}$ (33)

Moreover,

  1. 1.

the minimum eigenvalue function of $\lambda_{1}(\mathbf{f}_{A})$ is non-
negative with a zero of order 2 in $\boldsymbol{\theta}_{0}=(0,0)$.
Furthermore, $e_{4}$ is the eigenvector of
$\mathbf{f}_{A}(\boldsymbol{\theta}_{0})$ associated with
$\lambda_{1}(\mathbf{f}_{A}(\boldsymbol{\theta}_{0}))=0$, where ${\rm
e}_{4}=[1,1,1,1]^{T}$.

  2. 2.

the minimal eigenvalue of $T_{\bf{n}}(\mathbf{f}_{A})$ (and consequently of
$A$) goes to zero as $(4n^{2})^{-1}$;

  3. 3.

the condition number $\kappa(T_{\bf{n}}(\mathbf{f}_{A}))$ of
$T_{\bf{n}}(\mathbf{f}_{A})$ (and consequently of $A$) is proportional to
$O(n^{2})$.

###### Proof.

The function $\mathbf{f}_{A}$ can be written in a more compact form as

$\mathbf{f}_{A}(\theta_{1},\theta_{2})=\frac{9}{3}I_{4}-\frac{1}{3}h(\theta_{2})\otimes
h(\theta_{1})=3I_{4}-\frac{1}{3}h(\theta_{2})\otimes h(\theta_{1}),$ (34)

where $h(\theta)=\begin{bmatrix}1&{\rm e}^{-\iota\theta}+1\\\ {\rm
e}^{\iota\theta}+1&1\end{bmatrix}.$ This implies that
$T_{\bf{n}}(\mathbf{f}_{A})$ is the SPD matrix

$T_{\bf{n}}(\mathbf{f}_{A})=3I_{4n^{2}}-\frac{1}{3}T_{{n}}(h(\theta_{2}))\otimes
T_{{n}}(h(\theta_{1})).$ (35)

Moreover, formula (34) and the properties of the tensor product imply that the
four eigenvalue functions $\lambda_{k}(\mathbf{f}_{A}),\,k=1,\dots,4$ of
$\mathbf{f}_{A}$ are given by the following combination of the two eigenvalue
functions of $h({\theta})$

$\left\\{3-\frac{1}{3}\lambda_{i}(h(\theta_{2}))\lambda_{j}(h(\theta_{1})\right\\}_{i,j=1}^{2}.$

The eigenvalue functions of $h({\theta})$ can be computed analytically and
they are given by

$\lambda_{1}(h(\theta))=-{\rm e}^{-\iota\frac{\theta}{2}}({\rm
e}^{\iota{\theta}}-{\rm
e}^{\iota\frac{\theta}{2}}+1);\quad\lambda_{2}(h(\theta))={\rm
e}^{-\iota\frac{\theta}{2}}({\rm e}^{\iota{\theta}}+{\rm
e}^{\iota\frac{\theta}{2}}+1).$

Then, using their combination, we derive the expressions of the four
eigenvalue functions of $\mathbf{f}_{A}(\theta_{1},\theta_{2})$ which are
exactly those given in formula (33). Moreover, by direct computation, we have

$\mathbf{f}_{A}(0,0)e_{4}=\begin{bmatrix}&8/3&-2/3&-2/3&-4/3\\\
&-2/3&8/3&-4/3&-2/3\\\ &-2/3&-4/3&8/3&-2/3\\\ &-4/3&-2/3&-2/3&8/3\\\
\end{bmatrix}e_{4}=0e_{4}.$

Then, $e_{4}$ is an eigenvector of $\mathbf{f}_{A}(\boldsymbol{\theta}_{0})$
associated with $0$. In addition, the minimal eigenvalue function is

$\lambda_{1}(\mathbf{f}_{A})=3-\frac{1}{3}{\rm
e}^{-\iota\left(\frac{\theta_{1}+\theta_{2}}{2}\right)}\left({\rm
e}^{\iota\theta_{1}}+{\rm e}^{\iota\frac{\theta_{1}}{2}}+1\right)\left({\rm
e}^{\iota\theta_{2}}+{\rm e}^{\iota\frac{\theta_{2}}{2}}+1\right),$

and it is straightforward to see that is a non-negative function over
$[-\pi,\pi]$ and it is such that

$\displaystyle\lambda_{1}(\mathbf{f}_{A}){|_{(0,0)}}$ $\displaystyle=0,$
$\displaystyle\frac{\partial\,\lambda_{1}(\mathbf{f}_{A})(\theta_{1},\theta_{2})}{\partial{\theta_{1}}}{\left|{}_{(0,0)}\right.}$
$\displaystyle=\frac{\partial\,\lambda_{1}(\mathbf{f}_{A})(\theta_{1},\theta_{2})}{\partial{\theta_{2}}}{\left|{}_{(0,0)}\right.}=0,$
$\displaystyle\frac{\partial^{2}\,\lambda_{1}(\mathbf{f}_{A})(\theta_{1},\theta_{2})}{\partial{\theta_{2}}\partial{\theta_{1}}}{\left|{}_{(0,0)}\right.}$
$\displaystyle=\frac{\partial^{2}\,\lambda_{1}(\mathbf{f}_{A})(\theta_{1},\theta_{2})}{\partial{\theta_{1}}\partial{\theta_{2}}}{\left|{}_{(0,0)}\right.}=0,$
$\displaystyle\frac{\partial^{2}\,\lambda_{1}(\mathbf{f}_{A})(\theta_{1},\theta_{2})}{\partial{\theta_{1}^{2}}}{\left|{}_{(0,0)}\right.}$
$\displaystyle=\frac{\partial^{2}\,\lambda_{1}(\mathbf{f}_{A})(\theta_{1},\theta_{2})}{\partial{\theta_{2}^{2}}}{\left|{}_{(0,0)}\right.}=\frac{1}{2}.$

Therefore, $\lambda_{1}(\mathbf{f}_{A})$ has a zero of order 2 in
$\boldsymbol{\theta}_{0}=(0,0)$. In light of the third item of Theorem 2, we
conclude that the minimal eigenvalue of $T_{\bf{n}}(\mathbf{f}_{A})$ goes to
zero as $(4n^{2})^{-1}$ and, by similitude, that of $A$ as well. Furthermore,

$M_{s}=\max_{(\theta_{1},\theta_{2})\in[-\pi,\pi]}\lambda_{4}(\mathbf{f}_{A})=\max_{(\theta_{1},\theta_{2})\in[-\pi,\pi]}\lambda_{3}(\mathbf{f}_{A})=4,$
(36)

and from Theorem 2 we have that the spectrum of $T_{\bf{n}}(\mathbf{f}_{A})$
is contained in $(0,4).$ Consequently, the condition number
$\kappa(T_{\bf{n}}(\mathbf{f}_{A}))$ of $T_{\bf{n}}(\mathbf{f}_{A})$ is

$\kappa(T_{\bf{n}}(\mathbf{f}_{A}))=\frac{\lambda_{\max}(T_{\bf{n}}(\mathbf{f}_{A}))}{\lambda_{\min}(T_{\bf{n}}(\mathbf{f}_{A}))}\leq\frac{4}{\lambda_{\min}(T_{\bf{n}}(\mathbf{f}_{A}))}\thickapprox
n^{2}.$

The similitude between $A$ and $T_{\bf{n}}(\mathbf{f}_{A})$ concludes the
proof of the Proposition. ∎

A graphical confirmation of the result of the Proposition 1 concerning the
behavior of the eigenvalue functions is given in Figure 1 where the plot of
$\lambda_{i}(\mathbf{f}_{A})$, $i=1,\dots,4,$ is shown.

Figure 1: Plot of the eigenvalues functions $\lambda_{l}(\mathbf{f}_{A})$,
$l=1,\ldots,4$ when $n=10$.

######  Remark 7.

Because of the spectral properties of the symbol $\mathbf{f}_{A}$ and
$T_{\bf{n}}(\mathbf{f}_{A})$, we need to choose

$\alpha<2\max_{j=1,\dots,4}\left\\{\left(\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A})\right)^{-1}\|\lambda_{j}\left(\mathbf{f}_{A}\right)\|_{\infty}\right\\}^{-1}=2\left(\,\frac{3}{8}\,4\right)^{-1}=\frac{4}{3}$

Then, in the following we will fix $\alpha$ equal to $\frac{2}{3}$ since is
the middle value in the range $(0,4/3].$

As suggested by Theorem 5, we need to construct

${P}_{\tilde{A}}=\begin{bmatrix}\mathcal{C}_{\mathbf{n}}(\mathbf{p}_{A_{x}})(K_{\mathbf{n}}^{T}\otimes
I_{s_{a}})&\\\
&\mathcal{C}_{\mathbf{n}}(\mathbf{p}_{A_{y}})(K_{\mathbf{n}}^{T}\otimes
I_{s_{a}})\end{bmatrix},\quad
P_{\hat{C}}=\mathcal{C}_{\mathbf{n}}(\mathbf{p}_{\hat{C}})(K_{\mathbf{n}}^{T}\otimes
I_{s_{c}})$ (37)

with $\mathbf{p}_{A_{x}}$, $\mathbf{p}_{A_{y}}$, and $\mathbf{p}_{\hat{C}}$
proper bi-variate matrix-valued trigonometric polynomials. Since
$\mathbf{f}_{A_{x}}$ and $\mathbf{f}_{A_{y}}$ have the same spectral
properties, we fix $\mathbf{p}_{A_{x}}=\mathbf{p}_{A_{y}}=:\mathbf{p}_{A}.$ In
the following we discuss the choice of the polynomial $\mathbf{p}_{A}$.

We recall that a classical approach is that of constructing $P_{{A}}$ as the
multilevel linear interpolation projector. That is,

${P}_{{A}}=\Pi_{4}^{T}\left(T_{\mathbf{n}}(\mathbf{p}^{(4)})\left(K_{\mathbf{n,k}}\otimes
I_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}4}}\right)\right)\Pi_{4},$
(38)

where $K_{\mathbf{n,k}}=K_{{n,k}}\otimes K_{{n,k}}$, $k=\frac{n-1}{2}$, $n$
odd, and

$K_{n,k}=\left[\begin{array}[]{cccccccc}0&&&&\\\ 1&&&&\\\ 0&0&&&\\\
\vdots&1&&&\\\ &0&&&\\\ &\vdots&&&\\\ &&&&\vdots\\\ &&&&1\\\ &&&&0\\\
\end{array}\right]_{n\times k}.$

Moreover, we choose as
$\mathbf{p}^{(4)}(\theta_{1},\theta_{2})=\mathbf{p}^{(2)}(\theta_{1})\otimes\mathbf{p}^{(2)}(\theta_{2})$,
that is the Kronecker product of the univariate trigonometric polynomial

$\mathbf{p}^{(2)}(\theta)=\begin{bmatrix}1+{\rm e}^{-\iota\theta}&{\rm
e}^{\iota\theta}+1\\\ 2{\rm e}^{-\iota\theta}&2\end{bmatrix}.$ (39)

Next theorem highlights the properties of the polynomial
$\mathbf{p}^{(2)}(\theta)$ and its multivariate version
$\mathbf{p}^{(4)}(\theta_{1},\theta_{2})$ in order to show that the grid
transfer operator ${P}_{{A}}$ fulfills the conditions for the convergence and
optimality of the TGM methods proposed in [8, Section 6], here conditions (8)
\- (10).

######  Theorem 6.

Let $\mathbf{p}^{(2)}(\theta)$ be the $2\times 2$ trigonometric polynomial
defined in (39) and
$\mathbf{p}^{(4)}(\theta_{1},\theta_{2})=\mathbf{p}^{(2)}(\theta_{1})\otimes\mathbf{p}^{(2)}(\theta_{2})$.
Consider $\mathbf{f}_{A}(\theta_{1},\theta_{2})$ defined in (30). Then,
$\mathbf{p}^{(4)}(\theta_{1},\theta_{2})=\mathbf{p}^{(2)}(\theta_{1})\otimes\mathbf{p}^{(2)}(\theta_{2})$
is such that:

  * 1.

$\sum_{\xi\in\Omega(\boldsymbol{\theta})}\mathbf{p}^{(4)}(\xi)^{H}\mathbf{p}^{(4)}(\xi)>0,\quad\forall\,\boldsymbol{\theta}\in[0,2\pi)^{{2}},$
(40)

where
$\Omega(\theta)=\left\\{({\theta_{1}},\theta_{2}),({\theta_{1}+\pi},\theta_{2}),({\theta_{1}},\theta_{2}+\pi),({\theta_{1}}+\pi,\theta_{2}+\pi)\right\\}$.

  * 2.

The function
$\mathbf{s}^{(4)}(\boldsymbol{\theta})=\mathbf{p}^{(4)}(\boldsymbol{\theta})\left(\sum_{\xi\in\Omega(\boldsymbol{\theta})}\mathbf{p}^{(4)}(\xi)^{H}\mathbf{p}^{(4)}(\xi)\right)^{-1}\mathbf{p}^{(4)}(\boldsymbol{\theta})^{H}$
is well-defined for all $\boldsymbol{\theta}\in[0,2\pi)^{{2}}$ and

$\mathbf{s}^{(4)}(\boldsymbol{\theta}_{0})e_{4}=e_{4},$ (41)

where $e_{4}=[1,1,1,1]^{T}$ and $\boldsymbol{\theta}_{0}=(0,0)$.

  * 3.

$\lim_{\boldsymbol{\theta}\rightarrow\boldsymbol{\theta}_{0}}(\lambda_{1}(\mathbf{f}_{A}(\boldsymbol{\theta})))^{-1}(1-\lambda_{1}(\mathbf{s}^{(4)}(\boldsymbol{\theta})))=c,$
(42)

where $c\in\mathbb{R}$ is a constant and $\lambda_{1}(\mathbf{f}_{A})$ and
$\lambda_{1}(\mathbf{s}^{(4)})$ are the minimal eigenvalue functions of
$\mathbf{f}_{A}$ and $\mathbf{s}^{(4)}$, respectively.

###### Proof.

Given $e_{2}=[1,1]^{T}$, the trigonometric polynomial $\mathbf{p}^{(2)}$ by
direct computation fulfils the following:

  1. 1)

$\mathbf{p}^{(2)}(0)\,{\rm e}_{2}=4\,{\rm e}_{2}$.

  2. 2)

$\mathbf{p}^{(2)}(\pi)\,{\rm e}_{2}=0\,{\rm e}_{2}$.

  3. 3)

$\mathbf{p}^{(2)}(0)^{H}\,{\rm e}_{2}=4\,{\rm e}_{2}$.

  4. 4)

$\mathbf{p}^{(2)}(\theta)^{H}\mathbf{p}^{(2)}(\theta)+\mathbf{p}^{(2)}(\theta+\pi)^{H}\mathbf{p}^{(2)}(\theta+\pi)=\begin{bmatrix}12&2{\rm
e}^{2\iota\theta}+2\\\ 2{\rm e}^{-2\iota\theta}+2&12\end{bmatrix}$.

Item $4)$ clearly implies that
$\mathbf{p}^{(2)}(\theta)^{H}\mathbf{p}^{(2)}(\theta)+\mathbf{p}^{(2)}(\theta+\pi)^{H}\mathbf{p}^{(2)}(\theta+\pi)$
is a positive definite matrix, $\forall\theta\in[0,2\pi)$. The items $1)-3)$,
together with [8, Lemma 4.3], imply that

$\mathbf{s}^{(2)}(0)e_{2}=e_{2},$

where
$\mathbf{s}^{(2)}=\mathbf{p}^{(2)}(\theta)\left(\mathbf{p}^{(2)}(\theta)^{H}\mathbf{p}^{(2)}(\theta)+\mathbf{p}^{(2)}(\theta+\pi)^{H}\mathbf{p}^{(2)}(\theta+\pi)\right)^{-1}\mathbf{p}^{(2)}(\theta)^{H}$,
is well-defined for all $\theta\in[0,2\pi)$. Because of the tensor structure
of $\mathbf{p}^{(4)}$, the result in [8, Lemma 6.2] ensures that it verifies
the analogous positivity condition given in (40). Consequently, the quantity

$\mathbf{s}^{(4)}(\boldsymbol{\theta})=\mathbf{p}^{(4)}(\boldsymbol{\theta})\left(\sum_{\xi\in\Omega(\boldsymbol{\theta})}\mathbf{p}^{(4)}(\xi)^{H}\mathbf{p}^{(4)}(\xi)\right)^{-1}\mathbf{p}^{(4)}(\boldsymbol{\theta})^{H}$

is well-defined for all
$\boldsymbol{\theta}=(\theta_{1},\theta_{2})\in[0,2\pi)^{{2}}$. In addition,
by [8, Lemma 6.2], we have
$\mathbf{s}^{(4)}(\theta_{1},\theta_{2})=\mathbf{s}^{(2)}(\theta_{1})\otimes\mathbf{s}^{(2)}(\theta_{2})$.
Writing $e_{4}$ as $e_{2}\otimes e_{2}$ and exploiting the properties of the
Kronecker product, we have

$\begin{split}&\mathbf{s}^{(4)}(\boldsymbol{\theta}_{0})e_{4}=(\mathbf{s}^{(2)}(0)\otimes\mathbf{s}^{(2)}(0))(e_{2}\otimes
e_{2})=\\\
&(\mathbf{s}^{(2)}(0)e_{2})\otimes(\mathbf{s}^{(2)}(0)e_{2})=e_{2}\otimes
e_{2}=e_{4}.\end{split}$

which concludes the proof of (41). The last thing to be proved is that

$\limsup_{\boldsymbol{\theta}\rightarrow\boldsymbol{\theta}_{0}}\lambda_{1}(\mathbf{f}_{A}(\boldsymbol{\theta}))^{-1}(1-\lambda_{1}(\mathbf{s}^{(4)}(\boldsymbol{\theta})))<+\infty.$
(43)

Note that $\mathbf{s}^{(2)}$ is an algebraic projector, that is, it can be
verified by direct computation that
$\mathbf{s}^{2}(\theta)-\mathbf{s}(\theta)=\mathbf{0}$. Since
$\lambda_{1}(\mathbf{s}(0))=1$ and the eigenvalue functions are continuous, we
have $\lambda_{1}(\mathbf{s}^{(2)}(\theta))\equiv 1$ and so
$\lambda_{1}(\mathbf{s}^{(4)})$ is identically 1 as well. Moreover, by (33),
the function $\lambda_{1}(\mathbf{f}_{A})$ is different from $0$ in a
neighbourhood of $\boldsymbol{\theta}_{0}$. Then, the bound in (42) holds. ∎

###  3.3 Analysis of $\hat{C}$ and choice of $P_{\hat{C}}$

In order to compute the sequence of grid transfer operators $P_{\hat{C}}$, we
study the structure of $\hat{C}$, which in the presented setting has a (bi-
level) scalar nature.

Then, the associated spectral symbol will be a (bi-variate) scalar-valued
function instead of matrix-valued. Multigrid methods for scalar non-negative
Toeplitz matrix sequences have been deeply investigated in [1, 2, 15, 16, 26],
so we will first draft the idea on which our choice of ${P}_{\hat{C}}$ is
based to see that it fits into a classical setting, referring in particular to
Lemma 4.3 in [26]. Then, we will exploit Remark 5 to theoretically validate
our choice.

Firstly, we recall that $\hat{C}$ is a non-negative matrix and
$\hat{C}=B(2\alpha
D_{\tilde{A}}^{-1}-\alpha^{2}D_{\tilde{A}}^{-1}\tilde{A}D_{\tilde{A}}^{-1})B^{T}$.
Because of the structure of $A$, $D_{\tilde{A}}=\frac{8}{3}I_{8n^{2}}$.
Moreover $\alpha=\frac{2}{3}$. Then, $\hat{C}$ is a low-rank correction of the
Toeplitz matrix $T_{\mathbf{n}}({f}_{\hat{C}})$ with

${f}_{\hat{C}}(\theta_{1},\theta_{2})=\left(\mathbf{f}_{B_{x}}\left(\frac{1}{2}I_{4}-\frac{1}{16}\mathbf{f}_{A_{x}}\right)\mathbf{f}^{H}_{B_{x}}+\mathbf{f}_{B_{y}}\left(\frac{1}{2}I_{4}-\frac{1}{16}\mathbf{f}_{A_{y}}\right)\mathbf{f}^{H}_{B_{y}}\right)(\theta_{1},\theta_{2}),$

see Remark 3 for more details on the construction of the function
${f}_{\hat{C}}$.

The function ${f}_{\hat{C}}$ is a bivariate scalar valued function such that
$\hat{a}_{0}({f}_{\hat{C}})=11/96$. Moreover, it is straightforward to show
that ${f}_{\hat{C}}$ has a zero of order 2 in the origin. For example we can
use the procedure analogous to the proof of Proposition 1 with the additional
simplification that ${f}_{\hat{C}}$ is scalar valued then its eigenvalue
function coincides with ${f}_{\hat{C}}$ itself. Figure 2 shows the plot of
${f}_{\hat{C}}$ on the grid $[0,\pi]^{2}$.

Figure 2: Plot of the function ${f}_{\hat{C}}$ on the grid $[0,\pi]^{2}$ with
$n=20$ points.

Then, we construct ${P}_{{C}}$ as

${P}_{{C}}=T_{\mathbf{n}}\left((2+2\cos\theta_{1})(2+2\cos\theta_{2})\right)K_{\mathbf{n,k}}.$
(44)

In order to be sure that with this choice of grid transfer operator assumption
2. of Theorem 5 is fulfilled, we recall Remark 5. Positivity condition (8) is
satisfied since the sum

$\displaystyle\left|(2+2\cos(\theta_{1}))(2+2\cos(\theta_{2}))\right|^{2}+\left|(2+2\cos(\theta_{1}+\pi))(2+2\cos(\theta_{2}))\right|^{2}+$
$\displaystyle+$
$\displaystyle\left|(2+2\cos(\theta_{1}))(2+2\cos(\theta_{2}+\pi))\right|^{2}+\left|(2+2\cos(\theta_{1}+\pi))(2+2\cos(\theta_{2}+\pi))\right|^{2}$

is different from 0 for all $(\theta_{1},\theta_{2})\in\mathbb{R}^{2}$.
Condition (27) is fulfilled since the functions in (32) vanish in $(0,0)$ with
a zero of order 1.

###  3.4 Choice of the smoothing parameter

In order to fulfill the assumption 3. of Theorem 5 in our setting we require
that

$\omega<2\min\left\\{2\alpha-\alpha^{2}\max_{j}\left\\{\frac{\|\mathbf{f}_{A_{x}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{x}})}},\frac{\|\mathbf{f}_{A_{y}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{y}})}}\right\\},\right.\\\
\left.\hat{a}_{0}({f}_{\hat{C}})\left\|\mathbf{f}_{C}+\mathbf{f}_{B_{x}}\mathbf{f}_{A_{x}}^{-1}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\mathbf{f}_{A_{y}}^{-1}\mathbf{f}_{B_{y}}^{H}\right\|_{\infty}^{-1}\right\\}.$

That is

$\displaystyle\omega$
$\displaystyle<2\min\left\\{2\left(\frac{2}{3}\right)-\left(\frac{2}{3}\right)^{2}\frac{4\cdot
3}{8},\right.\left.\frac{11}{96}\left\|\mathbf{f}_{B_{x}}\mathbf{f}_{A_{x}}^{-1}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\mathbf{f}_{A_{y}}^{-1}\mathbf{f}_{B_{y}}^{H}\right\|_{\infty}^{-1}\right\\}$
(45)
$\displaystyle=2\min\left\\{\frac{2}{3},\frac{11}{24}\right\\}=\frac{11}{12}.$
(46)

In the first equality above we are using the fact that the function
$f_{S}=\mathbf{f}_{B_{x}}\mathbf{f}_{A_{x}}^{-1}\mathbf{f}_{B_{x}}^{H}+\mathbf{f}_{B_{y}}\mathbf{f}_{A_{y}}^{-1}\mathbf{f}_{B_{y}}^{H}$
is a scalar valued function whose maximum is $\frac{1}{4}$ corresponding to
the point $(0,\pi)$ (equiv. $(\pi,0)$). Figure 3 shows the plot of ${f}_{S}$
on a uniform grid over $[0,\pi]^{2}$.

Figure 3: Plot of the function ${f}_{S}$ on a uniform grid over $[0,\pi]^{2}$
with $n=20$ points.

Therefore, Theorem 5 ensures that the TGM applied to the system having
$\mathcal{\hat{A}}$ as coefficient matrix converges with the choice
$\omega\in(0,11/12)$.

##  4 Extension to MGM

In Section 2 we presented a TGM procedure for solving a system of the form
(12) consisting in the transformation (21) and the application of the
projecting strategy given by $\mathcal{P}$. In this section we propose a
multigrid strategy with more than two grids. Moreover, in order to analyse the
multigrid method for multiple grids, we need to study the problem at the
coarser levels.

###  4.1 A 3-step procedure

Firstly, we define the CGC introducing a procedure consisting of 3 steps
(3SP):

  * (i)

transform the system into a new one by making use of the two invertible
matrices $\mathcal{L}$ and $\mathcal{U}$.

  * (ii)

Apply the projecting strategy given by $\mathcal{P}$.

  * (iii)

“Symmetrize” the problem with a proper variable transformation.

First of all, note that after step one the transformed system is not
symmetric. However, the block $\hat{C}=B(2\alpha
D_{\tilde{A}}^{-1}-\alpha^{2}D_{\tilde{A}}^{-1}\tilde{A}D_{\tilde{A}}^{-1})B^{T}$
is a SPD matrix as long as the parameter $\alpha$ is such that
$\alpha<2(\lambda_{\rm max}(D_{\tilde{A}}^{-1}A))^{-1}$. Equivalently, the
matrix in block position $(1,1)$ is SPD, since it is given by $\tilde{A}$.

After the application of the grid transfer operator

$\mathcal{P}=\begin{bmatrix}P_{\tilde{A}}&\\\ &P_{\hat{C}}\end{bmatrix},$

the matrix at the coarse level has the form

$\begin{bmatrix}P_{\tilde{A}}^{T}\tilde{A}P_{\tilde{A}}&P_{\tilde{A}}^{T}(I_{(d-1)n^{2}}-\alpha\tilde{A}D_{\tilde{A}}^{-1})B^{T}P_{\hat{C}}\\\
-P_{\hat{C}}^{T}B(I_{(d-1)n^{2}}-\alpha
D_{\tilde{A}}^{-1}\tilde{A})P_{\tilde{A}}&P_{\hat{C}}^{T}CP_{\hat{C}}+P_{\hat{C}}^{T}B(2\alpha
D_{\tilde{A}}^{-1}-\alpha^{2}D_{\tilde{A}}^{-1}\tilde{A}D_{\tilde{A}}^{-1})B^{T}P_{\hat{C}}\end{bmatrix},$
(47)

which is a matrix of a form which resembles that of the system (12), apart
from the signs of the blocks in position $(2,1)$ and $(2,2)$. This explains
why the third step of our procedure consists in the variable transformation
obtained multiplying both (47) and the right hand side on the left by the
matrix

$\begin{bmatrix}I_{2s_{a}n^{2}}&\\\ &-I_{s_{c}n^{2}}\end{bmatrix}.$ (48)

This last step restores the symmetry of the matrix and gives to (47) the same
block form of the original coefficient matrix so we can proceed recursively,
defining a V-cycle procedure.

Precisely, the base of the recursion is the original system written as

$\mathcal{A}_{N_{0}}\\{0\\}=\left[\begin{array}[]{ccc}\tilde{A}\\{0\\}&B\\{0\\}^{T}\\\
B\\{0\\}&-C\\{0\\}\end{array}\right],$ (49)

with $\tilde{A}\\{0\\}\in\mathbb{R}^{s_{a}n_{0}^{2}\times s_{a}n_{0}^{2}}$
SPD, $C\\{0\\}\in\mathbb{R}^{s_{c}n_{0}^{2}\times s_{c}n_{0}^{2}}$ non-
negative definite, and $B\\{0\\}\in\mathbb{R}^{s_{c}n_{0}^{2}\times
2s_{a}n_{0}^{2}}$.

At level $\ell$, we have
$P_{\tilde{A}\\{\ell\\}}\in\mathbb{R}^{2s_{a}n_{\ell}^{2}\times
2s_{a}n_{\ell+1}^{2}}$ and
$P_{\hat{C}\\{\ell\\}}\in\mathbb{R}^{s_{c}n_{\ell}^{2}\times
s_{c}n_{\ell+1}^{2}}$ the prolongation operators chosen for solving
efficiently the scalar systems with coefficient matrix $\tilde{A}\\{\ell\\}$
and
$\hat{C}\\{\ell\\}=C\\{\ell\\}+B\\{\ell\\}(2\alpha_{\ell}D_{\tilde{A}\\{\ell\\}}^{-1}-\alpha_{\ell}^{2}D_{\tilde{A}\\{\ell\\}}^{-1}\tilde{A}\\{\ell\\}D_{\tilde{A}\\{\ell\\}}^{-1})B\\{\ell\\}^{T}$,
respectively, where $\alpha_{\ell}$ is such that
$\alpha_{\ell}<2\lambda_{\max}\left(D_{\tilde{A}\\{\ell\\}}^{-1}\tilde{A}\right)^{-1}$.
Then, the inductive step of our recursive definition is given by

$\mathcal{A}_{N_{\ell+1}}\\{\ell+1\\}=\left[\begin{array}[]{ccc}\tilde{A}\\{\ell+1\\}&B\\{\ell+1\\}^{T}\\\
B\\{\ell+1\\}&-\hat{C}\\{\ell+1\\}\end{array}\right].$ (50)

where

$\tilde{A}\\{\ell+1\\}=P_{\tilde{A}\\{\ell\\}}^{T}\tilde{A}\\{\ell\\}P_{\tilde{A}\\{\ell\\}},$
(51)
$B\\{\ell+1\\}=P_{\hat{C}\\{\ell\\}}^{T}B\\{\ell\\}(I_{2s_{a}n_{\ell}^{2}}-\alpha_{\ell}D_{\tilde{A}\\{\ell\\}}^{-1}\tilde{A}\\{\ell\\})P_{\tilde{A}\\{\ell\\}},$
(52)
$\hat{C}\\{\ell+1\\}=\left(P_{\hat{C}\\{\ell\\}}^{T}\hat{C}\\{\ell\\}P_{\hat{C}\\{\ell\\}}+P_{\hat{C}}^{T}B\\{\ell\\}(2\alpha_{\ell}D_{\tilde{A}\\{\ell\\}}^{-1}-\alpha_{\ell}^{2}D_{\tilde{A}\\{\ell\\}}^{-1}\tilde{A}\\{\ell\\}D_{\tilde{A}\\{\ell\\}}^{-1})B\\{\ell\\}^{T}P_{\hat{C}\\{\ell\\}}\right).$
(53)

For the considerations that we made when we described our 3SP, for each $\ell$
we have that $\tilde{A}\\{\ell+1\\}\in\mathbb{R}^{2s_{a}n_{\ell+1}^{2}\times
2s_{a}n_{\ell+1}^{2}}$ is SPD,
$C\\{{\ell+1}\\}\in\mathbb{R}^{s_{c}n_{\ell+1}^{2}\times s_{c}n_{\ell+1}^{2}}$
is non-negative definite, and
$B\\{{\ell+1}\\}\in\mathbb{R}^{s_{c}n_{\ell+1}^{2}\times 2s_{a}n_{\ell}^{2}}$.

This procedure can be seen globally as one block multigrid algorithm with
restriction and prolongation operators $R\\{\ell\\}$ and $P\\{\ell\\}$.

######  Proposition 2.

Consider the three steps procedure described by items (i)-(ii). At level
$\ell$, define the prolongation and restriction operators as

$R\\{\ell\\}=\begin{bmatrix}I_{2s_{a}n^{2}\\{\ell\\}}&\\\
\alpha\\{\ell\\}P_{\hat{C}\\{\ell\\}}^{T}B\\{\ell-1\\}D_{\tilde{A}\\{\ell\\}}^{-1}P_{\tilde{A}\\{\ell\\}}&-I_{s_{c}n^{2}\\{\ell\\}}\end{bmatrix}\cdot\\\
\cdot\begin{bmatrix}I_{2s_{a}n^{2}\\{\ell\\}}&\\\
&-I_{s_{c}n^{2}\\{\ell\\}}&\\\
\end{bmatrix}\begin{bmatrix}P^{T}_{\tilde{A}\\{\ell\\}}&\\\
&P^{T}_{\tilde{A}\\{\ell\\}}\end{bmatrix}$ (54)
$P\\{\ell\\}=\begin{bmatrix}P^{T}_{\tilde{A}\\{\ell\\}}&\\\
&P^{T}_{\hat{C}\\{\ell\\}}\end{bmatrix}\begin{bmatrix}I_{2s_{a}n^{2}\\{\ell\\}}&-\alpha\\{\ell\\}P_{\tilde{A}\\{\ell\\}}^{T}D_{\tilde{A}\\{\ell\\}}^{-1}B^{T}\\{\ell-1\\}P_{\tilde{A}\\{\ell\\}}\\\
&I_{s_{c}n^{2}\\{\ell\\}}\end{bmatrix}$ (55)

After $\ell+1$ iterations of 3SP, the resulting coefficient matrix
$\mathcal{A}_{N_{\ell+1}}\\{\ell+1\\}$ is given by

$\mathcal{A}_{N_{\ell+1}}\\{\ell+1\\}=R\\{\ell\\}\begin{bmatrix}\tilde{A}\\{\ell\\}&B^{T}\\{\ell\\}\\\
-B\\{\ell\\}&\hat{C}\\{\ell\\}\end{bmatrix}P\\{\ell\\}.$ (56)

###### Proof.

The assertion follows immediately by direct computation. ∎

###  4.2 Smoother and grid transfer operator at coarser levels

In order to develop a multigrid method as described in the previous
subsection, we need to choose a sequence of smoothing parameters
$\omega_{\ell}$ and to construct two sequences of grid transfer operators
whose elements are $P_{\tilde{A}\\{\ell\\}}$ and $P_{\hat{C}\\{\ell\\}}$,
associated respectively with ${\tilde{A}\\{\ell\\}}$ and
${\hat{C}\\{\ell\\}}$.

In the circulant setting, the matrices ${\tilde{A}\\{\ell\\}}$ and
${\hat{C}\\{\ell\\}}$ are still circulant matrices and we denote by
$\mathbf{f}_{\tilde{A}\\{\ell\\}}$ and $\mathbf{f}_{\hat{C}\\{\ell\\}}$ the
respective generating functions. In the Toeplitz setting, the Toeplitz
structure is lost at coarser levels, in particular for the block
${\hat{C}\\{\ell\\}}$. However, ${\hat{C}\\{\ell\\}}$ can still be associated
with the symbol $\mathbf{f}_{\hat{C}\\{\ell\\}}$, computed like in the
circulant case. This is possible because the corrections at the corners have a
rank which remains bounded as $\ell$ increases, as a consequence of the bound
on the matrix bandwidth on which we will comment in the next section.

In principle, our strategy at level $\ell$ is to choose the smoothing
parameter $\omega_{\ell}$ and the grid transfer operators
$P_{\tilde{A}\\{\ell\\}}$ and $P_{\hat{C}\\{\ell\\}}$ such that the TGM
applied to the system with coefficient matrix
$\mathcal{A}_{N_{\ell}}\\{\ell\\}$ converges. The convergence of the “coarser
TGMs” is a necessary condition for the convergence of full multigrid methods
such as the W-cycle and the V-cycle, see [2] for a discussion on the topic.

Combining all the previous considerations, for constructing an MGM for
$\mathcal{A}$ in (2) we aim at fulfilling the assumptions of Theorem 5 applied
to the matrix $\mathcal{A}_{N_{\ell}}\\{\ell\\}$ for all $\ell$. A rigorous
theoretical analysis of the properties of the symbols
$\mathbf{f}_{\tilde{A}\\{\ell\\}}$ and $\mathbf{f}_{\hat{C}\\{\ell\\}}$ in the
block multilevel setting would require a significant amount of computations,
undermining the readability of the paper, so we omit it. Concerning
$\mathbf{f}_{\tilde{A}\\{\ell\\}}$, the reader can refer to [8]. In the
applications, it is usually possible to choose at the finest level the
trigonometric polynomials $\mathbf{p}_{A_{x}}$, $\mathbf{p}_{A_{y}}$, and
$\mathbf{p}_{\hat{C}}$ in such a way that the same choice is suitable for all
levels, that is, the properties of $\mathbf{f}_{\tilde{A}\\{\ell\\}}$ and
$\mathbf{f}_{\hat{C}\\{\ell\\}}$ which are relevant for the TGM convergence
remain unchanged level after level. For instance, this is the case for the
problem and TGM analysed in Section 3, as we will remark in Section 5.

Regarding the choice of $\omega_{\ell}$, considering equation (22) we derive
the following condition:

$\omega_{\ell}<2\min\left\\{2\alpha_{\ell}-\alpha_{\ell}^{2}\max_{j}\left\\{\frac{\|\mathbf{f}_{A_{x}\\{\ell\\}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{x}\\{\ell\\}})}},\frac{\|\mathbf{f}_{A_{y}\\{\ell\\}}\|_{\infty}}{{\hat{a}_{0}^{(j,j)}(\mathbf{f}_{A_{y}\\{\ell\\}})}}\right\\},\right.\\\
\left.\max_{j}\hat{a}_{0}^{(j,j)}(\mathbf{f}_{\hat{C}\\{\ell\\}})\left\|\mathbf{f}_{C\\{\ell\\}}+\mathbf{f}_{B_{x}\\{\ell\\}}\mathbf{f}_{A_{x}\\{\ell\\}}^{-1}\mathbf{f}_{B_{x}\\{\ell\\}}^{H}+\mathbf{f}_{B_{y}\\{\ell\\}}\mathbf{f}_{A_{y}\\{\ell\\}}^{-1}\mathbf{f}_{B_{y}\\{\ell\\}}^{H}\right\|_{\infty}^{-1}\right\\},$

which suggests how the smoothing parameter should be chosen at each level for
obtaining a convergent W-cycle.

###  4.3 Bandwidth at coarser levels

When the generating function is a (multivariate matrix-valued) trigonometric
polynomial, the matrix-vector product with the associated circulant or
Toeplitz matrix has a computational cost linear in the matrix size. If this is
the case for the blocks $A,B$, and $C$, the matrix-vector product with matrix
$\mathcal{A}$ in (18) is still linear in the matrix size. Therefore, in order
to bound the computational cost of an iteration of the multigrid method, it is
fundamental to make sure that the degree of the generating trigonometric
polynomials is bounded as the level increases. In [1] and [7] the authors
prove the latter in the unilevel scalar case, we report the results below.

######  Lemma 2.

[1, Proposition 2] Let ${A\\{\ell\\}}$ be defined in (51), with
$P_{A{\\{\ell\\}}}=\mathcal{C}_{n_{\ell}}(p_{A})K_{n_{\ell}}^{T}$, where
$f_{A\\{0\\}}$ and $p_{A}$ are trigonometric polynomials of degree
$z_{A\\{0\\}}$ and $q_{A}$, respectively. Then $f_{A\\{\ell\\}}$ is a
trigonometric polynomial of degree $z_{A\\{\ell\\}}$ such that:

  1. 1.

$z_{A\\{\ell\\}}\leq\max(z_{A\\{0\\}},2q_{A})$ for all $\ell$;

  2. 2.

$z_{A\\{\ell\\}}\leq 2q_{A}$ for $\ell$ large enough.

######  Lemma 3.

[7] Consider the matrices
${B}\\{\ell\\}=\mathcal{C}_{n_{\ell}}(f_{{B}\\{\ell\\}})$ and
${C}\\{\ell\\}=\mathcal{C}_{n_{\ell}}(f_{{C}\\{\ell\\}})$ defined by formulas
(52)-(53), with
$P_{A{\\{\ell\\}}}=\mathcal{C}_{n_{\ell}}(p_{A})K_{n_{\ell}}^{T}$ and
$P_{\hat{C}{\\{\ell\\}}}=\mathcal{C}_{n_{\ell}}(p_{\hat{C}})K_{n_{\ell}}^{T}$.
Let $z_{B\\{\ell\\}}$ and $z_{C\\{\ell\\}}$ be the degrees of
$f_{{B}\\{\ell\\}}$ and $f_{{C}\\{\ell\\}}$, respectively. Let
$q=\max\\{q_{A},q_{\hat{C}}\\},$ where $q_{A}$ and $q_{\hat{C}}$ are the
polynomial degrees associated with $p_{A}$ and $p_{\hat{C}}$, respectively.

Then, for $\ell$ large enough, it holds

  1. 1.

$z_{B\\{\ell\\}}\leq\max(2z_{B\\{0\\}},4q),$

  2. 2.

$z_{C\\{\ell\\}}\leq\max(4z_{B\\{0\\}},6q,2z_{C\\{0\\}}).$

Seeing a matrix-valued trigonometric polynomial as a matrix whose elements are
trigonometric polynomials, it is straightforward to extend the latter results
to block circulant matrices. Indeed, an element of a generating matrix-valued
trigonometric polynomial $\mathbf{f}_{{A}\\{\ell\\}}$,
$\mathbf{f}_{{B}\\{\ell\\}}$ or $\mathbf{f}_{{C}\\{\ell\\}}$ at coarser levels
is the sum of trigonometric polynomials whose degree is bounded by Lemmas 2–3.

##  5 Numerical Examples

In this section we present the numerical confirmation of the efficiency of the
proposed block multigrid method when applied on the linear system (29).

We exploit the theoretical results presented in Section 2 and Subsection 3.2
to choose the smoothing parameter $\omega$ and the grid transfer operator
$\mathcal{P}=\begin{bmatrix}P_{\tilde{A}}&\\\ &P_{\hat{C}}\end{bmatrix}$.
Precisely, $P_{\tilde{A}}$ and $P_{C}$ are defined by formulae (37) and (44),
respectively. The choice of the Jacobi smoothing parameter $\omega$, which
completes the the multigrid procedure, is related with the range of
admissibility in (45), that is $\omega\in(0,11/12)$.

All the tests are performed using MATLAB 2021a and the error equation at the
coarsest level is solved with the MATLAB backslash function. In all the
experiments we use the standard stopping criterion
$\frac{\|r^{(k)}\|_{2}}{\|\mathbf{b}\|_{2}}<\epsilon$, where
$r^{(k)}=\mathbf{b}-{\mathcal{M}}\mathbf{x}^{(k)}$ and $\epsilon=10^{-6}$.
Concerning $\mathbf{b}$, we consider the case where the true solution
$\mathbf{x}$ of the linear system ${\mathcal{M}}\mathbf{x}=\mathbf{b}$ is a
uniform sampling of $\sin(4t)+\cos(6t)+1$ on $[0,\pi]$ and we compute the
right-hand side $\mathbf{b}$ as $\mathbf{b}={\mathcal{M}}\mathbf{x}$.

Since the grid-transfer operator $P_{{A}}$ has a geometrical meaning—that is,
it represents standard bilinear interpolation—we take the partial dimension
$n$ of the form $2^{t}+1$, $t=5,\dots,8$.

We first consider a TGM consisting of only 1 post smoothing step of damped
Jacobi. We test the efficiency of the TGM with three different values,
$\omega=2/5,3/5,4/5$. The iterations shown in Table 1 remain almost constant
as the matrix size increases, confirming the theoretical optimal convergence
rate.

Table 1:  Two-Grid iterations for the matrix ${\mathcal{A}}_{N}$ with $\epsilon=10^{-6}$ and $\omega_{post}=[2/5,3/5,4/5]$. $N=9\cdot(2^{t}+1)^{2}$ |  $\\#$ Iterations
---|---
$t$ | $\omega_{post}=2/5$ | $\omega_{post}=3/5$ | $\omega_{post}=4/5$
5 | 38 | 23 | 17
6 | 36 | 22 | 16
7 | 35 | 20 | 15
8 | 35 | 19 | 15

As a second experiment we consider different number of steps of pre/post
smoother. In order to damp the error both in the middle and in the high
frequencies, we take a different parameter for the pre-smoother and the post-
smoother (when both are present). In particular, we choose $\omega_{\rm
post}=4/5$ and $\omega_{\rm pre}=3/5$, since they give the fastest convergence
over the considered sample of admissible values. The convergence rate in all
these cases is optimal (see Table 2) since the presence of the pre smoothing
or of additional smoothing steps can accelerate the convergence of the TGM.

Table 2:  Two-Grid iterations for the matrix ${\mathcal{A}}_{N}$ with $\epsilon=10^{-6}$ and with different number of pre/post smoothing steps. If both present, $\omega_{\rm post}=4/5$ and $\omega_{\rm pre}=3/5$. $N=9\cdot(2^{t}+1)^{2}$ |  $\\#$ Iterations
---|---
$t$ | TGM$(0,1)$ | TGM$(1,0)$ | TGM$(1,1)$ | TGM$(2,2)$
5 | 17 | 20 | 16 | 15
6 | 16 | 19 | 16 | 15
7 | 15 | 18 | 15 | 14
8 | 15 | 17 | 14 | 13

In Section 4 we described the structure of the blocks at coarser levels and
how the multigrid ingredients should be chosen when considering more than two
grids. In this case, our choice of trigonometric polynomials
$\mathbf{p}_{A_{x}}$, $\mathbf{p}_{A_{y}}$, and $\mathbf{p}_{\hat{C}}$ ensures
that the spectral properties of $\tilde{A}\\{\ell\\}$ and $\hat{C}\\{\ell\\}$
remain unchanged at coarser levels, as we verified numerically using the
associated symbols. Regarding the choice of $\omega_{\ell}$, we numerically
verified that choosing at all levels the smoothing parameter which is suitable
for the first level is a computationally cheaper and more robust choice in
order to accelerate the convergence of the V-cycle. This approach was already
taken in [7], where a comparison with an adaptive choice of $\omega_{\ell}$
was provided. In Table 3 we numerically prove the independence of the
convergence rate with respect to the matrix size also for the W-cycle and
V-cycle methods using 2 steps of pre/post smoother, with the choice
$\omega_{\rm post}=4/5$ and $\omega_{\rm pre}=3/5$ at all levels.

Table 3:  Comparison of TGM, W-cycle and V-cycle iterations for the matrix ${\mathcal{A}}_{N}$ with 2 iterations of pre/post smoother of Jacobi and the choice $\omega_{\rm post}=4/5$ and $\omega_{\rm pre}=3/5$, $\epsilon=10^{-6}$. $N=9\cdot(2^{t}+1)^{2}$ |  $\\#$ Iterations
---|---
$t$ | TGM | W-cycle | V-cycle
5 | 15 | 15 | 15
6 | 15 | 15 | 15
7 | 14 | 14 | 15
8 | 13 | 13 | 16

In many works [3, 13, 4, 11] problems with such saddle-point structure have
been treated, exploiting particular preconditioned iterative strategies. In
particular, a common choice is that of constructing preconditioning for two-
by-two block systems involving the Schur complement, whose computation
requires the solution of the linear system with coefficient matrix the
positive definite part ${\tilde{A}}$ of ${\mathcal{A}}_{N}$. Then, the
theoretical findings of Subsection 3.2 on $\tilde{A}$ and on the grid transfer
operator $P_{\tilde{A}\\{\ell\\}}$ can be exploited also for strategies
different from the proposed ones. Moreover, Subsection 3.2 suggests valid
tools that can be used for constructing grid transfer operators, since it
shows how to proceed for verifying the conditions which lead to convergence
and optimality. Furthermore, another approach to solve saddle-point problem of
the form (12) consists of multigrid methods in which the smoother takes into
account a special coupling, represented by the off-diagonal blocks in (12).
Then, as a last test we compare our results with the standard Vanka smoother
[28]. The analysis of block smoothers like the one mentioned can be found in
[17, 18, 21]. In Table 4 we compare the iterations of the V-cycle method using
the two different smoothers: column I shows the iterations for the system
$\mathcal{A}_{N}$ considering 2 pre/post smooting steps of Vanka method. The
iterations of the second column are related to the V-cycle method applied to
system $\mathcal{LAU}$, whose efficiency has been already discussed. Even if
the Vanka smoother performs well in combination with the linear interpolation
projecting strategy, we stress that the computational cost of a single Vanka-
based MGM iteration is higher than the one of an iteration of our procedure.

Table 4:  V-cycle iterations with $\epsilon=10^{-6}$, 2 pre 2 post smoothing steps $N=9\cdot(2^{t}+1)^{2}$ |  $\\#$ Iterations of V-cycle
---|---
$t$ | with Vanka for $\mathcal{A}_{N}$ | with Jacobi for $\hat{\mathcal{A}}_{N}$
5 | 10 | 15
6 | 12 | 15
7 | 12 | 15
8 | 11 | 16

##  6 Conclusions

In this paper we studied a symbol-based convergence analysis for problems that
have a saddle-point form with block Toeplitz structure. From theoretical point
of view, we extended the sufficient conditions for the TGM convergence and the
choice of smoothing parameters to the case where the associated generating
function $\mathbf{f}$ is a matrix-valued function. Moreover, we demonstrated
that a recursive procedure is possible, preserving the saddle-point structure
at the coarser levels, after proper symmetrization. In order to test the
efficiency of the theory we considered the saddle-point system involved in the
numerical solution of the Stokes problem in 2D discretized with FEM.

In particular, we exploited the properties of the generating functions of the
blocks and how these can be used to develop a block multigrid method which is
convergent and optimal. Numerically, the resulting methods show the expected
convergence behavior, confirming the validity of our analysis.

The promising results suggest further investigation into the efficiency of the
procedure for other practical applications with challenging problems. The
analysis is noteworthy as it provides with reasonable effort a theoretical
guarantee of convergence by solely examining the symbols of the matrix
sequences involved.

## Acknowledgments

The work of the second, third, and fourth authors is partly supported by
“Gruppo Nazionale per il Calcolo Scientifico” (GNCS-INdAM).

## References

  * [1] A. Aricò and M. Donatelli.  A v-cycle multigrid for multilevel matrix algebras: proof of optimality.  Numer. Math., 105(4):511–547, 2007.
  * [2] A. Aricò, M. Donatelli, and S. Serra-Capizzano.  V-cycle optimal convergence for certain (multilevel) structured linear systems.  SIAM J. Matrix Anal. Appl., 26(1):186–214, 2004.
  * [3] O. Axelsson, S. Farouq, and M. Neytcheva.  Comparison of preconditioned Krylov subspace iteration methods for PDE-constrained optimization problems: Poisson and convection-diffusion control.  Numer. Algorithms, 73(3):631–663, 2016.
  * [4] O. Axelsson and M. Neytcheva.  Eigenvalue estimates for preconditioned saddle point matrices.  Numer. Linear Algebra Appl., 13(4):339–360, 2006.
  * [5] M. Benzi, G. H. Golub, and J. Liesen.  Numerical solution of saddle point problems.  Acta Numer., 14:1–137, 2005.
  * [6] M. Benzi and V. Simoncini.  On the eigenvalues of a class of saddle point matrices.  Numer. Math., 103(2):173–196, 2006.
  * [7] M. Bolten, M. Donatelli, P. Ferrari, and I. Furci.  Symbol based convergence analysis in multigrid methods for saddle point problems (under revision) arXiv 2203.05959.
  * [8] M. Bolten, M. Donatelli, P. Ferrari, and I. Furci.  A symbol-based analysis for multigrid methods for block-circulant and block-Toeplitz systems.  SIAM J. Matrix Anal. Appl., 43(1):405–438, 2022.
  * [9] F. Brezzi and M. Fortin.  Mixed and hybrid finite element methods, volume 15 of Springer Series in Computational Mathematics.  Springer-Verlag, New York, 1991.
  * [10] P. Davis.  Circulant Matrices.  J. Wiley and Sons, New York, 1979.
  * [11] M. Donatelli, A. Dorostkar, M. Mazza, M. Neytcheva, and S. Serra-Capizzano.  Function-based block multigrid strategy for a two-dimensional linear elasticity-type problem.  Comput. Math. Appl., 74(5):1015–1028, 2017.
  * [12] M. Donatelli, P. Ferrari, I. Furci, D. Sesana, and S. Serra-Capizzano.  Multigrid methods for block-circulant and block-Toeplitz large linear systems: Algorithmic proposals and two-grid optimality analysis.  Numer. Linear Algebra Appl., e2356, 2020.
  * [13] A. Dorostkar, M. Neytcheva, and S. Serra-Capizzano.  Spectral analysis of coupled PDEs and of their Schur complements via generalized locally Toeplitz sequences in 2D.  Comput. Methods Appl. Mech. Engrg., 309:74–105, 2016.
  * [14] F. Durastante and I. Furci.  Spectral analysis of saddle-point matrices from optimization problems with elliptic PDE constraints.  Electron. J. Linear Algebra, 36:773–798, 2020.
  * [15] G. Fiorentino and S. Serra.  Multigrid methods for Toeplitz matrices.  Calcolo, 28(3-4):283–305 (1992), 1991.
  * [16] G. Fiorentino and S. Serra.  Multigrid methods for symmetric positive definite block Toeplitz matrices with nonnegative generating functions.  SIAM J. Sci. Comput., 17(5):1068–1081 (1996), 1996.
  * [17] Y. He and S. P. MacLachlan.  Local Fourier analysis of block-structured multigrid relaxation schemes for the Stokes equations.  Numer. Linear Algebra Appl., 25(3):e2147, 28, 2018.
  * [18] S. P. MacLachlan and C. W. Oosterlee.  Local Fourier analysis for multigrid with overlapping smoothers applied to systems of PDEs.  Numer. Linear Algebra Appl., 18(4):751–774, 2011.
  * [19] Y. Notay.  A new algebraic multigrid approach for Stokes problems.  Numer. Math., 132(1):51–84, 2016.
  * [20] D. Noutsos, S. Serra-Capizzano, and P. Vassalos.  Matrix algebra preconditioners for multilevel Toeplitz systems do not insure optimal convergence rate.  Theoret. Comput. Sci., 315(2):557–579, 2004.
  * [21] C. Rodrigo, F. J. Gaspar, and F. J. Lisbona.  On a local Fourier analysis for overlapping block smoothers on triangular grids.  Appl. Numer. Math., 105:96–111, 2016.
  * [22] J. W. Ruge and K. Stüben.  Algebraic multigrid.  In Multigrid methods, volume 3 of Frontiers Appl. Math., pages 73–130. SIAM, Philadelphia, PA, 1987.
  * [23] S. Serra.  Asymptotic results on the spectra of block Toeplitz preconditioned matrices.  SIAM J. Matrix Anal. Appl., 20(1):31–44, 1999.
  * [24] S. Serra-Capizzano.  Asymptotic results on the spectra of block Toeplitz preconditioned matrices.  SIAM Journal on Matrix Analysis and Applications, 20(1):31–44, 1998\.
  * [25] S. Serra-Capizzano.  Spectral and computational analysis of block Toeplitz matrices with nonnegative definite generating functions.  BIT, 39:152–175, 1999.
  * [26] S. Serra-Capizzano and C. Tablino-Possio.  Multigrid methods for multilevel circulant matrices.  SIAM J. Sci. Comput., 26(1):55–85, 2004.
  * [27] D. Sesana and V. Simoncini.  Spectral analysis of inexact constraint preconditioning for symmetric saddle point matrices.  Linear Algebra Appl., 438(6):2683–2700, 2013.
  * [28] S. P. Vanka.  Block-implicit multigrid solution of Navier-Stokes equations in primitive variables.  J. Comput. Phys., 65(1):138–158, 1986.