# Convergence analysis of GMRES applied to Helmholtz problems near resonances

Victorita Dolean, Pierre Marchand, Axel Modave, Timothée Raynaud

► To cite this version:

Victorita Dolean, Pierre Marchand, Axel Modave, Timothée Raynaud. Convergence analysis of GMRES applied to Helmholtz problems near resonances. 2025. hal-05078654

HAL Id: hal-05078654

https://hal.science/hal-05078654v1

Preprint submitted on 22 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# CONVERGENCE ANALYSIS OF GMRES APPLIED TO HELMHOLTZ PROBLEMS NEAR RESONANCES*

V. DOLEAN†, P. MARCHAND‡, A. MODAVE†, AND T. RAYNAUD†

**Abstract.** In this work we study how the convergence rate of GMRES is influenced by the properties of linear systems arising from Helmholtz problems near resonances or quasi-resonances. We extend an existing convergence bound to demonstrate that the approximation of small eigenvalues by harmonic Ritz values plays a key role in convergence behavior. Next, we analyze the impact of deflation using carefully selected vectors and combine this with a Complex Shifted Laplacian preconditioner. Finally, we apply these tools to two numerical examples near (quasi-)resonant frequencies, using them to explain how the convergence rate evolves.

**Key words.** Helmholtz equation, resonances, finite elements, GMRES, convergence, harmonic Ritz values, deflation, preconditioning

MSC codes. 65F10, 65N22, 65N30

1. Introduction. The solution of time-harmonic wave propagation problems is of paramount importance in many fields of science and engineering, including acoustics, electromagnetism, and structural mechanics. One of the most fundamental equation for modeling scalar waves is the Helmholtz equation,

$$ (1.1) \qquad -\Delta u - k^2 u = f, $$

where $k \in \mathbb{R}_{>0}$ is the wavenumber, $u$ is an unknown field, and $f$ is a given source. This equation is completed with boundary conditions. Here, we are interested in the solution of (1.1) when $k^2$ is close to an eigenvalue of the Laplacian. Discretizing the resulting problem with the finite element method (FEM) leads to a linear system $Au = b$, where $A \in \mathbb{C}^{N \times N}$ is a sparse non-singular matrix which may be non-symmetric or non-Hermitian, $u \in \mathbb{C}^N$ is unknown, $b \in \mathbb{C}^N$ is given. To ensure solution accuracy, the number of degrees of freedom per wavelength must be sufficiently high, as it must increase faster than the wavenumber to mitigate the pollution effect [1]. Consequently, for large wavenumber problems, the resulting linear system can be extremely large.

Direct solvers, with an algorithmic complexity of $O(N^3)$ in a worst-case scenario, are impractical for large systems [26, Chapter 2]. Instead, iterative methods such as Krylov subspace methods are commonly used. A few matrix-vector multiplications per iteration are required, giving a computational cost of $O(N^2l)$, where $l$ is the number of iterations [26, Chapter 7]. Iterative methods are also easier to parallelize, which is advantageous for solving large-scale problems. The key challenge is to control the number of iterations $l$ needed to achieve a given tolerance [30].

For non-normal matrices, the Generalized Minimal Residual (GMRES) method is a Krylov method of choice [32]. Although GMRES convergence has been widely studied and several bounds have been established, there is no universally applicable estimate for all classes of matrices [15]. In particular, GMRES convergence in the case of discretized Helmholtz problems for different case scenarios remains an open problem due to its highly indefinite or non-self-adjoint nature, as seen in previous studies [18, 24, 20]. In addition, when the wavenumber is close to a resonance value, the number of GMRES iterations required to achieve a given tolerance can significantly increase. Due to the slow convergence of GMRES for these problems, acceleration techniques are usually required [30, 37]. A well suited

<footnote>*Distributed under Creative Commons CC BY 4.0 license.</footnote>

<footnote>**Funding:** This work was funded in part by the ANR JCJC project WavesDG (research grant ANR-21-CE46-0010) and by the Agence de l'Innovation de Défense [AID] through Centre Interdisciplinaire d'Etudes pour la Défense et la Sécurité [CIEDS] (project 2022 ElectroMath).</footnote>

<footnote>† Department of Mathematics and Computer Science, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands (v.dolean.maini@tue.nl).</footnote>

<footnote>‡POEMS, CNRS, Inria, ENSTA, Institut Polytechnique de Paris, 91120 Palaiseau, France (pierre.marc-hand@inria.fr, axel.modave@ensta.fr, timothee.raynaud@ensta.fr).</footnote>
preconditioner for Helmholtz problems is the Complex Shifted Laplacian (CSL) [28, 29, 22]. Domain decomposition [14] and deflation [33, 23] are also used as methods of acceleration.

In this work, we focus on GMRES convergence for a class of quasi-resonant Helmholtz problems. Resonances are eigenvalues of the Laplace's operator in a given domain, where the Helmholtz problem (1.1) becomes ill-posed. Quasi-resonances refer to an increasing sequence of wavenumbers where the norm of the inverse Helmholtz operator grows rapidly [21]. While previous studies have examined the impact of quasi-resonances on GMRES convergence [24], further research is needed for a finer understanding in order to enhance the robustness of iterative solvers in these cases.

The main aim of this paper is to investigate the relationships between GMRES convergence rates and the properties of linear systems associated with Helmholtz problems near resonances. More specifically our main contributions are:

* to generalize a GMRES convergence bound based on harmonic Ritz values;

* to show that the approximation of every eigenvalue or group of eigenvalues by harmonic Ritz values affects convergence;

* to apply deflation techniques with a suitable choice of eigenvectors adapted to resonance problems to accelerate GMRES convergence and combine them with a CSL based preconditioner;

* to understand the numerical behavior of GMRES for different cases scenarios and acceleration methods.

The remainder of the paper is structured as follows. Section 2 introduces the GMRES method and reviews key convergence properties. We then extend a GMRES convergence bound based on harmonic Ritz values, showing that when these values approach a subset of eigenvalues closely enough, the corresponding eigenvalues cease to hinder convergence. In Section 3, we present the continuous Helmholtz problem and the settings in which resonances and quasi-resonances arise. To efficiently solve the discretized problem obtained via the finite element method (FEM), we introduce acceleration techniques tailored to its specific challenges. In particular, we employ a Complex Shifted Laplacian (CSL) preconditioner suited to Helmholtz problems, and we motivate the use of deflation strategies by the presence of resonances. Finally, in Sections 4 and 5, we investigate numerical benchmarks near resonant regimes. These case studies illustrate how the proposed tools influence GMRES convergence and allow us to assess the effectiveness and robustness of the acceleration techniques.

**2. GMRES.** The Generalized Minimal Residual (GMRES) method [32] is a popular iterative algorithm for solving general linear systems. In this section, we introduce its main components and discuss a standard convergence result. We then propose new convergence results based on harmonic Ritz (HR) values, which allow a better interpretation of the GMRES convergence.

**2.1. Key ideas.** The GMRES method belongs to the class of Krylov subspace iterative methods, leveraging the Krylov subspace defined as

$$K_l(\mathbf{A}, \mathbf{r}_0) := \operatorname{span}\{\mathbf{r}_0, \mathbf{A}\mathbf{r}_0, \mathbf{A}^2\mathbf{r}_0, \dots, \mathbf{A}^{l-1}\mathbf{r}_0\}, \quad \text{for } l < N,$$

where $\mathbf{r}_0 := \mathbf{b} - \mathbf{A}\mathbf{x}_0$ is the initial residual and $\mathbf{x}_0$ is an initial guess. At each iteration $l < N$, the approximate solution $\mathbf{x}_l$ belongs to $\mathbf{x}_0 + K_l(\mathbf{A}, \mathbf{r}_0)$, and it is computed such that the residual $\mathbf{r}_l = \mathbf{b} - \mathbf{A}\mathbf{x}_l$ has minimum norm, i.e.

$$ (2.1) \qquad \mathbf{x}_l := \arg\min_{\mathbf{x} \in \mathbf{x}_0 + K_l(\mathbf{A}, \mathbf{r}_0)} \|\mathbf{b} - \mathbf{A}\mathbf{x}\|_2. $$

A critical component of GMRES is the Arnoldi process, which generates an orthonormal basis for the Krylov subspace. The basis vectors are the columns of a matrix $\mathbf{U}_l \in \mathbb{C}^{N \times l}$. This basis is used to compute the restriction of $\mathbf{A}$ to the Krylov subspace, $\mathbf{H}_l = \mathbf{U}_l^* \mathbf{A} \mathbf{U}_l \in \mathbb{C}^{l \times l}$, which is an upper Hessenberg matrix. Each iteration of the GMRES algorithm can be summarized as follows: perform the Arnoldi process to compute $\mathbf{U}_l$ and $\mathbf{H}_l$, then solve the least-squares
<figure><img src="image_2.png" /><figcaption>Fig. 2.1: GMRES convergence history for a cavity problem close to a resonance.</figcaption></figure>

problem min$_{y}$ ||$\beta$e₁ − Hᵢy||₂ with $\beta = \|\mathbf{r}_0\|_2$ and e₁ the first canonical vector, and finally form the approximate solution xᵢ = x₀ + Uᵢy, see e.g. [34].

**2.2. Convergence results.** Several convergence bounds have been developed to analyze the behavior of GMRES, often by expressing the residual vector in terms of polynomials; see, for example, [15]. The residual at iteration *l* satisfies the following relation:

$$
(2.2) \qquad \|\mathbf{r}_l\|_2 = \min_{q_l \in \mathcal{P}_l^1} \|q_l(\mathbf{A})\mathbf{r}_0\|_2,
$$

where $\mathcal{P}_l^1$ is the set of polynomials of degree at most $l$ with the constraint $q_l(0) = 1$. Through-out this paper, we denote by $p_l \in \mathcal{P}_l^1$ the polynomial that achieves this minimum, such that $\mathbf{r}_l = p_l(\mathbf{A})\mathbf{r}_0$.

GMRES convergence can be understood in terms of the eigenvalue distribution of **A** and the conditioning of its eigenvectors. If **A** is diagonalizable, we can write **A** = **V****Δ**V**⁻¹. For any polynomial $q_l \in \mathcal{P}_l^1$, it follows that:

$$
\|q_l(\mathbf{A})\mathbf{r}_0\|_2 = \|\mathbf{V}q_l(\mathbf{\Lambda})\mathbf{V}^{-1}\mathbf{r}_0\|_2 \leq \kappa_2(\mathbf{V}) \|q_l(\mathbf{\Lambda})\|_2 \|\mathbf{r}_0\|_2,
$$

where $\kappa_2(\mathbf{V}) := \|\mathbf{V}\|_2 \|\mathbf{V}^{-1}\|_2$ is the condition number of the eigenvector matrix. This leads to the classical convergence bound [32]:

$$
(2.3) \qquad \frac{\|\mathbf{r}_l\|_2}{\|\mathbf{r}_0\|_2} \le \kappa_2(\mathbf{V}) \min_{q_l \in \mathcal{P}_l^1} \max_{\lambda \in \sigma(\mathbf{A})} |q_l(\lambda)|,
$$

where $\sigma(\mathbf{A})$ denotes the set of eigenvalues of $\mathbf{A}$. If $\mathbf{A}$ is normal, then $\kappa_2(\mathbf{V}) = 1$, and the convergence depends solely on the eigenvalue distribution. In contrast, for non-normal matrices, $\kappa_2(\mathbf{V})$ can be large, making the bound in (2.3) less predictive of actual convergence behavior [15].

Alternative convergence bounds relate to other properties of **A**, such as its field of values [3, 16, 12] or its pseudospectrum [34]. While these bounds have distinct advantages, none of them provides an accurate or systematic prediction of GMRES convergence. Moreover, they describe a linear asymptotic convergence rate [15], which fails to capture specific convergence behaviors, such as stagnation phases followed by rapid decays. Such patterns can be observed in the GMRES convergence history, and can occur close to resonances (see Figure 2.1 and Sections 4 and 5). To better understand this phenomenon, it is crucial to highlight the role of small eigenvalues related to these resonances. This need motivates our interest in harmonic Ritz values. These results also highlight how both spectral and non-normality properties of **A** shape the convergence of GMRES.
**2.3. Harmonic Ritz values.** Harmonic Ritz (HR) values provide estimates of the eigenvalues of $\mathbf{A}$; see, for example, [34, Chapter 26]. They are defined as follows:

**DEFINITION 2.1 (Harmonic Ritz (HR) values).** At iteration $l > 0$, the harmonic Ritz values $\{\nu_j^{(l)}\}_{j=1}^l$ are the roots of the minimizing polynomial $p_l \in \mathcal{P}_l^1$ from (2.2).

By construction, the number of HR values increases with each GMRES iteration. Since $p_l(0) = 1$, none of the HR values can be zero. As a result, the polynomial $p_l$ can be factored as:

$$ (2.4) \qquad p_l(z) = \prod_{j=1}^{l} \left( 1 - \frac{z}{\nu_j^{(l)}} \right). $$

The HR values satisfy a useful lower bound:

**PROPOSITION 2.2 (Lower bound [10, Theorem 4]).** Assume $\mathbf{A}$ is non-singular. At iteration $l > 0$, the HR values satisfy

$$ |\nu_j^{(l)}| \ge s_{\min}(\mathbf{A}), \quad \forall j = 1, \dots, l, $$

where $s_{\min}(\mathbf{A})$ is the smallest singular value of $\mathbf{A}$.

As iterations progress, the HR values tend to approach the eigenvalues of $\mathbf{A}$ increasingly well. They converge to the full spectrum after $N$ iterations, or earlier in cases of “lucky breakdown” [34, Chapter 26]. Proposition 2.2 guarantees that these approximations remain bounded away from zero.

**Remark 2.3 (Ritz values).** The Ritz values are the eigenvalues of the upper Hessenberg matrix $\mathbf{H}_l$, and they coincide with the roots of the minimizing polynomial in the Full Orthogonalization Method (FOM). Like the HR values, Ritz values serve as approximations of the eigenvalues of $\mathbf{A}$. In [36], Ritz values are used to analyze the convergence of the conjugate gradient method for symmetric positive definite systems. This approach is extended in [13] to study the convergence of GMRES for general non-normal matrices. However, because Ritz values are associated with FOM rather than GMRES, the analysis must relate the residual norms of the two methods. Specifically, the identity $\|\mathbf{r}_l^{\text{GMRES}}\|_2 = C_l \|\mathbf{r}_l^{\text{FOM}}\|_2$ is used, where the constant $C_l$ depends on the iteration. As a result, the convergence bounds derived from Ritz values involve the factor $C_l$, making them generally looser than those obtained directly using HR values.

**2.4. Convergence results with HR values.** While Ritz values offer useful spectral insights, their indirect connection to GMRES residuals often leads to looser bounds. In contrast, convergence estimates based on HR values—being directly tied to the GMRES residual polynomial—can yield sharper and more informative results. Cao [10] proposed such a bound by linking GMRES convergence to the approximation of the smallest eigenvalue by an HR value. Building on this idea, we introduce two generalizations that consider the approximation of an entire subset of eigenvalues by corresponding HR values, rather than just a single one. This leads to a natural partition of the spectrum into two groups: the eigenvalues closely approximated by HR values, and the remaining ones. This strategy is reminiscent of spectral clustering techniques that separate eigenvalue distributions into dominant clusters and peripheral outliers [9, 24].

To formalize this spectral separation and derive quantitative bounds, we now consider the application of GMRES to a linear system of the form $\mathbf{A}\mathbf{x} = \mathbf{b}$, with $\mathbf{A} \in \mathbb{C}^{N \times N}$ and $\mathbf{b} \in \mathbb{C}^N$. We denote the residual vector at iteration $l$ by $\mathbf{r}_l$. The convergence results presented in this section rely on the following assumption.

**ASSUMPTION 2.4.** The matrix $\mathbf{A}$ is non-singular and diagonalizable.

We aim to analyze the convergence between two iterations, $l$ and $l+m$, assuming $0 < l < l+m < N$. We begin with a preliminary result that sets the stage for the generalized bounds.
LEMMA 2.5. Suppose that Assumption 2.4 holds. Let $J \in [1, l]$, and consider a set $N_J^{(l)}$ of J HR values at iteration $l$, along with a set $\Lambda_J$ of J eigenvalues of $\mathbf{A}$. Then, for any $m > 0$ such that $l + m < N$, the following bound holds:

$$ \frac{\|\mathbf{r}_{l+m}\|_2}{\|\mathbf{r}_l\|_2} \leq \min_{q_m \in \mathcal{P}_m^1} \|q_m(\mathbf{A}) s_J^l(\mathbf{A})\|_2, $$

where

$$ (2.5) \qquad s_J^l(z) := \frac{\prod_{\lambda_j \in \Lambda_J} \left(1 - \frac{z}{\lambda_j}\right)}{\prod_{\nu_j^{(l)} \in N_J^{(l)}} \left(1 - \frac{z}{\nu_j^{(l)}}\right)}. $$

Proof. Let $p_l \in \mathcal{P}_l^1$ denote the GMRES minimizing polynomial at iteration $l$, so that $\mathbf{r}_l = p_l(\mathbf{A})\mathbf{r}_0$. For any polynomial $q_m \in \mathcal{P}_m^1$, define a new polynomial of degree $l+m$ by:

$$ \tilde{p}_{l+m}(z) := q_m(z) s_J^l(z) p_l(z). $$

Since the denominator of $s_J^l$ consists of HR values—i.e., roots of $p_l$ (2.4)—these cancel out, ensuring that $\tilde{p}_{l+m}(z)$ is a well-defined polynomial of degree at most $l+m$. Furthermore, $\tilde{p}_{l+m}(0) = q_m(0) \cdot s_J^l(0) \cdot p_l(0) = 1$, so $\tilde{p}_{l+m} \in \mathcal{P}_{l+m}^1$.

From the GMRES optimality property (2.2) at iteration $l+m$, we have:

$$ \|\mathbf{r}_{l+m}\|_2 \leq \min_{q_{l+m} \in \mathcal{P}_{l+m}^1} \|q_{l+m}(\mathbf{A})\mathbf{r}_0\|_2. $$

In particular, using $\tilde{p}_{l+m}$ yields:

$$ \begin{aligned} \|\mathbf{r}_{l+m}\|_2 &\leq \|\tilde{p}_{l+m}(\mathbf{A})\mathbf{r}_0\|_2 \\ &= \|q_m(\mathbf{A}) s_J^l(\mathbf{A}) p_l(\mathbf{A})\mathbf{r}_0\|_2 \\ &= \|q_m(\mathbf{A}) s_J^l(\mathbf{A}) \mathbf{r}_l\|_2 \\ &\leq \|q_m(\mathbf{A}) s_J^l(\mathbf{A})\|_2 \|\mathbf{r}_l\|_2. \end{aligned} $$

Taking the minimum over all $q_m \in \mathcal{P}_m^1$ completes the proof. $\square$

Lemma 2.5 provides a general convergence estimate based on the positions of the HR values relative to the eigenvalues. To make this bound more explicit and interpretable, we now use the eigendecomposition of $\mathbf{A}$ to express the norm $q_m(\mathbf{A})s_J^l(\mathbf{A})$ in spectral terms. This leads to our first main result, which quantifies the effect of the remaining eigenvalues—those not approximated by HR values—on the convergence behavior.

THEOREM 2.6 (First convergence bound). If the assumptions of Lemma 2.5 hold, then, for iterations $l$ and $l+m$, we have

$$ \frac{\|\mathbf{r}_{l+m}\|_2}{\|\mathbf{r}_l\|_2} \leq \left( \sum_{\lambda_i \in \Lambda_J^c} \kappa(\lambda_i) \right) \max_{\lambda_i \in \Lambda_J^c} |s_J^l(\lambda_i)| \min_{q_m \in \mathcal{P}_m^1} \max_{\lambda_i \in \Lambda_J^c} |q_m(\lambda_i)|, $$

with $\Lambda_J^c := \sigma(\mathbf{A}) \setminus \Lambda_J$ and $\kappa(\lambda_i) := \|\hat{\mathbf{v}}_i\|_2 \|\mathbf{v}_i\|_2$, where $\hat{\mathbf{v}}_i$ and $\mathbf{v}_i \in \mathbb{C}^n$ are the left and right eigenvectors associated to $\lambda_i \in \sigma(\mathbf{A})$, with $\hat{\mathbf{v}}_i^* \mathbf{v}_i = 1$ (see [34, Chapter 52]).

Proof. Since $\mathbf{A}$ is diagonalizable, it admits the decomposition $\mathbf{A} = \mathbf{V}\mathbf{D}\mathbf{V}^{-1}$, where $\mathbf{V} \in \mathbb{C}^{N \times N}$ is the matrix of right eigenvectors, and $\mathbf{D} = \text{diag}(\lambda_1, \dots, \lambda_N)$ is the diagonal matrix of eigenvalues. We denote by $\mathbf{v}_i$ the $i$-th column of $\mathbf{V}$ (the right eigenvector associated with $\lambda_i$), and by $\hat{\mathbf{v}}_i^*$ the $i$-th row of $\mathbf{V}^{-1}$, corresponding to the conjugate transpose of the left eigenvector associated with $\lambda_i$.
By using the eigenvalue decomposition from Lemma 2.5, we obtain

$$ \frac{\|\mathbf{r}_{l+m}\|_2}{\|\mathbf{r}_l\|_2} \le \min_{q_m \in \mathcal{P}_m^1} \left\| \mathbf{V} q_m(\mathbf{D}) s_J^l(\mathbf{D}) \mathbf{V}^{-1} \right\|_2 \\ \le \min_{q_m \in \mathcal{P}_m^1} \left\| \sum_{\lambda_i \in \sigma(\mathbf{A})} q_m(\lambda_i) s_J^l(\lambda_i) \mathbf{v}_i \hat{\mathbf{v}}_i^* \right\|_2. $$

Because $s_J^l(\lambda_j) = 0$ for all $\lambda_j \in \Lambda_J$, the sum can be restricted to $\lambda_i \in \Lambda_J^c$, leading to

$$ \frac{\|\mathbf{r}_{l+m}\|_2}{\|\mathbf{r}_l\|_2} \le \min_{q_m \in \mathcal{P}_m^1} \left\| \sum_{\lambda_i \in \Lambda_J^c} q_m(\lambda_i) s_J^l(\lambda_i) \mathbf{v}_i \hat{\mathbf{v}}_i^* \right\|_2 \\ \le \min_{q_m \in \mathcal{P}_m^1} \sum_{\lambda_i \in \Lambda_J^c} |q_m(\lambda_i)| |s_J^l(\lambda_i)| \|\mathbf{v}_i\|_2 \|\hat{\mathbf{v}}_i^*\|_2 \\ \le \min_{q_m \in \mathcal{P}_m^1} \sum_{\lambda_i \in \Lambda_J^c} |q_m(\lambda_i)| |s_J^l(\lambda_i)| \kappa(\lambda_i) \\ \le \left( \sum_{\lambda_i \in \Lambda_J^c} \kappa(\lambda_i) \right) \max_{\lambda_i \in \Lambda_J^c} |s_J^l(\lambda_i)| \min_{q_m \in \mathcal{P}_m^1} \max_{\lambda_i \in \Lambda_J^c} |q_m(\lambda_i)|. \quad \square $$

*Interpretation.* If, at iteration $l$, the HR values in $N_J^{(l)}$ closely approach the eigenvalues in $\Lambda_J$, then the ratio $\max_{\lambda_i \in \Lambda_J^c} |s_J^l(\lambda_i)|$ approaches 1. In this case, the convergence rate for subsequent iterations is effectively governed only by the remaining eigenvalues in $\Lambda_J^c$. This implies that GMRES behaves as though the eigenvalues in $\Lambda_J$ have been deflated—i.e., they no longer hinder convergence.

An important advantage of expressing convergence using the ratio $\|\mathbf{r}_{l+m}\|_2 / \|\mathbf{r}_l\|_2$ is that it allows for a stage-by-stage analysis of the convergence process. This formulation captures the evolution of the convergence rate between any two iterations and helps interpret the typically nonlinear behavior of GMRES convergence.

**Remark 2.7** (Link with Cao [10]). The convergence bound proposed in [10, Theorem 6] is given by:

$$ \frac{\|\mathbf{r}_{l+m}\|_2}{\|\mathbf{r}_l\|_2} \le \kappa_2(\mathbf{V}) \max_{\lambda_i \neq \lambda_1} \left| \frac{1 - \frac{\lambda_i}{\lambda_1}}{1 - \frac{\lambda_i}{\nu_1^{(l)}}} \right| \min_{q_m \in \mathcal{P}_m^1} \max_{\lambda_i \neq \lambda_1} |q_m(\lambda_i)|, $$

where $\lambda_1$ is the smallest eigenvalue of $\mathbf{A}$, and $\nu_1^{(l)}$ is an HR value at iteration $l$. This bound focuses on the approximation of a single eigenvalue and includes the global condition number $\kappa_2(\mathbf{V})$, which reflects the influence of all eigenvectors and the overall non-normality of $\mathbf{A}$. In contrast, the bound in Theorem 2.6 generalizes this idea by allowing for the approximation of a subset of eigenvalues and excludes the contribution of well-approximated ones. As a result, it provides a more localized and consistent characterization of GMRES convergence behavior.

The bound in Theorem 2.6 provides a convergence estimate based on eigenvalue separation using the eigendecomposition of $\mathbf{A}$. However, when $\mathbf{A}$ is non-normal or nearly defective, its eigenvalues alone may not fully capture the convergence behavior of GMRES. In such cases, a more robust analysis can be achieved using the concept of the pseudospectrum, which accounts for the sensitivity of the spectrum to perturbations and highlights the influence of non-normality. To establish the second main result, which leverages pseudospectral properties, we first introduce a technical lemma involving spectral projectors. This result allows us to isolate the effect of different parts of the spectrum through contour integrals and is crucial in the derivation of the forthcoming bound.
LEMMA 2.8 (Spectral projectors). Assume that the hypotheses of Lemma 2.5 hold. Let $\Gamma_J$ (resp. $\Gamma_J^c$) be a Jordan curve enclosing the eigenvalues in $\Lambda_J$ (resp. $\Lambda_J^c := \sigma(\mathbf{A}) \setminus \Lambda_J$) of $\mathbf{A}$ and no other eigenvalues. The spectral projectors associated with $\Lambda_J$ and $\Lambda_J^c$ are defined by

$$ \mathbf{P}_J := \frac{-1}{2i\pi} \int_{\Gamma_J} (\mathbf{A} - z\mathbf{I})^{-1} dz \quad \text{and} \quad \mathbf{P}_J^c := \frac{-1}{2i\pi} \int_{\Gamma_J^c} (\mathbf{A} - z\mathbf{I})^{-1} dz, $$

respectively. Then, at iteration l we have

$$ s_J^l(\mathbf{A}) = \mathbf{P}_J^c s_J^l(\mathbf{A}). $$

*Proof.* We begin by recalling standard properties of spectral projectors (see [31, Theorem 3.3, Lemma 3.1]):

(i) $\mathbf{P}_J^2 = \mathbf{P}_J$ and $\mathbf{P}_J^c = \mathbf{P}_J^c$ (idempotence),

(ii) $\mathbf{P}_J + \mathbf{P}_J^c = \mathbf{I}$ (completeness),

(iii) $\mathbf{P}_J \mathbf{P}_J^c = \mathbf{P}_J^c \mathbf{P}_J = 0$ (orthogonality),

(iv) $\mathbf{P}_J \mathbf{A} = \mathbf{A} \mathbf{P}_J$ and $\mathbf{P}_J^c \mathbf{A} = \mathbf{A} \mathbf{P}_J^c$ (commutativity with $\mathbf{A}$),

(v) $\ker [\prod_{\lambda_j \in \Lambda_J} (\mathbf{I} - \mathbf{A}/\lambda_j)] = \text{ran}(\mathbf{P}_J)$.

Let $q(\mathbf{A})$ be any matrix polynomial. Using property (ii), we can write:

$$ q(\mathbf{A}) = (\mathbf{P}_J + \mathbf{P}_J^c) q(\mathbf{A}) = \mathbf{P}_J q(\mathbf{A}) + \mathbf{P}_J^c q(\mathbf{A}). $$

Applying this decomposition to $s_J^l(\mathbf{A})$, we obtain:

$$ \begin{aligned} s_J^l(\mathbf{A}) &= \mathbf{P}_J^c s_J^l(\mathbf{A}) + \mathbf{P}_J s_J^l(\mathbf{A}) \\ &= \mathbf{P}_J^c \left[ \prod_{\lambda_j \in \Lambda_J} \left( \mathbf{I} - \frac{\mathbf{A}}{\lambda_j} \right) \right] \left[ \prod_{\nu_j^{(l)} \in N_J^{(l)}} \left( \mathbf{I} - \frac{\mathbf{A}}{\nu_j^{(l)}} \right)^{-1} \right] \\ &\quad + \mathbf{P}_J \left[ \prod_{\lambda_j \in \Lambda_J} \left( \mathbf{I} - \frac{\mathbf{A}}{\lambda_j} \right) \right] \left[ \prod_{\nu_j^{(l)} \in N_J^{(l)}} \left( \mathbf{I} - \frac{\mathbf{A}}{\nu_j^{(l)}} \right)^{-1} \right]. \end{aligned} $$

We now analyze the second term. The factor $\prod_{\lambda_j \in \Lambda_J} (\mathbf{I} - \mathbf{A}/\lambda_j)$ acts as a zero operator on $\text{ran}(\mathbf{P}_J)$ due to the property (v) of spectral projectors. Since $\mathbf{P}_J \mathbf{A} = \mathbf{A} \mathbf{P}_J$, all matrix functions of $\mathbf{A}$ also commute with $\mathbf{P}_J$. Therefore, the entire second term vanishes:

$$ \mathbf{P}_J s_J^l(\mathbf{A}) = 0, $$

and the conclusion follows.

THEOREM 2.9 (Second Convergence Bound). Let the assumptions of Lemma 2.5 hold. Let $\Gamma_J^c$ be a Jordan curve that encloses the eigenvalues in the set $\Lambda_J^c$ and no others. Suppose there exists a constant $\varepsilon_J^c > 0$ such that

$$ \|\mathbf{A} - z\mathbf{I}\|_2^{-1} \le \frac{1}{\varepsilon_J^c}, \quad \forall z \in \Gamma_J^c. $$

Then, for iteration indices $l$ and $l+m$, the residuals satisfy:

$$ \frac{\|\mathbf{r}_{l+m}\|_2}{\|\mathbf{r}_l\|_2} \le \frac{\mathcal{L}(\Gamma_J^c)}{2\pi\varepsilon_J^c} \cdot \max_{z \in \Gamma_J^c} |s_J^l(z)| \cdot \min_{q_m \in \mathcal{P}_m^1} \max_{z \in \Gamma_J^c} |q_m(z)|, $$

where $\mathcal{L}(\Gamma_J^c)$ denotes the length of the contour $\Gamma_J^c$.
*Proof.* From Lemma 2.5 and Lemma 2.8, we have:

$$ \frac{\|\mathbf{r}_{l+m}\|_2}{\|\mathbf{r}_l\|_2} \leq \min_{q_m \in \mathcal{P}_m^1} \|\mathbf{P}_J^c q_m(\mathbf{A}) s_J^l(\mathbf{A})\|_2. $$

Using the integral representation of the spectral projector $\mathbf{P}_J^c$, we obtain:

$$ \frac{\|\mathbf{r}_{l+m}\|_2}{\|\mathbf{r}_l\|_2} \leq \min_{q_m \in \mathcal{P}_m^1} \frac{1}{2\pi} \int_{\Gamma_J^c} |q_m(z)| \cdot |s_J^l(z)| \cdot \|(\mathbf{A} - z\mathbf{I})^{-1}\|_2 \, dz. $$

Bounding the integrand by its maximum on $\Gamma_J^c$, we find:

$$ \frac{\|\mathbf{r}_{l+m}\|_2}{\|\mathbf{r}_l\|_2} \leq \frac{\mathcal{L}(\Gamma_J^c)}{2\pi} \cdot \max_{z \in \Gamma_J^c} |s_J^l(z)| \cdot \min_{q_m \in \mathcal{P}_m^1} \max_{z \in \Gamma_J^c} |q_m(z)| \cdot \max_{z \in \Gamma_J^c} \|(\mathbf{A} - z\mathbf{I})^{-1}\|_2. $$

Applying the resolvent bound yields the desired result:

$$ \frac{\|\mathbf{r}_{l+m}\|_2}{\|\mathbf{r}_l\|_2} \leq \frac{\mathcal{L}(\Gamma_J^c)}{2\pi\varepsilon_J^c} \cdot \max_{z \in \Gamma_J^c} |s_J^l(z)| \cdot \min_{q_m \in \mathcal{P}_m^1} \max_{z \in \Gamma_J^c} |q_m(z)|. \quad \square $$

*Interpretation.* In Theorem 2.9, the convergence bound involves a maximum over the contour $\Gamma_J^c$, which is a continuous path in the complex plane. The choice of this contour is restricted by the requirement that the resolvent norm $\|(\mathbf{A} - z\mathbf{I})^{-1}\|_2$ remains bounded above by $1/\varepsilon_J^c$ for all $z \in \Gamma_J^c$. In cases where $\mathbf{A}$ is highly non-normal, especially when eigenvalues are clustered around a non-normal eigenvalue, ensuring this bound may force $\varepsilon_J^c$ to be very small.

As a result, even if an HR value is close to such an eigenvalue, its influence can still be felt implicitly in the convergence bound via the smallness of $\varepsilon_J^c$. This highlights the key difference between this *pseudospectrum-based* bound (Theorem 2.9) and the earlier *spectrum-based* bound (Theorem 2.6): the latter is completely unaffected by the eigenvalues that are already well-approximated by HR values, whereas the former can still be influenced by them indirectly through the geometry and norm behavior on the contour.

**3. Helmholtz problem and acceleration techniques.** In this section, we introduce the Helmholtz boundary value problem, the linear system arising from its finite element discretization, and techniques to accelerate the iterative solution process. We also discuss the difficulties caused by resonances and quasi-resonances occurring near the wavenumber, which can significantly impair convergence.

**3.1. Boundary value problem.** Let $\Omega \subset \mathbb{R}^d$ be a bounded domain with Lipschitz boundary $\Gamma$, which is partitioned into three disjoint parts: $\Gamma = \Gamma_D \cup \Gamma_N \cup \Gamma_R$, corresponding to Dirichlet, Neumann, and Robin boundary segments, respectively. The Helmholtz problem is defined as follows:

$$ (3.1) \qquad \left\{ \begin{array}{ll} -\Delta u - k^2 u = f & \text{in } \Omega, \\ u = g_D & \text{on } \Gamma_D, \\ \partial_n u = g_N & \text{on } \Gamma_N, \\ \partial_n u + iku = 0 & \text{on } \Gamma_R, \end{array} \right. $$

with source term $f$ and boundary data $g_D$, $g_N$, and $g_R$.

A variational formulation is obtained by multiplying the Helmholtz equation by a test function $v \in V_0 := \{v \in H^1(\Omega) \mid v = 0 \text{ on } \Gamma_D\}$ and integrating by parts. The solution space is $V := \{v \in H^1(\Omega) \mid v = g_D \text{ on } \Gamma_D\}$. Taking into account the boundary conditions yields the weak form:

$$ \text{Find } u \in V \text{ such that } a(u, v) = \ell(v) \text{ for all } v \in V_0. $$
Here, the sesquilinear form $a$ and linear form $\ell$ are defined by

$$a(u, v) = \int_{\Omega} (\nabla u \cdot \nabla \bar{v} - k^2 u \bar{v}) d\mathbf{x} + ik \int_{\Gamma_R} u \bar{v} ds, \quad \ell(v) = \int_{\Omega} f \bar{v} d\mathbf{x} + \int_{\Gamma_N} g_N \bar{v} ds.$$

If $\Gamma_R = \emptyset$, the problem may be ill-posed for certain wavenumbers $k$, known as resonances, which are eigenvalues of the Laplace operator on $\Omega$.

In unbounded domains, similar phenomena called quasi-resonances can occur. The Helmholtz equation is then supplemented with the Sommerfeld radiation condition at infinity, and solutions are sought in $H_{\text{loc}}^1(\Omega)$, the space of functions locally in $H^1(\Omega)$. Quasimodes associated with these quasi-resonances are defined below.

**DEFINITION 3.1** (Quasimodes [24, Definition 1.1]). A family of Dirichlet (or Neumann) quasimodes with quality $\epsilon(k)$ is a sequence

$$\{(u_j, k_j)\}_{j=1}^{\infty} \subset H_{\text{loc}}^1(\Omega) \times \mathbb{R},$$

where $u_j = 0$ (or $\partial_n u_j = 0$) on $\Gamma$, such that $k_j \to \infty$ as $j \to \infty$, and there exists a compact set $K \subset \Omega$ such that, for all $j$,

$$\text{supp}(u_j) \subset K, \quad \| -\Delta u_j - k_j^2 u_j \|_{L^2(\Omega)} \le \epsilon(k_j), \quad \text{and} \quad \| u_j \|_{L^2(\Omega)} = 1.$$

In practice, numerical simulations are carried out on truncated computational domains equipped with artificial boundary conditions—such as absorbing boundary conditions [17, 27] or perfectly matched layers (PMLs) [6, 7, 5, 8, 4, 25, 35]. In this work, we employ PMLs (see Section 5). Even with domain truncation, quasimodes can still arise, although the overall problem remains well-posed.

As the wavenumber $k$ approaches a resonance or quasi-resonance, the conditioning of the problem deteriorates, making the resulting linear system significantly harder to solve.

**3.2. Discretization.** To numerically solve the Helmholtz problem introduced above, we apply a polynomial finite element discretization on a conforming mesh of the domain. Each mesh element has a characteristic size $h$. The discrete solution $u_h$ is sought in a finite-dimensional subspace $V_h \subset V$ and satisfies the variational formulation

$$a(u_h, v_h) = \ell(v_h), \quad \forall v_h \in V_{h,0} := V_h \cap V_0.$$

The solution $u_h$ can be expanded as

$$u_h(\mathbf{x}) = \sum_{j=1}^{N} u_j \phi_j(\mathbf{x}) + \sum_{j=N+1}^{N+N_D} \hat{u}_{j-N} \phi_j(\mathbf{x}),$$

where $\{\phi_j\}_{j=1}^{N+N_D}$ is a basis of $V_h$, with $\{\phi_j\}_{j=1}^N$ forming a basis for $V_{h,0}$. The coefficients $\{u_j\}$ are unknowns, while $\{\hat{u}_j\}$ encode the Dirichlet data on $\Gamma_D$. Substituting this representation into the variational formulation and choosing the test functions from $V_{h,0}$ leads to the linear system

$$(3.2) \qquad | \text{Find } \mathbf{u} \in \mathbb{C}^N \text{ such that } \mathbf{A}\mathbf{u} = \mathbf{b},$$

with system matrix and right-hand side defined by

$$A_{i,j} = a(\phi_j, \phi_i), \quad b_i = \ell(\phi_i) - \sum_{j=N+1}^{N+N_D} a(\phi_j, \phi_i) \hat{u}_{j-N}, \quad \text{for } i,j = 1, \dots, N.$$

The properties of the discrete system mirror those of the continuous problem:

* The matrix $\mathbf{A}$ is non-singular, since $k$ is assumed not to correspond to a resonance.
* The sesquilinear form is typically not coercive. As a result, **A** is often highly *indefinite*, and this indefiniteness increases with *k*.

* The matrix **A** may also be *non-normal*, meaning its eigenvectors are not mutually orthogonal.

These features make the Helmholtz linear system particularly challenging for iterative solvers. At high frequencies or near (quasi-)resonances, these difficulties become more pronounced: [21] shows that quasimodes give rise to small eigenvalues. In the boundary element setting, [24] analyzes how GMRES convergence deteriorates near quasimodes. Here, with a finite element discretization, we likewise observe the emergence of small eigenvalues and corresponding GMRES stagnation.

**3.3. Acceleration techniques and preconditioning.** In this section, we describe two techniques to accelerate the convergence of GMRES for solving the discrete problem (3.2): the deflation technique and the CSL preconditioning technique.

**Deflation.** The principle of the deflation technique is to remove some well-chosen eigenvalues from the spectrum of **A** by using projection operators, see e.g. [33]. We assume that we have an exact or approximate expression of $n_{def} \ll N$ eigenvectors, which we store in the columns of a matrix $\mathbf{Z} \in \mathbb{C}^{N \times n_{def}}$.

**DEFINITION 3.2.** (*Deflation matrices*) Let $\mathbf{Z} \in \mathbb{C}^{N \times n_{def}}$ be a full rank matrix such that $\ker(\mathbf{Z}^*) \cap \text{range}(\mathbf{A}\mathbf{Z}) = \{0\}$, and let

$$ \mathbf{E} := \mathbf{Z}^* \mathbf{A} \mathbf{Z} \in \mathbb{C}^{n_{def} \times n_{def}} \quad \text{and} \quad \mathbf{Q} := \mathbf{Z} \mathbf{E}^{-1} \mathbf{Z}^*. $$

The deflation matrices are defined as

$$ \mathbf{P}_{\text{def}} := \mathbf{I} - \mathbf{A}\mathbf{Q} \quad \text{and} \quad \mathbf{Q}_{\text{def}} := \mathbf{I} - \mathbf{Q}\mathbf{A}. $$

Matrix **E** is the restriction of **A** to the deflated subspace. It is non-singular and a priori easy to invert, since $n_{def}$ is assumed to be small.

The deflation matrices are projectors onto the complement of the deflated subspace and satisfy $\mathbf{P}_{\text{def}}^2 = \mathbf{P}_{\text{def}}$, $\mathbf{Q}_{\text{def}}^2 = \mathbf{Q}_{\text{def}}$ and $\mathbf{P}_{\text{def}}\mathbf{A} = \mathbf{A}\mathbf{Q}_{\text{def}}$, see [33, Lemma 2]. Applying $\mathbf{P}_{\text{def}}$ or $\mathbf{Q}_{\text{def}}$ to **A** removes the components associated with the deflated eigenvectors. The deflated problem is [23, Lemma 3.2]

$$ (3.3) \qquad | \text{Find } \tilde{\mathbf{u}} \in \mathbb{C}^N \text{ such that } (\mathbf{P}_{\text{def}}\mathbf{A})\tilde{\mathbf{u}} = (\mathbf{P}_{\text{def}}\mathbf{b}), $$

and the solution of the original problem (3.2) is obtained by using

$$ \mathbf{u} = \mathbf{Q}\mathbf{b} + \mathbf{Q}_{\text{def}}\tilde{\mathbf{u}}. $$

This can be shown by noting that **u** can be decomposed as

$$ (3.4) \qquad \mathbf{u} = (\mathbf{I} - \mathbf{Q}_{\text{def}}) \mathbf{u} + \mathbf{Q}_{\text{def}} \mathbf{u} = \mathbf{Q}\mathbf{b} + \mathbf{Q}_{\text{def}} \mathbf{u}. $$

Using the property $\mathbf{P}_{\text{def}}\mathbf{A} = \mathbf{A}\mathbf{Q}_{\text{def}}$, (3.3) gives

$$ \mathbf{A}\mathbf{Q}_{\text{def}}\tilde{\mathbf{u}} = \mathbf{P}_{\text{def}}\mathbf{A}\tilde{\mathbf{u}} = \mathbf{P}_{\text{def}}\mathbf{b} = \mathbf{P}_{\text{def}}\mathbf{A}\mathbf{u} = \mathbf{A}\mathbf{Q}_{\text{def}}\mathbf{u}. $$

Since **A** is non-singular, one has $\mathbf{Q}_{\text{def}}\tilde{\mathbf{u}} = \mathbf{Q}_{\text{def}}\mathbf{u}$.

Because $\mathbf{P}_{\text{def}}$ projects out the deflated subspace, the matrix $\mathbf{P}_{\text{def}}\mathbf{A}$ is singular, with zero as an eigenvalue of multiplicity *s*. But GMRES applied to this system still converges to a solution if the following condition is satisfied [33, Theorem 1]:

$$ (3.5) \qquad \ker(\mathbf{P}_{\text{def}}\mathbf{A}) \cap \text{range}(\mathbf{P}_{\text{def}}\mathbf{A}) = \{0\}. $$

If the columns of $\mathbf{Z}$ are exact eigenvectors of **A**, then the hypothesis of **Definition 3.2** holds. In this case, the condition (3.5) is automatically satisfied, and GMRES is guaranteed to converge to the solution [33, Lemma 2].
Although the deflated system resembles a system preconditioned on the left by $\mathbf{P}_{\text{def}}$, GMRES applied to (3.3) still minimizes the residual norm of the original problem. This follows from the decomposition (3.4) and the definition of $\mathbf{P}_{\text{def}}$; see [19, Theorem 3.1].

Each GMRES iteration becomes more expensive due to the application of $\mathbf{P}_{\text{def}}$. Nevertheless, if $n_{\text{def}}$ is small, computing $\mathbf{E}^{-1}$ is inexpensive. In addition, noting that $\mathbf{P}_{\text{def}}\mathbf{A} = \mathbf{A} - (\mathbf{A}\mathbf{Z})(\mathbf{E}^{-1}\mathbf{Z}^*\mathbf{A})$, the $N \times n_{\text{def}}$ matrix $\mathbf{A}\mathbf{Z}$ and the $n_{\text{def}} \times N$ matrix $\mathbf{E}^{-1}\mathbf{Z}^*\mathbf{A}$ can be stored, and then applying $\mathbf{P}_{\text{def}}\mathbf{A}$ requires only $4Nn_{\text{def}} - N - n_{\text{def}}$ additional operations compared to applying $\mathbf{A}$.

Having described deflation as a technique to accelerate GMRES by removing selected eigenvalues, we now turn to another complementary approach: preconditioning. In particular, we study the Complex Shifted Laplacian (CSL) preconditioner, which aims to reduce iteration counts by modifying the spectral properties of the system. We then explore how deflation and CSL preconditioning can be combined to further improve solver performance.

**CSL preconditioning.** We consider the CSL preconditioner as a representative example, following [22, 29]. The preconditioner is based on solving a modified Helmholtz problem with a complex-shifted wavenumber.

$$ \int_{\Omega} \nabla u \cdot \nabla v \, d\Omega - (k^2 + i\epsilon) \int_{\Omega} uv \, d\Omega = \int_{\Omega} fv \, d\Omega, \quad \forall v \in H^1(\Omega), \epsilon > 0. $$

As suggested in [22], we use $\epsilon = k$. Discretizing this problem yields the matrix $\mathbf{A}_{\epsilon}$, and the original system is then right-preconditioned as:

$$ (3.6) \qquad | \text{Find } \mathbf{x} \in \mathbb{C}^N \text{ such that } (\mathbf{A}\mathbf{A}_{\epsilon}^{-1})\mathbf{x} = \mathbf{b}, $$

with the solution recovered as $\mathbf{u} = \mathbf{A}_{\epsilon}^{-1}\mathbf{x}$.

The CSL preconditioner aims to reduce iteration counts without introducing excessive computational cost. Its assembly is inexpensive, since the mass and stiffness matrices are reused from the original Helmholtz discretization. However, applying $\mathbf{A}_{\epsilon}^{-1}$ can be as costly as solving the original system. To address this, we use an incomplete LU (ILU) factorization as an efficient approximation. Each GMRES iteration is then more expensive due to the cost of applying this approximate inverse.

**Deflation and preconditioning.** The deflation and CSL preconditioning techniques can be combined to enhance solver robustness and convergence. Substituting the preconditioned formulation (3.6) into the deflated system (3.3) leads to the doubly modified system:

$$ (3.7) \qquad | \text{Find } \tilde{\mathbf{x}} \in \mathbb{C}^N \text{ such that } (\mathbf{P}_{\text{def}}\mathbf{A}\mathbf{A}_{\epsilon}^{-1})\tilde{\mathbf{x}} = \mathbf{P}_{\text{def}}\mathbf{b}. $$

The final solution is reconstructed using

$$ \mathbf{u} = \mathbf{Q}\mathbf{b} + \mathbf{Q}_{\text{def}}\mathbf{A}_{\epsilon}^{-1}\tilde{\mathbf{x}}. $$

If the conditions

$$ \ker(\mathbf{Z}^*) \cap \text{range}(\mathbf{A}\mathbf{Z}) = \{0\} \quad \text{and} \quad \ker(\mathbf{Z}^*) \cap \text{range}(\mathbf{A}_{\epsilon}^{-1}\mathbf{Z}) = \{0\}, $$

hold, then GMRES applied to this system converges for all initial guesses [33, Theorem 2]. In practice, $\mathbf{A}_{\epsilon}^{-1}$ is approximated by its ILU factorization. Each GMRES iteration then involves the additional cost of deflation mentioned in Subsection 3.3 plus the cost of applying the preconditioner.

Alternatively, the deflation and preconditioning ideas can be combined additively, yielding the so-called *additive coarse grid correction* method [2]. This leads to the modified system:

$$ (3.8) \qquad | \text{Find } \mathbf{x} \in \mathbb{C}^N \text{ such that } \mathbf{A} (\mathbf{A}_{\epsilon}^{-1} + \mathbf{Q}) \mathbf{x} = \mathbf{b}, $$
<figure><img src="image_3.png" /><figcaption>Fig. 4.1: Cavity benchmark. Relative L²-error and number of GMRES iterations to reach the tolerance 10⁻⁶ on the relative residual as a function of k. The vertical lines correspond to resonance wavenumbers.</figcaption></figure>

with the solution given by $\mathbf{u} = (\mathbf{A}_{\varepsilon}^{-1} + \mathbf{Q})\mathbf{x}$. This approach has the same additional cost per iteration of (3.7), but it does not require the initial computation of $\mathbf{AZ}$ and $\mathbf{Z}^*\mathbf{A}$. However, since the numerical results obtained with this method are nearly identical to those with deflation, and because deflation enjoys a more developed theoretical foundation, we do not pursue this approach further in this work.

**4. Numerical experiments: a cavity problem with resonances.** In this section and the next, we present numerical experiments designed to illustrate the specific challenges of solving Helmholtz problems. These benchmarks serve to validate the interpretation tools and acceleration techniques discussed previously. The first set of experiments focuses on a simple cavity problem with a wavenumber chosen close to resonance, a setting known to lead to slow convergence.

**Description of the benchmark.** We solve the Helmholtz equation on the unit square domain $\Omega = (0,1)^2$, subject to homogeneous Dirichlet boundary conditions on $\partial\Omega$ and a constant source term equal to one in the interior. The problem admits a unique solution provided that the square of the wavenumber does not coincide with an eigenvalue of the Laplacian operator.

For this configuration, the Laplacian eigenvalues and eigenfunctions are given by

$$k_{n,m}^2 = \pi^2(n^2 + m^2), \quad u_{n,m}(x,y) = \sin(n\pi x)\sin(m\pi y), \quad n, m > 0.$$

After discretization, the system matrix $\mathbf{A}$ is real, symmetric, and therefore normal. Its eigenvalues are real, and it admits an orthonormal basis of eigenvectors.

All experiments are conducted using a Matlab implementation. We use $P_2$ (quadratic) finite elements on a structured mesh with characteristic size $h = 1/32$, leading to approximately $N_{\text{dof}} = 5,000$ degrees of freedom. The mesh resolution is chosen so that the number of degrees of freedom per wavelength satisfies $D_{\lambda} = (\sqrt{N_{\text{dof}} - 1})/\lambda \ge 30$, where $\lambda = k/2\pi$ denotes the wavelength.

**Influence of the wavenumber on the relative error and the number of iterations.** We now analyze the behavior of the solution and solver convergence as a function of the wavenumber $k$ for the cavity benchmark introduced above.

In Figure 4.1, we plot the relative error and the number of GMRES iterations required to reduce the relative residual below 10⁻⁶, as functions of the wavenumber $k$. The selected
<figure><img src="image_4.png" /><figcaption>Fig. 4.2: Cavity benchmark. Convergence analysis for $k = \hat{k} = 3.01\sqrt{2\pi}$ near the resonance wavenumber $k_{3,3}$. The history of the relative residual over the iterations is shown at the top, and the trajectory of the HR values is shown at the bottom. The horizontal green lines correspond to the eigenvalues of $\mathbf{A}$.</figcaption></figure>

range of $k$ includes two resonance values: $k_{3,3} = 3\sqrt{2\pi} \approx 13.33$ and $k_{4,1} = k_{1,4} \approx 12.95$, corresponding to eigenmodes of the Laplacian.

A pronounced peak in the relative error appears near $k_{3,3}$, while no such peak is observed near $k_{4,1}$ (see Figure 4.1a). This difference is explained by the structure of the source term: it excites only those modes for which both indices $n$ and $m$ are odd. As a result, modes such as $u_{4,1}$ and $u_{1,4}$, which are not excited, contribute negligibly to the solution, and no significant error is observed even when $k$ approaches their corresponding eigenvalues.

By contrast, the number of GMRES iterations is sensitive to the proximity of any eigenmode, regardless of whether it is excited by the source term. This is evident in Figure 4.1b, where two iteration peaks appear near both resonance values. This suggests that the convergence of GMRES deteriorates when the wavenumber is close to an eigenvalue, a behavior we explore in more detail in the following discussion.

**Influence of the wavenumber on the residual history.** To better understand the convergence behavior near resonance, we now examine the residual history of GMRES without acceleration for the wavenumber $\hat{k} = 3.01\sqrt{2\pi} \approx 13.37$, which lies close to the resonance $k_{3,3}$. At this frequency, one eigenvalue of the system matrix $\mathbf{A}$ lies close to the origin, and convergence issues are expected. As shown in Figure 4.2 (top), the residual exhibits several distinct plateau phases followed by sharp drops (highlighted by black circles labeled 1, 2, and 3). This pattern can be interpreted using the HR values and the convergence bound in Theorem 2.6.

Figure 4.2 (bottom) tracks the evolution of the HR values near zero throughout the GMRES iterations. These values, which are real in this case, converge toward the eigenvalues
of **A** (indicated by green lines), as expected. The first plateau in the residual history aligns with the moment when an HR value begins to approach the smallest eigenvalue $\lambda_1$ (circle 1 in both figures). To analyze this behavior, we apply Theorem 2.6 with $\Lambda_J^c = \sigma(\mathbf{A}) \setminus \{\lambda_1\}$ and examine the quantity

$$ (4.1) \qquad \max_{\lambda_i \in \Lambda_J^c} |s_J^l(\lambda_i)| = \max_{\lambda_i \in \Lambda_J^c} \left| \frac{(1 - \lambda_i / \lambda_1)}{(1 - \lambda_i / \nu_1^{(l)})} \right| = \max_{\lambda_i \in \Lambda_J^c} \left| \frac{\nu_1^{(l)}}{\lambda_1} \right| \left| \frac{\lambda_1 - \lambda_i}{\nu_1^{(l)} - \lambda_i} \right|, $$

which is a factor in the upper bound of the residual.

We identify three distinct phases in the evolution of the GMRES residual, labeled I, II, and III in both plots of Figure 4.2. Each phase reflects a different stage in the approximation of the smallest eigenvalue $\lambda_1$ by an HR value $\nu_1^{(l)}$, and corresponds to a different regime in the theoretical bound from Theorem 2.6.

* **Phase I: HR values are far from** $\lambda_1$. At this stage, $\nu_1^{(l)}$ is not yet close to $\lambda_1$, and both terms in (4.1) are large. The bound from Theorem 2.6 is then loose, and cannot provide an accurate description of the convergence rate.

* **Phase II: HR value approaches** $\lambda_1$. As the iteration progresses, one HR value $\nu_1^{(l)}$ begins to converge towards $\lambda_1$. However, the term $|\nu_1^{(l)} / \lambda_1|$ remains large, and the convergence bound still does not fully capture the residual behavior. This corresponds to the plateau in the residual curve, where the influence of $\lambda_1$ is strong but not yet neutralized.

* **Phase III: HR value accurately captures** $\lambda_1$. Once $\nu_1^{(l)} \approx \lambda_1$, we have $|s_J^l(\lambda_i)| \approx 1$ for all $\lambda_i \in \Lambda_J^c$. The impact of $\lambda_1$ is effectively removed from the bound, and the convergence behavior improves significantly. The residual drops, and the theoretical prediction becomes tight.

This analysis reveals a direct link between the evolution of HR values and the shape of the GMRES residual curve. Each plateau reflects the presence of a small eigenvalue not yet well captured by the Krylov subspace. Once a small eigenvalue is well approximated, its effect is factored out of the residual bound, and the solver progresses more efficiently. The same mechanism explains the subsequent plateaus and drops: each corresponds to another small eigenvalue of **A** being successively captured by HR values. The more small eigenvalues are well approximated, the smaller the quantity in the convergence bound of Theorem 2.6, and the more precise the theoretical prediction becomes.

In summary, GMRES convergence slows down when a small eigenvalue of **A** is in the process of being resolved by the Krylov subspace. Once the associated HR value converges, the eigenvalue's influence fades, and the residual decays more rapidly, as if the eigenvalue were no longer present.

**Influence of the deflation.** Having established that the plateaus in the residual history correspond to small eigenvalues of **A**, we now investigate the effect of deflating the associated eigenvectors on GMRES convergence. Recall that the eigenvectors of **A** approximate those of the Laplace operator with Dirichlet boundary conditions.

We begin by deflating the eigenvector associated with the smallest eigenvalue of **A**—corresponding to the mode $u_{3,3}$, illustrated in Figure 4.3. In preliminary tests (not shown), we compared two approaches: deflating the discrete eigenvector of **A** directly, and deflating the finite element projection of the continuous mode $u_{3,3}$. Both yielded nearly identical results on fine meshes. In the remainder of the analysis, we focus on the deflation of projected modes $u_{m,n}$.

Figure 4.4a shows the effect of deflating $u_{3,3}$ on the number of GMRES iterations as a function of the wavenumber $k$. Comparing the baseline case without deflation (green ◦) and the deflated case (green ×), we see that the peak near $k_{3,3}$ disappears when $u_{3,3}$ is deflated, while the peak near $k_{4,1} = k_{1,4}$ remains unaffected. Notably, the iteration count is reduced across the full range of $k$, even for values not directly associated with the deflated mode.

To better understand this effect, we examine the residual history for $k = \hat{k}$, a value close
<table><caption>Table 4.1: Cavity benchmark. Mode numbers and negative eigenvalues of $\mathbf{A}$ for $k = \hat{k}$.</caption><thead><tr><th>(n, m)</th><th>(3, 3)</th><th>(1, 4)<br>(4, 1)</th><th>(2, 3)<br>(3, 2)</th><th>(1, 3)<br>(3, 1)</th><th>(2, 2)</th><th>(1, 2)<br>(2, 1)</th><th>(1, 1)</th></tr></thead><tbody><tr><th>λ</th><td>-0.0010</td><td>-0.0092</td><td>-0.0425</td><td>-0.0676</td><td>-0.0846</td><td>-0.1101</td><td>-0.1359</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Fig. 4.3: Cavity benchmark. Mode (3,3) of the Laplacian.</figcaption></figure>

to $k_{3,3}$. As shown in Figure 4.4b, deflating $u_{3,3}$ eliminates the corresponding plateau in the residual curve. The remainder of the curve remains nearly unchanged, as deflation modifies only the contribution of the targeted eigenvector and does not affect the rest of the spectrum.

Finally, we test deflation of the negative part of the spectrum to remove the indefinite part of the problem. All the eleven negative eigenvalues correspond to Laplace modes whose wavenumbers are below $\hat{k}$ (Table 4.1). The resulting residual history (green □) shows that all plateaus present in the non-deflated case have disappeared, confirming the strong connection between small eigenvalues and the observed stagnation behavior.

**Influence of the preconditioning combined with the deflation.** Building on the previous observations, we now examine the effect of combining ILU preconditioning with deflation. Figure 4.4a shows that, for all wavenumbers, the use of ILU($\mathbf{A}_\epsilon$) as a preconditioner consistently reduces the number of GMRES iterations compared to the unpreconditioned case. When deflation of $u_{3,3}$ is added to the preconditioning strategy (pink ×), the iteration count decreases even further, and the small peak near $k_{3,3}$—previously associated with resonance—is eliminated.

This trend is confirmed by the residual histories in Figure 4.4c. While ILU preconditioning alone significantly shortens the plateaus and accelerates convergence after each stagnation phase, some residual plateaus still remain. However, deflating the eigenvector associated with the resonant mode $u_{3,3}$ removes the first plateau entirely. Moreover, when deflation is extended to the eleven eigenvectors corresponding to the first resonance modes, all plateaus vanish, and the residual decays rapidly in a nearly linear fashion.

**Influence of the restarted version of GMRES.** Having explored the effects of preconditioning and deflation in full GMRES, we now turn to a more practical variant commonly used in large-scale applications: restarted GMRES.

In the standard implementation of GMRES, the cost per iteration grows with the number of iterations, as each step requires storing and orthogonalizing against an expanding Krylov basis. To limit memory and computational overhead, GMRES is often restarted every $m$ iterations, retaining only the last Krylov vector to start the next cycle. This variant, known as GMRES($m$), reduces per-iteration cost but sacrifices theoretical guarantees of convergence.

Figure 4.4b and Figure 4.4c show the residual histories for GMRES($m$) with a restart parameter $m = 25$, in both unpreconditioned and ILU-preconditioned settings (blue and
<figure><img src="image_6.png" /><figcaption>Fig. 4.4: Cavity benchmark. Number of GMRES iteration to reach the relative residual $10^{-6}$, and residual history for $k = \hat{k}$ in different configurations: with/without deflation, with/without preconditioning, with/without restart.</figcaption></figure>
yellow curves, respectively). In both cases, and without deflation, convergence stagnates rapidly and the residual plateaus—even when preconditioning is applied. We observed similar behavior with larger restart values (e.g., $m = 80$, not shown), indicating that the issue is robust.

However, when deflation is introduced, GMRES($m$) recovers its ability to converge. In both the unpreconditioned and preconditioned cases, the residual drops more steadily, and the stagnation phases are significantly reduced or eliminated. This improvement can be explained as follows: in the restarted setting, the Krylov subspace may not have enough iterations to accurately approximate the small eigenvalues responsible for stagnation. Deflation explicitly removes these difficult modes, making GMRES($m$) more effective. This highlights an important practical benefit: deflation enhances the robustness of restarted GMRES, particularly near resonant wavenumbers where convergence is otherwise problematic.

**5. Numerical experiments: scattering problem with quasi-resonances.** We now turn to a more challenging configuration: the diffraction of a plane wave by an open cavity. Unlike the previous benchmark, where the system matrix was real, symmetric, and associated with exact resonance modes, the matrix $A$ here is complex and non-normal. Although exact resonances are absent, the problem features *quasimodes*, which can lead to similar convergence difficulties for iterative solvers. As before, we use numerical experiments to study GMRES convergence and the effect of acceleration techniques.

**Description of the benchmark.** The problem consists of the scattering of a time-harmonic plane wave $u_{\text{inc}}(x, y) = e^{ik(\cos(\theta)x + \sin(\theta)y)}$ by an obstacle $\mathcal{O}$ in free space $\mathbb{R}^2$, where $\theta$ is the incident angle. The boundary of the obstacle is assigned either a Dirichlet condition (Dirichlet case) or a Neumann condition (Neumann case).

To model the unbounded domain, we truncate the computational domain to a square region $\Omega_{\text{dom}} = (-L, L)^2 \setminus \mathcal{O}$ and surround it with perfectly matched layers (PMLs) of thickness $L_{\text{pml}}$, denoted $\Omega_{\text{pml}}$ [6, 7]. The PMLs absorb outgoing waves and eliminate reflections from the artificial boundary, on which a homogeneous Dirichlet condition is imposed. The PML formulation introduces complex coordinate stretchings:

$$ \gamma_x(x) := 1 + \frac{i\sigma_x(x)}{k}, \quad \gamma_y(y) := 1 + \frac{i\sigma_y(y)}{k}, $$

with absorption profiles:

$$ \sigma_x(x) := \frac{1_{\text{pml}}(x)}{L_{\text{pml}} - |x| + L}, \quad \sigma_y(y) := \frac{1_{\text{pml}}(y)}{L_{\text{pml}} - |y| + L}. $$

The scattered field $u$ then satisfies:

$$ \left\{ \begin{array}{ll} \displaystyle \frac{\partial}{\partial x} \left( \frac{\gamma_y}{\gamma_x} \frac{\partial u}{\partial x} \right) + \frac{\partial}{\partial y} \left( \frac{\gamma_x}{\gamma_y} \frac{\partial u}{\partial y} \right) + \gamma_x \gamma_y k^2 u = 0 & \text{in } \Omega_{\text{dom}} \cup \Omega_{\text{pml}}, \\ u = -u_{\text{inc}} & \text{or } \partial_n u = -\partial_n u_{\text{inc}} \quad \text{on } \Gamma_{\text{obs}}, \\ u = 0 & \text{on } \Gamma_{\text{ext}}. \end{array} \right. $$

The obstacle $\mathcal{O}$ is an open rectangular cavity of length $L_{\mathcal{O}}$ and opening width $l_{\mathcal{O}}$ (see Figure 5.1). In contrast to closed domains, this setup does not admit true resonance modes but supports *quasimodes*, i.e., modes that approximately satisfy the homogeneous equation and are spatially localized.

In the Dirichlet case, *quasimodes* are known to accumulate near wavenumbers $\pi n/l_{\mathcal{O}}$ for $n \in \mathbb{N}$ [11, Theorem 5.1]. However, numerical simulations indicate the presence of additional *quasimodes* at other wavenumbers. These are often close to the resonance frequencies of a closed rectangular cavity with the same dimensions. In that case, the Dirichlet resonance wavenumbers are given by

$$ k_{n,m} = \pi \sqrt{\frac{m^2}{L_{\mathcal{O}}^2} + \frac{n^2}{l_{\mathcal{O}}^2}}, \quad m, n \in \mathbb{N}. $$
<figure><img src="image_7.png" /><figcaption>Fig. 5.1: Scattering benchmark. Computational domain with PML.</figcaption></figure>

For the Neumann case (used in our experiments), the expected quasimodes are close to the resonance frequencies of a closed rectangular cavity with Neumann boundary conditions on all sides except the left edge (Dirichlet), given by

$$k_{n,m} = \pi \sqrt{\frac{(m + 1/2)^2}{L_O^2} + \frac{n^2}{l_O^2}}, \quad m, n > 0.$$

In both cases, the quasimodes inside the cavity resemble the eigenmodes of the corresponding closed configuration.

After discretization using finite elements, the problem leads to a linear system $Au = b$. The resulting matrix $A$ is complex, non-normal, and has a spectrum that includes complex-valued eigenvalues.

All numerical results presented in this section correspond to the Neumann case. The Dirichlet case shows similar behavior and is omitted for brevity. We use the following parameters $L_O = 1.3$, $l_O = 0.4$, $\theta = \frac{4\pi}{10}$, $h = \frac{1}{20}$, $P_3$ elements, $\text{tol} = 10^{-6}$, $D_\lambda \ge 30$. This configuration leads to approximately 13,500 degrees of freedom.

**Influence of the wavenumber on the number of iterations.** We now investigate how the wavenumber $k$ affects the convergence of GMRES for the open cavity problem. Figure 5.2a shows the number of GMRES iterations required to reach a residual tolerance of $10^{-6}$ as a function of $k \in [23.5, 24.5]$ (green ○).

As in the closed cavity benchmark, we observe peaks in the iteration count near the expected wavenumbers of quasimodes, specifically around $k_{0,3} \approx 23.591$, $k_{1,3}$, and $k_{2,3}$. However, the peaks here are less pronounced than in the cavity benchmark, reflecting the fact that quasimodes—unlike true resonance modes—do not make the problem singular, but still degrade the performance of iterative solvers. Among the peaks, the one near $k_{0,3}$ is the most prominent. To better understand the behavior of GMRES in this regime, we focus in the following on the wavenumber $\hat{k} := 23.591$.

**Influence of the wavenumber on the residual history.** We now examine the convergence behavior near the wavenumber $\hat{k} = 23.591$, where the influence of a quasimode was observed in the iteration count.

Figure 5.2b shows the GMRES residual history without any acceleration technique. Compared to the cavity problem, convergence is significantly slower. Only a single plateau is observed, between iterations 800 and 1050, after which the residual begins to decrease more rapidly. As in the previous benchmark, this stagnation can be explained using HR values and Theorem 2.6.

Since the matrix $A$ has a complex spectrum, we examine the distribution of HR values in the complex plane. Figure 5.3 (left) compares the positions of HR values at iteration 850—during the plateau—and iteration 1100—just after the plateau ends. At iteration 850,
<figure><img src="image_8.png" /><figcaption>Fig. 5.2: Scattering benchmark. Number of GMRES iteration to reach the relative residual $10^{-6}$, and residual history for $k = \hat{k}$ in different configurations: with/without deflation, with/without preconditioning, with/without restart.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 5.3: Scattering benchmark. HR values and deflated spectrum.</figcaption></figure>

we observe that an HR value is approaching the smallest eigenvalue of **A**, located in the lower left region of the complex plane. By iteration 1100, this eigenvalue is well approached, marking the end of the plateau and the onset of faster convergence.

Although the spectrum is now complex, the convergence behavior remains consistent with the interpretation provided by Theorem 2.6. In this benchmark, the slower convergence of HR values towards the smallest eigenvalue explains the extended duration of the plateau phase.

**Influence of the deflation.** To mitigate the slow convergence near the quasimode, we now assess the impact of deflation on the solver's performance.

Although no analytical expression is available for the quasimodes in this configuration, it is possible to construct approximate deflation vectors based on physical intuition. We use eigenvectors of the corresponding closed cavity problem—with Dirichlet conditions on the left edge and Neumann conditions on the remaining edges—extended by zero outside the cavity.

Figure 5.3 (right) displays the spectrum of the deflated matrix $\mathbf{P}_{\text{def}}\mathbf{A}$. The eigenvalue associated with the quasimode is reduced to zero, confirming that the deflated vector effectively captures the problematic mode.

The influence of this deflation on convergence is shown in Figure 5.2a (green ×). As expected, the peak near $k_{0,3}$—previously linked to the quasimode—is entirely removed. Moreover, the number of GMRES iterations is slightly reduced across other values of $k$, even away from the targeted quasimode.

To explore this further, we compare the residual histories with and without deflation at $k = \hat{k}$, near $k_{0,3}$ (Figure 5.2b). The plateau observed in the unaccelerated case disappears with deflation. However, between iterations 600 and 800, the residual with deflation is briefly higher than in the non-deflated case. This is likely due to the fact that the deflated vector only approximates the true quasimode.

We also test the deflation of multiple approximate quasimodes (green □), but the improvement in convergence is modest, suggesting diminishing returns beyond the first key mode.

**Influence of preconditioning combined with deflation.** We now assess how ILU preconditioning interacts with deflation in improving GMRES convergence.

As shown in Figure 5.2a, applying the ILU($\mathbf{A}_{\varepsilon}$) preconditioner systematically reduces
the number of iterations required to reach the convergence tolerance, across all wavenumbers. Although the benefit of combining preconditioning with deflation is not visually striking at the scale of the plot, the effect is measurable. At $k = \hat{k}$, for example, the number of iterations drops from 350 (preconditioning only) to 319 when deflation is added—an improvement of approximately 9%.

The residual histories in Figure 5.2c provide a more detailed view. With preconditioning alone (pink ○), the plateau associated with the quasimode persists, though it is shorter and followed by a sharp decline in the residual. When deflation of the approximate quasimode is added, the plateau is eliminated entirely, and convergence becomes both faster and smoother.

These results confirm that $ILU(A_{\varepsilon})$ is effective in mitigating spectral clustering effects and improving convergence. Nonetheless, the combination of preconditioning and deflation remains the most effective strategy, even when the deflated vectors are only approximations of the true quasimodes.

**Influence of the restarted version of GMRES.** We now turn to the restarted GMRES($m$), to evaluate its behavior in this more challenging setting. As with the cavity problem, restarting is used to reduce the memory footprint and per-iteration cost. Figure 5.2b (for $m = 25$) and Figure 5.2c (for $m = 250$) display the residual histories with and without preconditioning.

Without deflation, we note that GMRES($m$) fails to converge—regardless of preconditioning—even when the restart value is increased (e.g., $m = 300$, not shown). This highlights the difficulty of resolving the quasimode and overcoming the effects of non-normality and indefiniteness with restarts alone.

In contrast to the cavity benchmark, deflation alone (blue × and □) is not sufficient to guarantee convergence within a reasonable number of iterations. This can be attributed to the compounded challenges of non-normality, indefiniteness, and the complexity of the quasimodal structure in this benchmark.

Only the combination of deflation and preconditioning (yellow ×) yields successful convergence, though a high restart value is required: for $m \le 250$, the tolerance is not reached within 2000 iterations. However, increasing the number of deflated vectors (yellow □) improves convergence significantly, making it possible to reduce the restart parameter while still achieving convergence (detailed results not shown here). These results underscore the importance of using both preconditioning and deflation in restarted settings, especially when dealing with non-normal systems influenced by quasimodes.

**6. Conclusion.** This study examined the convergence behavior of GMRES for the finite element solution of Helmholtz problems, emphasizing the impact of resonance and quasimodes.

Novel convergence bounds based on harmonic Ritz values have been proposed for general systems. These bounds are able to describe nonlinear residual histories, offering a detailed explanation of residual plateaus. Through numerical experiments on both closed and open cavity benchmarks, we observed that convergence stalls align with wavenumbers associated with resonances or quasimodes. These stalls were accurately predicted by the evolution of harmonic Ritz values.

The main observations of the numerical studies can be summarized as follows:

* GMRES convergence is highly slowed and plateau because of small eigenvalues due to (quasi-)resonances.

* Deflation of the eigenvectors corresponding to (quasi-)resonances, even when approximated, removes plateaus and accelerates convergence. This approach is robust, but requires knowledge of the eigenvectors and slightly increases the cost per iteration, which depends on the number of the deflated eigenvectors.

* $ILU(A_{\varepsilon})$ preconditioning enhanced performance further by addressing non-normality and indefiniteness, particularly in open domain problems, at the cost of an additional matrix product. Although deflation and preconditioning contributed individually, their combination proved most effective, ensuring fast and stable convergence across
all regimes tested.

* Restarted GMRES cannot achieve convergence near (quasi-)resonances, and only a combination of deflation and ILU($A_{\epsilon}$) preconditioning restored it in the most difficult cases.

These results underscore the value of combining spectral analysis with tailored numerical techniques to design robust solvers for high-frequency wave propagation problems. Future work includes extending these strategies to 3D configurations and exploring more scalable preconditioning techniques based on domain decomposition.

REFERENCES

[1] I. M. BABUŠKA AND S. A. SAUTER, *Is the Pollution Effect of the FEM Avoidable for the Helmholtz Equation Considering High Wave Numbers?*, SIAM Journal on Numerical Analysis, 34 (1997), pp. 2392–2423, https://doi.org/10.1137/s0036142994269186.

[2] P. BASTIAN, G. WITTUM, AND W. HACKBUSCH, *Additive and multiplicative multi-grid — A comparison*, Computing, 60 (1998), pp. 345–364, https://doi.org/10.1007/bf02684380.

[3] B. BECKERMANN, S. A. GOREINOV, AND E. E. TYRTYSHNIKOV, *Some Remarks on the Elman Estimate for GMRES*, SIAM Journal on Matrix Analysis and Applications, 27 (2005), pp. 772–778, https://doi.org/10.1137/040618849.

[4] J.-P. BERENGER, *A perfectly matched layer for the absorption of electromagnetic waves*, Journal of Computational Physics, 114 (1994), pp. 185–200, https://doi.org/10.1006/jcph.1994.1159.

[5] H. BÉRIOT AND A. MODAVE, *An automatic perfectly matched layer for acoustic finite element simulations in convex domains of general shape*, International Journal for Numerical Methods in Engineering, (2020), https://doi.org/10.1002/nme.6560.

[6] A. BERMÜDEZ, L. HERVELLA-NIETO, A. PRIETO, AND R. RODRÍGUEZ, *An exact bounded PML for the Helmholtz equation*, Comptes Rendus Mathematique, 339 (2004), pp. 803–808, https://doi.org/10.1016/j.crma.2004.10.006.

[7] A. BERMÜDEZ, L. HERVELLA-NIETO, A. PRIETO, AND R. RODRÍGUEZ, *An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems*, Journal of Computational Physics, 223 (2007), pp. 469–488, https://doi.org/10.1016/j.jcp.2006.09.018.

[8] E. BÉCACHE, A. S. BONNET-BEN DHIA, AND G. LEGENDRE, *Perfectly Matched Layers for Time-Harmonic Acoustics in the Presence of a Uniform Flow*, SIAM Journal on Numerical Analysis, 44 (2006), pp. 1191–1217, https://doi.org/10.1137/040617741.

[9] S. L. CAMPBELL, I. C. F. IPSEN, C. T. KELLEY, AND C. D. MEYER, *GMRES and the minimal polynomial*, BIT Numerical Mathematics, 36 (1996), pp. 664–675, https://doi.org/10.1007/bf01733786.

[10] Z.-H. CAO, *A note on the convergence behavior of GMRES*, Applied Numerical Mathematics, 25 (1997), pp. 13–20, https://doi.org/10.1016/s0168-9274(97)00039-1.

[11] S. N. CHANDLER-WILDE, I. G. GRAHAM, S. LANGDON, AND M. LINDNER, *Condition number estimates for combined potential boundary integral operators in acoustic scattering*, Journal of Integral Equations and Applications, 21 (2009), https://doi.org/10.1216/jie-2009-21-2-229.

[12] M. CROUZEIX AND C. PALENCIA, *The Numerical Range is a (1 + √2)-Spectral Set*, SIAM Journal on Matrix Analysis and Applications, 38 (2017), pp. 649–655, https://doi.org/10.1137/17m1116672.

[13] H. V. DER VORST AND C. VUIK, *The superlinear convergence behaviour of GMRES*, Journal of Computational and Applied Mathematics, 48 (1993), pp. 327–341, https://doi.org/10.1016/0377-0427(93)90028-a.

[14] V. DOLEAN, P. JOLIVET, AND F. NATAF, *An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation*, Society for Industrial and Applied Mathematics, Nov. 2015, https://doi.org/10.1137/1.9781611974065.

[15] M. EMBREE, *How Descriptive are GMRES Convergence Bounds?*, 2022, https://doi.org/10.48550/ARXIV.2209.01231.

[16] M. EMBREE, *Extending Elman's Bound for GMRES*, 2023, https://doi.org/10.48550/ARXIV.2312.15022.

[17] B. ENGQUIST AND A. MAJDA, *Absorbing boundary conditions for numerical simulation of waves*, Proceedings of the National Academy of Sciences, 74 (1977), pp. 1765–1766, https://doi.org/10.1073/pnas.74.5.1765.

[18] Y. A. ERLANGGA, *Advances in Iterative Methods and Preconditioners for the Helmholtz Equation*, Archives of Computational Methods in Engineering, 15 (2007), pp. 37–66, https://doi.org/10.1007/s11831-007-9013-7.

[19] Y. A. ERLANGGA AND R. NABBEN, *Deflation and Balancing Preconditioners for Krylov Subspace Methods Applied to Nonsymmetric Matrices*, SIAM Journal on Matrix Analysis and Applications, 30 (2008), pp. 684–699, https://doi.org/10.1137/060678257.

[20] O. G. ERNST AND M. J. GANDER, *Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods*, Springer Berlin Heidelberg, Aug. 2012, pp. 325–363, https://doi.org/10.1007/978-3-642-22061-6_10.

[21] J. GALKOWSKI, P. MARCHAND, AND E. A. SPENCE, *Eigenvalues of the Truncated Helmholtz Solution Operator*, 2011, https://doi.org/10.1093/soc/dab037.
erator under Strong Trapping, SIAM Journal on Mathematical Analysis, 53 (2021), pp. 6724-6770,
https://doi.org/10.1137/21m1399658.

[22] M. J. GANDER, I. G. GRAHAM, AND E. A. SPENCE, *Applying GMRES to the Helmholtz equation with shifted Laplacian preconditioning: what is the largest shift for which wavenumber-independent convergence is guaranteed?*, Numerische Mathematik, 131 (2015), pp. 567-614, https://doi.org/10.1007/s00211-015-0700-2.

[23] L. GARCÍA RAMOS, R. KEHL, AND R. NABBEN, *Projections, Deflation, and Multigrid for Nonsymmetric Matrices*, SIAM Journal on Matrix Analysis and Applications, 41 (2020), pp. 83-105, https://doi.org/10.1137/18m1180268.

[24] P. MARCHAND, J. GALKOWSKI, E. A. SPENCE, AND A. SPENCE, *Applying GMRES to the Helmholtz equation with strong trapping: how does the number of iterations depend on the frequency?*, Advances in Computational Mathematics, 48 (2022), https://doi.org/10.1007/s10444-022-09931-9.

[25] P. MARCHNER, H. BÉRIOT, X. ANTOINE, AND C. GEUZAINE, *Stable Perfectly Matched Layers with Lorentz transformation for the convected Helmholtz equation*, Journal of Computational Physics, 433 (2021), p. 110180, https://doi.org/10.1016/j.jcp.2021.110180.

[26] G. A. MEURANT, *Computer solution of large linear systems*, vol. 28 of Studies in mathematics and its applications, Elsevier, 1999.

[27] A. MODAVE, C. GEUZAINE, AND X. ANTOINE, *Corner treatments for high-order local absorbing boundary conditions in high-frequency acoustic scattering*, Journal of Computational Physics, 401 (2020), p. 109029, https://doi.org/10.1016/j.jcp.2019.109029.

[28] L. G. RAMOS AND R. NABBEN, *A two-level shifted Laplace preconditioner for Helmholtz problems: Field-of-values analysis and wavenumber-independent convergence*, 2020, https://doi.org/10.48550/ARXIV.2006.08750.

[29] L. G. RAMOS, O. SÊTE, AND R. NABBEN, *Preconditioning the Helmholtz equation with the shifted Laplacian and Faber polynomials*, ETNA - Electronic Transactions on Numerical Analysis, 54 (2021), pp. 534-557, https://doi.org/10.1553/etna_vol54s534.

[30] Y. SAAD, *Iterative Methods for Sparse Linear Systems*, Society for Industrial and Applied Mathematics, 2003, https://doi.org/10.1137/1.9780898718003.

[31] Y. SAAD, *Numerical Methods for Large Eigenvalue Problems*, Society for Industrial and Applied Mathematics, 2011, https://doi.org/10.1137/1.9781611970739.

[32] Y. SAAD AND M. H. SCHULTZ, *GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems*, SIAM Journal on Scientific and Statistical Computing, 7 (1986), pp. 856-869, https://doi.org/10.1137/0907058.

[33] N. SPILLANE AND D. B. SZYLD, *New Convergence Analysis of GMRES with Weighted Norms, Preconditioning and Deflation, Leading to a New Deflation Space*, 2023, https://doi.org/10.48550/ARXIV.2312.13625.

[34] L. N. TREFETHEN AND M. EMBREE, *Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators*, Princeton University Press, 2005, https://doi.org/10.1515/9780691213101.

[35] E. TURKEL AND A. YEFET, *Absorbing PML boundary layers for wave-like equations*, Applied Numerical Mathematics, 27 (1998), pp. 533-557, https://doi.org/10.1016/s0168-9274(98)00026-9.

[36] A. VAN DER SLUIS AND H. A. VAN DER VORST, *The rate of convergence of Conjugate Gradients*, Numerische Mathematik, 48 (1986), pp. 543-560, https://doi.org/10.1007/bf01389450.

[37] C. VUIK, *Krylov Subspace Solvers and Preconditioners*, ESAIM: Proceedings and Surveys, 63 (2018), pp. 1-43, https://doi.org/10.1051/proc/201863001.