# Symbol based convergence analysis in multigrid methods for saddle point
problems

Matthias Bolten  Marco Donatelli  Paola Ferrari  Isabella Furci  furci@uni-
wuppertal.de School of Mathematics and Natural Sciences, University of
Wuppertal, Wuppertal, Germany  Department of Science and high Technology,
University of Insubria, Como, Italy

###### Abstract

Saddle point problems arise in a variety of applications, e.g., when solving
the Stokes equations. They can be formulated such that the system matrix is
symmetric, but indefinite, so the variational convergence theory that is
usually used to prove multigrid convergence cannot be applied. In a 2016 paper
in Numerische Mathematik Notay has presented a different algebraic approach
that analyzes properly preconditioned saddle point problems, proving
convergence of the Two-Grid method.

In the present paper we analyze saddle point problems where the blocks are
circulant within this framework. We are able to derive sufficient conditions
for convergence and provide optimal parameters for the preconditioning of the
saddle point problem and for the point smoother that is used. The analysis is
based on the generating symbols of the circulant blocks. Further, we show that
the structure can be kept on the coarse level, allowing for a recursive
application of the approach in a W- or V-cycle and proving the “level
independency” property. Numerical results demonstrate the efficiency of the
proposed method in the circulant and the Toeplitz case.

###### keywords:

Multigrid methods, saddle-point systems, spectral symbol, Toeplitz-like
matrices

††journal: Linear Algebra and its Applications

##  1 Introduction

Saddle point linear systems arise in different cases. One of the most
important examples is the discretization of the Stokes equations that are
given by

$\displaystyle\xi\mathbf{u}-\nu\Delta\mathbf{u}+\nabla p$
$\displaystyle=\mathbf{f},\quad\text{in $\Omega$},$
$\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0,\quad\text{in
$\Omega$},$

where $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ and suitable boundary conditions
are imposed. Here, $\mathbf{u}$ represents the velocity and $p$ the pressure.
They give rise to linear systems

$\mathcal{A}x=b,$ (1)

with

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

where $A$ is symmetric positive definite, $C$ is symmetric nonnegative
definite and $B$ has full rank. The iterative solution of saddle point
problems has been studied extensively, for an introductory overview we refer
to [4]. Here we focus on multigrid methods. Different approaches to solve
saddle point problems using multigrid exist. In most cases more powerful
smoothers are used to take into account the special coupling, represented by
the off-diagonal blocks in (2). This includes the Braess-Sarazin smoother [9],
the Uzawa smoother [18] and the Vanka smoother [29], where the latter one is
probably the most widely used. Usually, these smoothers are applied in
geometric multigrid methods. As such they are often analyzed using local
Fourier analysis (LFA). An introduction to LFA can be found in [30] and the
analysis of block smoothers like the ones mentioned, e.g., in [14, 17, 24].

Instead of altering the smoother Notay has recently presented a method that
applies a point-wise smoother and a coarse grid correction to (1) after
preconditioning from the left and the right using lower and upper triangular
matrices, respectively [22].

In this paper we consider the case where the matrices $A,B,C$ in (2) are
$n\times n$ circulant matrices. Multigrid for circulant matrices are well-
understood and have been analyzed in [28] and the optimality of the V-cycle
has been shown in [2]. Its relation to LFA is described in [11]. Systems of
PDEs, like the Stokes equations, yield block matrices. Multigrid for Toeplitz
matrices with small blocks has been studied in [16] and recently in [6, 12].
These analyses share that they are based on a variational principle, that
cannot be applied here, as the considered systems do not induce a scalar
product. The analysis of numerical methods for structured matrices not only is
of interest when the problem to solve is posed in a rectangular domain and
possesses constant coefficients. The results also carry over to the case of
non-constant coefficients by means of generalized locally Toeplitz (GLT)
sequences [26, 27], and the developed methods can also be used on more complex
domains, e.g., by using fictitious domain techniques or suitable
discretizations as in [1, 5, 7].

Analyzing the system matrix in the case of circulant blocks $A,B,C$ within the
framework presented in [22] we are able to derive sufficient conditions based
on the symbols of the matrices such that the requirements presented there are
fulfilled. The symbol-based analysis allows to choose optimal parameters for
the left and right preconditioners defined in [22] and for the damped Jacobi
methods used as smoothing procedure. Further, the analysis motivates the
choice of the projector in the case of $C$ approaching the zero matrix. For
the multigrid case, we propose a strategy that keeps the same structure on the
coarse level, allowing for a recursive application in a W- or V-cycle. For
this strategy we are able to show that the degree of the polynomial that
represents the generating symbol is bounded, i.e., the bandwidth of the
matrices on the coarse levels is bounded, as well. Moreover, we prove the
“level independency” property, i.e., the two-grid optimality at a generic
level of the multigrid method, which ensures a robust W-cycle method. Finally,
numerical tests demonstrate the efficiency of the proposed method and the
validity of the theoretical analysis.

This paper is organized as follows. Section 2 defines the notation used in the
paper, in particular concerning the symbol of circulant matrices. Section 3 is
devoted to recall the main results on the convergence of multigrid methods
both for circulant matrices and for the Stokes problem. In particular, Section
3.2 contains the main results on the convergence of the Two-Grid method (TGM)
for matrices of the form (2). Section 4 is dedicated to the theoretical
analysis of the TGM convergence in terms of the generating functions and in
Section 4.1 we present a key example showing the numerical efficiency of the
derived convergence results. Section 5 extends the TGM convergence analysis to
multigrid methods providing a strategy to preserve the same structure of the
coefficient matrices at the coarser levels and proving the “level
independency” property. The numerical results in Section 6 confirm the linear
convergence of the W-cycle, while in Section 7 the proposed multigrid method
is applied to Toeplitz matrices obtaining the same optimal convergence
behavior also with the V-cycle method. Some conclusions and future research
lines are drawn in Section 8.

##  2 Notation and definitions

For $X\in\mathbb{C}^{n\times}$, we denote by $\Lambda(X)$ the set of all the
eigenvalues of $X$ and with $\rho(X)$ its spectral radius. If $X$ and $Y$ are
Hermitian matrices, then the notation $X\leq Y$ (resp. $X<Y$ ) means that
$Y-X$ is a nonnegative definite (resp. positive definite) matrix. Moreover, we
numerate such eigenvalues adopting the following notation

$\lambda_{\min}(X)=\lambda_{1}\leq\lambda_{2}\leq\dots\leq\lambda_{n}=\lambda_{\max}(X).$

If $X$ is a Hermitian positive definite (HPD) matrix, then
$\|{{\mathbf{v}}}\|_{X}={\|}X^{1/2}{{\mathbf{v}}}{\|}_{2}$ (resp.
$\|Y\|_{X}=\|X^{1/2}YX^{-1/2}\|_{2}$) denotes the Euclidean norm weighted by
$X$ on $\mathbb{C}^{n}$ (resp. on $\mathbb{C}^{n\times n}$).

We denote by $I_{n}$ and $O_{n}$ the $n\times n$ identity matrix and the
matrix of all zeros respectively. Moreover, ${\mathbf{e}}_{n}$ and
${\mathbf{0}}_{n}$ are respectively the vectors of length $n$ of all ones and
zeros. When the dimension is clear from the context, we omit the subscript
$n$. Given a matrix $X$, we denote by $D_{X}=\operatorname*{diag}(X)$ the
diagonal matrix having ones as elements on the main diagonal of $X$.

Defining the $n$ equispaced grid points

$\theta_{j,n}=\frac{(j-1)2\pi}{n},\qquad j=1,\ldots,n,$

for the interval $[0,2\pi)$, the matrix $\mathbb{F}_{n}$ is the so called
Fourier matrix of order $n$ given by

$\displaystyle(\mathbb{F}_{n})_{i,j}=\frac{1}{\sqrt{n}}{\rm
e}^{\hat{\imath}(i-1)\theta_{j,n}},\quad i,j=1,\ldots,n.$ (3)

######  Definition 1.

Let the Fourier coefficients of a given function $f\in L^{1}(-\pi,\pi)$ be

$\displaystyle\hat{a}_{j}(f):=\frac{1}{2\pi}\int_{Q}f(\theta){\rm
e}^{-\hat{\imath}j\theta}d\theta\in\mathbb{C},\qquad\hat{\imath}^{2}=-1,\,j\in\mathbb{Z}.$
(4)

Then, the ${n}$th circulant matrix $\mathcal{C}_{n}(f)$ associated with $f$ is
given by

$\mathcal{C}_{n}(f)=\sum_{|j|<n}\hat{a}_{j}(f)Z_{n}^{j}=\mathbb{F}_{n}D_{n}(f)\mathbb{F}_{n}^{H},$
(5)

where $Z_{n}^{(j)}$ is the $n\times n$ matrix whose $(i,k)$ entry equals 1 if
$\mathrm{mod}(i-k,n)=1$ and zero otherwise. Moreover,

$D_{n}(f)=\operatorname*{diag}\left(s_{n}(f(\theta_{j,n}))\right),\quad
j=1,\ldots,n,$

where $s_{n}(f(\theta))$ is the $n$th Fourier sum of $f$ given by

$s_{n}(f({\theta}))=\sum_{k=1-n}^{n-1}\hat{a}_{k}(f){\rm e}^{\iota k\theta}.$

The set $\\{\mathcal{C}_{n}(f)\\}_{n\in\mathbb{N}}$ is called the family of
circulant matrices generated by $f$, that in turn is referred to as the
generating function or the symbol of
$\\{\mathcal{C}_{n}(f)\\}_{n\in\mathbb{N}}$. Note that
$\rho(\mathcal{C}_{n}(f))\leq\|f\|_{\infty}$ for all $n\in\mathbb{N}$. The set
of circulant matrices of the same size $n$ defines a _matrix algebra_ since it
is closed by sum, product, and inversion. From a computational point of view,
allocating only the vector of the Fourier coefficients $\hat{a}_{j}(f)$, all
the computations (matrix-vector product, inversion, etc.) involving the matrix
$\mathcal{C}_{n}(f)$ can be computed using the FFT.

######  Remark 1.

If ${f}$ is a trigonometric polynomial of fixed degree less than $n$, the
entries of $D_{n}(f)$ are the eigenvalues of $\mathcal{C}_{n}(f)$, explicitly
given by sampling $f$ using the grid $\theta_{j,n}$:

$\displaystyle\lambda_{j}(\mathcal{C}_{n}(f))$
$\displaystyle=f\left(\theta_{j,n}\right),\quad j=1,\ldots,n,$ $\displaystyle
D_{n}(f)$
$\displaystyle=\operatorname*{diag}\left(f\left(\theta_{j,n}\right)\right),\quad
j=1,\ldots,n.$

Given a circulant matrix $A\in\mathbb{C}^{n\times n}$, we denote by $f_{A}$
its symbol such that $A=\mathcal{C}_{n}(f_{A})$. Moreover, if $f_{A}$ is a
trigonometric polynomial, its degree is denoted by $z_{A}$ and $A$ is a band
matrix with bandwidth $z_{A}+1$.

##  3 Multigrid methods

This section collects relevant results concerning the convergence theory of
algebraic multigrid methods. We first recall the approximation property and
then its equivalent condition in terms of the symbol of circulant matrices.
Next we report the proposal in [22] for the Stokes problem that will be
exploited in the case of circulant matrices in Section 4.

In the general case, we are interested in solving a linear system
$X_{n}x_{n}=b_{n}$ where $X_{n}$ is HPD. Assume $k<n$ and define a full-rank
rectangular matrix $P_{n,k}\in\mathbb{C}^{n\times k}$, which is used as a grid
transfer operator to reduce the problem size.

###  3.1 TGM and approximation property

A TGM combines smoothing iterations with a coarse grid correction, which
requires the solution of the error equation on a subspace of reduced
dimension. In this paper, we consider only a post smoother, which consists in
a single step of the damped Jacobi method with iteration matrix $I_{n}-\omega
D_{X_{n}}^{-1}X_{n}.$

The global iteration matrix of TGM is given by

${\rm TGM}(X_{n},P_{n,k},\omega)=(I_{n}-\omega
D_{X_{n}}^{-1}X_{n})\left[I_{n}-P_{n,k}\left(P_{n,k}^{H}X_{n}P_{n,k}\right)^{-1}P_{n,k}^{H}X_{n}\right].$

The convergence results focus on the choice of $P_{n,k}$ and $\omega$ such
that the spectral radius of ${\rm TGM}(X_{n},P_{n,k},\omega)$ is strictly less
than 1.

Following the Ruge and Stüben approach [25], later generalized by Notay in
[20, 22, 23], we introduce the so-called approximation property.

######  Definition 2.

Let $X_{n}$ be an HPD $n\times n$ matrix. Let $P_{n,k}$ a full-rank $n\times
k$ matrix, $k<n$. If there exists a constant
${\kappa}({X_{n},P_{n,k}})\in\mathbb{R}$ such that

$\min_{\mathbf{v}\in\mathbb{C}^{k}}\|\mathbf{u}-P_{n,k}\mathbf{v}\|^{2}_{D_{X_{n}}}\leq{\kappa}({X_{n}},P_{n,k})\|\mathbf{u}\|^{2}_{X_{n}},\quad\forall\,\mathbf{u}\in\mathbb{C}^{n},$
(6)

then the pair $(X_{n},P_{n,k})$ is said to fulfil the _approximation property_
and the constant ${\kappa}({X_{n}},P_{n,k})$ is an associated approximation
property constant.

Starting from the study in [13], many results have been given on the choice of
the prolongation and restriction operators for Toeplitz and circulant systems
[10, 15].

Let $n$ be even, the common approach for a circulant matrix
$\mathcal{C}_{n}(f)$, where $f$ is a nonegative trigonometric polynomial,
consists in choosing the grid transfer operator

$P_{\mathcal{C}_{n}(f)}=\mathcal{C}_{n}(p)K_{n}^{T}\in\mathbb{C}^{n\times\frac{n}{2}},$

where the trigonometric polynomial $p$ is chosen according to Lemma 1 and the
matrix

$K_{n}=\left[\begin{array}[]{cccccccc}1&0&&&&&\\\ &&1&0&&&&\\\
&&&&\ddots&\ddots&&\\\ &&&&&&1&0\end{array}\right]_{\frac{n}{2}\times n}.$ (7)

is the down-sampling operator.

Using the Galerkin approach, the classical TGM convergence theorem for
circulant matrices was proved in [28] for the approximation property as
formulated in [25]. Here we prove that the same conditions satisfy the
approximation property according to Definition 2.

######  Lemma 1.

Let $A=\mathcal{C}_{n}(f)$, with $f$ being a nonnegative trigonometric
polynomial such that $f(\theta_{0})=0$ and $f(\theta)>0$ for all
$\theta\in[0,2\pi)$. Let $P_{A}=\mathcal{C}_{n}(p)K_{n}^{T}$, with $p$
satisfying:

  1. 1.

$|p|^{2}(\theta)+|p|^{2}(\theta+\pi)>0\quad\forall\theta\in[0,2\pi),$

  2. 2.

$\limsup_{\theta\to\theta_{0}}\frac{|p|^{2}(\theta+\pi)}{f(\theta)}<c$.

Then the pair $(A,P_{A})$ fulfills the approximation property in equation (6)
with

${\kappa}({A},P_{A})=2\hat{a}_{0}(f)\left\|\frac{|p|^{2}(\theta+\pi)}{f(\theta)}\right\|_{\infty}\left\|\frac{1}{|p|^{2}(\theta)+|p|^{2}(\theta+\pi)}\right\|_{\infty}.$
(8)

###### Proof.

Fixing ${\mathbf{v}}=(P_{A}^{H}P_{A})^{-1}P_{A}^{H}{\mathbf{u}}$, the
condition (6) is implied by

$\|\mathbf{u}-P_{n,k}\mathbf{v}\|^{2}_{D_{A}}\leq{\kappa}(A,P_{A})\|\mathbf{u}\|^{2}_{A},\quad\forall\,\mathbf{u}\in\mathbb{C}^{n},$

which is equivalent to the matrix inequality

$\hat{a}_{0}(f)(I-P_{A}(P_{A}^{H}P_{A})^{-1}P_{A}^{H})\leq{\kappa}(A,P_{A})A$

since $D_{A_{n}}=\hat{a}_{0}(f)I_{n}$. By performing a block diagonalization
of all the involved matrices (see [28]), to have (6), it is enough to prove

${\kappa}({A},P_{A})\begin{bmatrix}&f(\theta)&0\\\
&0&f(\theta+\pi)\end{bmatrix}\geq\frac{\hat{a}_{0}(f)}{|p|^{2}(\theta)+|p|^{2}(\theta+\pi)}\begin{bmatrix}&|p|^{2}(\theta+\pi)&-p(\theta)p(\theta+\pi)\\\
&-p(\theta)p(\theta+\pi)&|p|^{2}(\theta)\end{bmatrix},$

for all grid point
$\theta\in\left\\{\theta_{j,n}=\frac{(j-1)2\pi}{n}\,\big{|}\;j=1,\dots,n\right\\}$.
In conclusion, ${\kappa}({A},P_{A})$, defined as in (8), satisfies

${\kappa}({A},P_{A})\geq\frac{\hat{a}_{0}(f)}{|p(\theta)|^{2}+|p(\theta+\pi)|^{2}}\left(\frac{|p|^{2}(\theta+\pi)}{f(\theta)}+\frac{|p|^{2}(\theta)}{f(\theta+\pi)}\right),\qquad\forall\theta\in[0,2\pi),$

and hence the inequality (6) is true. ∎

######  Remark 2.

Let $A=\mathcal{C}_{n}(f)$, if $f$ vanishes at a grid point
$\theta_{0}=\theta_{j,n}=\frac{(j-1)2\pi}{n},j=1,\dots,n$, then $A$ is
singular and the approximation property in Definition 2 cannot be applied.
Nevertheless, Lemma 1 still holds applying a small rank one correction to $A$
obtaining a HPD matrix as in [3], or the system matrix is “naturally” singular
in the sense of [21, section 3.2] and the kernel is in the range of the
prolongation.

###  3.2 TGM for the Stokes problem

Considering the Stokes problem, the TGM proposed in [22] converges thanks to
the following result.

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

Let

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

be a matrix such that $A$ is an $n\times n$ HPD matrix and $C$ is an $m\times
m$ nonnegative definite matrix. Assume that $B$ has rank $m\leq n$ or that $C$
is positive definite on the null space of $B^{T}$.

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_{m}\end{bmatrix},\qquad\mathcal{U}=\begin{bmatrix}I_{n}&-\alpha
D_{A}^{-1}B^{T}\\\ &I_{m}\end{bmatrix},$ (10)

and

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

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

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

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
(6).

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),$
(13)

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}})}.$

The transformed linear system with coefficient matrix $\hat{\mathcal{A}}$
defined in (10) allows to study separately the approximation property for the
two matrices $A$ and $\hat{C}$, which is HPD thanks to the choice of $\alpha$.
All $\alpha\in\left(0,2/\lambda_{\max}(D_{A}^{-1}A)\right)$ ensure the
convergence of TGM, but $\hat{\gamma}_{A}$, and hence $\rho({\rm
TGM}(\mathcal{\hat{A}},\mathcal{P},\omega))$, depends on $\alpha$. A good
compromise is to choose $\alpha\approx\lambda_{\max}(D_{A}^{-1}A)^{-1}$ which
minimizes $\hat{\gamma}_{A}$. On the other hand, for fixed $\alpha$, the
relaxation parameter $\omega$ of the post smoother should be chosen in order
to minimize the bound on $\rho({\rm
TGM}(\mathcal{\hat{A}},\mathcal{P},\omega))$, i.e., the maximum in inequality
(13).

###  3.3 Multigrid methods

The TGM is useful for practical and preliminary convergence analysis of a
multigrid method, but in practical applications, even the coarser problem is
too large to be solved directly. Hence, a simple strategy is to apply
recursively the same algorithm at the coarser error equation instead of
solving it directly. Such a recursive application, until a small size problem
is obtained, is known as V-cycle. A more robust multigrid method can be
obtained concatenating two recursive calls resulting in the so-called W-cycle.

Since the solutions of the error equations at the coarser levels are only
approximated, the TGM convergence study is necessary but not sufficient to
have a robust multigrid method. A more robust result, known as _level
independence_ , is obtained by applying a TGM at a generic level of the
multigrid hierarchy. This allows to obtain a linear convergence rate for the
W-cycle but it is still not enough for the V-cycle [3, 25]. Therefore further
algebraic tools have been introduced for proving the V-cycle optimality, see
e.g. [19, 25].

For circulant matrices, using the Galerkin approach, the following lemma
states that the circulant structure is preserved at the coarser levels with a
symbol depending explicitly on the restriction
$K_{n}\mathcal{C}_{n}(p_{1})^{H}$ and the prolongation
$\mathcal{C}_{n}(p_{2})K_{n}^{T}$.

######  Lemma 2.

[11, Proposition 6] Let $f$ be a trigonometric polynomial. Let $K_{n}$ be
defined as in formula (7). Then the matrix
$\left(\mathcal{C}_{n}(p_{1})K_{n}^{T}\right)^{H}\mathcal{C}_{n}(f)\left(\mathcal{C}_{n}(p_{2})K_{n}^{T}\right)\in\mathbb{C}^{k\times
k}$, $k=\frac{n}{2}$, coincides with $\mathcal{C}_{k}(\hat{f})$ where

$\displaystyle\hat{f}(\theta)=\frac{1}{2}\left(\overline{p}_{1}fp_{2}\left(\frac{\theta}{2}\right)+\overline{p}_{1}fp_{2}\left(\frac{\theta}{2}+\pi\right)\right).$
(14)

Moreover, if $\theta_{0}\in[0,2\pi)$ is a zero of $f$ and the projectors
$p_{1}$ and $p_{2}$ satisfy the conditions 1 and 2 in Lemma 1, then
$\hat{\theta}_{0}=2\theta_{0}\,{\rm mod}\,2\pi$ is a zero of $\hat{f}$ and the
two zeros have the same order.

The properties of the coarser symbols are crucial to study the multigrid
convergence. Indeed, thanks to Lemma 2, it is possible to prove the level
independence under the same TGM assumptions in Lemma 1. While for the V-cycle
optimality, the condition 2 on the projector has to be replaced with the
stronger condition
$\limsup_{\theta\to\theta_{0}}\frac{p(\theta+\pi)}{f(\theta)}<c$ (see [3]).

For saddle point problems, at least to our knowledge, there are no results in
the literature that extend the general algebraic TGM analysis in Theorem 1 to
multigrid methods. Therefore, in Section 5, we will provide a multigrid
analysis in the case of circulant blocks.

##  4 TGM convergence for saddle point matrices with circulant blocks

In this section, we prove how to take advantage of the result in Theorem 1 in
the case where the blocks in the saddle-point problem (1)-(2) involve
circulant matrices. In this case, the following preliminary lemma is useful to
obtain an HPD matrix $\hat{C}$ in (11).

######  Lemma 3.

Let $A=\mathcal{C}_{n}(f)$ with $f\geq 0$. If
$0<\alpha<\frac{2\hat{a}_{0}(f)}{\|f\|_{\infty}}$, then the matrix $2\alpha
D^{-1}_{A}-\alpha^{2}D^{-1}_{A}AD^{-1}_{A}$ is HPD.

###### Proof.

Since $A=\mathcal{C}_{n}(f)$ with $f\geq 0$, then
$D^{-1}_{A}=\frac{1}{\hat{a}_{0}(f)}I_{n}$ with $\hat{a}_{0}(f)>0$. Moreover,
the matrix $2\alpha D^{-1}_{A}-\alpha^{2}D^{-1}_{A}AD^{-1}_{A}$ is circulant
and generated by the real function
$g(\theta)=\frac{2\alpha}{\hat{a}_{0}(f)}-\frac{\alpha^{2}}{\hat{a}^{2}_{0}(f)}f(\theta)$,
which is positive for all $\theta$ if
$\alpha\in\big{(}0,\frac{2\hat{a}_{0}(f)}{f(\theta)}\big{)}$. ∎

The following theorem provides a deeper analysis of Theorem 1 using the symbol
analysis of circulant matrices. In particular, it allows to define optimal
projectors and smoothers as numerically confirmed in the example in Subsection
4.1. This is the first step towards a multigrid analysis provided in Section
5.

######  Theorem 2.

Let us define the matrix

$\mathcal{A}=\begin{bmatrix}A&B^{T}\\\ B&-C\end{bmatrix}_{2n\times 2n},$ (15)

where

  * 1.

$A=\mathcal{C}_{n}(f_{A})$ with $f_{A}$ trigonometric polynomial such that
$f_{A}(\theta_{0})=0$ and $f_{A}(\theta)>0$ for all
$\theta\in[0,2\pi]\setminus\\{\theta_{0}\\}$;

  * 2.

$C=\mathcal{C}_{n}(f_{C})$ with $f_{C}$ nonnegative trigonometric polynomial;

  * 3.

$B=\mathcal{C}_{n}(f_{B})$ with $f_{B}$ trigonometric polynomial such that
$|f_{B}|^{2}({\theta}_{0})=0$ and $|f_{B}|^{2}(\theta)>0$ for all
$\theta\in[0,2\pi]\setminus\\{\theta_{0}\\}$, and
$\limsup_{\theta\to\theta_{0}}\frac{|f_{B}|^{2}(\theta)}{f_{A}(\theta)}<\infty$.

Let $\alpha$ be a positive number such that
$\alpha<\frac{2\hat{a}_{0}(f_{A})}{\|f_{A}\|_{\infty}}$ and define
$\mathcal{\hat{A}}=\mathcal{LAU}$ such as in Theorem 1. Let $\hat{C}$ be
defined as in (11), then $\hat{C}=\mathcal{C}_{n}(f_{\hat{C}})$ with

$f_{\hat{C}}(\theta)=f_{C}(\theta)+\frac{\alpha|f_{B}|^{2}(\theta)}{\hat{a}_{0}(f_{A})}\left(2-\frac{\alpha}{\hat{a}_{0}(f_{A})}f_{A}(\theta)\right).$
(16)

Let $\mathcal{P}$ be defined as in (12), where
$P_{A}=\mathcal{C}_{n}(p_{A})K_{n}^{T}$ and
$P_{\hat{C}}=\mathcal{C}_{n}(p_{\hat{C}})K_{n}^{T}$, with $p_{A}$ and
$p_{\hat{C}}$ trigonometric polynomials.

Consider a TGM associated with one iteration of damped Jacobi as postsmoothing
with relaxation parameter $\omega$. If

  1. 1.

$|p_{A}|^{2}(\theta)+|p_{A}|^{2}(\theta+\pi)>0$ and
$|p_{\hat{C}}|^{2}(\theta)+|p_{\hat{C}}|^{2}(\theta+\pi)>0$ for all
$\,\theta\in[0,2\pi];$

  2. 2.

$\lim\sup_{\theta\to\theta_{0}}\frac{|p_{A}|^{2}(\theta+\pi)}{f_{A}(\theta)}<\infty$
and
$\lim\sup_{\theta\to{\theta_{0}}}\frac{|p_{\hat{C}}|^{2}(\theta+\pi)}{f_{\hat{C}}(\theta)}<\infty;$

  3. 3.

$\omega<2\min\left(2\alpha-\frac{\alpha^{2}}{\hat{a}_{0}(f_{A})}\|f_{A}\|_{\infty}\,,\,\hat{a}_{0}(f_{\hat{C}})\left\|f_{C}+\frac{|f_{B}|^{2}}{f_{A}}\right\|_{\infty}^{-1}\right)$.

Then,

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

###### Proof.

For the sake of simplicity, let us assume that $\theta_{0}$ is different from
each grid point $\theta_{i,n}$, $i=1,\dots,n$ (for the case of $\theta_{0}$
equal to a grid point see Remark 3), then $A$ is HPD. Analogously, $C$ is
nonnegative definite and $B$ is a full-rank matrix. Since
$\alpha<\frac{2\hat{a}_{0}(f_{A})}{\|f_{A}\|_{\infty}}$, thanks to Lemma 3 and
the nonnegativity of $C$, we have that $\hat{C}$ is HPD. Hence, thanks to
hypotheses 1., 2., and Lemma 1, the pairs $(A,P_{A})$ and
$({\hat{C}},P_{\hat{C}})$ fulfil the approximation property (6). Therefore, we
can apply Theorem 1, which implies that the TGM applied to $\mathcal{\hat{A}}$
satisfies the inequality (13), i.e.,

$\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).$

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}_{A}-1$ and
$\omega\hat{\gamma}_{\hat{C}}-1$ in the maximum, we estimate

$\hat{\gamma}_{A}=\left(\alpha\left(2-\alpha\lambda_{\max}(D_{A}^{-1}A)\right)\right)^{-1}\leq\left(2\alpha-\frac{\alpha^{2}}{\hat{a}_{0}(f_{A})}\|f_{A}\|_{\infty}\right)^{-1},$
(17)
$\hat{\gamma}_{\hat{C}}=\lambda_{\max}\left(D_{\hat{C}}^{-1}(C+BA^{-1}B^{T})\right)\leq\frac{\left\|f_{C}+\frac{|f_{B}|^{2}}{f_{A}}\right\|_{\infty}}{\hat{a}_{0}(f_{\hat{C}})},$
(18)

where in the latter inequality we are using
$\limsup_{\theta\to\theta_{0}}\frac{|f_{B}|^{2}(\theta)}{f_{A}(\theta)}<\infty$
and hence the generating function $f_{C}+\frac{|f_{B}|^{2}}{f_{A}}$ belongs to
$L^{1}([-\pi,\pi])$. Hence, using the majorization of $\omega$ in the
hypothesis 1., it holds

$\omega\hat{\gamma}_{A}-1<2\left(2\alpha-\frac{\alpha^{2}}{\hat{a}_{0}(f_{A})}\|f_{A}\|_{\infty}\right)\left(2\alpha-\frac{\alpha^{2}}{\hat{a}_{0}(f_{A})}\|f_{A}\|_{\infty}\right)^{-1}-1=1,$
$\omega\hat{\gamma}_{\hat{C}}-1<2\hat{a}_{0}(f_{\hat{C}})\left\|f_{C}+\frac{|f_{B}|^{2}}{f_{A}}\right\|_{\infty}^{-1}\frac{1}{\hat{a}_{0}(f_{\hat{C}})}\left\|f_{C}+\frac{|f_{B}|^{2}}{f_{A}}\right\|_{\infty}-1=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

$\omega<\left(2\alpha-\frac{\alpha^{2}}{\hat{a}_{0}(f_{A})}\|f_{A}\|_{\infty}\right)+\hat{a}_{0}(f_{\hat{C}})\left\|f_{C}+\frac{|f_{B}|^{2}}{f_{A}}\right\|_{\infty}^{-1}.$
(19)

Moreover, using the estimations (17) and (18), we have

$\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-\frac{\alpha^{2}}{\hat{a}_{0}(f_{A})}\|f_{A}\|_{\infty}\right)^{-1}\frac{1}{\hat{a}_{0}(f_{\hat{C}})}\left\|f_{C}+\frac{|f_{B}|^{2}}{f_{A}}\right\|_{\infty}}{\left(2\alpha-\frac{\alpha^{2}}{\hat{a}_{0}(f_{A})}\|f_{A}\|_{\infty}\right)^{-1}+\frac{1}{\hat{a}_{0}(f_{\hat{C}})}\left\|f_{C}+\frac{|f_{B}|^{2}}{f_{A}}\right\|_{\infty}}\right),$
(20)

where we used the fact that the function $(x,y)\mapsto 2\frac{xy}{x+y}$ is
increasing in $(0,+\infty)\times(0,+\infty)$. Then, combining (19) and (20),
we obtain $\omega\tilde{\gamma}<2$. ∎

######  Remark 3.

Similar to Lemma 1, the proof of the previous theorem requires that the
symbols $f_{A}$ and $f_{B}$ do not vanish at a grid point. Otherwise, the
result still holds applying one of the two techniques in Remark 2.

###  4.1 Example: 1D elasticity problem

In this subsection we want to show the numerical efficiency of the TGM
convergence results in Theorem 2 when applied to linear systems stemming from
the finite difference approximation of a one dimensional elasticity problem.

Consider the coupled system of one-dimensional scalar equations

$\begin{cases}-u^{\prime\prime}-v^{\prime}&=g_{1}(x),\\\ u^{\prime}-\rho
v&=g_{2}(x),\end{cases}$ (21)

with $x\in\Omega=[0,1]$, $\rho>0$, and periodic boundary conditions.
Discretizing the problem using standard finite difference methods with
stepsize $h=1/(n+1)$ and scaling by the diagonal matrix $\mathcal{D}^{(1)}$,
we obtain

$\mathcal{A}=\mathcal{D}^{(1)}\begin{bmatrix}A&hB^{T}\\\
hB&-h^{2}C\end{bmatrix}\mathcal{D}^{(1)}=\begin{bmatrix}A&B^{T}\\\
B&-C\end{bmatrix},\qquad\mathcal{D}^{(1)}=\begin{bmatrix}I&O\\\
O&\frac{1}{h}I\\\ \end{bmatrix},$ (22)

where $A,$ $C$, and $B$ are circulant matrices defined by the symbols

$\displaystyle f_{A}(\theta)=2-2\cos(\theta),\quad f_{B}(\theta)=1-{\rm
e}^{\hat{\imath}\theta},\quad f_{C}(\theta)=\frac{2\rho}{3}(2+\cos(\theta)).$
(23)

We prove that these symbols fulfil the hypothesis of Theorem 2. In particular,
$f_{C}>0$ and $|f_{B}|^{2}=f_{A}$, which vanishes in $\theta_{0}=0$ and is
positive in $(0,2\pi)$.

Since $\hat{a}_{0}(f_{A})=2$ and $\|f_{A}\|_{\infty}=4$, we have
$\frac{2\hat{a}_{0}(f_{A})}{\|f_{A}\|_{\infty}}=1$ and hence we choose
$\alpha=\frac{1}{2}$, which is the middle point of the admissible interval
$(0,1)$.

According to (16), the symbol of $\hat{C}$ is

$\displaystyle
f_{\hat{C}}(\theta)=\frac{2\rho}{3}(2+\cos(\theta))+\frac{\alpha(2-2\cos(\theta))}{2}\left(2-\frac{\alpha}{2}(2-2\cos(\theta)\right),$

that for $\alpha=\frac{1}{2}$ is

$f_{\hat{C}}(\theta)=\frac{32\rho+15}{24}+\frac{4\rho-3}{6}\cos\theta-\frac{1}{8}\cos
2\theta.$ (24)

Concerning the grid transfer operators, the symbol

$p_{A}(\theta)=\sqrt{2}(1+\cos(\theta))$

fulfills the assumptions 1. and 2. of Theorem 2 because

$\displaystyle|p_{A}|^{2}(\theta)+|p_{A}|^{2}(\theta+\pi)$
$\displaystyle=8+8\cos^{2}(\theta)>0,\qquad\forall\theta\in[0,2\pi],$
$\displaystyle\limsup_{\theta\to
0}\frac{|p_{A}|^{2}(\theta+\pi)}{f_{A}(\theta)}$
$\displaystyle=\limsup_{\theta\to
0}\frac{|2-2\cos(\theta)|^{2}}{2-2\cos(\theta)}=0.$

On the other hand, $p_{\hat{C}}$ could be chosen as the trivial downsampling
operator since $f_{\hat{C}}>0$ for all $\rho$, but in practice

$\limsup_{\theta\to 0}\frac{1}{f_{\hat{C}}(\theta)}=\frac{1}{2\rho}.$ (25)

Nevertheless, when $\rho$ approaches zero the previous limit goes to infinity
and hence, when $\rho$ is small, it could be useful to choose $p_{\hat{C}}$
such that $p_{\hat{C}}(\pi)=0$. Therefore, we choose

$p_{\hat{C}}(\theta)=p_{A}(\theta)=\sqrt{2}(1+\cos(\theta))$

such that

$\limsup_{\theta\to
0}\frac{|p_{\hat{C}}|^{2}(\theta+\pi)}{f_{\hat{C}}(\theta)}=0,\qquad\forall\rho>0.$

To fulfill the assumption 3. of Theorem 2 we require that

$\displaystyle\omega$
$\displaystyle<2\min\left(2\alpha-\frac{\alpha^{2}}{\hat{a}_{0}(f_{A})}\|f_{A}\|_{\infty}\,,\,\hat{a}_{0}(f_{\hat{C}})\left\|f_{C}+\frac{|f_{B}|^{2}}{f_{A}}\right\|_{\infty}^{-1}\right).$

For $f_{\hat{C}}$ defined as in (24), it holds

$\hat{a}_{0}(f_{\hat{C}})\left\|f_{C}+\frac{|f_{B}|^{2}}{f_{A}}\right\|_{\infty}^{-1}=\,\frac{32\rho+15}{48\rho+24}\,.$

and hence, since $\alpha=\frac{1}{2}$, we have

$\displaystyle\omega$
$\displaystyle<2\min\left(\frac{1}{2}\,,\,\frac{32\rho+15}{48\rho+24}\right)=1.$

Therefore, Theorem 2 ensures that the TGM applied to the system having
$\mathcal{\hat{A}}$ as coefficient matrix converges, but it remains to
estimate the best $\omega\in(0,1)$. This can be done minimizing the upper
bound in (13). Note that such $\omega_{\rm opt}$ could be different from the
value that minimizes $\rho({\rm TGM}(\mathcal{\hat{A}},\mathcal{P},\omega))$.
Nevertheless, the numerical results confirm that it gives the minimum number
of iterations to convergence.

We already provided upper bounds for the quantities $\hat{\gamma}_{A},$
$\hat{\gamma}_{\hat{C}}$, and $\tilde{\gamma}$ in the proof of Theorem 2.
Precisely,

$\displaystyle\hat{\gamma}_{A}\leq\left(2\alpha-\frac{\alpha^{2}}{\hat{a}_{0}(f_{A})}\|f_{A}\|_{\infty}\right)^{-1}=2,$
$\displaystyle\hat{\gamma}_{\hat{C}}\leq\frac{1}{\hat{a}_{0}(f_{\hat{C}})}\left\|f_{C}+\frac{|f_{B}|^{2}}{f_{A}}\right\|_{\infty}=\frac{48\rho+24}{32\rho+15}=\frac{48}{31},$

and $\tilde{\gamma}\leq 96/55.$

In the remaining part of the section we fix $\rho=1/2$ in order to show how to
make the explicit computation of the optimal parameter $\omega$. Firstly, we
focus on $\kappa(A,P_{A})$, $\kappa(\hat{C},P_{\hat{C}})$ and $\tilde{k}$. In
particular, we use the fact that the pairs
$\left(\mathcal{C}_{n}(f_{A}),\mathcal{C}_{n}(p_{A})K_{n}^{T}\right)$ and
$\left(\mathcal{C}_{n}(f_{\hat{C}}),\mathcal{C}_{n}(p_{\hat{C}})K_{n}^{T}\right)$
satisfy the approximation property (6) and, by Lemma 1,

$\displaystyle{\kappa}(\mathcal{C}_{n}(f_{A}),\mathcal{C}_{n}(p_{A})K_{n}^{T})$
$\displaystyle\leq
2\hat{a}_{0}(f_{A})\left\|\frac{|p_{A}|^{2}(\theta+\pi)}{f_{A}(\theta)}\right\|_{\infty}\cdot\left\|\frac{1}{|p_{A}|^{2}(\theta)+|p_{A}|^{2}(\theta+\pi)}\right\|_{\infty}$
$\displaystyle=4\cdot\left\|\frac{(2+2\cos(\theta+\pi))^{2}}{2(2-2\cos(\theta))}\right\|_{\infty}\cdot\left\|\frac{2}{(2+2\cos(\theta))^{2}+(2-2\cos(\theta))^{2}}\right\|_{\infty}$
$\displaystyle=4\cdot\left\|\frac{(2-2\cos(\theta))}{2}\right\|_{\infty}\cdot\left\|\frac{2}{12+4\cos(2\theta)}\right\|_{\infty}=4\cdot
2\cdot\frac{1}{4}=2.$

Since $\hat{a}_{0}(f_{\hat{C}})=31/24$,

$\displaystyle{\kappa}(\mathcal{C}_{n}(f_{\hat{C}}),\mathcal{C}_{n}(p_{\hat{C}})K_{n}^{T})$
$\displaystyle\leq
2\hat{a}_{0}(f_{\hat{C}})\left\|\frac{|p_{\hat{C}}|^{2}(\theta+\pi)}{f_{\hat{C}}(\theta)}\right\|_{\infty}\cdot\left\|\frac{1}{|p_{\hat{C}}|^{2}(\theta)+|p_{\hat{C}}|^{2}(\theta+\pi)}\right\|_{\infty}$
$\displaystyle=\frac{31}{24}\cdot\left\|\frac{(2-2\cos\theta)^{2}}{\frac{31}{24}-\frac{1}{6}\cos(\theta)-\frac{1}{8}\cos(2\theta)}\right\|_{\infty}\cdot\frac{1}{4}=\frac{31}{8}$

and $\tilde{k}\leq 124/47$.

Then we need to minimize the function

$\mu:\omega\mapsto\max\left(1-\frac{\omega}{2},1-\frac{31}{8}\omega,2\omega-1,\frac{48}{31}\omega-1,\sqrt{\frac{1128}{1705}\omega^{2}-\frac{47}{62}\omega+1}\right)$
(26)

in order to choose the relaxation parameter $\omega$ such that the TGM
convergence is as fastest as possible. Figure 1 depicts the five functions in
(26) and the minimizer of $\mu$, denoted as $\omega_{\rm opt}$ and computed as
the minimum of the parabola
$\frac{1128}{1705}\omega^{2}-\frac{47}{62}\omega+1$. This parameter
$\omega_{\rm opt}=\frac{55}{96}$ gives the upper bound

$\rho({\rm TGM}(\mathcal{\hat{A}},\mathcal{P},\omega_{\rm
opt}))\leq\sqrt{\frac{1128}{1705}\omega_{\rm opt}^{2}-\frac{47}{62}\omega_{\rm
opt}+1}\approx 0.8848$

and hence the TGM has a linear convergence since the bound of the spectral
radius does not depend on the matrix size.

Figure 1: Plot of the functions in (26).

##  5 Multigrid analysis

In this section, we extend the previous TGM convergence to the multigrid
method. In particular, for each recursion level of the multigrid method,
firstly we apply a diagonal scaling to preserve the same structure of the
coefficient matrices, then we prove the TGM convergence (level independence)
and the band structure of the involved matrices keeping a linear cost of the
matrix-vector product. It follows that the W-cycle has a constant convergence
rate thanks to an automatic estimation of the smoothing parameters.

For the TGM defined in Theorem 1, after the projection by $\mathcal{P}$ the
coarser matrix has the same structure of $\mathcal{A}$ in (9) except for the
sign of the last block row. Therefore, changing this sign by a left diagonal
scaling we can apply recursively the TGM. In detail, fix the finer level
matrix as

$\mathcal{A}_{2n_{0}}\\{0\\}=\begin{bmatrix}{A}\\{0\\}&B\\{0\\}^{T}\\\
B\\{0\\}&-C\\{0\\}\end{bmatrix}=\\\ \left[\begin{array}[]{ccc}{A}&B^{T}\\\
B&-C\end{array}\right],$ (27)

with $A,B,C\in\mathbb{C}^{n_{0}\times n_{0}}$ for $n_{0}=2^{\beta}$,
$\beta\in\mathbb{N}$, and $\hat{C}\\{0\\}=\hat{C}$. For each level $\ell\geq
0$, let

$\mathcal{A}_{2n_{\ell}}\\{\ell\\}=\begin{bmatrix}{A}\\{\ell\\}&B\\{\ell\\}^{T}\\\
B\\{\ell\\}&-{C}\\{\ell\\}\end{bmatrix}$

and compute the transformation (10), i.e.,

$\mathcal{\hat{A}}_{2n_{\ell}}\\{\ell\\}=\mathcal{L}\\{\ell\\}\mathcal{A}_{2n_{\ell}}\\{\ell\\}\mathcal{U}\\{\ell\\},$
(28)

where

$\mathcal{L}\\{\ell\\}=\begin{bmatrix}I_{n_{\ell}}\\\
\alpha_{\ell}B\\{\ell\\}D_{A\\{\ell\\}}^{-1}&-I_{n_{\ell}}\end{bmatrix}_{2n_{\ell}\times
2n_{\ell}},\quad\mathcal{U}\\{\ell\\}=\begin{bmatrix}I_{n_{\ell}}&-\alpha_{\ell}D^{-1}_{{A\\{\ell\\}}}B\\{\ell\\}^{T}\\\
&I_{n_{\ell}}\end{bmatrix}_{2n_{\ell}\times 2n_{\ell}}.$

Then we define recursively the sequence of matrices

$\displaystyle\mathcal{A}_{2n_{\ell+1}}\\{\ell+1\\}$
$\displaystyle=\begin{bmatrix}I_{n}&\\\
&-I_{n}\end{bmatrix}\begin{bmatrix}{P}_{{A}\\{\ell\\}}&\\\
&-{P}_{\hat{C}\\{\ell\\}}\end{bmatrix}^{T}\mathcal{\hat{A}}_{2n_{\ell}}\\{\ell\\}\begin{bmatrix}{P}_{{A}\\{\ell\\}}&\\\
&-{P}_{\hat{C}\\{\ell\\}}\end{bmatrix},$ (29)

where

$\begin{split}{A}\\{\ell+1\\}&={P}_{{A}\\{\ell\\}}^{T}{A}\\{\ell\\}{P}_{{A}\\{\ell\\}},\\\
B\\{\ell+1\\}&={P}_{\hat{C}\\{\ell\\}}^{T}B\\{\ell\\}(I_{n_{\ell}}-\alpha_{\ell}D_{{A}\\{\ell\\}}^{-1}{A}\\{\ell\\}){P}_{{A}\\{\ell\\}},\\\
{C}\\{\ell+1\\}&={P}_{\hat{C}\\{\ell\\}}^{T}\hat{C}\\{\ell\\}{P}_{\hat{C}\\{\ell\\}},\\\
\hat{C}\\{\ell\\}&={C}\\{\ell\\}+B\\{\ell\\}(2\alpha_{\ell}D_{{A}\\{\ell\\}}^{-1}-\alpha_{\ell}^{2}D_{{A}\\{\ell\\}}^{-1}{A}\\{\ell\\}D_{{A}\\{\ell\\}}^{-1})B\\{\ell\\}^{T}.\end{split}$
(30)

In formula (30) the matrices ${P}_{{A}\\{\ell\\}}\in\mathbb{R}^{n_{\ell}\times
n_{\ell+1}}$ and ${P}_{\hat{C}\\{\ell\\}}\in\mathbb{R}^{n_{\ell}\times
n_{\ell+1}}$, with $n_{\ell+1}=n_{\ell}/2$, are the prolongation operators
chosen for solving efficiently the scalar systems with coefficient matrix
${A}\\{\ell\\}$ and $\hat{C}\\{\ell\\}$, respectively, where $\alpha_{\ell}$
is such that

$\alpha_{\ell}=\|D_{{A}\\{\ell\\}}^{-1}A\\{\ell\\}\|_{2}^{-1}.$ (31)

For the sake of simplicity, we prove the level independence of the multigrid
procedure described above involving the matrices
$\hat{\mathcal{A}}_{2n_{\ell}}\\{\ell\\}$ defined by formulas (27)-(30) in the
case of symbols vanishing in $\theta_{0}=0$, which arises from the
discretization of boundary values problems. For the more general case
$\theta_{0}\neq 0$ we can apply a change of variable shifting the zero like in
Remark 7 in [2].

######  Lemma 4.

Consider the matrices
$\mathcal{A}_{2n_{\ell}}\\{\ell\\}\in\mathbb{C}^{2n_{\ell}\times 2n_{\ell}}$
defined by formulae (27)-(30). Define for all $\ell$ the grid transfer
operators

$\mathcal{P}\\{\ell\\}=\begin{bmatrix}P_{A{\\{\ell\\}}}&\\\
&P_{{\hat{C}\\{\ell\\}}}\end{bmatrix},\qquad
P_{A{\\{\ell\\}}}=\mathcal{C}_{n_{\ell}}(p_{A})K_{n_{\ell}}^{T},\qquad
P_{\hat{C}\\{\ell\\}}=\mathcal{C}_{n_{\ell}}(p_{\hat{C}})K_{n_{\ell}}^{T}.$

Suppose that $f_{A\\{\ell\\}}$, $f_{C\\{\ell\\}}$, $f_{B\\{\ell\\}}$,
$\alpha_{\ell}$, $p_{A}$ and $p_{\hat{C}}$ fulfil the hypotheses of Theorem 2
with $\theta_{0}=0$. Moreover, assume that

$|p_{A}|^{2}(\theta)+|p_{\hat{C}}|^{2}(\theta+\pi)>0,\qquad\forall\,\theta\in[0,2\pi].$

Then $f_{A\\{\ell+1\\}}(0)=f_{B\\{\ell+1\\}}(0)=0$,
$f_{A\\{\ell+1\\}}(\theta)>0$ and $f_{B\\{\ell+1\\}}(\theta)\neq 0$ for all
$\theta\in(0,2\pi)$, and

$\limsup_{\theta\rightarrow
0}\frac{f_{A\\{\ell\\}}(\theta)}{f_{A\\{\ell+1\\}}(\theta)}=c_{1},\qquad\limsup_{\theta\rightarrow
0}\frac{\left|f_{B\\{\ell+1\\}}(\theta)\right|^{2}}{\left|f_{B\\{\ell\\}}(\theta)\right|^{2}}=c_{2},\qquad
0<c_{1},c_{2}<\infty.$ (32)

Moreover, $f_{{C}\\{\ell+1\\}}(\theta)\geq 0$ for all $\theta\in[0,2\pi]$ and

$\limsup_{\theta\rightarrow
0}\frac{f_{\hat{C}\\{\ell\\}}(\theta)}{f_{\hat{C}\\{\ell+1\\}}(\theta)}<\infty.$
(33)

###### Proof.

The assertion follows directly from Lemmas 8 and 9 in A. ∎

######  Theorem 3.

Consider the matrices
$\mathcal{A}_{2n_{\ell}}\\{\ell\\}\in\mathbb{C}^{2n_{\ell}\times 2n_{\ell}}$
defined by formulae (27)-(30). Define for all $\ell$ the grid transfer
operators

$\mathcal{P}\\{\ell\\}=\begin{bmatrix}P_{A{\\{\ell\\}}}&\\\
&P_{{\hat{C}\\{\ell\\}}}\end{bmatrix},\qquad
P_{A{\\{\ell\\}}}=\mathcal{C}_{n_{\ell}}(p_{A})K_{n_{\ell}}^{T},\qquad
P_{\hat{C}\\{\ell\\}}=\mathcal{C}_{n_{\ell}}(p_{\hat{C}})K_{n_{\ell}}^{T}.$
(34)

Suppose that $f_{A\\{0\\}}$, $f_{C\\{0\\}}$, $f_{B\\{0\\}}$, $\alpha_{0}$,
$p_{A}$ and $p_{\hat{C}}$ fulfil the hypotheses of Theorem 2 with
$\theta_{0}=0$. Moreover, assume that

$|p_{A}|^{2}(\theta)+|p_{\hat{C}}|^{2}(\theta+\pi)>0,\qquad\forall\theta\in[0,2\pi].$

For each $\ell\geq 0$, consider a TGM associated with one iteration of damped
Jacobi as postsmoothing with relaxation parameter

$\omega_{\ell}<2\min\left(2\alpha_{\ell}-\frac{\alpha_{\ell}^{2}}{\hat{a}_{0}(f_{A\\{\ell\\}})}\|f_{A\\{\ell\\}}\|_{\infty}\,,\,\hat{a}_{0}(f_{\hat{C}\\{\ell\\}})\left\|f_{{C}\\{\ell\\}}+\frac{|f_{B\\{\ell\\}}|^{2}}{f_{{A}\\{\ell\\}}}\right\|_{\infty}^{-1}\right).$
(35)

Then, the TGM iteration matrix involving the matrices
$\hat{\mathcal{A}}_{2n_{\ell}}\\{\ell\\}$ defined in (28) is such that

$\rho\left(TGM\left(\hat{\mathcal{A}}_{2n_{\ell}}\\{\ell\\},{\mathcal{P}}\\{\ell\\},\omega_{\ell}\right)\right)<1.$

###### Proof.

The algebra structure of circulant matrices and Lemma 2 imply that the
matrices $A\\{\ell\\}$, $C\\{\ell\\}$ and $B\\{\ell\\}$ defined by formulae
(27)-(30) are circulant matrices themselves. We prove by induction that
$f_{A\\{\ell\\}}$, $f_{C\\{\ell\\}}$, $f_{B\\{\ell\\}}$, $\alpha_{\ell}$,
$p_{A}$ and $p_{\hat{C}}$ fulfil the hypotheses of Theorem 2 for all $\ell$
with $\theta_{0}=0$.

For $\ell=0$, the hypotheses of Theorem 2 are fulfilled by assumption.

Then, we suppose that the hypotheses of Theorem 2 are fulfilled for $\ell$ and
we prove them for $\ell+1$. By Lemma 4 we have

  * 1.

$f_{A\\{\ell+1\\}}(0)=0$, $f_{A\\{\ell+1\\}}(\theta)>0$ for all
$\theta\in(0,2\pi)$;

  * 2.

$f_{{C}\\{\ell+1\\}}(\theta)\geq 0$ for all $\theta\in[0,2\pi]$;

  * 3.

$f_{B\\{\ell+1\\}}(0)=0$, $f_{B\\{\ell+1\\}}(\theta)\neq 0$ for all
$\theta\in(0,2\pi)$;

  * 4.

$\limsup_{\theta\rightarrow
0}\frac{\left|f_{B\\{\ell+1\\}}(\theta)\right|^{2}}{f_{A\\{\ell+1\\}}(\theta)}=\limsup_{\theta\rightarrow
0}\frac{\left|f_{B\\{\ell+1\\}}(\theta)\right|^{2}}{\left|f_{B\\{\ell\\}}(\theta)\right|^{2}}\frac{\left|f_{B\\{\ell\\}}(\theta)\right|^{2}}{f_{A\\{\ell\\}}(\theta)}\frac{f_{A\\{\ell\\}}(\theta)}{f_{A\\{\ell+1\\}}(\theta)}<\infty,$

combining the induction hypothesis with the two limits in (32).

Concerning assumption 2. of Theorem 2, for $A\\{\ell+1\\}$ we use the
induction hypothesis $\limsup_{\theta\rightarrow
0}\frac{\left|p_{A}(\theta+\pi)\right|^{2}}{f_{A\\{\ell\\}}(\theta)}<\infty$
and (32) so that we can write

$\limsup_{\theta\rightarrow
0}\frac{\left|p_{A}(\theta+\pi)\right|^{2}}{f_{A\\{\ell+1\\}}(\theta)}=\limsup_{\theta\rightarrow
0}\frac{\left|p_{A}(\theta+\pi)\right|^{2}}{f_{A\\{\ell\\}}(\theta)}\frac{f_{A\\{\ell\\}}(\theta)}{f_{A\\{\ell+1\\}}(\theta)}<\infty.$
(36)

Concerning $\hat{C}\\{\ell+1\\}$, we write

$\limsup_{\theta\to{0}}\frac{|p_{\hat{C}}|^{2}(\theta+\pi)}{f_{\hat{C}\\{\ell+1\\}}(\theta)}=\limsup_{\theta\to{0}}\frac{|p_{\hat{C}}|^{2}(\theta+\pi)}{f_{\hat{C}\\{\ell\\}}(\theta)}\frac{f_{\hat{C}\\{\ell\\}}(\theta)}{f_{\hat{C}\\{\ell+1\\}}(\theta)}<\infty,$

where, for bounding the two terms, we used the induction hypothesis and
formula (33), respectively.

In conclusion, we proved that $f_{A\\{\ell\\}}$, $f_{C\\{\ell\\}}$,
$f_{B\\{\ell\\}}$, $\alpha_{\ell}$, $p_{A}$ and $p_{\hat{C}}$ fulfil the
hypotheses of Theorem 2 for all $\ell$ with $\theta_{0}=0$. Hence, Theorem 2
guarantees that
$\rho\left(TGM\left(\hat{\mathcal{A}}\\{\ell\\},{\mathcal{P}}\\{\ell\\},\omega_{\ell}\right)\right)<1$,
where $\omega_{\ell}$ is chosen according to (35), and the proof is complete.
∎

When the matrices $A,B$, and $C$ are banded, then the matrix-vector product
with matrix $\mathcal{A}$ in (2) has a computational cost linear in $n$.
Therefore, we would like to preserve the band structure of each block at the
coarser levels such that each iteration of the V-cycle has a computational
cost proportional to $n$. This property is a consequence of the following
lemmas that state that, at every level of the multigrid procedure, the
generating functions of each block are trigonometric polynomials of degree
lower than a constant independent of $n$.

######  Lemma 5.

[2] Let $g$ be a trigonometric polynomial of degree $z_{g}$. Define

$\psi(g)(\theta)=\frac{1}{2}\left[g\left(\frac{\theta}{2}\right)+g\left(\frac{\theta}{2}+\pi\right)\right].$
(37)

Then, $\psi(g)$ is a trigonometric polynomial of degree at most
$\left\lfloor{\frac{z_{g}}{2}}\right\rfloor$.

This lemma implies that, for the classical multigrid method, the bandwidth of
the coefficient matrix at the coarser levels becomes equal to the double of
the bandwidth of the grid transfer operator, even when the coefficient matrix
at the finer level has a large bandwidth.

To clearly distinguish the bandwidth of the grid transfer operators with
respect to the coefficient matrices, we denote by $q_{A}=z_{P_{A}}$ and
$q_{\hat{C}}=z_{P_{\hat{C}}}$ the degrees of the trigonometric polynomials
$P_{A}$ and $P_{\hat{C}}$, respectively.

######  Lemma 6.

[2, Proposition 2] Let ${A\\{\ell\\}}$ be defined in (30), 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.

A similar result holds for the matrix (2) when $A$, $B$, and $C$ are band
matrices.

######  Lemma 7.

Consider the matrices $\mathcal{A}\\{\ell\\}\in\mathbb{C}^{2n_{\ell}\times
2n_{\ell}}$ defined by formulas (27)-(30), with $P_{A{\\{\ell\\}}}$ and
$P_{\hat{C}{\\{\ell\\}}}$ defined as in (34). 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\\}}).$

###### Proof.

To prove item 1. we consider the function $\psi$ defined in equation (37).
Exploiting the structure of $B\\{\ell+1\\}$ in (30), we have that the
associated generating function is

$f_{B\\{\ell+1\\}}=\psi(g_{B\\{\ell\\}}),\qquad
g_{B\\{\ell\\}}(\theta)=\overline{p}_{\hat{C}}p_{A}f_{B\\{\ell\\}}(\theta)\left(1-\frac{\alpha_{\ell}}{\hat{a}_{0}(f_{A\\{\ell\\}})}f_{A\\{\ell\\}}(\theta)\right).$

Therefore, by Lemma 5, we have

$z_{B\\{\ell+1\\}}\leq\left\lfloor{\frac{2q+z_{B\\{\ell\\}}+z_{A\\{\ell\\}}}{2}}\right\rfloor$

and, for $\ell$ large enough, item 2. of Lemma 6 and $q_{A}\leq q$ give

$z_{B\\{\ell+1\\}}\leq
2q+\left\lfloor{\frac{z_{B\\{\ell\\}}}{2}}\right\rfloor.$ (38)

We can now prove item 1. by induction over $\ell$. For $\ell=0$ is trivial.
For the induction step, inserting the induction assumption
$z_{B\\{\ell\\}}\leq\max(2z_{B\\{0\\}},4q)$ in (38), we have

$z_{B\\{\ell+1\\}}\leq
2q+\left\lfloor{\frac{\max(2z_{B\\{0\\}},4q)}{2}}\right\rfloor.$

Distinguishing the two cases in the maximum:

  * 1.

Case $\max(2z_{B\\{0\\}},4q)=4q$ implies $z_{B\\{\ell+1\\}}\leq
4q=\max(2z_{B\\{0\\}},4q)$.

  * 2.

Case $\max(2z_{B\\{0\\}},4q)=2z_{B\\{0\\}}$ implies $z_{B\\{\ell+1\\}}\leq
2q+z_{B\\{0\\}}\leq 2\max(2q,z_{B\\{0\\}})$.

Therefore, the item 1. follows in both cases.

We now prove item 2. exploiting the structure of ${C}\\{\ell+1\\}$ in (30).
Since $f_{C\\{\ell+1\\}}=\psi(|p_{\hat{C}}|^{2}f_{\hat{C}\\{\ell\\}})$,
recalling the definition $\hat{C}\\{\ell\\}$ in (30) and applying Lemma 5, we
have

$z_{C\\{\ell+1\\}}<\left\lfloor\frac{2q_{\hat{C}}+\max(z_{C\\{\ell\\}},2z_{B\\{\ell\\}}+z_{A\\{\ell\\}})}{2}\right\rfloor.$
(39)

Moreover, for $\ell$ large enough, from item 1 follows that
$z_{B\\{\ell\\}}\leq\max(2z_{B\\{0\\}},4q)$ while item 2. of Lemma 6 implies
$z_{A\\{\ell\\}}\leq 2q_{A}\leq 2q$. Therefore, inserting these two
majorizations in (39), we have

$z_{C\\{\ell+1\\}}<q+\left\lfloor\frac{\max(z_{C\\{\ell\\}},\max(4z_{B\\{0\\}},8q)+2q)}{2}\right\rfloor.$
(40)

We can now prove item 2. by induction on $\ell$. For $\ell=0$ is trivial. For
the induction step, inserting the induction assumption
$z_{C\\{\ell\\}}\leq\max(4z_{B\\{0\\}},6q,2z_{C\\{0\\}})$ in (40), we have

$\displaystyle z_{C\\{\ell+1\\}}$
$\displaystyle<q+\left\lfloor\frac{\max\left(\max(4z_{B\\{0\\}},6q,2z_{C\\{0\\}}),\,\max(4z_{B\\{0\\}},8q)+2q\right)}{2}\right\rfloor$
$\displaystyle=q+\max\left(z_{C\\{0\\}},2z_{B\\{0\\}}+q,5q\right).$

Distinguishing the three cases in the maximum
$\xi=\max\left(z_{C\\{0\\}},2z_{B\\{0\\}}+q,5q\right)$, we have:

  * 1.

Case $\xi=z_{C\\{0\\}}$ implies $z_{C\\{\ell+1\\}}\leq q+z_{C\\{0\\}}\leq
2\max(q,z_{C\\{0\\}})$.

  * 2.

Case $\xi=2z_{B\\{0\\}}+q$ implies $z_{C\\{\ell+1\\}}\leq 2q+2z_{B\\{0\\}}\leq
2\max(2q,2z_{B\\{0\\}})$.

  * 3.

Case $\xi=5q$ implies $z_{C\\{\ell+1\\}}\leq 6q$.

Therefore, combining the three cases, we obtain
$z_{C\\{\ell+1\\}}\leq\max(4z_{B\\{0\\}},6q,2z_{C\\{0\\}})$ which is item 2. ∎

##  6 Numerical results

The present section is devoted to the numerical validation of the theoretical
results presented in Sections 4 and 5. In all the experiments, we use the
standard stopping criterion
$\frac{\|r^{(k)}\|_{2}}{\|\mathbf{b}\|_{2}}<\epsilon$, where
$r^{(k)}=\mathbf{b}-{\mathcal{LAU}}\mathbf{x}^{(k)}$ and $\epsilon=10^{-6}$.
The true solution $\mathbf{x}$ of the linear system
${\mathcal{LAU}}\mathbf{x}=\mathbf{b}$ is a uniform sampling of $\sin(t)$ on
$[0,\pi]$ and we consider the right-hand side $\mathbf{b}$ defined as
$\mathbf{b}={\mathcal{LAU}}\mathbf{x}$, which automatically ensures that the
requirements of Remark 3 are satisfied. We take the null initial guess. All
the tests are performed using MATLAB 2021a and the error equation at the
coarsest level is solved with the MATLAB backslash function after the proper
projection of the residual into the range of the coefficient matrix.

Firstly, we investigate the numerical behavior of the TGM applied to the
elasticity problem described in Subsection 4.1 with $\rho=1/2$. In Table 1 we
test the efficiency of the TGM using 1 step of damped Jacobi as post smoother
with four different values of $\omega$. The results show that the number of
iterations needed for reaching the tolerance $\epsilon$ remains constant or
decreases as the matrix size increases, confirming the theoretical optimal
convergence rate. Moreover, the value $\omega_{\rm opt}$ obtained by
minimizing the spectral radius estimate in (26) proves to be the one
associated to the minimum number of iterations, when compared to the other
values belonging to a uniform sampling in the admissible interval $(0,1)$.

Table 1:  Two-Grid iterations with different values of the Jacobi relaxation parameter $\omega$ and stopping tolerance $\epsilon=10^{-6}$.  $N=2\cdot(2^{t})$ |  $\\#$ Iterations
---|---
$t$ | $\omega=1/4$ | $\omega=1/2$ | $\omega_{\rm opt}=55/96$ | $\omega=3/4$
9 | 34 | 14 | 12 | 15
10 | 33 | 14 | 12 | 15
11 | 32 | 14 | 11 | 14
12 | 30 | 13 | 11 | 14
13 | 29 | 13 | 11 | 13
14 | 28 | 12 | 10 | 13

As a second experiment, we numerically validate the convergence results of
Theorem 3 applying the W-cycle strategy analyzed in Section 5. In particular,
the matrices at the finest level
$A_{n_{0}}\\{0\\}=\mathcal{C}_{n_{0}}(f_{A\\{0\\}})$,
$C_{n_{0}}\\{0\\}=\mathcal{C}_{n_{0}}(f_{C\\{0\\}})$,
$B_{n_{0}}\\{0\\}=\mathcal{C}_{n_{0}}(f_{B\\{0\\}})$ are the circulant
matrices given in equations (23), with $n_{0}$ equal to $n=2^{t}$. If we
choose $p_{\hat{C}}(\theta)=p_{A}(\theta)=\sqrt{2}(1+\cos(\theta))$, the
hypotheses the Theorem 3 are fulfilled as long as we make a proper choice of
the relaxation parameters $\omega_{\ell}$. In this way, the thesis of Theorem
3 suggests the convergence and optimality of the W-cycle procedure.

Concerning the choice of $\omega_{\ell}$, we propose an adaptive strategy
approximating at each level $\ell$ the quantities

$2\alpha_{\ell}-\frac{\alpha_{\ell}^{2}}{\hat{a}_{0}(f_{A\\{\ell\\}})}\|f_{A\\{\ell\\}}\|_{\infty}\,,\quad\hat{a}_{0}(f_{\hat{C}\\{\ell\\}})\left\|f_{{C}\\{\ell\\}}+\frac{|f_{B\\{\ell\\}}|^{2}}{f_{{A}\\{\ell\\}}}\right\|_{\infty}^{-1}.$
(41)

In particular, we exploit the fact that for each $\ell$ we have
$f_{A\\{\ell\\}}(\theta)=2-2\cos(\theta)$ and hence
${\hat{a}_{0}(f_{A\\{\ell\\}})}=2$ and, according to formula (31),
$\alpha_{\ell}=1/2$. The invariance of the generating functions
$f_{A\\{\ell\\}}$ at coarser levels with the choice
$p_{A}(\theta)={\sqrt{2}}(1+\cos(\theta))$ can be derived from Lemma 2 by
direct computation.

For the second quantity in (41), we read the Fourier coefficients of
$f_{{B}\\{\ell\\}}$ and $f_{{C}\\{\ell\\}}$ from the first rows and columns of
${{B}\\{\ell\\}}$ and ${{C}\\{\ell\\}}$ and, in particular, we compute the
Fourier coefficient $\hat{a}_{0}(f_{\hat{C}\\{\ell\\}})$. Moreover, we
approximate the
$\|f_{{C}\\{\ell\\}}+\frac{|f_{B\\{\ell\\}}|^{2}}{f_{{A}\\{\ell\\}}}\|_{\infty}$
choosing the maximum value among its uniform sampling over $[0,\pi]$ with
step-size $h=1/100$.

Then, according to (35), we set at each level

$\omega_{\ell}=\min\left(2\alpha_{\ell}-\frac{\alpha_{\ell}^{2}}{\hat{a}_{0}(f_{A\\{\ell\\}})}\|f_{A\\{\ell\\}}\|_{\infty}\,,\,\hat{a}_{0}(f_{\hat{C}\\{\ell\\}})\left\|f_{{C}\\{\ell\\}}+\frac{|f_{B\\{\ell\\}}|^{2}}{f_{{A}\\{\ell\\}}}\right\|_{\infty}^{-1}\right),$
(42)

which is the middle point of the interval of admissible values, inspired by
the fact that at the finest level the $\omega_{\rm opt}=55/96$ is close to the
middle point of the interval $(0,1)$.

In Table 2 we report the number of iterations needed by the W-cycle method for
reaching the convergence with tolerance $\epsilon$, comparing the previous
adaptive choice of the parameters $\omega_{\ell}$ with the fixed choice of
$\omega_{\ell}=1/2$. We observe that the adaptive choice of the relaxation
parameter at each level permits to obtain a number of W-cycle iterations which
remains constant, confirming the theoretical optimal convergence rate. For
this example, the choice of the fixed $\omega_{\ell}=1/2$ for the W-cycle is
also valid and it provides the same behavior in terms of iterations. Indeed,
the presence of the positive definite mass term $C$ implies that the choice of
$\omega_{\ell}$ depends only on the quantities involving $A\\{\ell\\}$ that do
not change level by level.

Table 2: W-cycle iterations with the adaptive choice of $\omega_{\ell}$ according to Theorem 3 and with the fixed choice $\omega_{\ell}=1/2$ for all $\ell$. $N=2\cdot(2^{t})$ |  $\\#$ Iterations
---|---
$t$ | adaptive $\omega_{\ell}$ | $\omega_{\ell}=1/2$
9 | 14 | 14
10 | 14 | 14
11 | 14 | 14
12 | 13 | 13
13 | 13 | 13
14 | 12 | 12

An analogous behavior, although a higher number of iterations, is obtained
using $p_{\hat{C}}\equiv 1$ in the construction of $P_{\hat{C}}$, see column 5
of Table 3. Indeed, as already mentioned in Subsection 4.1, this choice of
$p_{\hat{C}}$ is still acceptable for reasonable large value of $\rho$. Taking
smaller values of $\rho$, for example $\rho=1/20,1/200$ corresponds to weaken
the positive definiteness of the $C\\{0\\}$ part. Then, formula (25)
highlights that in this case is crucial to choose
$p_{\hat{C}}(\theta)={\sqrt{2}}(1+\cos(\theta))$ in order to numerically
satisfy item 2 of Theorem 2. This can be numerically confirmed comparing
columns 2-4 with columns 5-7 of Table 3, where we observe a dramatic
increasing of the number of iteration needed to achieve the convergence of the
W-cycle passing from $p_{\hat{C}}(\theta)={\sqrt{2}}(1+\cos(\theta))$ to
$p_{\hat{C}}(\theta)=1$.

Table 3: W-cycle iterations with the adaptive choice of $\omega_{\ell}$ according to Theorem 3 for different values of $\rho$ and two different projectors by $p_{\hat{C}}(\theta)={\sqrt{2}}(1+\cos(\theta))$ and $p_{\hat{C}}(\theta)=1$.  $N=2\cdot(2^{t})$ | $p_{\hat{C}}(\theta)={\sqrt{2}}(1+\cos(\theta))$ | $p_{\hat{C}}(\theta)=1$
---|---|---
$t$ | $\rho=1/2$ | $\rho=1/20$ | $\rho=1/200$ | $\rho=1/2$ | $\rho=1/20$ | $\rho=1/200$
9 | 14 | 17 | 18 | 24 | 105 | 817
10 | 14 | 16 | 18 | 24 | 107 | 842
11 | 14 | 16 | 18 | 24 | 107 | 866
12 | 13 | 16 | 18 | 24 | 107 | 890
13 | 13 | 15 | 17 | 24 | 108 | 912
14 | 12 | 15 | 17 | 24 | 108 | 932

##  7 Saddle point matrices with Toeplitz blocks

In the present section, we discuss the applicability of our multigrid method.
For completeness, we report the definition of Toeplitz matrix generated by a
function.

######  Definition 3.

The Toeplitz matrix associated with $f\in L^{1}(-\pi,\pi)$ is the matrix of
order $n$ given by

$\displaystyle T_{n}(f)=\sum_{|j|<n}\hat{a}_{j}(f)J_{n}^{(j)},$

where $\hat{a}_{j}(f)$ are the Fourier coefficients of $f$ defined in (4) and
$J_{n}^{(j)}$ is the matrix of order $n$ whose $(i,k)$ entry equals $1$ if
$i-k=j$ and zero otherwise.

###  7.1 Multigrid methods for Toeplitz matrices

The theoretical results of Sections 4 and 5 can be naturally extended to
systems in saddle-point form where the matrices $A$, $B$, $C$ belong to matrix
algebras different from the circulant one. Indeed, as stressed in [2, 3], it
is sufficient to substitute the Fourier transform in (5) with the proper
unitary transform and to choose the corresponding grid points. In particular,
we first assume that $A=\tau_{n}(f_{A})$, $BB^{H}=\tau_{n}(|f_{B}|^{2})$,
$C=\tau_{n}(f_{C})$ belong to the $\tau$ algebra [8]. In this case the
structure of the grid transfer operators has to slightly change with respect
to (7) in order to preserve the algebra structure at coarser levels and the
dimension of the problem should be $2n$ with odd

$n=2^{t}-1.$

Then, the $(2^{t-1}-1)\times(2^{t}-1)$ cutting matrix $K_{n}$ takes the form

$K_{n}=\left[\begin{array}[]{ccccccccc}0&1&0&&&&&\\\ &&0&1&0&&&&\\\
&&&&\ddots&\ddots&\ddots&&\\\ &&&&&&0&1&0\end{array}\right].$

Consequently, the grid transfer operators are chosen of the form
$\tau_{n}(p)K_{n}^{T}$. Note that if $f_{A}$, $|f_{B}|^{2}$, and $f_{C}$ are
trigonometric polynomials of degree at most one, the related Toeplitz matrices
belong to the $\tau$-algebra. If, for instance, the degree $\delta$ of $p_{A}$
is greater than 1, then $T_{n}(p_{A})$ is a Toeplitz matrix which differs from
the $\tau$-matrix $\tau_{n}(p_{A})$ for a low rank correction of rank at most
$2(\delta-1)$. Therefore, the associated grid transfer operator

$P_{A}=T_{n}(p_{A})K_{n}^{T}$ (43)

should be adapted whenever $\delta>1$ to preserve the Toeplitz structure at
the coarser levels (see [3]).

###  7.2 Numerical results for the elasticity problem with Toeplitz blocks

In the present subsection, we consider the problem (21) with Dirichlet
boundary conditions. In this case, the analogous finite difference
discretization with stepsize $h=1/(n+1)$ leads to a linear system with
coefficient matrix of the form

$\mathcal{A}=\mathcal{D}^{(1)}\begin{bmatrix}A&hB^{T}\\\
hB&-h^{2}C\end{bmatrix}\mathcal{D}^{(1)}=\begin{bmatrix}A&B^{T}\\\
B&-C\end{bmatrix},\qquad\mathcal{D}^{(1)}=\begin{bmatrix}I&O\\\
O&\frac{1}{h}I\\\ \end{bmatrix},$ (44)

where $A,$ $C$, and $B$ are the Toeplitz matrices

$\displaystyle A=T_{n}(2-2\cos(\theta)),\quad B=T_{n}(1-{\rm
e}^{\hat{\imath}\theta}),\quad
C=T_{n}\left(\frac{1}{3}(2+\cos(\theta))\right).$ (45)

We construct the grid transfer operators according to the $\tau$ algebra
requirements, as described in Subsection 7.1. Exploiting the analysis for the
circulant case in Subsection 4.1, we firstly construct a TGM for
$\mathcal{LAU}$ in the Toeplitz setting with a basic scheme, using only one
step of damped Jacobi post-smoothing with relaxation parameter $\omega=55/96$.
Table 4 and Figure 2 show that the linear convergence independent of the
matrix size of the TGM is preserved also in the Toeplitz case.

Table 4: Toeplitz case: TGM iterations with $\omega=55/96$.  $t$ |  $\\#$ Iterations
---|---
9 | 13
10 | 12
11 | 12
12 | 11
13 | 11
14 | 11

Figure 2: Convergence rate for TGM method corresponding to $t=14$ of Table 4.

We conclude the section showing the efficiency of the proposed method in the
Toeplitz case with more grids. For the definition of all the objects at
coarser levels, we follow the recursive procedure described at the beginning
of Section 5. Here the size of the matrix is $N=2\cdot(2^{t}-1)$. The grid
transfer operators are associated to the trigonometric polynomials
$p_{\hat{C}}(\theta)=p_{A}(\theta)={\sqrt{2}}(1+1\cos(\theta))$ as in Section
6, while $P_{\hat{C}}$ and $P_{A}$ are defined according to Section 7.1.
Finally, the relaxation parameter $\omega_{\ell}$ at each level is set
according to formula (42).

The possible extension of Theorem 3 to the $\tau$ algebra suggests that the
convergence and optimality of the W-cycle method can be achieved also with
Toeplitz blocks, even though in our case low-rank corrections are present.
Moreover, a linear convergence independent of the matrix size is obtained for
both the W-cycle and V-cycle methods (see Table 5) with a behavior similar to
that of TGM, shown in Table 4. The linear convergence is shown in Figure 3.

Table 5: Toeplitz case: W-cycle and V-cycle iterations with the adaptive choice of $\omega_{\ell}$. $t$ | W-cycle | V-cycle
---|---|---
9 | 13 | 14
10 | 12 | 14
11 | 12 | 14
12 | 11 | 13
13 | 11 | 13
14 | 11 | 13

Figure 3: Convergence rate for W-cycle and V-cycle methods corresponding to
$t=14$ of Table 5.

##  8 Conclusions and future work

Using the results from [22] we have provided sufficient conditions for the
convergence of the TGM in the case of circulant blocks. Further, we have shown
that the structure is kept on coarser levels and that the convergence rate is
bounded independent from the level. While these are important results for the
convergence of the W-cycle method, the convergence of the V-cycle remains
open.

Using the analysis based on the generating symbols of the circulant blocks, we
have been further able to provide optimal choices for the parameters $\alpha$
and $\omega$ in the left and right preconditioning and in the smoothing,
respectively.

Numerically the resulting methods show the expected convergence behavior,
demonstrating the validity of our analysis. Further, the W- and the V-cycle
converge, as well.

In the future, we will extend the analysis to the multilevel case, where the
generating symbols of $A,B$, and $C$ are multi-variate functions. This is of
importance for applications that usually are posed in 2D or 3D, resulting in
2- or 3-level circulant or Toeplitz matrices. In this case we will have to
consider different sizes of $A$ and $C$ and thus rectangular matrices $B$, as
well.

##  Appendix A Level Independency Proofs

######  Lemma 8.

Consider the matrices
$\mathcal{A}_{2n_{\ell}}\\{\ell\\}\in\mathbb{C}^{2n_{\ell}\times 2n_{\ell}}$
defined by formulae (27)-(30). Define for all $\ell$ the grid transfer
operators

$\mathcal{P}\\{\ell\\}=\begin{bmatrix}P_{A{\\{\ell\\}}}&\\\
&P_{{\hat{C}\\{\ell\\}}}\end{bmatrix},\qquad
P_{A{\\{\ell\\}}}=\mathcal{C}_{n_{\ell}}(p_{A})K_{n_{\ell}}^{T},\qquad
P_{\hat{C}\\{\ell\\}}=\mathcal{C}_{n_{\ell}}(p_{\hat{C}})K_{n_{\ell}}^{T}.$

Suppose that $f_{A\\{\ell\\}}$, $f_{C\\{\ell\\}}$, $f_{B\\{\ell\\}}$,
$\alpha_{\ell}$, $p_{A}$ and $p_{\hat{C}}$ fulfil the hypotheses of Theorem 2
with $\theta_{0}=0$. Moreover, assume that

$|p_{A}|^{2}(\theta)+|p_{\hat{C}}|^{2}(\theta+\pi)>0,\qquad\forall\theta\in[0,2\pi].$

Then $f_{B\\{\ell+1\\}}(0)=0$, $f_{B\\{\ell+1\\}}(\theta)\neq 0$ for all
$\theta\in(0,2\pi)$ and

$\limsup_{\theta\rightarrow
0}\frac{\left|f_{B\\{\ell+1\\}}(\theta)\right|^{2}}{\left|f_{B\\{\ell\\}}(\theta)\right|^{2}}=c,\quad
0<c<\infty.$ (46)

###### Proof.

By Lemma 2, for each $\ell$ we have

$2f_{B\\{\ell+1\\}}(\theta)=\left(\overline{p_{\hat{C}}}f_{B\\{\ell\\}}\left(1-\frac{\alpha_{\ell}}{\hat{a}_{0}\left(f_{A\\{\ell\\}}\right)}f_{A\\{\ell\\}}\right)p_{A}\right)\left(\frac{\theta}{2}\right)+\left(\overline{p_{\hat{C}}}f_{B\\{\ell\\}}\left(1-\frac{\alpha_{\ell}}{\hat{a}_{0}\left(f_{A\\{\ell\\}}\right)}f_{A\\{\ell\\}}\right)p_{A}\right)\left(\frac{\theta}{2}+\pi\right).$

Hence, the function $\left|f_{B\\{\ell+1\\}}(\theta)\right|^{2}$ can be
written as the sum of three terms:

$\begin{split}{\left|f_{B\\{\ell+1\\}}(\theta)\right|^{2}}&=\frac{1}{4}\left[\left|f_{B\\{\ell\\}}\left(\frac{\theta}{2}\right)\right|^{2}\left|1-\frac{\alpha_{\ell}}{\hat{a}_{0}\left(f_{A\\{\ell\\}}\right)}f_{A\\{\ell\\}}\left(\frac{\theta}{2}\right)\right|^{2}\left|p_{A}\left(\frac{\theta}{2}\right)\right|^{2}\left|p_{\hat{C}}\left(\frac{\theta}{2}\right)\right|^{2}+\right.\\\
&+\left|f_{B\\{\ell\\}}\left(\frac{\theta}{2}+\pi\right)\right|^{2}\left|1-\frac{\alpha_{\ell}}{\hat{a}_{0}\left(f_{A\\{\ell\\}}\right)}f_{A\\{\ell\\}}\left(\frac{\theta}{2}+\pi\right)\right|^{2}\left|p_{A}\left(\frac{\theta}{2}+\pi\right)\right|^{2}\left|p_{\hat{C}}\left(\frac{\theta}{2}+\pi\right)\right|^{2}+\\\
&+\left.2\mathfrak{Re}\left(\overline{f_{B\\{\ell\\}}\left(\frac{\theta}{2}\right)}f_{B\\{\ell\\}}\left(\frac{\theta}{2}+\pi\right)\left|1-\frac{\alpha_{\ell}}{\hat{a}_{0}\left(f_{A\\{\ell\\}}\right)}f_{A\\{\ell\\}}\left(\frac{\theta}{2}\right)\right|^{2}\overline{p_{\hat{C}}\left(\frac{\theta}{2}\right)}p_{\hat{C}}\left(\frac{\theta}{2}+\pi\right)\overline{p_{A}\left(\frac{\theta}{2}\right)}p_{A}\left(\frac{\theta}{2}+\pi\right)\right)\right].\end{split}$
(47)

By assumption we have $f_{B\\{\ell\\}}(0)=0$, $f_{B\\{\ell\\}}(\theta)\neq 0$
for all $\theta\in(0,2\pi)$, $|p_{A}|^{2}(\theta)+|p_{A}|^{2}(\theta+\pi)>0$,
$|p_{\hat{C}}|^{2}(\theta)+|p_{\hat{C}}|^{2}(\theta+\pi)>0$ and
$|p_{A}|^{2}(\theta)+|p_{\hat{C}}|^{2}(\theta+\pi)>0$ for all
$\,\theta\in[0,2\pi]$ and $p_{A}(\pi)=0$. These assumptions guarantee that for
all $\theta\in(0,2\pi)$ at least one term of the sum in (47) is different from
0. Hence, it is straightforward to check that $f_{B\\{\ell+1\\}}(0)=0$ and
$f_{B\\{\ell+1\\}}(\theta)>0$ for all $\theta\in(0,2\pi)$.

If we prove that the limit involving each one of the three terms in (47)
divided by $\left|f_{B\\{\ell\\}}(\theta)\right|^{2}$ is finite, then we can
exploit the algebraic properties of limits and conclude that the $\limsup$ in
(46) is finite. We have to focus on the following quantities:

  1. 1.

For the first term, the only quantity that should be checked is
$\frac{\left|f_{B\\{\ell\\}}\left(\frac{\theta}{2}\right)\right|^{2}}{\left|f_{B\\{\ell\\}}(\theta)\right|^{2}}.$
Indeed, the other terms are proportional to a positive quantity
$G_{1}(\theta)$ such that $0<c_{1}<G_{1}(\theta)<c_{2}<\infty$ in a
neighborhood of 0. Consequently, we can write

$\limsup_{\theta\rightarrow
0}\frac{\left|f_{B\\{\ell\\}}\left(\frac{\theta}{2}\right)\right|^{2}}{\left|f_{B\\{\ell\\}}(\theta)\right|^{2}}G_{1}(\theta)=G_{1}(0)<\infty.$

  2. 2.

For the second term, the quantity that should be checked is the ratio
$\limsup_{\theta\rightarrow
0}\frac{\left|p_{A}\left(\frac{\theta}{2}+\pi\right)\right|^{2}}{\left|f_{B\\{\ell\\}}(\theta)\right|^{2}}$,
since the other terms are proportional to a positive quantity $G_{2}(\theta)$
with $0<c_{1}<G_{2}(\theta)<c_{2}<\infty$ in a neighbourhood of 0.
Consequently, we can write

$\limsup_{\theta\rightarrow
0}\frac{\left|p_{A}\left(\frac{\theta}{2}+\pi\right)\right|^{2}}{\left|f_{B\\{\ell\\}}(\theta)\right|^{2}}G_{2}(\theta)=\limsup_{\theta\rightarrow
0}\frac{\left|p_{A}\left(\frac{\theta}{2}+\pi\right)\right|^{2}}{f_{A\\{\ell\\}}(\theta)}\frac{f_{A\\{\ell\\}}(\theta)}{\left|f_{B\\{\ell\\}}(\theta)\right|^{2}}G_{2}(\theta)$
(48)

and the finiteness of this $\limsup$ is implied by the hypotheses of Theorem
2.

  3. 3.

For the third term, we bound it from above with its modulus and then it
suffices to prove that

$\limsup_{\theta\rightarrow
0}\frac{\left|f_{B\\{\ell\\}}\left(\frac{\theta}{2}\right)\right|}{\left|f_{B\\{\ell\\}}(\theta)\right|}\frac{\left|p_{A}(\frac{\theta}{2}+\pi)\right|}{\left|f_{B\\{\ell\\}}(\theta)\right|}G_{3}(\theta)=G_{3}(0)\limsup_{\theta\rightarrow
0}\frac{\left|p_{A}(\frac{\theta}{2}+\pi)\right|}{\left|f_{B\\{\ell\\}}(\theta)\right|}$

with $0<c_{1}<G_{3}(\theta)<c_{2}<\infty$ in a neighbourhood of 0. The
previous limit is finite thanks to (48).

We conclude the proof by highlighting that the limit in Item 1. is different
from 0, which implies that the $\limsup$ in (46) is different from 0. ∎

A further result can be proven for the symbols $f_{A\\{\ell\\}}$,
$f_{{C}\\{\ell\\}}$ and $f_{\hat{C}\\{\ell\\}}$.

######  Lemma 9.

Consider the matrices
$\mathcal{A}_{2n_{\ell}}\\{\ell\\}\in\mathbb{C}^{2n_{\ell}\times 2n_{\ell}}$
defined by formulae (27)-(30). Define for all $\ell$ the grid transfer
operators

$\mathcal{P}\\{\ell\\}=\begin{bmatrix}P_{A{\\{\ell\\}}}&\\\
&P_{{\hat{C}\\{\ell\\}}}\end{bmatrix},\qquad
P_{A{\\{\ell\\}}}=\mathcal{C}_{n_{\ell}}(p_{A})K_{n_{\ell}}^{T},\qquad
P_{\hat{C}\\{\ell\\}}=\mathcal{C}_{n_{\ell}}(p_{\hat{C}})K_{n_{\ell}}^{T}.$

Suppose that $f_{A\\{\ell\\}}$, $f_{C\\{\ell\\}}$, $f_{B\\{\ell\\}}$,
$\alpha_{\ell}$, $p_{A}$ and $p_{\hat{C}}$ fulfil the hypotheses of Theorem 2
with $\theta_{0}=0$. Moreover, assume that

$|p_{A}|^{2}(\theta)+|p_{\hat{C}}|^{2}(\theta+\pi)>0,\qquad\forall\theta\in[0,2\pi].$

Then $f_{A\\{\ell+1\\}}(0)=0$, $f_{A\\{\ell+1\\}}(\theta)>0$ for all
$\theta\in(0,2\pi)$ and

$\limsup_{\theta\rightarrow
0}\frac{f_{A\\{\ell\\}}(\theta)}{f_{A\\{\ell+1\\}}(\theta)}=c,\quad
0<c<\infty.$ (49)

Moreover, the symbol $f_{{C}\\{\ell+1\\}}(\theta)$ associated with
${C}\\{\ell+1\\}$ is non-negative and

$\limsup_{\theta\rightarrow
0}\frac{f_{\hat{C}\\{\ell\\}}(\theta)}{f_{\hat{C}\\{\ell+1\\}}(\theta)}<\infty.$
(50)

###### Proof.

By Lemma 2, we have

$f_{A\\{\ell+1\\}}(\theta)=\frac{1}{2}\left(|p_{A}|^{2}f_{A\\{\ell\\}}\left(\frac{\theta}{2}\right)+|p_{A}|^{2}f_{A\\{\ell\\}}\left(\frac{\theta}{2}+\pi\right)\right).$

Moreover, by assumption we have $f_{A\\{\ell\\}}(0)=0$,
$f_{A\\{\ell\\}}(\theta)>0$ for all $\theta\in(0,2\pi)$,
$|p_{A}|^{2}(\theta)+|p_{A}|^{2}(\theta+\pi)>0$ for all $\,\theta\in[0,2\pi]$
and $p_{A}(\pi)=0$. Hence, as a consequence of Lemma 2, we have
$f_{A\\{\ell+1\\}}(0)=0$, $f_{A\\{\ell+1\\}}(\theta)>0$ for all
$\theta\in(0,2\pi)$, and

$\limsup_{\theta\rightarrow
0}\frac{f_{A\\{\ell\\}}(\theta)}{f_{A\\{\ell+1\\}}(\theta)}=c,\quad
0<c<\infty.$

Concerning $f_{\hat{C}\\{\ell+1\\}}$, we exploit the two definitions

$f_{\hat{C}\\{\ell\\}}(\theta)=f_{{C}\\{\ell\\}}(\theta)+\frac{\alpha_{\ell}\left|f_{{B}\\{\ell\\}}\right|^{2}(\theta)}{\hat{a}_{0}\left(f_{{A}\\{\ell\\}}\right)}\left(2-\frac{\alpha_{\ell}}{\hat{a}_{0}\left(f_{{A}\\{\ell\\}}\right)}f_{{A}\\{\ell\\}}(\theta)\right)=f_{{C}\\{\ell\\}}(\theta)+g_{1}(\theta)\left|f_{{B}\\{\ell\\}}(\theta)\right|^{2}$

and

$f_{\hat{C}\\{\ell+1\\}}(\theta)=f_{{C}\\{\ell+1\\}}(\theta)+\frac{\alpha_{\ell+1}\left|f_{{B}\\{\ell+1\\}}\right|^{2}(\theta)}{\hat{a}_{0}\left(f_{{A}\\{\ell+1\\}}\right)}\left(2-\frac{\alpha_{\ell+1}}{\hat{a}_{0}\left(f_{{A}\\{\ell+1\\}}\right)}f_{{A}\\{\ell+1\\}}(\theta)\right)=f_{{C}\\{\ell+1\\}}(\theta)+g_{2}(\theta)\left|f_{{B}\\{\ell+1\\}}(\theta)\right|^{2}.$
(51)

The functions $g_{1}(\theta)$ and $g_{2}(\theta)$ are such that
$0<c_{1}<g_{1}(\theta),g_{2}(\theta)<c_{2}<\infty$ in a neighborhood of 0
since $f_{{A}\\{\ell\\}}$ and $f_{{A}\\{\ell+1\\}}$ are non-negative and the
$\alpha_{\ell}$ are chosen such that the function that maps $\theta$ into
$2-\frac{\alpha_{\ell}}{\hat{a}_{0}\left(f_{{A}\\{\ell\\}}\right)}f_{{A}\\{\ell\\}}(\theta)$
is strictly positive for all $\ell$.

Note that $f_{{C}\\{\ell\\}}(\theta)\geq 0$ by assumption, which implies that
$f_{\hat{C}\\{\ell\\}}(\theta)$ is non-negative. From Lemma 2

$f_{{C}\\{\ell+1\\}}(\theta)=\frac{1}{2}\left(|p_{\hat{C}}|^{2}f_{\hat{C}\\{\ell\\}}\left(\frac{\theta}{2}\right)+|p_{\hat{C}}|^{2}f_{\hat{C}\\{\ell\\}}\left(\frac{\theta}{2}+\pi\right)\right)$

which is a non-negative trigonometric polynomial.

In order to prove the bound in (50) we exploit the expression of
$\hat{C}\\{\ell\\}$ and $\hat{C}\\{\ell+1\\}$. In particular,

$\displaystyle\limsup_{\theta\rightarrow
0}\frac{f_{\hat{C}\\{\ell\\}}(\theta)}{f_{\hat{C}\\{\ell+1\\}}(\theta)}$
$\displaystyle=\limsup_{\theta\rightarrow
0}\frac{f_{{C}\\{\ell\\}}(\theta)+g_{1}(\theta)\left|f_{{B}\\{\ell\\}}(\theta)\right|^{2}}{f_{{C}\\{\ell+1\\}}(\theta)+g_{2}(\theta)\left|f_{{B}\\{\ell+1\\}}(\theta)\right|^{2}}$
$\displaystyle=\limsup_{\theta\rightarrow
0}\frac{f_{{C}\\{\ell\\}}(\theta)}{f_{{C}\\{\ell+1\\}}(\theta)+g_{2}(\theta)\left|f_{{B}\\{\ell+1\\}}(\theta)\right|^{2}}+\limsup_{\theta\rightarrow
0}\frac{g_{1}(\theta)\left|f_{{B}\\{\ell\\}}(\theta)\right|^{2}}{f_{{C}\\{\ell+1\\}}(\theta)+g_{2}(\theta)\left|f_{{B}\\{\ell+1\\}}(\theta)\right|^{2}}.$

Now, we prove that both limits are finite. Concerning the first term, we have

$\displaystyle\limsup_{\theta\rightarrow
0}\frac{f_{{C}\\{\ell\\}}(\theta)}{f_{{C}\\{\ell+1\\}}(\theta)+g_{2}(\theta)\left|f_{{B}\\{\ell+1\\}}(\theta)\right|^{2}}$
$\displaystyle\leq\limsup_{\theta\rightarrow
0}\frac{f_{{C}\\{\ell\\}}(\theta)}{f_{{C}\\{\ell+1\\}}(\theta)}$
$\displaystyle\leq\limsup_{\theta\rightarrow
0}\frac{2f_{{C}\\{\ell\\}}(\theta)}{|p_{\hat{C}}|^{2}f_{{C}\\{\ell\\}}\left(\frac{\theta}{2}\right)+|p_{\hat{C}}|^{2}f_{{C}\\{\ell\\}}\left(\frac{\theta}{2}+\pi\right)+g_{3}(\theta)}<\infty$

where
$g_{3}(\theta)=|p_{\hat{C}}|^{2}g_{1}|f_{{B}\\{\ell\\}}|^{2}\left(\frac{\theta}{2}\right)+|p_{\hat{C}}|^{2}g_{1}|f_{{B}\\{\ell\\}}|^{2}\left(\frac{\theta}{2}+\pi\right)$
is a non-negative function and the latter bound is given by the result of
Lemma 2.

Finally, concerning the second term, the bound

$\displaystyle\limsup_{\theta\rightarrow
0}\frac{g_{1}(\theta)\left|f_{{B}\\{\ell\\}}(\theta)\right|^{2}}{f_{{C}\\{\ell+1\\}}(\theta)+g_{2}(\theta)\left|f_{{B}\\{\ell+1\\}}(\theta)\right|^{2}}\leq\limsup_{\theta\rightarrow
0}\frac{g_{1}(\theta)\left|f_{{B}\\{\ell\\}}(\theta)\right|^{2}}{g_{2}(\theta)\left|f_{{B}\\{\ell+1\\}}(\theta)\right|^{2}}<\infty$

follows by the fact that $f_{\hat{C}\\{\ell+1\\}}$ is non-negative and by
equation (46) in Lemma 4. ∎

## Acknowledgements

The work of Marco Donatelli, Paola Ferrari, Isabella Furci is partially
supported by Gruppo Nazionale per il Calcolo Scientifico (GNCS-INdAM).
Moreover, the work of Isabella Furci was also supported by the Young
Investigator Training Program 2020 (YITP 2019) promoted by ACRI.

## References

  * [1] D. J. Agress and P. Q. Guidotti.  The smooth extension embedding method.  SIAM J. Sci. Comput., 43(1):A446–A471, 2021.
  * [2] A. Aricò and M. Donatelli.  A v-cycle multigrid for multilevel matrix algebras: proof of optimality.  Numer. Math., 105(4):511–547, 2007.
  * [3] 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.
  * [4] M. Benzi, G. H. Golub, and J. Liesen.  Numerical solution of saddle point problems.  Acta Numer., 14:1–137, 2005.
  * [5] M. Bolten, E. De Sturler, and C. Hahn.  Krylov subspace recycling for evolving structures.  Comput. Methods Appl. Mech. Engrg., (accepted), 2021.
  * [6] 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., accepted, 2022.
  * [7] M. Bolten and C. Hahn.  Using composite finite elements for shape optimization with a stochastic objective functional.  In I. Faragó, F. Izsák, and P. L. Simon, editors, Progress in Industrial Mathematics at ECMI 2018, volume 30 of Mathematics in Industry, pages 515–520, Cham, 2019. Springer.
  * [8] E. Bozzo and C. Di Fiore.  On the use of certain matrix algebras associated with discrete trigonometric transforms in matrix displacement decomposition.  SIAM J. Matrix Anal. Appl., 16(1):312–326, 1995.
  * [9] D. Braess and R. Sarazin.  An efficient smoother for the Stokes problem.  volume 23, pages 3–19. 1997.  Multilevel methods (Oberwolfach, 1995).
  * [10] R. H. Chan, Q.-S. Chang, and H.-W. Sun.  Multigrid method for ill-conditioned symmetric toeplitz systems.  SIAM J. Sci. Comput., 19(2):516–529, 1998.
  * [11] M. Donatelli.  An algebraic generalization of local fourier analysis for grid transfer operators in multigrid based on toeplitz matrices.  Numer. Linear Algebra Appl., 17(2-3):179–197, 2010.
  * [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] G. Fiorentino and S. Serra.  Multigrid methods for Toeplitz matrices.  Calcolo, 28(3-4):283–305 (1992), 1991.
  * [14] 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.
  * [15] T. Huckle and J. Staudacher.  Multigrid preconditioning and Toeplitz matrices.  Electron. Trans. Numer. Anal., 13:81–105, 2002.
  * [16] T. Huckle and J. Staudacher.  Multigrid methods for block Toeplitz matrices with small size blocks.  BIT, 46(1):61–83, 2006.
  * [17] 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.
  * [18] J.-F. Maitre, F. Musy, and P. Nigon.  A fast solver for the Stokes equations using multigrid with a Uzawa smoother.  In Advances in multigrid methods (Oberwolfach, 1984), volume 11 of Notes Numer. Fluid Mech., pages 77–83. Friedr. Vieweg, Braunschweig, 1985.
  * [19] A. Napov and Y. Notay.  Smoothing factor, order of prolongation and actual multigrid convergence.  Numer. Math., 118(3):457–483, 2011.
  * [20] Y. Notay.  Algebraic theory of two-grid methods.  Numer. Math. Theory Methods Appl., 8(2):168–198, 2015.
  * [21] Y. Notay.  Algebraic two-level convergence theory for singular systems.  SIAM J. Matrix Anal. Appl., 37(4):1419–1439, 2016.
  * [22] Y. Notay.  A new algebraic multigrid approach for Stokes problems.  Numer. Math., 132(1):51–84, 2016.
  * [23] Y. Notay.  Algebraic multigrid for Stokes equations.  SIAM J. Sci. Comput., 39(5):S88–S111, 2017.
  * [24] 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.
  * [25] 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.
  * [26] S. Serra-Capizzano.  Generalized locally Toeplitz sequences: spectral analysis and applications to discretized partial differential equations.  Linear Algebra Appl., 366:371–402, 2003.  Special issue on structured matrices: analysis, algorithms and applications (Cortona, 2000).
  * [27] S. Serra-Capizzano.  The GLT class as a generalized Fourier analysis and applications.  Linear Algebra Appl., 419(1):180–233, 2006.
  * [28] S. Serra-Capizzano and C. Tablino-Possio.  Multigrid methods for multilevel circulant matrices.  SIAM J. Sci. Comput., 26(1):55–85, 2004.
  * [29] S. P. Vanka.  Block-implicit multigrid solution of Navier-Stokes equations in primitive variables.  J. Comput. Phys., 65(1):138–158, 1986.
  * [30] R. Wienands and W. Joppich.  Practical Fourier analysis for multigrid methods, volume 4 of Numerical Insights.  Chapman & Hall/CRC, Boca Raton, FL, 2005.  With 1 CD-ROM (Windows and UNIX).