remarkRemark hypothesisHypothesis claimClaim Stand-alone Multigrid for
HelmholtzVandana Dwarka, Cornelis Vuik

# Stand-alone multigrid for Helmholtz revisited: towards convergence using
standard components

Vandana Dwarka Department of Applied Mathematics, Delft University of
Technology, Delft, the Netherlands<EMAIL_ADDRESS>Cornelis Vuik
Department of Applied Mathematics, Delft University of Technology, Delft, the
Netherlands (, homepage.tudelft.nl/d2b4e/<EMAIL_ADDRESS>
###### Abstract

Getting standard multigrid to work efficiently for the high-frequency
Helmholtz equation has been an open problem in applied mathematics for years.
Much effort has been dedicated to finding solution methods which can use
multigrid components to obtain solvers with a linear time complexity. In this
work we present one among the first stand-alone multigrid solvers for the 2D
Helmholtz equation using both a constant and non-constant wavenumber model
problem. We use standard smoothing techniques and do not impose any
restrictions on the number of grid points per wavelength on the coarse-grid.
As a result we are able to obtain a full V- and W-cycle algorithm. The key
features of the algorithm are the use of higher-order inter-grid transfer
operators combined with a complex constant in the coarsening process. Using
weighted-Jacobi smoothing, we obtain a solver which is $h-$independent and
scales linearly with the wavenumber $k$. Numerical results using 1 to 5
GMRES(3) smoothing steps approach $k-$ and $h-$ independent convergence, when
combined with the higher-order inter-grid transfer operators and a small or
even zero complex shift. The proposed algorithm provides an important step
towards the perpetuating branch of research in finding scalable solvers for
challenging wave propagation problems.

###### keywords:

Helmholtz Equation, Multigrid, Jacobi, Indefinite Matrix, CSLP, Convergence

65N55, 65F10, 65F15

##  1 Introduction

For over three decades, researchers have studied ways to design convergent and
efficient multigrid schemes for the indefinite Helmholtz equation, which to
this day remains an open problem in applied mathematics. The Helmholtz
equation finds its application in many engineering practices, ranging from
geoscience to plasma physics and quantum mechanics. Solving the Helmholtz
equation numerically is a challenging task. Multigrid methods are widely known
to be efficient and scalable numerical solvers for elliptic PDEs. It provided
a state-of-the-art numerical approach, where the amount of computational work
scales with the number of unknowns. In essence, the method combines solutions
on coarser grids with relaxation (smoothing) techniques to arrive at a
solution on the fine grid.

While multigrid works well for the positive-definite variant of the Helmholtz
equation or when the wavenumber is small, the method generally fails to
converge as a stand-alone solver for the indefinite version. This break-down
can be attributed to two things. First of all, the coarse-grid needs to be
fine enough in order to resolve the waves. A difference between the wavenumber
on the fine vs. coarse-grid(s) can cause a phase error leading to
instabilities [17]. Moreover, eigenvalues approaching zero also hamper the
convergence as they lead to near-singularity within the coarse-grid operator
[17]. Several strategies have been proposed to study the convergence behavior
of such indefinite systems [2, 1, 34, 20]. This boils down to prove
convergence of the two-grid method for the positive definite variant (see
[20], chapter 11.6). Another approach lies in trying to find upper and lower
bounds to the field of values of the two-grid operator [28, 26, 27, 29]. This
requires that the Hermitian part of the linear system
$\frac{1}{2}\left(A+A^{T}\right)$ is positive definite, which is not the case
for the indefinite Helmholtz equation. While we can still use the condition
that the spectral radius of the two-grid iteration matrix should be smaller
than one, the computation of the eigenvalues becomes difficult in case of
Sommerfeld boundary conditions. Moreover, it has been noted that for nonnormal
(complex) systems in particular, the Helmholtz equation being an example of
such, assessing convergence solely based on the spectral radii can be
misleading as it only rigorously describes asymptotic convergence [29]. Bounds
in terms of norms would be more appropriate. However, these often go hand in
hand with assumptions which are difficult to check and contain unknown
constants, which reduces their practical applicability.

For positive-definite and complex Hermitian problems, Geometric and/or
Algebraic Multigrid are powerful solvers which provide solutions in linear
time complexity by constructing coarser linear systems based on the
information provided by the fine-grid linear system. For complex Hermitian
matrices, the coarse linear systems can be build using the Galerkin condition,
which has been studied in [22, 25]. If the use of a simple smoother by means
of a stationary iterative method such as weighted-Jacobi or Gauss-Seidel is
preferred, the coarse grid correction scheme should balance the simplicity of
the smoother by accurately transferring smooth errors from the coarse grid to
the fine grid [25]. When the linear system is complex symmetric, linear
interpolation can be used to construct the prolongation operator. The
restriction operator can then be chosen as the transpose of the prolongation
operator. In fact, for real matrices optimality studies have been conducted in
order to construct the interpolation and prolongation operators [4, 5, 19, 3].
In [32], a similar optimality condition from [19, 3] has been extended to the
Hermitian case. These studies show that the optimal prolongation and
restriction operator are based on the eigenvectors of the system matrix and
several approximations are given which could function as a satisfactory proxy
in the case of symmetric positive definite (SPD) and Hermitian systems
respectively.

Including a complex shift, leads to the Complex Shifted Laplacian (CSL), which
can be solved using multigrid as a stand-alone solver. Consequently, for the
original indefinite Helmholtz problem, multigrid applied to the CSL can be
used as a preconditioner,  together with an iterative Krylov subspace method,
such as GMRES. Here, the original indefinite Helmholtz problem is solved by
approximately inverting the CSL using geometric multigrid [16, 15, 14, 8].
While this works well for homogeneous problems, it is not suitable for
heterogeneous problems with high contrasts [33, 15].

When it comes to using multigrid as a stand-alone solver for the original
Helmholtz operator, several studies have focused on getting multigrid to
comply with the underlying characteristics of the indefinite nonsymmetric
system. One study uses a smoothed aggregation approach, where an effective
prolongation operator is constructed by minimizing the energy norm defined
using the normal equations $A^{*}A$. This is achieved by running the CGNR
algorithm to solve the equation $A^{*}AP=0$, where $P$ is the tentative
prolongation operator [30, 31]. Two other methods are the wave-ray approach
[24] and the first-order least squares (FOLS) approach [23]. Both compute
multiple coarse spaces to approximate the plane waves in the near null space.
In both works, the development of more efficient approaches for such linear
systems has been stressed. Moreover, for variable wave numbers, it remains
challenging to choose an effective prolongation and restriction operator,
especially if the aim is to keep the algorithm as light as possible.

Another difficulty arises in picking a suitable smoother, as for Helmholtz
problems, the (weighted)-Jacobi method diverges. While the latter has natural
parallelization properties and is often preferred for SPD problems, no clear
remedy exists to apply it suitably for Helmholtz problems at large wavenumbers
and tailored adaptions are needed [21, 18, 16]. Essentially, this relates to
the indefiniteness of the underlying operators, leading to negative
eigenvalues which ensure that the spectral radius of the weighted-Jacobi
iteration matrix is always greater than 1. One important recent study finds
that replacing the weighted-Jacobi smoother with a few GMRES iterations could
lead to improved convergence [13]. In [17, 18], it was shown that for the one-
dimensional indefinite Helmholtz equation, convergence can be achieved by
considering a two-step Jacobi smoother. Scheduling the smoothing steps on all
levels as a function of the wave number provides a uniform error reduction in
the high-frequency modes. More recently a two-dimensional dispersion
correction has been developed for the original Helmholtz operator when using a
constant wave number. Here, a larger number of smoothing steps combined with
at least 4.5 grid points per wavelength at the coarsest level leads to a
convergent multigrid scheme. [9].

A full multigrid hierarchy, without any restrictions on the grid resolution at
the coarsest level, is given in [10]. Here, a re-discretization technique is
used where a level-dependent complex shift is incorporated, except the finest
level. The authors show that the altered two-grid correction scheme, combined
with 3 GMRES iterations as a smoother, remains robust and does not lead to
agitation of smooth modes. The method works well on structured grids and in
cases where the behavior of the wave number on the fine and coarse(r) grids
can be defined a-priori. Complex geometries, high-contrast and irregularly
varying wave number would pose difficulties in constructing the coarser linear
systems along the multilevel hierarchy, when re-discretization is used
combined with a level-dependent parameter.

In this work we aim to solve the complex indefinite Helmholtz problem using
constant and variable wave numbers with a stand-alone multigrid method using
Galerkin coarsening. For the baseline convergence study, we impose two
requirements: the method should converge with a small number of smoothing
steps and the intergrid transfer operators should be level-independent and
easy to construct and implement.
Our contribution to this long standing field of research can be summarized
into the following points:

  * •

To the best of our knowledge, we are the first to present a fully convergent
multigrid solver for the two-dimensional indefinite nonsymmetric Helmholtz
equation using classical multigrid components for elliptic problems (V- and
W-cycle hierarchy, damped Jacobi smoothing and no restriction on the
(coarsest) levels.

  * •

We provide a novel generalized proof to assess two-grid convergence for any
(highly) indefinite system. Our proof is not limited to the Helmholtz equation
and does not require the calculation of eigenvalues.

  * •

The novelty in our theory also lies in the fact that we do not require
assumptions on normality of the operator.

  * •

For Helmholtz problems in particular, we provide a optimality condition which
can easily be verified without computing eigenvalues.

  * •

The results are robust to large wavenumbers and non-constant wavenumbers
representing high-contrast heterogeneities.

In particular, we show that the use of coarsening on the CSL combined with a
higher-order prolongation and restriction operator leads to a convergent
scheme, even when using the classical divergent weighted-Jacobi smoother.

Several remarks are in place. First of all, in this work we solely consider
getting a convergent solver rather than a solver with wavenumber independent
convergence. This is a topic for future research due to the sheer difficulty
in getting a convergent stand-alone multigrid solver in the first place.
Second of all, given that multigrid generally already diverges for one-
dimensional Helmholtz problems with small wavenumbers, this research will, in
accordance with the literature, focus on one- and two-dimensional model
problems. The study of three-dimensional problems is a topic for future
research as it requires a thorough understanding of the convergence behavior
in one- and two-dimensions first.

The structure of this paper is as follows. In Section 2 we introduce the model
problems which will be the focus of this work. In Section 3 we give an
introduction to the multigrid scheme and its parameters, to continue with the
analysis of the convergence of multigrid methods for Helmholtz problems in
Section 4. Numerical experiments are presented in Section 5. Our conclusions
are summarized in Section 6.

##  2 Problem Description

We continue by defining the model problems which will be studied in this
paper. We focus on the 2D constant wavenumber model using $k>0$, which we call
MP 2-A.

$\displaystyle-\Delta u(x,y)-k^{2}u(x,y)\,u(x,y)$ $\displaystyle=$
$\displaystyle\delta{(x-\frac{1}{2},y-\frac{1}{2})},\quad{(x,y)}\in\Omega=[0,1]^{2},$
$\displaystyle\left(\frac{\partial{}}{\partial{{\mathbf{n}}}}-ik\right)u({x,y})$
$\displaystyle=$ $\displaystyle 0,\quad(x,y)\in\partial{\Omega}.$

The final test problem uses a non-constant wavenumber $k(x,y)$. This gives us
MP 2-B as defined below

$\displaystyle-\Delta u(x,y)-{k(x,y)}^{2}u(x,y)\,u(x,y)$ $\displaystyle=$
$\displaystyle\delta{(x-\frac{1}{2},y-\frac{1}{2})},\quad{(x,y)}\in\Omega=[0,1]^{2},$
$\displaystyle\left(\frac{\partial{}}{\partial{{\mathbf{n}}}}-ik(x,y)\right)u({x,y})$
$\displaystyle=$ $\displaystyle 0,\quad(x,y)\in\partial{\Omega},$

where $\chi(x,y)$ is a random real function in the range of $[0,1]$ and
$k_{1},k_{2}$ are positive real numbers. Note that in this case $k$ varies
exactly between $k_{1}$ and $k_{2}$. We use second-order finite differences to
discretize the model problems and define the step-size $h=\frac{1}{n}$, where
$n$ is chosen such that for each $k_{1},k_{2}$ we have $kh=0.625$. This is
equivalent to having 10 grid points per wavelength.

##  3 Multigrid Methods

Multigrid method as a stand-alone solver combines two ingredients in a
recursive way: a coarse grid correction and a smoother. If we let $P$ and
$P^{\prime}$ denote the matrix variant of the prolongation and restriction
operators, then we define the coarse grid correction matrix $CGC$ as

$\displaystyle CGC=I-P^{\prime}{A_{c}}^{-1}PA,\quad A_{c}=PAP^{\prime}.$

Note that we here denote the coarse grid version of $A$ by $A_{c}$ instead of
$A_{H}$ as the matrix $A$ is now complex due to the Sommerfeld boundary
conditions and we reserve the $H-$script to denote the conjugate transpose.

As for the smoother, we work with the standard weighted-Jacobi smoother, which
is defined as

$\displaystyle S_{JAC}=I-X^{-1}A,\quad X=\omega\Lambda_{A},$

where $\Lambda_{A}$ denotes the matrix containing the diagonal elements of our
matrix $A$. In general, for symmetric elliptic problems such as the Poisson
equation, the optimal relaxation parameter is set at $\omega=\frac{3}{2}$. For
the indefinite nonsymmetric Helmholtz problem, more tailored values of
relaxation parameters should be used, especially outside the context of using
multigrid as a preconditioner [18, 21].

To develop the theory, we start with the assumption that we are only using a
post-smoother. In this case, the two-grid iteration matrix is given by

(1) $\displaystyle
T_{0}=\left(I-P^{\prime}{A_{c}}^{-1}PA\right)\left(I-X^{-1}A\right).$

Generally the multigrid iteration scheme converges if the spectral radius of
the two-grid iteration matrix is less than 1. Thus, the straightforward way to
prove convergence is by showing that the spectral radius of the two-grid
iteration matrix is strictly less than zero. If Dirichlet conditions are used,
this naturally boils down to using Rigorous Fourier Analysis (RFA) to find
analytical expressions for the eigenvalues of $T_{0}$. However, once we use
Sommerfeld boundary conditions, this is no longer an option and we have to
find another way to assess the boundedness of the spectral radius,  which will
be the focus of the next section, where we develop a general theoretical
framework to assess the convergence, irrespective of the boundary conditions
used.

###  3.1 Convergence

For convergence of the two-grid method we require that
$\left\|T_{0}\right\|_{2}<1$ independent of $h$. We can write $T_{0}$ as
$T_{0}=I-DA$, where $D$ represents an approximation to $A^{-1}$. This is shown
in Lemma 3.1. We use this to show that if the two-grid operator $T_{0}$ can be
written in this form and the two-norm is bounded by one, then the multigrid
method will converge.

######  Lemma 3.1.

Convergence - I Let $A\in\mathbb{C}^{n\times n}$ be a non-singular and non-
defective matrix, and let the two-grid operator be given by
$T_{0}=\left(I-P^{\prime}{A_{c}}^{-1}PA\right)\left(I-X^{-1}A\right)$. Then,
$T_{0}$ can be written as $T_{0}=I-DA$, with

(2) $\displaystyle
D=X^{-1}+P^{\prime}{A_{c}}^{-1}P-P^{\prime}{A_{c}}^{-1}PX^{-1}$

Moreover we have ${T_{0}}^{H}=I-A^{H}D^{H}.$

######  Proof 3.2.

Expanding $T_{0}$, and factoring in terms of $A$ we have

(3) $\displaystyle T_{0}$
$\displaystyle=I-X^{-1}A-P^{\prime}{A_{c}}^{-1}PA+P^{\prime}{A_{c}}^{-1}PX^{-1}A,$
(4)
$\displaystyle=I-\left(X^{-1}+P^{\prime}{A_{c}}^{-1}P-P^{\prime}{A_{c}}^{-1}PX^{-1}\right)A,$
(5) $\displaystyle=I-DA,$

where $D=X^{-1}+P^{\prime}{A_{c}}^{-1}P-P^{\prime}{A_{c}}^{-1}PX^{-1}$. We
similarly obtain

(6) $\displaystyle{T_{0}}^{H}$
$\displaystyle={\left(\left(I-P^{\prime}{A_{c}}^{-1}PA\right)\left(I-X^{-1}A\right)\right)}^{H},$
(7)
$\displaystyle={\left(I-X^{-1}A\right)}^{H}{\left(I-P^{\prime}{A_{c}}^{-1}PA\right)}^{H},$
(8)
$\displaystyle=\left(I-A^{H}X^{-1}\right)\left(I-A^{H}P^{\prime}{\left({A_{c}}^{H}\right)}^{-1}P\right),$
(9)
$\displaystyle=I-A^{H}X^{-1}-A^{H}P^{\prime}{\left({A_{c}}^{H}\right)}^{-1}P+A^{H}X^{-1}A^{H}P^{\prime}{\left({A_{c}}^{H}\right)}^{-1}P,$
(10)
$\displaystyle=I-A^{H}\left(X^{-1}+P^{\prime}{\left({A_{c}}^{H}\right)}^{-1}P-X^{-1}A^{H}P^{\prime}{\left({A_{c}}^{H}\right)}^{-1}P\right),$
(11) $\displaystyle=I-A^{H}D^{H},$

where
$D^{H}=X^{-1}+P^{\prime}{\left({A_{c}}^{H}\right)}^{-1}P-X^{-1}P^{\prime}{\left({A_{c}}^{H}\right)}^{-1}P$.

It has been shown that convergence of the two-grid operator implies multigrid
convergence, see [20]. The proof is based on an induction argument, where for
the $W-$cycle, it is shown that the multigrid iteration matrix can be written
as the two-grid iteration matrix plus a perturbation term. The perturbation
term is bounded if the prolongation operator is a bounded linear operator and
the smoother can be bounded by a constant. Note that a convergent smoothing
scheme is not required. In fact, the smoothing property also holds for non-
convergent iterations [20]. It is only required that the error is smoothed up
to a certain $\nu$, where $\nu$ is the number of smoothing steps. In other
words, there exists a $\nu$ such that the error is reduced. It is not
imperative that the smoothing property holds for $\nu=\infty$. We now work
towards conditions for which $\left\|T_{0}\right\|_{2}<1$ holds.

######  Theorem 3.3 (Hermitian Positive Definiteness (HPD) implies
Convergence).

Let $A\in\mathbb{C}^{n\times n}$ be a non-singular and non-defective matrix.
Let the two-grid operator be given by
$T_{0}=\left(I-P^{\prime}{A_{c}}^{-1}PA\right)\left(I-X^{-1}A\right)$, where
$X$ is the $\omega-$Jacobi smoother, $P,P^{\prime}\in\mathbb{R}^{n\times m}$
are the prolongation and restriction operator and $A_{c}=P^{\prime}AP$. Let

(12) $\displaystyle\Gamma=A^{H}D^{H}+DA-A^{H}D^{H}DA,$

then ${T_{0}}^{H}{T_{0}}=I-\Gamma$. If $\Gamma$ is positive definite, then
${\lambda_{j}}\left(\Gamma\right)\in(0,2)$ and $\left\|T_{0}\right\|_{2}<1$
independent of $h$.

######  Proof 3.4.

Using equation Lemma 3.1, equation Eq. 5 and Eq. 11, we can write

(13) $\displaystyle{T_{0}}^{H}{T_{0}}$
$\displaystyle=\left(I-A^{H}D^{H}\right)\left(I-DA\right),$ (14)
$\displaystyle=I-\left(DA+A^{H}D^{H}-A^{H}D^{H}DA\right),$ (15)
$\displaystyle=I-\Gamma.$

We can show a stronger result than just positive definiteness since $\Gamma$
is also Hermitian. To see this, note

(16) $\displaystyle{\Gamma}^{H}$
$\displaystyle={\left(A^{H}D^{H}+DA-A^{H}D^{H}DA\right)}^{H},$ (17)
$\displaystyle={\left(DA+A^{H}D^{H}\right)}^{H}-{\left(A^{H}D^{H}DA\right)}^{H},$
(18)
$\displaystyle={\left(DA\right)}^{H}+{\left(A^{H}D^{H}\right)}^{H}-{\left(A^{H}D^{H}DA\right)}^{H},$
(19) $\displaystyle=A^{H}D^{H}+DA-A^{H}D^{H}DA=\Gamma.$

Given that $\Gamma$ is HPD, the singular value decomposition (SVD) and
eigendecomposition coincide, so we can find a unitary matrix
$U\in\mathbb{C}^{n\times n}$ such that $\Gamma=U\Sigma_{\Gamma}U^{H}$, where
$\Sigma$ is the diagonal matrix containing the singular resp. eigenvalues
$(\sigma_{j})$, where we have
$\sigma_{1}\geq\sigma_{2}\geq\sigma_{3}\ldots\geq\sigma_{n}>0$. We have

(20) $\displaystyle\Gamma$ $\displaystyle=A^{H}D^{H}+DA-A^{H}D^{H}DA,$

which is equivalent to

(21) $\displaystyle\Sigma_{\Gamma}$
$\displaystyle=UA^{H}D^{H}U^{H}+UDAU^{H}-UA^{H}D^{H}DAU^{H}.$

By our assumption on $\Gamma$ being positive definite, we have that for all
$x\in\mathbb{C}^{n}\setminus\\{0\\}$,

(22) $\displaystyle 0<x^{H}\Gamma x$ $\displaystyle\Leftrightarrow
0<x^{H}A^{H}D^{H}x+x^{H}DAx-x^{H}A^{H}D^{H}DAx.$

Similarly, by using Eq. 21, for all $x\in\mathbb{C}^{n}\setminus\\{0\\}$,

(23) $\displaystyle 0<x^{H}\Sigma_{\Gamma}x$ $\displaystyle\Leftrightarrow
0<x^{H}UA^{H}D^{H}U^{H}x+x^{H}UDAU^{H}x-x^{H}UA^{H}D^{H}DAU^{H}x,$ (24)
$\displaystyle\Leftrightarrow
0<x^{H}UA^{H}D^{H}DAU^{H}x<x^{H}UA^{H}D^{H}U^{H}x+x^{H}UDAU^{H}x.$

Using the unitary invariance of the SVD, Eq. 24 leads to

(25) $\displaystyle 0$
$\displaystyle<x^{H}UA^{H}D^{H}DAU^{H}x<x^{H}UA^{H}D^{H}U^{H}x+x^{H}UDAU^{H}x,$
(26) $\displaystyle\Rightarrow
0<{\sigma_{j}(DA)}^{2}<{\sigma_{j}(A^{H}D^{H}+DA)}\leq 2{\sigma_{j}(DA)},$
(27) $\displaystyle\Rightarrow\sigma_{j}(DA)<2.$

This gives an upperbound for $\sigma_{j}(DA)$. We now develop the upper bound
for the largest singular value of $\Gamma$, which is denoted by
$\sigma_{1}(\Gamma)$. We know that $0\leq\sigma_{j}(DA)<2$. The case where
$0\leq\sigma_{j}(DA)<1$ is trivial as that automatically renders
$\left\|T_{0}\right\|_{2}<1$. We focus on the case where
$1\leq\sigma_{j}(DA)<2$. Recall, we assumed $\Gamma$ is HPD, but $DA$ does not
necessarily have to be. As a result, the singular value decomposition for
$\Gamma$ with unitary matrix $U$ might not coincide with the singular value
decomposition for $DA$.

Suppose, we can write $DA$ as $DA=W\Gamma_{DA}V^{H}$, where $WV^{H}\neq I$. We
then still have

(28) $\displaystyle A^{H}D^{H}DA$
$\displaystyle=V{\Gamma_{DA}}^{H}W^{H}W\Gamma_{DA}V^{H},$ (29) $\displaystyle
V^{H}A^{H}D^{H}DAV$
$\displaystyle={\Gamma_{DA}}^{H}\Gamma_{DA}={\left(\Gamma_{DA}\right)}^{2}.$

Moreover, we also have for $v_{j}\in V$

(30) $\displaystyle DAv_{j}$ $\displaystyle=\sigma_{j}(DA)w_{j}\mbox{ and }$
(31) $\displaystyle\left(DAv_{j}\right)^{H}$
$\displaystyle={v_{j}}^{H}A^{H}D^{H}={w_{j}}^{H}\sigma_{j}(DA).$

For $WV^{H}\neq I$ we will always have that for some index $j$ and $v_{j}\in
V$ and $w_{j}\in W$

(32) $\displaystyle{v_{j}}^{H}w_{j}+{w_{j}}^{H}v_{j}<2,$

where we used that the complex parts cancel out. We therefore have

(33) $\displaystyle
v^{H}DAv+v^{H}A^{H}D^{H}v={v_{j}}^{H}\sigma_{j}(DA)w_{j}+{w_{j}}^{H}\sigma_{j}(DA)v_{j}<2\sigma_{j}(DA).$

Taking $x=v_{j}$ in Eq. 22 and using Eq. 33, $x^{H}\Gamma x$ can be bounded by

(34) $\displaystyle 0$ $\displaystyle<{v_{j}}^{H}\Gamma
v_{j}={v_{j}}^{H}A^{H}D^{H}v_{j}+{v_{j}}^{H}DAv_{j}-{v_{j}}^{H}A^{H}D^{H}DAv_{j},$
(35)
$\displaystyle=\sigma_{j}(DA)\left({v_{j}}^{H}w_{j}+{w_{j}}^{H}v_{j}\right)-{\sigma_{j}(DA)}^{2},$
(36)
$\displaystyle<2\sigma_{j}(DA)-{\sigma_{j}(DA)}^{2}=\sigma_{j}(DA)\left(2-\sigma_{j}(DA)\right),$
(37) $\displaystyle<\sigma_{j}(DA)<2.$

As $\Gamma$ is positive definite, the eigenvalues and singular values coincide
and we have

(38) $\displaystyle\lambda_{\min}(\Gamma)>0\mbox{ and
}\lambda_{\max}(\Gamma)<2.$

We are now ready to bound the two-grid operator $T_{0}$ using the 2-norm over
$\mathbb{C}^{n\times n}$

$\displaystyle\|{T_{0}}\|_{2}$
$\displaystyle=\sup_{x\in\mathbb{C}^{n},\|x\|_{2}=1}{\|\left({T_{0}}^{H}{T_{0}}\right)x\|_{2}}=\sqrt{\rho\left({T_{0}}^{H}{T_{0}}\right)},$
(39)
$\displaystyle=\sqrt{\rho\left(I-\left(A^{H}D^{H}+DA-A^{H}D^{H}DA\right)\right)},$
(40) $\displaystyle=\sqrt{\max_{1\leq j\leq
n}|\lambda_{j}\left(I-\left(A^{H}D^{H}+DA-A^{H}D^{H}DA\right)\right)|},$ (41)
$\displaystyle\leq\sqrt{|1-\min_{1\leq j\leq
n}\lambda_{j}\left(A^{H}D^{H}+DA-A^{H}D^{H}DA\right)|},$ $\displaystyle<1,$

where we used that trivially
$0<\min\lambda_{j}\left(A^{H}D^{H}+DA-A^{H}D^{H}DA\right)$ by positive
definiteness of $\Gamma$. As regards, the largest eigenvalue
$\max\lambda_{j}\left(\Gamma\right)<2$, we obtain

$\displaystyle\|{T_{0}}\|_{2}$
$\displaystyle=\sup_{x\in\mathbb{C}^{n},\|x\|_{2}=1}{\|\left({T_{0}}^{H}{T_{0}}\right)x\|_{2}}=\sqrt{\rho\left({T_{0}}^{H}{T_{0}}\right)},$
(42)
$\displaystyle=\sqrt{\rho\left(I-\left(A^{H}D^{H}+DA-A^{H}D^{H}DA\right)\right)},$
(43) $\displaystyle=\sqrt{\max_{1\leq j\leq
n}|\lambda_{j}\left(I-\left(A^{H}D^{H}+DA-A^{H}D^{H}DA\right)\right)|},$ (44)
$\displaystyle\leq\sqrt{|1-\max_{1\leq j\leq
n}\lambda_{j}\left(A^{H}D^{H}+DA-A^{H}D^{H}DA\right)|},$ (45)
$\displaystyle<\sqrt{|1-2|}=1,\mbox{ using<EMAIL_ADDRESS>
Thus, if $\Gamma$ is positive definite then $\lambda_{j}(\Gamma)\in(0,2)$ and
consequently $\lambda_{j}(T_{0})\in(0,1)$.

Theorem 3.3 shows that if $\Gamma$ is positive definite, then the eigenvalues
of $T_{0}$ will lie in the interval $(0,1)$ which automatically leads to a
bound for the 2-norm of the two-grid operator. Note that this is independent
of $h$, as no constant appears in Eq. 38 which depends on $h$. Now that we
know that a sufficient conditions for convergence is to have $\Gamma$ be
positive definite, we need a way to practically assess this. One approach is
by simplifying the expression for $\Gamma$ and finding another sufficient
condition for positive definiteness, which is more easy to analyze. This will
be constructed in Corollary 3.6. We again start by defining the simplified
operator in Definition 3.5.

######  Definition 3.5 (Simplified operator).

Let $A\in\mathbb{C}^{n\times n}$ be a non-singular and non-defective matrix.
$P,P^{\prime}\in\mathbb{R}^{n\times m}$ are real valued prolongation and
restriction operators and $A_{c}=P^{\prime}AP$. The $\omega-$Jacobi smoother
is given by $X=\omega\Lambda_{A}$, where $\Lambda$ denotes the diagonal matrix
containing the diagonal entries of $A$. Then we define
$\tilde{D}=X^{-1}+P^{\prime}{A_{c}}^{-1}P$ such that
$\tilde{\Gamma}=A^{H}\tilde{D}^{H}+\tilde{D}A-A^{H}\tilde{D}^{H}\tilde{D}A$.

Note that in the definition given above, $\tilde{D}$ does not contain the term
$P^{\prime}{A_{c}}^{-1}PX^{-1}A$, which is part of the original operator $D$.

######  Corollary 3.6.

Convergence - II Let $A\in\mathbb{C}^{n\times n}$ be a non-singular and non-
defective matrix. Let $\tilde{\Gamma}$ be defined as in Definition 3.5. If
$\tilde{\Gamma}$ is positive definite, then $\Gamma$ from Theorem 3.3 is
positive definite.

######  Proof 3.7.

To see this, we write $\Gamma$ as

(46) $\displaystyle\Gamma$ $\displaystyle=A^{H}D^{H}+DA-A^{H}D^{H}DA,$ (47)
$\displaystyle=\left(A^{H}\tilde{D}^{H}+\tilde{D}A-A^{H}\tilde{D}^{H}\tilde{D}A\right)+{A^{H}}{X^{-1}P^{\prime}{{\left({A_{c}}^{H}\right)}^{-1}}P}{P^{\prime}{A_{c}}^{-1}PX^{-1}}A$
(48)
$\displaystyle-A^{H}\tilde{D}^{H}P^{\prime}{A_{c}}^{-1}PX^{-1}A-{A^{H}}X^{-1}P^{\prime}{{\left({A_{c}}^{H}\right)}^{-1}}P\tilde{D}A.$

We assume $\tilde{\Gamma}$ is positive definite and HPD respectively. Now,
assuming the opposite, i.e. $\Gamma$ is not HPD implies that $\exists
x\in\mathbb{C}^{n}\setminus\ \\{0\\}$ such that

$\displaystyle
x^{H}\tilde{\Gamma}x+{x^{H}}{A^{H}}{X^{-1}P^{\prime}{{\left({A_{c}}^{H}\right)}^{-1}}P}{P^{\prime}{A_{c}}^{-1}PX^{-1}}Ax\hskip
5.69054pt(\ast)$
$\displaystyle-{x^{H}}A^{H}\tilde{D}^{H}P^{\prime}{A_{c}}^{-1}PX^{-1}Ax-{x^{H}}{A^{H}}X^{-1}P^{\prime}{{\left({A_{c}}^{H}\right)}^{-1}}P\tilde{D}Ax\leq
0.\hskip 5.69054pt(\ast\ast)$

By assumption and the fact that the second term in $(\ast)$ contains a
quadratic form, $(\ast)$ is always positive for any
$x\in\mathbb{C}^{n}\setminus\ \\{0\\}$. But if $\Gamma$ is not positive
definite then using $(\ast\ast)$ we must also have

(49) $\displaystyle
0<x^{H}\tilde{\Gamma}x+{x^{H}}{A^{H}}{X^{-1}P^{\prime}{{\left({A_{c}}^{H}\right)}^{-1}}P}{P^{\prime}{A_{c}}^{-1}PX^{-1}}Ax,$
(50)
$\displaystyle\leq{x^{H}}A^{H}\tilde{D}^{H}P^{\prime}{A_{c}}^{-1}PX^{-1}Ax+{x^{H}}{A^{H}}X^{-1}P^{\prime}{{\left({A_{c}}^{H}\right)}^{-1}}P\tilde{D}Ax$

We show that $P^{\prime}{A_{c}}^{-1}PX^{-1}A$ has singular values equal to
zero. Observe that the rank of $P^{\prime}{A_{c}}^{-1}PX^{-1}A$ is
$\frac{n}{2}$. Suppose we can write the singular value decomposition as
$\hat{\Gamma}=W^{H}\left(P^{\prime}{A_{c}}^{-1}PX^{-1}A\right)V$. In this case
we have $V$ and $W$ to represent two unitary complex matrices as we have not
assumed that $P^{\prime}{A_{c}}^{-1}PX^{-1}A$ is HPD. We know that
$P^{\prime}{A_{c}}^{-1}P$ is singular and half of all singular values are
equal to zero. We thus need to show that $P^{\prime}{A_{c}}^{-1}PX^{-1}A$ has
zero singular values as well. In order to do this, we use an analogous version
of Ostrowski’s theorem. Here instead of using eigenvalues, we use singular
values. According to the theorem, for $K\in\mathbb{C}^{n\times n}$ and
$Y\in\mathbb{C}^{n\times n}$ with $Y$ is non-singular, we have

(51) $\displaystyle\sigma_{j}(KY)\leq\sigma_{1}(Y)\sigma_{j}(A),\quad
j=1,2,\ldots n,$

where $\sigma_{1}\geq\sigma_{2}\geq\ldots\sigma_{n}$. Applying Eq. 51 to our
case by taking $K=P^{\prime}{A_{c}}^{-1}P$ and $Y=X^{-1}A$, we obtain

(52)
$\displaystyle\sigma_{j}\left(P^{\prime}{A_{c}}^{-1}PX^{-1}A\right)\leq\sigma_{1}(X^{-1}A)\sigma_{j}(P^{\prime}{A_{c}}^{-1}P),\quad
j=1,2,\ldots n.$

We know that $\sigma_{1}(X^{-1}A)>0$ given that both $X$ and $A$ are non-
singular. However, we know that $P^{\prime}{A_{c}}^{-1}P$ is singular and has
zero singular values. Thus there exists indices $\tilde{j}$ such that

(53)
$\displaystyle\sigma_{\tilde{j}}\left(P^{\prime}{A_{c}}^{-1}PX^{-1}A\right)\leq
0,\quad\tilde{j}\in\\{1,2,\ldots n\\}.$

Trivially $0\leq\sigma_{\tilde{j}}\left(P^{\prime}{A_{c}}^{-1}PX^{-1}A\right)$
is a lower bound for the smallest singular value. Combining this with Eq. 53
leads to

(54) $\displaystyle
0\leq\sigma_{\tilde{j}}\left(P^{\prime}{A_{c}}^{-1}PX^{-1}A\right)\leq
0,\quad\exists\tilde{j}\in\\{1,2,\ldots n\\},$

and so by the interlacing theorem, we must have that there exist at least one
singular value at an index $\tilde{j}$ such that
$\sigma_{\tilde{j}}\left(P^{\prime}{A_{c}}^{-1}PX^{-1}A\right)=0$. Naturally,
the eigenvalue of $P^{\prime}{A_{c}}^{-1}PX^{-1}A$ corresponding to the index
$\tilde{j}$ must be zero as well, given that the singular value is zero.

We know that for any complex matrix, there exists a unitary matrix such that
we can write it in the Schur decomposition. Without the assumption of
normality, it solves a general eigenvalue problem of the form $Kv=\lambda Bv$
for some matrix $B$.

We take $x=v_{\tilde{j}}$, corresponding to the zero eigenvalue, which we
denote by $\lambda_{\tilde{j}}^{0}$ such that

(55) $\displaystyle
P^{\prime}{A_{c}}^{-1}PX^{-1}Av_{\tilde{j}}=\lambda_{\tilde{j}}(P^{\prime}{A_{c}}^{-1}PX^{-1}A)Bv_{\tilde{j}}=0Bv_{\tilde{j}}.$

Substituting Eq. 55 into Eq. 50 gives

(56) $\displaystyle 0$
$\displaystyle<{v_{\tilde{j}}}^{H}\tilde{\Gamma}v_{\tilde{j}}+{{v_{\tilde{j}}}^{H}}{A^{H}}{X^{-1}P^{\prime}{{\left({A_{c}}^{H}\right)}^{-1}}P}{P^{\prime}{A_{c}}^{-1}PX^{-1}}Av_{\tilde{j}},$
(57)
$\displaystyle\leq{{v_{\tilde{j}}}^{H}}A^{H}\tilde{D}^{H}P^{\prime}{A_{c}}^{-1}PX^{-1}Av_{\tilde{j}}+{{v_{\tilde{j}}}^{H}}{A^{H}}X^{-1}P^{\prime}{{\left({A_{c}}^{H}\right)}^{-1}}P\tilde{D}Av_{\tilde{j}},$
(58)
$\displaystyle=\lambda_{\tilde{j}}^{0}{{v_{\tilde{j}}}^{H}}A^{H}\tilde{D}^{H}B{v_{\tilde{j}}}+\lambda_{\tilde{j}}^{0}{{v_{\tilde{j}}}^{H}}B^{H}\tilde{D}A{v_{\tilde{j}}}=0.$

This gives a contradiction as $\tilde{\Gamma}$ is assumed to be positive
definite and must always be larger than zero for any
$x\in\mathbb{C}^{n}\setminus\ \\{0\\}$. We thus have that if $\tilde{\Gamma}$
is positive definite, then $x^{H}\Gamma x>0\mbox{ }\forall
x\in\mathbb{C}^{n}\setminus\ \\{0\\}$, i.e. $\Gamma$ must be positive definite
as well.

Using the above, we now have that it is sufficient to check if
$\tilde{\Gamma}$ is positive definite. We therefore test for convergence by
verifying whether $\tilde{\Gamma}$ is positive definite. We use the Cholesky
decomposition to determine whether $\tilde{\Gamma}$ is positive definite, as
any HPD system matrix $B$ can be written as $B=LL^{\ast}$, where
$L\in\mathbb{C}^{n\times n}$ is a lower-triangular matrix. Another method to
check for positive definiteness without having to compute the Cholesky factors
is by checking all of the following

  1. 1.

$b_{ii}>0\forall i$

  2. 2.

$b_{ii}+b_{jj}>2\left|\Re\left[b_{ij}\right]\right|\mbox{ for }i\neq j$

  3. 3.

The element with the largest modulus lies on the diagonal

  4. 4.

$\det(B)>0$

The results from this section show that if we have an HPD ${\Gamma}$ matrix,
then the two-grid method converges. This can also be tested for a simplified
operator $\tilde{\Gamma}$. Note that so far, all theoretical considerations
are for general matrices. Next, we will apply this theory to the Helmholtz
operator.

##  4 Convergence for indefinite nonsymmetric Helmholtz operators

The previous section has stipulated properties which we can check in order to
find a convergent configuration for the indefinite nonsymmetric Helmholtz
equation. For this we need to consider optimality conditions for the smoother
and coarse grid correction scheme.

Starting with the smoother, we will the use classical $\omega-$Jacobi smoother
due to its natural parallelizability. Even though the smoother diverges for
Helmholtz problems, we will show that combined with CSL-coarsening and higher-
order prolongation and restriction, the $\tilde{\Gamma}$ and consequently
$\Gamma$ operator becomes HPD, leading to a convergent solver on a complete
V-cycle hiearchy. In the remainder of this paper, unless stated otherwise, we
will use aggregate values from the literature presented for the two-step
Jacobi as our $\omega$ in the weighted-Jacobi method, which boils down to
keeping $\omega=4.5$. Our own experimental results in Section 4.2 are in line
with the aggregate values from the literature, confirming relaxation
parameters which are generally lower than the underrelaxation parameters
reported for symmetric elliptic problems, such as Poisson’s equation.

We proceed by evaluating these properties for the Helmholtz operator. We use
MP2-A with constant wavenumber $k$. We numerically investigate the changes we
can make to the design of the coarse-grid system in order to achieve
convergence. In this work, we consider two options; we either use higher-order
interpolation schemes to construct the prolongation and restriction operator
and/or we use the CSL, which will be denoted by $C$, as the system on which we
apply the coarsening operations. We use the higher-order interpolation scheme
based on the quadratic rational Beziér curve from [11, 12]. There, the scheme
was used to construct the deflation vectors for a two-level and multi-level
deflation preconditioner, which showed wavenumber independent convergence for
the two-level scheme and close to wavenumber independent convergence for the
multi-level scheme. The main motivation for studying these type of higher-
order schemes is that they have been shown to reduce the projection error
[11]. Using this scheme, the prolongation operator acts on a grid function as
follows

(59)
$P\left[u_{2h}\right]_{i}=\left\\{\begin{matrix}\begin{matrix}\frac{1}{8}\left(\left[u_{2h}\right]_{\left(i-2\right)/2}+6\left[u_{2h}\right]_{\left(i\right)/2}+\left[u_{2h}\right]_{\left(i+2\right)/2}\right)&\mbox{if
$i$ is even,}\\\
\frac{1}{2}\left(\left[u_{2h}\right]_{\left(i-1\right)/2}+\left[u_{2h}\right]_{\left(i+1\right)/2}\right)&\mbox{if
$i$ is odd,}\end{matrix}\end{matrix}\right\\},$

for $i=1,\dots,N$ and for $i=1,\dots,\frac{N}{2}$, where $N$ denotes the size
of the fine-level linear system. Note that the second line ($i$ = odd) is in
fact the linear interpolation scheme. Thus, we now have a combination of both
a higher-order and linear interpolation scheme for the nodes $i$ within the
numerical domain.

The inclusion of the complex shift was first applied to multigrid schemes in
[10]. There, the main motivation relies on spectral analysis of the two-grid
operator and the observation that the coarse-grid eigenvalues can never
approach zero due to the complex part. We proceed by a different yet general
theoretical argument, which is independent of the boundary conditions.
Moreover, we include a constant complex shift compared to a level-dependent
one and also evaluate the effect of a classic smoother compared to the
polynomial GMRES(3) smoother.

In Table 1, we show that the use of the CSL in combination with a higher-order
interpolation scheme in fact leads to an HPD system. Note that while our
conditions do not require the actual computation of the spectral norm, we will
do so for sake of illustration and completeness.

| Linear | Beziér
---|---|---
$k$ | $A$ | $C$ | $A$ | $C$
$5$ | $\times{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>| $\times{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>| $\checkmark{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>| $\checkmark{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>
$10$ | $\times{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>| $\times{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>| $\times{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>| $\checkmark{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>
$20$ | $\times{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>| $\times{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>| $\times{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>| $\checkmark{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>
$30$ | $\times{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>| $\times{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>| $\times{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>| $\checkmark{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>

Table 1: At the right of each entry, the spectral norm of the two-grid
operator is given. Linear uses linear interpolation to construct
$P,P^{\prime}$. Beziér uses rational quadratic Beziér interpolation. $A$
represents $A_{c}=P^{\prime}AP$. $C$ represents $A_{c}=P^{\prime}CP$, where
$C$ denotes the CSL with complex shift $\beta_{2}=0.7$. In all cases, one
post-smoothing step is used. Left of each entry, $\checkmark$ denotes that
$\tilde{\Gamma}$ is HPD and $\times$ denotes it is not.

Several interesting observations can be made. First of all, in almost call
cases when using the classical smoothing operator, a combination of both
higher-order interpolation and CSL for coarsening is needed to obtain an HPD
matrix and consequently a spectral norm which is bounded by one. Thus, using
either level independent CSL-coarsening or higher-order prolongatrion and
restriction, will not solve the problem.

###  4.1 Optimality

Apart from having a condition to check for convergence, we would like to have
a measure of how fast we can expect the two-grid method to converge. By having
a positive definite $\Gamma$, we know that the spectral norm is bounded from
above by one. To obtain a more qualitative estimate, we study the condition
number of $\Gamma$ and $\tilde{\Gamma}$. Using the condition numbers we obtain
a sharper bound in Theorem 3.3. This entails that the smallest eigenvalue of
$\tilde{\Gamma}$ will be a lower bound to the smallest eigenvalue of $\Gamma$.

######  Corollary 4.1.

Convergence Optimality Condition - I Let $A\in\mathbb{C}^{n\times n}$ be a
non-singular and non-defective matrix and let
$\tilde{D}=X^{-1}+P^{\prime}{A_{c}}^{-1}P$. If
$\tilde{\Gamma}=A^{H}\tilde{D}^{H}+\tilde{D}A-A^{H}D^{H}DA$ is positive
definite and ${\|{{\Gamma}}\|}\leq{\|\tilde{\Gamma}\|}$ together with
$\kappa(\tilde{\Gamma})\leq\kappa(\Gamma)$ then
$\lambda_{\min}(\tilde{\Gamma})$ is given and bounded by

(60)
$\displaystyle\frac{{\|{\tilde{\Gamma}}\|}}{\kappa(\tilde{\Gamma})}\leq\lambda_{\min}({\Gamma}).$

######  Proof 4.2.

Note that $\tilde{\Gamma}$ and $\Gamma$ are HPD and so all eigenvalues are
real, so we can consider ${\|{{\Gamma}}^{-1}\|}$ and
${\|{\tilde{\Gamma}}^{-1}\|}$. Moreover, $\forall\mbox{ matrices
}B\in\mathbb{C}^{n\times n},\rho(B)=\max_{1\leq j\leq
n}|\lambda_{j}(B)|\leq{\|B\|}$ in any associative norm. If
$\kappa(\tilde{\Gamma})\leq\kappa(\Gamma)$, then we obtain

(61)
$\displaystyle{\|{{\Gamma}}\|}\leq{\|\tilde{\Gamma}\|}\Rightarrow\frac{1}{{\|{\tilde{\Gamma}}\|}}\leq\frac{1}{{\|{\Gamma}\|}}\Rightarrow\frac{\kappa(\tilde{\Gamma})}{{\|{\tilde{\Gamma}}\|}}\leq\frac{\kappa(\Gamma)}{{\|{\Gamma}\|}}\Rightarrow\|\tilde{\Gamma}^{-1}\|\leq\|{\Gamma}^{-1}\|.$

Using that the systems are HPD, we have

(62) $\displaystyle\rho(\Gamma^{-1})=\max_{1\leq j\leq
n}|\lambda_{j}(\Gamma^{-1})|=\min_{1\leq j\leq
n}{|\lambda_{j}(\Gamma)|}^{-1}\leq\|{\Gamma}^{-1}\|.$

Note that this gives the reciprocal of the smallest eigenvalue of $\Gamma$. To
get the actual smallest eigenvalue of $\Gamma$, we have to take the
reciprocal, which provides us with inequality Eq. 60 and a lower bound to
$\lambda_{\min}({\Gamma})$.

We can use Corollary 4.1 to construct a sharper bound in Theorem 3.3.

######  Corollary 4.3.

Convergence Optimality Condition - II Let $A\in\mathbb{C}^{n\times n}$ be a
non-singular and non-defective indefinite matrix. Let $T_{0}$ be such that we
can write $T_{0}$ as $T_{0}=I-DA$ with
$D=X^{-1}+P^{\prime}{A_{c}}^{-1}P-P^{\prime}{A_{c}}^{-1}PX^{-1}$. Let
$\Gamma=A^{H}{D}^{H}+DA-A^{H}{D}^{H}{D}A$. If
$\tilde{\Gamma}=A^{H}\tilde{D}^{H}+\tilde{D}A-A^{H}\tilde{D}^{H}\tilde{D}A$
with $\tilde{D}=X^{-1}+P^{\prime}{A_{c}}^{-1}P$, is positive definite and
${\|\tilde{\Gamma}\|}\leq{\|{{\Gamma}}\|}$ together with
$\kappa(\Gamma)\leq\kappa(\tilde{\Gamma})$, then

$\displaystyle\|{T_{0}}\|_{2}<\sqrt{\left|1-\frac{\|\tilde{\Gamma}\|_{\text{}}}{\kappa_{\text{}}(\tilde{\Gamma})}\right|}<1.$

######  Proof 4.4.

Using Theorem 3.3 we can write ${T_{0}}^{H}T_{0}=I-\Gamma$. Using Eq. 41 from
Theorem 3.3 and substituting Eq. 60 from Corollary 4.3 have

$\displaystyle\|{T_{0}}\|_{2}$
$\displaystyle=\sup_{x\in\mathbb{C}^{n},\|x\|_{2}=1}{\|\left({T_{0}}^{H}{T_{0}}\right)x\|_{2}}=\sqrt{\rho\left({T_{0}}^{H}{T_{0}}\right)},$
(63) $\displaystyle=\sqrt{|1-\min\lambda_{j}\left(\Gamma\right)|},$ (64)
$\displaystyle<\sqrt{|1-\min\lambda_{j}\left(\tilde{\Gamma}\right)|},$ (65)
$\displaystyle<\sqrt{\left|1-\frac{\|\tilde{\Gamma}\|_{\text{}}}{\kappa_{\text{}}(\tilde{\Gamma})}\right|}<1.$

###  4.2 Optimal convergence for Helmholtz

From Section 4, we have observed that we need a clever combination of
coarsening using the CSL, higher-order interpolation schemes and the right
amount of smoothing steps in order to obtain HPD systems. The results from
Table 1, however, indicate that convergence, in terms of total number of
iterations, will be slow as the spectral norm is still close to one.

| $\omega=1.5$ | $\omega=2$ | $\omega=2.5$ | $\omega=4.5$ | $\omega=7$
---|---|---|---|---|---
$k$ | $\nu=1$ | $\nu=2$ | $\nu=1$ | $\nu=2$ | $\nu=1$ | $\nu=2$ | $\nu=1$ | $\nu=2$ | $\nu=1$ | $\nu=2$
$5$ | $0.373$ | $0.413$ | $0.273$ | $0.405$ | $0.206$ | $0.389$ | $0.088$ | $0.200$ | $0.031$ | $0.123$
$10$ | $0.137$ | $0.140$ | $0.128$ | $0.139$ | $0.112$ | $0.137$ | $0.065$ | $0.116$ | $0.028$ | $0.081$
$20$ | $0.030$ | $0.028$ | $0.029$ | $0.029$ | $0.028$ | $0.030$ | $0.022$ | $0.028$ | $0.012$ | $0.025$
$30$ | $0.011$ | $0.009$ | $0.011$ | $0.010$ | $0.010$ | $0.011$ | $0.008$ | $0.011$ | $0.001$ | $0.009$

Table 2: We report the value of
$\|\tilde{\Gamma}\|\kappa(\tilde{\Gamma})^{-1}$ in the $p=1$ norm, with
$\tilde{\Gamma}$ from Corollary 4.1. $\nu$ denotes the number of
$\omega$-Jacobi smoothing steps.

The results from Table 2 reveal that a lower $\omega$ is in favor when we use
one post-smoothing steps and the wavenumber is small. For $\omega=1.5$, we
observe that increasing $k$ and $\nu$ leads to an decrease in
$\|\tilde{\Gamma}\|\kappa(\tilde{\Gamma})^{-1}$. This means that as we perform
more smoothing steps, the two-grid spectral norm moves closer and closer to
one. For convergence, we require the opposite: the addition of a few more
smoothing steps should lead to an increase rather than a decrease in
$\|\tilde{\Gamma}\|\kappa(\tilde{\Gamma})^{-1}$. Based on these results, the
optimal estimate for $\omega$ lies between $[2,4.5]$ for the reported values
of $k$.  These values are different from the general consensus on the optimal
relaxation parameter for poisson-type of problems. While we focus on the
classical smoother in this work, our reported values do seem to coincide with
underrelaxation parameters observed in other studies where a two-step weighted
Jacobi smoother is analyzed [21, 18]

##  5 Numerical Experiments

In this section we examine the convergence behavior of the multigrid solver.
We use the model problems from Section 2. Unless stated otherwise, we use 10
grid points per wavelength, which is equivalent to $kh=0.625$. All experiments
are implemented sequentially on a Dell laptop using 8GB RAM and a i7-8665U
processor. In all experiments, we set the relative tolerance to $10^{-5}$. We
allow for coarsening until the dimension of the underlying linear system
$N_{c}$ is smaller than 10. The maximum size of the linear system on the
coarsest grid is therefore $10\times 10$. We use the $\omega-$Jacobi smoother
with the approximate optimal relaxation parameter from Section 4.2, i.e.
$\omega=4.5$ as we want to test for large $k$. Moreover, we perform the
coarse-grid corrections using the higher-order prolongation and restriction
operator, together with the CSL. For the latter, the complex shift is set at
$\beta_{2}=0.7$ unless stated otherwise.

###  5.1 2D Constant $k$

For the constant wavenumber problem (MP 2-A), we first start by confirming
that we have $h-$independent convergence. Recall that the constant from
Theorem 3.3 does not depend on $h$, which for the proposed setup should lead
to $h-$independent convergence. Thus, while the convergence bounds can be
shown to be independent of $h$, the practical convergence speed may still
depend on $h$.

####  5.1.1 $h-$independence

We report the results for various $k$ using $h=2^{-p}$, with $p=5$ up to
$p=9$.

| $k=15$ | $k=30$
---|---|---
$h$ | $\nu=1$ | $\nu=2$ | $\nu=4$ | $\nu=1$ | $\nu=2$ | $\nu=4$
$2^{-5}$ | 45 | 24 | 18 | $\odot$ | $\odot$ | $\odot$
$2^{-6}$ | 34 | 22 | 18 | 66 | 37 | 28
$2^{-7}$ | 36 | 22 | 18 | 52 | 33 | 27
$2^{-8}$ | 40 | 24 | 18 | 54 | 34 | 27
$2^{-9}$ | 42 | 23 | 18 | 58 | 36 | 27
Table 3: Number of V-cycles for $k=15$ and $k=30$. $\nu$ denotes the number of
$\omega$-Jacobi smoothing steps using $\omega=4.5$ $\odot$ denotes the case
where $kh>0.625$, which is excluded.

We observe that the solver indeed performs independently of $h$. In fact, it
appears that with 4 smoothing steps, the solver converges in about
$\tilde{C}k$ iterations, where $\tilde{C}\approx 1$ is a constant. The latter
holds irrespective of the problem size.

####  5.1.2 General results

Now that we have established, both theoretically and numerically, that the
convergence is $h-$independent once we take 4 smoothing steps, we revert back
to letting $kh=0.625$. Note that in order to obtain improved accuracy, $kh$
can be decreased without affecting the convergence behavior due to the
$h-$independence. In Table 4 we report the results for MP 2-A using a constant
wavenumber and Sommerfeld boundary conditions.

Table 4: Number of V- $(\gamma=1)$ and W-cycles $(\gamma=2)$ for constant $k$
(MP 2-A) using tol. $10^{-5}$. $\nu$ denotes the number of $\omega$-Jacobi
smoothing steps. $N_{D}$ is the size of the coarsest system.

| $k=50$ | $k=100$ | $k=150$ | $k=200$ | $k=250$
---|---|---|---|---|---
| $N=6724$ | $N=26244$ | $N=57600$ | $N=102400$ | $N=160000$
| $N_{D}=8$ | $N_{D}=8$ | $N_{D}=4$ | $N_{D}=8$ | $N_{D}=4$
$\gamma$ | $1$ | $2$ | $1$ | $2$ | $1$ | $2$ | $1$ | $2$ | $1$ | $2$
$\nu=4$ | $58$ | $58$ | $104$ | $108$ | $155$ | $159$ | $209$ | $213$ | $267$ | $271$
$\nu=5$ | $58$ | $58$ | $104$ | $104$ | $150$ | $166$ | $194$ | $229$ | $238$ | $287$
$\nu=6$ | $55$ | $58$ | $99$ | $102$ | $139$ | $167$ | $183$ | $222$ | $226$ | $283$
$\nu=7$ | $53$ | $60$ | $97$ | $101$ | $136$ | $163$ | $179$ | $219$ | $221$ | $280$
$\nu=8$ | $53$ | $60$ | $95$ | $104$ | $131$ | $161$ | $178$ | $212$ | $218$ | $277$

From Table 4 we observe that for the V-cycle the convergence improves if we
use more smoothing steps. Note that these additional smoothing steps are
computationally cheap, as we use the $\omega-$Jacobi smoother. In general, we
observe that we again need approximately $\tilde{C}k$ iterations with
$\tilde{C}\approx 1$, in order to reach convergence. While the convergence
behavior is promising, the results do provide some further insights into the
behavior of the smoother. For example, we observe that moving from a V-cycle
to a W-cycle does not seem to improve the performance. In fact, we observe
that we need more iterations instead of less. One potential explanation could
be that we know that in general the $\omega-$Jacobi smoother diverges for
Helmholtz-type of problems. While this in itself does not have to lead to
divergence (the smoothing property is still satisfied for a divergent
smoother, see [20]), the use of the W-cycle requires more smoothing iterations
than the V-cycle. At some point we expect the efficiency of these smoothing
steps to reduce, which could be reflected in the higher number of iterations.

Figure 1: Real part of the 2D solution for MP2-A $k=100$ (left) and $k=150$
(right)

###  5.2 2D results non-constant $k(x,y)$

In this subsection we report on the numerical results for the non-constant
wavenumber model problem. We distinguish two cases; a medium-varying and a
highly-varying wavenumber profile problem. In this respect the profile
contains either smooth or sharp changes between a range of wavenumbers. For an
illustration of the wavenumber profile, see Fig. 2 for the smooth changing
wavenumber profile and Fig. 3 for the profile where the wavenumber changes
sharply. In all cases we make sure that the largest value of $k(x,y)$ still
obeys the rule of thumb $kh\approx 0.625$.

####  5.2.1 Smooth Changes

Fig. 2 illustrates how the wavenumber varies between 10 and 75 for the medium-
varying problem. Note that the transition from a low to high wavenumber goes
gradually.

Figure 2: Left: real part of the 2D numerical solution. Right: $k(x,y)$

Table 5 contains the number of V-cycles needed to reach convergence for
$k(x,y)$ varying between $(10,50)$ and $(10,75)$ respectively. We again
observe that the number of iterations scales linearly with $k$. The iteration
count follows approximately $\tilde{C}k$, with $\tilde{C}\approx 1$, where $k$
in this case denotes the largest wavenumber. In this event, the W-cycle does
lead to a smaller number of iterations, compared to the case where we have a
constant wavenumber.

Table 5: Number of V- $(\gamma=1)$ and W-cycles $(\gamma=2)$ for MP 2-B
(medium variation). $\nu$ denotes the number of $\omega$-Jacobi smoothing
steps.

| $(k_{1},k_{2})=(10,50)$ | $(k_{1},k_{2})=(10,75)$
---|---|---
$\gamma$ | $1$ | $2$ | $1$ | $2$
$\nu=4$ | $65$ | $60$ | $90$ | $88$
$\nu=5$ | $62$ | $59$ | $86$ | $86$
$\nu=6$ | $61$ | $58$ | $85$ | $85$
$\nu=7$ | $60$ | $57$ | $84$ | $84$
$\nu=8$ | $59$ | $57$ | $83$ | $83$

####  5.2.2 Sharp Changes

In this subsection we report on the results for the highly-varying profile
model problem, containing sharp changes between the wavenumbers $k=10$ and
$k=75$. Fig. 3 illustrates how the wavenumber varies between 10 and 75 for the
medium-varying problem. In this case, we let $k(x,y)$ vary randomly, which
gives a wavenumber which varies highly across the entire numerical domain.

Figure 3: Left: real part of the 2D numerical solution. Right: $k(x,y)$

Table 6 contains the number of V-cycles needed to reach convergence for
$k(x,y)$ varying between $(10,50)$ and $(10,75)$ respectively. We immediately
observe that the number of iterations goes up by a factor of approximately 1.5
once we allow for the sharp and rapid changes in the wavenumber. For example
for the medium-varying case containing more smooth transition between the
wavenumbers, using 8 smoothing steps, we had 59 and 57 iterations respectively
for the V- and W-cycle. In this case, we have 94 iterations for both the V-
and W-cycle. Another interesting observation is that for the W-cycle, $k(x,y)$
varying between 10 and 75 leads to a larger number of iterations when $\nu>4$.
This could again be an indication of the smoother losing some of its efficacy
when we allow for more smoothing operations using the W-cycle. Also for this
model problem, we conclude that as the wavenumber grows, the V-cycle is
preferred over the W-cycle for this choice of smoothing scheme. In the
reported cases, the method converges in approximately $\leavevmode\nobreak\
1.5k$ iterations, where $k$ is the largest admissible wavenumber.

| $(k_{1},k_{2})=(10,50)$ | $(k_{1},k_{2})=(10,75)$
---|---|---
$\gamma$ | $1$ | $2$ | $1$ | $2$
$\nu=4$ | $102$ | $96$ | $111$ | $107$
$\nu=5$ | $97$ | $95$ | $103$ | $105$
$\nu=6$ | $95$ | $95$ | $101$ | $104$
$\nu=7$ | $94$ | $94$ | $102$ | $104$
$\nu=8$ | $94$ | $94$ | $102$ | $104$

Table 6: Number of V- $(\gamma=1)$ and W-cycles $(\gamma=2)$ for MP 2-B (high
variation). $\nu$ denotes the number of $\omega$-Jacobi smoothing steps.

###  5.3 GMRES-smoothing

In this section we explore how the solver performs when we use GMRES(3)
smoothing. In [10], where the authors propose coarsening on the CSL instead of
the original Helmholtz operator, GMRES(3) is used as the only smoother. The
motivation for this lies in the fact that indeed the $\omega-$Jacobi diverges
when combined with the original and standard configuration of the multigrid
method. By using a few GMRES iterations, the smoother acts as a polynomial
smoother. Already in [13], the use of GMRES(3) was proposed as a smoother.
There, the authors manually optimized the smoothing schedule, where the number
of smoothing steps ranges between 13 and 39 smoothing steps for a constant $k$
ranging from $4\pi$ to $32\pi$. More generally in the literature, on each
level where the Jacobi smoother becomes unstable, the smoother is replaced by
the Krylov iterations ranging between 5-40 iterations per smoothing step [13,
6, 7]. In this work we let the number of smoothing steps using GMRES(3) vary
between 1 and 5. One advantage of this approach is that by allowing for only 3
GMRES iterations per smoothing step, the cost of applying GMRES does not
increase significantly with respect to memory and computation time. We perform
the same numerical experiments using the model problems mentioned previously.

####  5.3.1 2D results constant $k$

In this subsection we report the results for the constant wavenumber model
problem (MP 2-A). As we are coarsening using the CSL, we want to distinguish
between the cases where the complex shift is large versus small. We therefore
start by keeping the shift $\beta_{2}=0.7$, similar to the one used in all
previous experiments. Results are reported in Table 7 for 1 up to 5 smoothing
steps. We immediately observe that the number of iterations for both the V-
and W-cycle are drastically reduced. If we compare these results to the ones
obtained when using the $\omega-$Jacobi smoother in Table 4, the reduction in
the number of iterations is approximately 2.5 times. For example, when using 5
smoothing steps within a V-cycle, the $\omega-$Jacobi smoother requires 238
iterations to reach convergence. When GMRES(3) is used as a smoother with a
similar amount of smoothing steps, 88 iterations are required to reach
convergence. Note that there appears to be no difference in the number of
iterations when using the V-cycle and W-cycle.

Table 7: Number of V- $(\gamma=1)$ and W-cycles $(\gamma=2)$ for constant $k$
(MP 2-A) using tol. $10^{-5}$. $\nu$ denotes the number of GMRES(3) smoothing
steps with $\beta_{2}=0.7$

| $k=50$ | $k=100$ | $k=150$ | $k=200$ | $k=250$
---|---|---|---|---|---
| $N=6724$ | $N=26244$ | $N=57600$ | $N=102400$ | $N=160000$
| $N_{D}=8$ | $N_{D}=8$ | $N_{D}=4$ | $N_{D}=8$ | $N_{D}=4$
$\gamma$ | $1$ | $2$ | $1$ | $2$ | $1$ | $2$ | $1$ | $2$ | $1$ | $2$
$\nu=1$ | $37$ | $36$ | $68$ | $67$ | $99$ | $98$ | $132$ | $131$ | $162$ | $161$
$\nu=2$ | $29$ | $29$ | $53$ | $53$ | $78$ | $78$ | $104$ | $104$ | $128$ | $128$
$\nu=3$ | $24$ | $24$ | $45$ | $45$ | $67$ | $67$ | $89$ | $89$ | $112$ | $112$
$\nu=4$ | $22$ | $22$ | $40$ | $40$ | $59$ | $59$ | $78$ | $78$ | $98$ | $98$
$\nu=5$ | $20$ | $20$ | $36$ | $36$ | $53$ | $53$ | $71$ | $71$ | $88$ | $88$

Next, we repeat the same numerical experiments, however this time we set the
complex shift to $\beta_{2}=k^{-1}$. Note that this is a very small shift, in
fact the CSL matrix starts being a close resemblance of the original Helmholtz
operator. Results are reported in Table 8. For all cases, we report that there
is a drastic improvement in the number of iterations. We obtain a solver which
is close to $k-$ independent when we use $\nu>2$ and $\gamma=2$. Unlike the
previous case where $\beta_{2}=0.7$ in Table 7, we now do observe a lower
number of iterations for the W-cycle.

In fact, additional numerical experiments confirm that the results are similar
if we let $\beta_{2}=0$, which results in the original Helmholtz operator.
Thus, when using GMRES(3) as a smoother, the original Helmholtz operator can
in fact be used instead of the CSL for coarsening within the multigrid
hierarchy. However, the fast convergence is only observed when combined with
the high-order prolongation and restriction operators.

Given that we obtained the best numerical results with $\beta_{2}=k^{-1}$
without any additional costs, we continue with this shift for the CSL in the
upcoming sections.

| $k=50$ | $k=100$ | $k=150$ | $k=200$ | $k=250$
---|---|---|---|---|---
| $N=6724$ | $N=26244$ | $N=57600$ | $N=102400$ | $N=160000$
| $N_{D}=8$ | $N_{D}=8$ | $N_{D}=4$ | $N_{D}=8$ | $N_{D}=4$
$\gamma$ | $1$ | $2$ | $1$ | $2$ | $1$ | $2$ | $1$ | $2$ | $1$ | $2$
$\nu=1$ | $14$ | $7$ | $24$ | $10$ | $39$ | $19$ | $51$ | $24$ | $64$ | $29$
$\nu=2$ | $8$ | $5$ | $13$ | $7$ | $22$ | $10$ | $28$ | $13$ | $34$ | $16$
$\nu=3$ | $6$ | $5$ | $10$ | $6$ | $16$ | $9$ | $20$ | $10$ | $24$ | $12$
$\nu=4$ | $6$ | $5$ | $8$ | $5$ | $12$ | $7$ | $15$ | $9$ | $18$ | $10$
$\nu=5$ | $5$ | $5$ | $7$ | $5$ | $11$ | $7$ | $13$ | $8$ | $15$ | $9$

Table 8: Number of V- $(\gamma=1)$ and W-cycles $(\gamma=2)$ for constant $k$
(MP 2-A) using tol. $10^{-5}$. $\nu$ denotes the number of GMRES(3) smoothing
steps with $\beta_{2}=k^{-1}$.

####  5.3.2 2D results non-constant $k(x,y)$

Finally, we investigate the convergence behavior using the GMRES(3) smoother
on the highly-varying problem (MP 2-B2). In Section 5.2, we observed that this
was the hardest problem to solve in terms of heterogeneous problem, due to the
wavenumber varying highly between $k=10$ and $k=75$ across the numerical
domain. Thus, we only test for this case in this subsection.

####  5.3.3 Sharp Changes

Results are reported in Table 9 and again indicate a drastic improvement
compared to the case where we use the $\omega-$Jacobi smoother. In fact, even
for this highly-varying model containing sharp changes between the wavenumbers
$k=10$ and $k=75$, where the wavenumber is allowed to vary randomly between 10
and 75 across the domain, we obtain a $k-$ and $h-$independent multigrid
solver for $\nu\geq 3$. Similar to the previous case, even when using the
original Helmholtz operator, we reach convergence when using GMRES(3) as a
smoother. However, the $k-$independence is only observed when the higher-order
prolongation and restriction scheme is used. These results are promising and
provide a solid framework for future research, given that the current industry
standard (CSL preconditioner with multigrid inversion) works well for
homogeneous problems, but is less suitable for heterogeneous problems with
sharp changes between wavenumbers due to the instability of the inexact
inversion using geometric multigrid [33, 15].

| $(k_{1},k_{2})=(10,50)$ | $(k_{1},k_{2})=(10,75)$
---|---|---
$\gamma$ | $1$ | $2$ | $1$ | $2$
$\nu=1$ | $28$ | $12$ | $31$ | $12$
$\nu=2$ | $16$ | $8$ | $17$ | $7$
$\nu=3$ | $12$ | $7$ | $12$ | $6$
$\nu=4$ | $10$ | $6$ | $10$ | $6$
$\nu=5$ | $9$ | $6$ | $9$ | $6$

Table 9: Number of V- $(\gamma=1)$ and W-cycles $(\gamma=2)$ for MP 2-B (high
variation). $\nu$ denotes the number of GMRES(3) smoothing steps and
$\beta_{2}={k_{\max}}^{-1}$.

##  6 Conclusion

In this paper, we present the first full multigrid hierarchy for both the V-
and W-cycles using classic multigrid components for the indefinite
nonsymmetric Helmholtz problem. Here, we emphasized the use of classical
components, such as the weighted Jacobi smoother. The resulting algorithm
shows $h-$independent convergence and thus adheres to the classic multigrid
features. Two novel and striking features should be mentioned. First of all,
no restriction is imposed on the number of grid points on the coarsest grid.
Second of all, no level-dependent parameters are needed.

Apart from the numerical results, we provide a new theory to assess
convergence of highly indefinite linear systems. Where most proofs require
some adjointness and/or symmetry assumptions on the underlying linear system,
we have constructed theoretical notions which do not require such assumptions.
In fact, we have found that the addition of a complex shift, solely for the
purpose of coarsening in the multigrid hierarchy, combined with higher-order
interpolation schemes for the inter-grid transfer operators leads to an HPD
system. The positive definiteness can be verified for medium-sized wavenumbers
by assessing properties of the matrix. As a result, it can be shown that the
convergence is independent of $h$. Thus, the computation of the spectral
radius can be circumvented, especially since the analytical eigenvalues can
only be computed in case of Dirichlet boundary conditions. We also study the
behavior of the solver once we use GMRES(3) as a smoother instead of
$\omega-$Jacobi. The convergence significantly improves and less smoothing
steps are needed. Second of all, if we use GMRES(3) as a smoother, we can keep
the complex shift very small without paying a penalty in terms of additional
costs. As a result, we obtain close to $k-$ and $h-$independent convergence,
both for the constant wavenumber model problem and the non-constant wavenumber
model problem. In fact, we can even use the original Helmholtz operator for
coarsening in this configuration without losing efficiency and effectiveness
as the GMRES iterations provide polynomial smoothing.

This work is one of the few to demonstrate, both theoretically and
numerically, the existence and potential of multigrid methods for such linear
systems. However, a lot of future work remains to be done. For example, future
research could focus on exploring different combinations of interpolation
schemes and smoothers and extending proofs. Also, instead of coarsening on the
CSL, other candidates in terms of ’near positive definite operators’ could
also be explored. We have shown convergence for $\omega-$Jacobi, mainly
because of its simplicity and parallelizability, but many other options remain
to be examined. Moreover, once a convergent solver has been established, the
next step in shifting the paradigm would be to study wavenumber independent
convergence.

## References

  * [1] J. H. Bramble, J. E. Pasciak, and J. Xu, The analysis of multigrid algorithms for nonsymmetric and indefinite elliptic problems, Mathematics of Computation, 51 (1988), pp. 389–414.
  * [2] A. Brandt and S. Ta’asan, Multigrid method for nearly singular and slightly indefinite problems, in Multigrid Methods II, Springer, 1986, pp. 99–121.
  * [3] J. Brannick, F. Cao, K. Kahl, R. D. Falgout, and X. Hu, Optimal interpolation and compatible relaxation in classical algebraic multigrid, SIAM Journal on Scientific Computing, 40 (2018), pp. A1473–A1493.
  * [4] M. Brezina, A. J. Cleary, R. D. Falgout, V. E. Henson, J. E. Jones, T. A. Manteuffel, S. F. McCormick, and J. W. Ruge, Algebraic multigrid based on element interpolation (amge), SIAM Journal on Scientific Computing, 22 (2001), pp. 1570–1592.
  * [5] M. Brezina, R. Falgout, S. MacLachlan, T. Manteuffel, S. McCormick, and J. Ruge, Adaptive smoothed aggregation ($\alpha$ sa) multigrid, SIAM review, 47 (2005), pp. 317–346.
  * [6] H. Chen, P. Lu, and X. Xu, A robust multilevel method for hybridizable discontinuous galerkin method for the helmholtz equation, Journal of Computational Physics, 264 (2014), pp. 133–151.
  * [7] H. Chen, H. Wu, and X. Xu, Multilevel preconditioner with stable coarse grid corrections for the helmholtz equation, SIAM Journal on Scientific Computing, 37 (2015), pp. A221–A244.
  * [8] P.-H. Cocquet and M. J. Gander, How large a shift is needed in the shifted helmholtz preconditioner for its effective inversion by multigrid?, SIAM Journal on Scientific Computing, 39 (2017), pp. A438–A478.
  * [9] P.-H. Cocquet, M. J. Gander, and X. Xiang, Closed form dispersion corrections including a real shifted wavenumber for finite difference discretizations of 2d constant coefficient helmholtz problems, SIAM Journal on Scientific Computing, 43 (2021), pp. A278–A308.
  * [10] S. Cools, B. Reps, and W. Vanroose, A new level-dependent coarse grid correction scheme for indefinite helmholtz problems, Numerical Linear Algebra with Applications, 21 (2014), pp. 513–533.
  * [11] V. Dwarka and C. Vuik, Scalable convergence using two-level deflation preconditioning for the helmholtz equation, SIAM Journal on Scientific Computing, 42 (2020), pp. A901–A928.
  * [12] V. Dwarka and C. Vuik, Scalable multi-level deflation preconditioning for highly indefinite time-harmonic waves, Journal of Computational Physics, 469 (2022), p. 111327.
  * [13] H. C. Elman, O. G. Ernst, and D. P. O’leary, A multigrid method enhanced by krylov subspace iteration for discrete helmholtz equations, SIAM Journal on scientific computing, 23 (2001), pp. 1291–1315.
  * [14] Y. Erlangga and E. Turkel, Iterative schemes for high order compact discretizations to the exterior helmholtz equation, ESAIM: Mathematical Modelling and Numerical Analysis, 46 (2012), pp. 647–660.
  * [15] Y. A. Erlangga and R. Nabben, On a multilevel krylov method for the helmholtz equation preconditioned by shifted laplacian, Electronic Transactions on Numerical Analysis, 31 (2008), p. 3.
  * [16] Y. A. Erlangga, C. W. Oosterlee, and C. Vuik, A novel multigrid based preconditioner for heterogeneous helmholtz problems, SIAM Journal on Scientific Computing, 27 (2006), pp. 1471–1492.
  * [17] O. G. Ernst and M. J. Gander, Why it is difficult to solve helmholtz problems with classical iterative methods, in Numerical analysis of multiscale problems, Springer, 2012, pp. 325–363.
  * [18] O. G. Ernst and M. J. Gander, Multigrid methods for helmholtz problems: a convergent scheme in 1d using standard components, Direct and Inverse Problems in Wave Propagation and Applications. De Gruyer, (2013), pp. 135–186.
  * [19] R. D. Falgout and P. S. Vassilevski, On generalizing the algebraic multigrid framework, SIAM Journal on Numerical Analysis, 42 (2004), pp. 1669–1693.
  * [20] W. Hackbusch, Multi-grid methods and applications, vol. 4, Springer Science & Business Media, 2013.
  * [21] G. R. Hadley, A complex jacobi iterative method for the indefinite helmholtz equation, Journal of computational physics, 203 (2005), pp. 358–370.
  * [22] D. Lahaye, H. De Gersem, S. Vandewalle, and K. Hameyer, Algebraic multigrid for complex symmetric systems, IEEE Transactions on Magnetics, 36 (2000), pp. 1535–1538.
  * [23] B. Lee, T. A. Manteuffel, S. F. McCormick, and J. Ruge, First-order system least-squares for the helmholtz equation, SIAM Journal on Scientific Computing, 21 (2000), pp. 1927–1949.
  * [24] I. Livshits and A. Brandt, Accuracy properties of the wave-ray multigrid algorithm for helmholtz equations, SIAM Journal on Scientific Computing, 28 (2006), pp. 1228–1251.
  * [25] S. P. MacLachlan and C. W. Oosterlee, Algebraic multigrid solvers for complex-valued matrices, SIAM Journal on scientific computing, 30 (2008), pp. 1548–1571.
  * [26] Y. Notay, Convergence analysis of perturbed two-grid and multigrid methods, SIAM journal on numerical analysis, 45 (2007), pp. 1035–1044.
  * [27] Y. Notay, An aggregation-based algebraic multigrid method, Electronic transactions on numerical analysis, 37 (2010), pp. 123–146.
  * [28] Y. Notay, Algebraic analysis of two-grid methods: The nonsymmetric case, Numerical Linear Algebra with Applications, 17 (2010), pp. 73–96.
  * [29] Y. Notay, Analysis of two-grid methods: The nonnormal case, Mathematics of Computation, 89 (2020), pp. 807–827.
  * [30] L. N. Olson and J. B. Schroder, Smoothed aggregation for helmholtz problems, Numerical Linear Algebra with Applications, 17 (2010), pp. 361–386.
  * [31] L. N. Olson, J. B. Schroder, and R. S. Tuminaro, A general interpolation strategy for algebraic multigrid using energy minimization, SIAM Journal on Scientific Computing, 33 (2011), pp. 966–991.
  * [32] L. G. Ramos and R. Nabben, On optimal algebraic multigrid methods, arXiv preprint arXiv:1906.01381, (2019).
  * [33] A. Sheikh, D. Lahaye, L. G. Ramos, R. Nabben, and C. Vuik, Accelerating the shifted laplace preconditioner for the helmholtz equation by multilevel deflation, Journal of Computational Physics, 322 (2016), pp. 473–490.
  * [34] J. Wang, Convergence analysis of multigrid algorithms for nonselfadjoint and indefinite elliptic problems, SIAM journal on numerical analysis, 30 (1993), pp. 275–285.