Multi-resolution Localized Orthogonal Decomposition for Helmholtz problems

We introduce a novel multi-resolution Localized Orthogonal Decomposition (LOD) for time-harmonic acoustic scattering problems that can be modeled by the Helmholtz equation. The method merges the concepts of LOD and operator-adapted wavelets (gamblets) and proves its applicability for a class of complex-valued, non-hermitian and indefinite problems. It computes hierarchical bases that block-diagonalize the Helmholtz operator and thereby decouples the discretization scales. Sparsity is preserved by a novel localization strategy that improves stability properties even in the elliptic case. We present a rigorous stability and a-priori error analysis of the proposed method for homogeneous media. In addition, we investigate the fast solvability of the blocks by a standard iterative method. A sequence of numerical experiments illustrates the sharpness of the theoretical findings and demonstrates the applicability to scattering problems in heterogeneous media.


Introduction
The concept of Localized Orthogonal Decomposition (LOD) has first been introduced in [MP14,HP13] for elliptic problems. It is an approximately orthogonal (w.r.t. the energy inner product) decomposition of the energy space into a finite-dimensional, problemadapted, mesh-based space with locally supported basis functions and an infinitedimensional remainder space. Using the finite-dimensional space as ansatz space for a Galerkin method (and allowing a moderate increase of the support of the basis functions), yields an optimally convergent multi-scale method that is capable of approximating solutions corresponding to arbitrarily rough L ∞ -coefficients, without any preasymptotic effects.
This method has also proved to posses stabilizing effects for high frequency wave scattering applications [GP15,Pet16,Pet17,LPS18]. For Helmholtz problems, the LOD method yields faithful approximations under the minimal resolution condition that the mesh size is coupled linearly to the wave number κ, and, if additionally, the support of the basis functions is allowed to increase logarithmically with κ. Its intrinsic ability to deal with high heterogeneity makes the LOD very appealing for such Helmholtz problems; see [BGP17,PV20,MV20].
Recently, the numerical homogenization approach has been extended to multi-resolution approximation schemes that allow for a decoupling into a hierarchy of discrete scales [Owh17]. This concept has been popularized under the name gamblets, and, up to now, has mostly been studied for the homogenization of elliptic operators [OZ17,OS19].
By bridging the gamblet approach with the well-developed LOD framework, we derive a multi-resolution (multi-level) LOD method that is applicable to a wide class of (complex-valued) non-hermitian and indefinite problems. As problem specific parameters make an unified (parameter-explicit) analysis nearly impossible, we henceforth consider the Helmholtz problem for demonstration purposes. Computing two sets of hierarchical bases for trial and test spaces, respectively, allows to block-diagonalize the operator. Note that, in contrast to the elliptic case, one hierarchical basis for both, trial and test space, does not suffice.
This paper features a κ-explicit stability and a-priori error analysis of the proposed multi-resolution method. For the elliptic case, gamblets suffer from severe instability, whenever the mesh size of the coarsest level is decreased; see also [Mai20]. This is especially problematic for Helmholtz problems, where the mesh size of the coarsest level is coupled to κ (for numerical illustrations, see Section 8.2). Using a novel basis construction enables us to cure this stability issue. The novel construction yields significant improvements also for single-level LOD methods as in [Mai20,Mai21] for elliptic problems. We prove optimal rates of convergence of the proposed method, if a moderate increase of the support of the basis functions is allowed.
The proposed multi-resolution method condenses the negative features of the linear system of equations arising from the discretization of the Helmholtz problem (indefinite, deteriorating condition for increasing κ) to the first (small) block. For the remaining blocks, we prove that they are coercive and can be solved within a fixed number of GMRES iterations, independent of levels, mesh sizes, and κ. This also inspires multi-level solver techniques for the Helmholtz equation, similarly as in [OS19] for elliptic problems.
For the sake of simplicity, the analysis covers only Helmholtz problems in homogeneous media. Nevertheless, the theoretical results carry over also to the case of heterogeneous media (for numerical experiments, see Section 8.3). However, the analytical well-posedness of the heterogeneous Helmholtz problem is problematic. In recent years, a lot progress has been made in this field; see [ST18,GPS19,GS19]. We believe that the results of the paper are also relevant for the simulation of other time-harmonic wave phenomena such as [BG16,GHV18,MV20].
The remaining part of the paper is outlined as follows. Section 2 defines the Helmholtz problem and recalls some of its fundamental properties. Section 3 introduces the hierarchy of meshes, the Haar basis and corresponding mesh-based operators that will be the basis for the derivation of a prototypical multi-resolution method in Section 4. Sections 5 and 6 will then turn this ideal approach into a feasible practical method including a rigorous stability and a-priori error analysis. In Section 7, we prove that all blocks besides the first one can be inverted easily. Finally, Section 8 illustrates the performance of the method in a sequence of numerical experiments.

Model Helmholtz problem
We consider the Helmholtz problem on a bounded Lipschitz domain D ⊂ R d , (d = 1, 2, 3), which is scaled to unit size with boundary conditions of Dirichlet, Neumann, and Robin type Here, κ ∈ (κ 0 , ∞) with κ 0 > 0 denotes the wave number, i is the imaginary unit and ν is the outer unit normal vector. The right-hand side f is assumed to be in L 2 (D) (the space of complex-valued square-integrable functions on D). It is also possible to use other types of boundary conditions like perfectly matched layers; see [CFGNT18]. We assume a decomposition of the boundary into closed components where the intersection of the interior of the components is supposed to be pairwise disjoint. We allow Γ D and Γ N to be empty, however, for Γ R , we suppose a positive surface measure, i.e., Given the Sobolev space W 1,2 (D) (the space of complex-valued square-integrable functions on D with square-integrable weak derivatives), we define the subspace which is endowed with the usual κ-dependent norm . The weak formulation of (2.1) seeks u ∈ V such that for all v ∈ V with the sesquilinear form a : Here, (u , v) L 2 (D) := D u · v dx ( · is the complex conjugation) denotes the inner product of scalar or vector-valued functions in L 2 (D). Similarly, (u , v) L 2 (Γ R ) := Γ R uv ds (integration with respect to the surface measure) is the inner product of the space L 2 (Γ R ). The sesquilinear form a is continuous, i.e., with a generic constant c > 0 independent of κ. In the following, we denote the real and the imaginary parts of a complex number by R and I, respectively. The well-posedness of weak formulation (2.3) depends on the shape of D and the choice of boundary conditions. The presence of Robin boundary conditions (c.f. (2.2)) ensures unique solvability, i.e., for all f ∈ L 2 (D), the solution u ∈ V satisfies For the stability and a-priori error analysis in Section 6, we assume a polynomial κdependence of c stab (c.f. Assumption 6.4), i.e., there exists a constant c > 0 and n ∈ N such that (2.6) c stab ≤ cκ n .
Clearly, condition (2.6) does not hold in general; see counter-examples with an at least exponential-in-κ growth of c stab for trapping domains [BCWG + 11]. Nevertheless, (2.6) can be proven under certain geometric assumptions. In [Mel95], condition (2.6) is proved with n = 0 for bounded star-shaped domains with smooth boundary or bounded convex domains. Among the known admissible setups is also the case of pure Robin boundary conditions (Γ R = ∂D) on Lipschitz domains [EM11]. Another possible example are truncated exterior Dirichlet problems, where the Sommerfeld radiation condition is approximated by truncating the (unbounded) exterior domain and applying Robin boundary conditions on the artificial boundary [Het07].
Remark 2.1 (Variable coefficients). In [ST18, GPS19, GS19] the well-posedness of the Helmholtz problem in heterogeneous media is analyzed; see (8.2) for the heterogeneous problem.

Preliminaries
This section lays the foundations for the construction of the multi-resolution LOD in the following sections.

Mesh hierarchy and Haar basis.
In what follows, we define a Haar basis given a sequence of nested meshes. The Haar construction is an essential ingredient for the construction of the proposed method as it induces its multi-resolution structure. Using the orthogonality properties of the Haar basis, a decoupling of discretization scales can be achieved.
Let T 1 denote a (coarse) Cartesian mesh of D with mesh size H 1 and let {T } =1,...,L , L ∈ N be a hierarchy of meshes obtained by successive uniform refinements of T 1 . The mesh size at level is given by H := 2 − +1 H 1 . The patch of a union of elements S in T is defined as In the setting of a Cartesian mesh hierarchy, the corresponding Haar wavelet basis can be given explicitly. Denoting with 1 S the characteristic function of the set S, we define The Haar basis can be written as H := =1,...,L H with H 1 := {|T | −1/2 1 T , T ∈ T 1 } and Here, a T,k denotes the k-th coordinate of the bottom, front, left corner of T ∈ T and 0 is the d-dimensional zero vector. It is straightforward that =1,...,L H is an othonormal basis of Q 0 (T L ). Henceforth, let φ ,j , j ∈ {1, ..., N } with N := #H denote the elements of H .
and a T,k as above (note that 1 0 b(x) dx = 1). We define the linear mapping Π : Q 0 (T ) → V uniquely by setting Π 1 T := b T for T ∈ T . It maps Q 0 (T )-functions to conforming bubble functions with the same element averages. The operator is the right inverse of Π , where Π : L 2 (D) → Q 0 (T ) denotes the L 2orthogonal projection onto the space of element-wise constants with respect to T . The operators Π and Π satisfy for all T ∈ T where c inv > 0 is the constant from the inverse inequality (see e.g. [DPE12, Lemma 1.44]).
Remark 3.1 (Construction of Haar basis and bubbles). The construction of the Haar basis and the bubbles is not restricted to Cartesian meshes. It can be performed for non-structured quadrilateral/hexahedral or simplicial meshes; see [Alp93,FP20]. If one considers tensor-product elements, higher order versions can be constructed; see [Alp93,Mai20].
Remark 3.2 (Tilde notation). Throughout the article, we write a b in short for a ≤ cb with a constant c independent of κ, , H , L, and m ( is defined analogously). The notation a b is used if a ≤ cb and b ≤ c a with constants c and c independent of κ, , H , L, and m.
3.3. V-stable projection. The V-stable projection introduced in the following, plays a crucial role in the construction of the localized ansatz spaces in Section 6. The construction of the operator P is based on the quasi-interpolation operator I , from which it inherits the V-stability. The quasi-interpolation is defined as I := E • Π , where E : Q 0 (T ) → Q 1 (T ) ∩ V (Q 1 (T ) is the space of element-wise tensor-product polynomials of degree up to one) denotes the nodal averaging operator, i.e., for an interior node z, let with T (z) := {T ∈ T : z ∈ T }. For boundary nodes, let (E v)(z) := 0. Since P should satisfy ker P = ker Π , we need to correct I v such that element averages of v on T are preserved: The following stability and approximation properties are shown in [AHP21] and are based on classical finite element theory [BS08,EG17]. For any v ∈ V (3.4)

Block-diagonal multi-resolution decomposition
In this section, we introduce a hierarchical multi-resolution Petrov-Galerkin method in order to discretize (2.3). In contrast to [Owh17,FP20], where only symmetric and elliptic problems are considered, we need two different sets of bases for trial and test spaces. This is due to the lack of hermiticity of the Helmholtz problem.
We propose a special choice of bases, which block-diagonalizes the (discretized) operator, see [FP20]. Each basis function is assigned to (exactly) one Haar basis function φ ,j . Since the basis functions, in general, have global support (c.f. Figure 4.1), the ideal method is not suited for practical purposes. A practical method is derived in Section 6 by localizing the basis functions to element patches. This procedure is justified by the exponential decay of the basis functions (Section 5).
4.1. Correction operators. The concept of correction operators is key in the LOD framework. In the multi-level context, one considers correction operators on every level projecting onto W := ker Π . The spaces W can be interpreted as (infinite-dimensional) subspaces of V consisting of functions oscillating on length scales smaller than H . The fine scale spaces are nested, i.e., For the correction operators to be well-defined, we need to assume a minimal resolution condition.
Assumption 4.1 (Resolution condition). Given the wave number κ, we assume that the mesh size H 1 of T 1 satisfies Due to the lack of hermiticity, we need to define for = 1, ..., L the two correction operators C and C * as The well-posedness of C and C * follows from the resolution condition.
Proof. a is coercive on W 1 × W 1 , since using Assumption 4.1. The norm equivalence of ∇· L 2 (D) and · V(D) follows from using again the same assumption.
Remark 4.3 (Well-posedness of correctors). The well-posedness of (4.1) can be concluded using the Lax-Milgram lemma using the coercivity of a on W 1 × W 1 shown in Lemma 4.2 and the continuity of a (2.4). In general, it holds that the complementary projection and the projection itself have the same operator norms (c.f. [Szy06]). Thus, we get boundedness of C , C * , 1 − C , and 1 − C * all with the same constants.

4.2.
Trial and test spaces. Using the correction operators, we can define the hierarchical multi-resolution trial and test spaces of the method. For levels = 1, ..., L, let The canonical bases of Φ and Ψ are given by {b ,j := (1 − C ) Π φ ,j , j = 1, ..., N } and {b * ,j := (1 − C * ) Π φ ,j , j = 1, ..., N }, respectively. Let us define the spaces U : Remark 4.4 (Conforming version of H ). Note that Π H is a conformal version of H . It holds that Π φ ,j and φ ,j have the same element-averages on T . The precise choice of Π is not important as long as its image is conforming, it preserves element-averages, and (3.2) holds.
Remark 4.5 (Relationship to [AHP21]). In [AHP21], the LOD method is introduced using a very general framework. So called 'quantities of interest' determine the precise method.
Using the element-averages on the mesh T L as quantities of interest, we can prove the relationships are the trial and test spaces from [AHP21], respectively. This relationship, in general, is not true for the practical method derived in Section 6. From Definition (4.1), one can deduce for k, ∈ {1, ..., L} with k ≤ the 'orthogonality' (a is not an inner-product) relations The single-level method discretizing (2.3) is based on the trial-test pairing ( U, V) and seeksũ ∈ U, such that for allṽ ∈ V Henceforth, we refer to (4.4) as ideal method, since this method only has theoretical purposes and cannot be applied in practice; see Section 6.
4.3. Stability and accuracy of ideal method. Since (4.4) coincides with the ideal method from [Pet17], also the stability result and convergence results are the same.
Lemma 4.6 (Stability of ideal method). Let Assumption 4.1 be satisfied. Then the trial space U and the test space V satisfy the inf-sup condition Proof. For the proof, see Appendix A.
Lemma 4.7 (Accuracy of ideal method). Let u ∈ V be the solution of (2.3) for the righthand side f . If Assumption 4.1 is satisfied, thenũ = (1 − C L )u ∈ U is the unique solution of (4.4), that is, the Petrov-Galerkin approximation of u in the subspace U with respect to the test space V. Moreover, there is a constant c independent of κ, H L and L such that Proof. For the proof, see Appendix A. 4.4. Block-diagonalization. The special choice of the hierarchy of bases (4.2) decouples the equations on every level; see also [Owh17,FP20].
Together with (4.3), this implies the assertion. If k < , one can similarly show that Ψ ⊂ W k , which again yields the assertion.
Corollary 4.9. Let Assumption 4.1 be satisfied, then for ∈ {2, ..., L} it holds that Proof. This is an immediate consequence of the proof of Lemma 4.8.
The following lemma shows that the sub-scale problems (4.5) are well-posed.
Lemma 4.11 (Stability of sub-scale problems for ≥ 2). Let Assumption 4.1 be satisfied. Then, for ≥ 2, the trial space Φ and the test space Ψ satisfy the inf-sup condition Using thatφ ∈ W 1 for ≥ 2, yields where we used Assumption 4.1 and Lemma 4.2. Sinceψ ∈ W 1 , we obtain Combining the previous estimates yields the assertion.
Remark 4.12 (Stability of sub-scale problem for = 1). Since Φ W 1 for = 1, the argument from Lemma 4.11 is not valid in this case. For = 1, one obtains a stability estimate similarly to Lemma 4.6 with inf-sup constant α

Localized numerical correctors
In this section, we decompose the correction operator into element-correctors, which decay exponentially fast away from the corresponding element (Lemma 5.2). The decay makes it possible to localize the element-corrector problems to element patches. Approximation properties are shown (Lemma 5.4).
5.1. Exponential decay of element-correctors. For ∈ {1, ..., L} and all T ∈ T , we define element-correctors Here, we used the restricted sesquilinear form a T defined as Remark 5.1 (Well-posedness of element-correctors). The element-correctors are well-posed due to the same arguments as in Remark 4.3. Similarly as before, we get boundedness of C ,T , C * ,T , 1 − C ,T , and 1 − C * ,T all with the same constants.
Next, we show that the moduli of the element-correctors decay exponentially fast.
Lemma 5.2 (Exponential decay of element-correctors). If Assumption 4.1 is satisfied, then there exists a constant 0 < β < 1 independent of κ, H and L such that for all ∈ {1, ..., L}, T ∈ T , v ∈ V and m ∈ N, the element correctors C ,T v satisfy An analogous decay result holds for C * ,T .
Proof. For the proof, see Appendix A.
Remark 5.3 (Rate of decay). From the proof of Lemma 5.2, it becomes clear that, if Assumption 4.1 is satisfied, the rate of decay β is independent of κ, H and L.

Localized correctors.
The exponential decay motivates the localization of the (global) problems (5.1) to patches, i.e., for ∈ {1, ..., L} and an oversampling parameter m ∈ N, we substitute the (global) ansatz space W by its localized counterparts We then define the localized element-correctors for all v ∈ V as the solution of and set C m, * ,T v := C m ,T v. The element-correctors C ,T and C * ,T are constructed such that their sum (over all T ∈ T ) equals the global correctors C and C * , respectively. Hence, the sum of the localized element-correctors should be good approximation to the actual corrector C and C * , respectively. The (exponential) approximation properties are shown in the next lemma.
where β is the constant from Lemma 5.2. An analogous result holds for C m, * − C * .
Proof. For the proof, see Appendix A.

Sparsification
In this section, we introduce localized variants of the hierarchical multi-resolution trial and test spaces from Section 4. The localization procedure results in a practical method with sparse block-diagonal system matrix; see Figure 6.1 for a comparison of the sparsity patterns for the ideal and localized methods.
6.1. Localized trial and test spaces. Here, we introduce the novel localization strategy that improves stability even for the elliptic case. Instead of just replacing the correctors in (4.2) by their localized counterparts, we additionally utilize the operator P in the definition of the localized trial and test spaces. For = 1, ..., L, let  Remark 6.1 (Effect of operator P ). Since P preserves element-averages, we have (1 − C l ) Πφ = (1 − C )P Π φ for all φ ∈ span{H } (analogously for C * ). Thus, (6.1) is consistent with (4.2) in the ideal (non-localized) case. The trial and test spaces (6.1) define a method that is uniformly stable in H . In contrast, omitting P in (6.1) would result in a method, where the oversampling parameter m had to be increased if H 1 was decreased; see [Owh17, OS19, Mai20].
Remark 6.2 (Optimality of exponential decay of basis functions). In the recent work [HP21], a multi-scale basis is constructed with super-exponential localization properties. While some theoretical aspects are still unclear, numerical experiments show the supremacy of the novel localization approach compared to existing approaches. Remark 6.3 (Exponential smallness of off-diagonal blocks). After localization, the operator is not exactly block-diagonal anymore. However, it can be shown that the off-diagonal blocks are exponentially small; see proof of Theorem 6.5. Thus, for sufficiently large m, this error does not impact stability and convergence results.
Let us propose the practical variant of the multi-level method (4.5), (4.6): First, on every level, solve Petrov-Galerkin problems, which seekφ m ∈ Φ m , such that for allψ m ∈ Ψ m Second, define the solution as the sum of the level contributions (6.4)ũ m :=φ m 1 + ... +φ m L . 6.2. Stability and accuracy of the practical multi-level method. Provided that the localized method is sufficiently close to the ideal method, the stability and convergence results of the ideal method (Lemma 4.6 and 4.7) carry over to the localized method. The closeness is ensured by a coupling of the oversampling parameter m to the stability constant c stab (κ). We assume a polynomial dependence of the stability constant on the wave number κ; for a discussion of this assumption, see Section 2.
Proof. Let ∈ {1, ..., L} be fixed and letφ m ∈ Φ m be arbitrary. Defineφ := (1−C )P φ m andψ m := (1 − C m, * )P ψ , whereψ ∈ Ψ is chosen such that (Lemma 4.11) We obtain The second term can be bounded by It remains to estimate the V-norm ofφ andψ against the V-norm of its localized counterparts and vice-versa. With (3.4) and Remark 4.3, we obtain where we used that c m β m can be bounded independently of m. The corresponding estimates forψ andψ m can be obtained using similar arguments. We conclude . In the last inequality, we employed the oversampling condition (6.5).
Remark 6.6 (Level oversampling parameter). In the proof of Theorem 6.5, it is possible to relax oversampling condition (6.5), choosing level-dependent oversampling parameters Proof. We first estimate using the triangle inequality withφ m ∈ Φ m solving (6.3). The first term can be estimated with Lemma 4.7. The functionsφ m are non-conforming, non-consistent Petrov-Galerkin approximations ofφ ∈ Φ solving (4.5). Using Strang's Lemma (c.f. [KA03]) yields For the first term, we chooseφ m := (1 − C m )P φ . Using Lemma 4.11, we get Letψ ∈ Ψ be arbitrary. Then, for the second term, we obtain after some algebraic manipulations The choiceψ := (1 − C * )P ψ m yields . After dividing by ψ m V(D) , the desired estimate follows. Combining both estimates, we obtain Using Lemma 4.11, Remark 4.12, and Theorem 6.5, we can derive the bounds c 1,m,κ c m κ 2n+2 for = 1 and c ,m,κ c m for ≥ 2. Summing over all levels proves the assertion with the constant c κ,m,L := c m (κ 2n+2 + L).
Remark 6.9 (L 2 -error estimate). In Theorem 6.7 an error estimate with respect to the norm · V(D) is stated. Defining e := u −ũ, an estimate in the weaker L 2 -norm can be proved for the ideal method where we used that e ∈ W L . This directly implies convergence of order 2 + s for the ideal method (4.4) provided that f ∈ H s (D). It is straightforward to extend this result to the localized multi-level method, i.e., Remark 6.10 (Comparison to classical Helmholtz analysis). The analysis of standard finite element discretizations for Helmholtz problems is based on an observation by Schatz [Sch74]. He observed that, if an indefinite sesquilinear form satisfies Gårding's inequality, the stability and quasi-optimality of a finite element discretization can be shown under the two assumptions: (i) approximation properties of the finite element test space and (ii) a sufficiently small mesh size.
In the original paper introducing the LOD for Helmholtz problems, Schatz's argument is employed for the error analysis of the localized method; see [Pet17,Theorem 5.5]. However, due to the sophisticated multi-resolution structure of the discrete problem, Schatz's argument could not be employed in the present work. Nevertheless, using Strang's lemma and the special properties of the multi-resolution ansatz spaces, we were able to prove similar error estimates also for the multi-resolution method with a slightly larger constant; see Theorem 6.7. 6.3. Practical aspects. The localized element-corrector problems (5.1) are still infinitedimensional problems and need to be discretized itself. For the ease of presentation, we restrict ourselves in this paper to the classical case of Q 1 finite elements. We emphasize that a large variety of discretization schemes can be applied, in particular, also hp adaptive methods.
The patch-local corrector problems (5.2) are discretized and solved on a fine mesh of the respective patch with mesh-size h satisfying that κ 2 h is sufficiently small. Under additional geometric requirements, this guarantees stability and quasi-uniformity of the Q 1 finite element discretization; see [Mel95]. Note that all element-corrector problems are independent and thus, can be solved in parallel. The number of fine-scale problems to be solved on level is independent of H and only depends on the oversampling parameter m. In Lemma 4.2, it was shown that the corrector problems are coercive. However, this is not a real practical benefit, since the spaces W are difficult to construct numerically. It is more practicable to derive a saddle point formulation of the problem with constraints enforcing that the solution is in W . The number of constraints (#c) is small and only depends polynomially on m.
The solution procedure proposed in [MP20] computes #c Helmholtz problems on every (reference) patch. Since the patches have a diameter at most of order mH 1 , the effective wave number of the patch problems is at most of order m (Assumption 4.1). For such Helmholtz problems, there exist effective preconditioners; see [GGS15,GZ19,RN20]. With the solution of the Helmholtz patch problems at hand, it is then possible to calculate the Schur complement explicitly. For a more detailed discussion of the numerical solution of the corrector problems, see [EHMP19,MP20].

Fast solvers
In this section, we propose a strategy for solving the level problems (6.3) efficiently. As Remark 6.6 already indicates, for ≥ 2, problems (6.3) are good-natured in the sense that it can be solved up to a given tolerance within a fixed number of GMRES steps (Theorem 7.1). However, for = 1, the κ-dependence in Remark 6.6 makes this impossible. Since this block is relatively small, we suggest a direct solver for the first block.
Recalling the choice of bases (6.2), we can rewrite (6.3) for = 1, ..., L into (decoupled) linear systems of equations 7.1. GMRES method. Since the linear systems of equations (7.1) are sparse but nonhermitian, our preferred iterative solver is the GMRES method; see [Saa03]. For general invertible matrices A ∈ C N ×N , its convergence can be expressed in terms of the field of values F(A) which is defined as with (·, ·) C N denoting the standard inner product on C N . Henceforth, let | · | 2 denote the norm induced by (·, ·) C N . Consider a right-hand side f ∈ C N and let x (k) denote the k-th iterate of the GMRES method applied to Ax = f . Then, the residuals r (k) := Ax (k) − f converge to zero with the rate This result can be found in [LT20].
7.2. Uniform number of GMRES iterations for ≥ 2. The following theorem relies on uniform upper and lower bounds of the fields of values. With these bounds, the convergence rate of the GMRES method can be uniformly bounded away from one.
Theorem 7.1 (Uniform convergence of GMRES applied to (7.1) for ≥ 2). Suppose that ≥ 2 and let Assumption 4.1 be satisfied as well as the oversampling condition Then the linear systems (7.1) can be solved up to a specified relative tolerance within a fixed number of iterations depending only on the tolerance.
Proof. Let ∈ {2, ..., L} be fixed. First, we consider the ideal (non-localized) case, with basis functions b ,j := (1 − C ) Π φ ,j and b * ,j := (1 − C * ) Π φ ,j . We derive upper and lower bounds for the field of values of A := a(b ,k , b * ,j ) j,k=1,...,N . For any ξ ∈ C N , we obtain the upper bound where we used (3.2), (4.3), as well as Assumption 4.1, Lemma 4.2 and Remark 4.3. The lower bound can be derived as follows = 1 H 2 |ξ| 2 2 using (3.1) and Lemma 4.2. Next, we consider the localized case. With (4.3), we can deduce the relation For the first term, we can use the bounds from the ideal case. The second term is exponentially small, since where we used

Numerical experiments
In this section, we present a sequence of numerical results. We begin with a convergence test, which confirms the convergence result from Theorem 6.7 numerically. Next, we demonstrate the improved stability properties (even for the elliptic case) by the novel construction (6.1); see Remark 6.1. Moreover, we demonstrate that the proposed method is also applicable to heterogeneous Helmholtz problems. A high-frequency scenario underlines the effectiveness of the proposed numerical method also for such regimes. Lastly, we consider a scattering scenario with relatively large wave number and evaluate the condition numbers of the respective levels. This confirms the statement of Theorem 7.1 numerically. 8.1. Convergence test. We consider D := (0, 1) 2 with Γ R := ∂D and the smooth righthand side (8.1) f (x 1 , x 2 ) = sin(πx 1 ) cos(πx 2 ).
The convergence test is performed for κ = 2 j , j = 0, ..., 3. For all κ, we fix a coarse (Cartesian) mesh with H 1 = 2 −j−1 and perform the numerical computations for mesh hierarchies with H L = 2 −j−1 , ..., 2 −7 . The corrector problems (5.2) are discretized using a fine mesh with mesh size h = 2 −9 . Errors are computed to the reference solution on the same fine mesh. In the convergence plots of Figure 8.1, one can clearly see the LOD-typical behavior [MP20] that the error first decreases with order 2 (f ∈ H 1 (D)), until the localization error dominates, then, the overall error stagnates. This is the expected outcome and is in line with Theorem 6.7. For increasing κ, one observes that for small oversampling parameters m, the error does not decrease, but stays of order O(1). In this case, the oversampling parameter does not fulfill condition (6.5) and convergence cannot be expected.
8.2. Improved stability by novel basis construction. A already mentioned in Remark 6.1, the novel basis construction (6.1) (henceforth referred to as stabilized gamblets) has superior stability properties compared to previously known constructions (henceforth referred to as normal gamblets); see [Owh17,OS19,Mai20]. Since this is true even for the elliptic case, we use the Poisson problem for demonstration purposes For the numerical experiment, we set D := (0, 1) 2 and use the source (8.1). For the sake of simplicity, we consider a single-level method (L = 1) for demonstrating the improved stability properties. For the multi-level method, the observed phenomena are qualitatively the same. We compute numerical approximations for H 1 = H L = 2 −3 , ..., 2 −7 and calculate correctors on a fine mesh with mesh size h = 2 −9 . Errors are calculated to the reference solution on the same fine mesh. In Figure 8.2a, one clearly observes the improved stability properties of stabilized gamblets (6.1) compared to normal gamblets. The error for normal gamblets increases as H 1 decreases. Thus, in order to preserve stability, the oversampling parameter m has to be increased as the mesh is refined. However, this is not true for the stabilized gamblets, where we have stability for a fixed oversampling parameter. The error even decreases further with a smaller order. This can be explained by the absence of the block-diagonalization error (c.f. Remark 6.3) in the single-level case.
For comparison, we also performed this numerical experiment for the Helmholtz problem with κ = 2 3 and Robin-type boundary conditions. Here, the exact same phenomena can be observed; see Figure 8.2b.
On each inclusion in D, the coefficient A is chosen to be constant with value uniformly distributed in [1,16]. Everywhere else, A is set to 1. For the numerical experiment, we used the realization of A shown in Figure 8.3. Since A| ∂D = 1, this choice is compatible with the Robin boundary conditions from (2.1). The right-hand side is chosen as with r = 0.125 and x 0 = (0.5, 0.5) T . In this example, we consider κ = 2 3 . Note that the effective wave number in the inclusions is actually smaller than κ. As fixed coarse mesh, we use a Cartesian mesh with H 1 = 2 −3 . The numerical computations are performed for hierarchy of meshes with H L = 2 −3 , ..., 2 −7 . The correctors are computed on a fine mesh with h = 2 −9 . Errors are computed to the reference solution on the same fine mesh. As norm for the errors, we choose In Figure 8.3, one observes a convergence behavior similar as in Figure 8.1. In contrast to the convergence test for κ = 2 3 in Figure 8.1, the oversampling parameter m = 2 already yields a good approximation in this example. This can be explained by the smaller effective wave number, caused by the coefficient A.
8.4. High-frequency. In this numerical experiment, we again choose D := (0, 1) 2 with pure Robin boundary conditions, i.e., Γ R := ∂D. We consider the wave number κ = 2 6 and use the right-hand side (8.3) with the same parameters as in the previous experiment. As fixed coarse mesh, we use a Cartesian mesh with H 1 = 2 −6 . The numerical computations are performed for hierarchy of meshes with H L = 2 −6 , ..., 2 −8 . The correctors are computed on a fine mesh with h = 2 −10 . Errors are computed to the reference solution on the same fine mesh.
In Figure 8.4 one can observe the convergence behavior as expected from Theorem 6.7, i.e., if the oversampling parameter m is chosen sufficiently large (i.e., m log(κ)), one has κ-independent convergence of the proposed multi-resolution method. This numerical example proves the effectiveness of the numerical method also for the high-frequency regime. We choose κ = 2 5 . The hierarchy used for the computations is specified in Figure 8.5. The correctors are computed on a fine mesh with h = 2 −10 . The oversampling parameter is chosen as m = 2. We use the GMRES method with restart after 50 iterations. The GMRES iteration terminates if a relative residual of 10 −6 is reached.  The table in Figure 8.5 clearly shows that uniform boundedness of the number of GM-RES iterations for ≥ 1. As expected, the GMRES iterations performs very poorly for = 1, which suggests to use a direct solver for the relatively small linear system of equations corresponding to = 1. The respective condition numbers, which are a good indicator for how fast iterative solvers converge, underline this observation.
If f ∈ H 1 (D), then (3.1) yields an additional order in H L . Dividing by e V(D) , the result follows for s ∈ {0, 1}. By arguments from interpolation theory (see e.g. [BS08]), the assertion can be concluded for any s ∈ [0, 1] with a non-explicit constant c.
Here, we applied (3.1) and (3.2) multiple times. The element-correctors satisfy the estimate ∇C ,T v L 2 (D) v V(T ) , since Using this and Lemma 4.2, Lemma A.1 and the finite overlap of the patches, we get after summing over all elements