A Super-Localized Generalized Finite Element Method

This paper presents a novel multi-scale method for elliptic partial differential equations with arbitrarily rough coefficients. In the spirit of numerical homogenization, the method constructs problem-adapted ansatz spaces with uniform algebraic approximation rates. Localized basis functions with the same super-exponential localization properties as the recently proposed Super-Localized Orthogonal Decomposition enable an efficient implementation. The method's basis stability is enforced using a partition of unity approach. A natural extension to higher order is presented, resulting in higher approximation rates and enhanced localization properties. We perform a rigorous a priori and a posteriori error analysis and confirm our theoretical findings in a series of numerical experiments. In particular, we demonstrate the method's applicability for challenging high-contrast channeled coefficients.


Introduction
We consider the numerical solution of a second-order linear elliptic partial differential equation with a strongly heterogeneous coefficient.The coefficient may be non-periodic, with oscillations appearing on several non-separated scales.For such coefficients, classical finite element methods based on problem-independent polynomial ansatz spaces typically yield unsatisfactory approximations, cf.[BO00].It is possible to overcome this issue by incorporating problem-specific information into the method's ansatz space, which is commonly known under the term numerical homogenization and has been an active research field throughout the past decades.For an overview on numerical homogenization, we refer to the recent textbooks [OS19,MP20] and the review article [AHP21].
Under minimal assumptions on the coefficient, numerical homogenization is able to achieve optimal orders of approximation without any pre-asymptotic effects.However, this is not possible without a computational overhead.Compared to classical finite element methods, it is either necessary to consider basis functions with an enlarged support, or to increase the number of basis functions per mesh entity.As prominent examples, we mention the Multiscale Spectral Generalized Finite Element Method (MS-GFEM) [BL11, EGH13, MSD22], the Adaptive Local Basis (AL-Basis) [GGS12,Wey17], the Localized Orthogonal Decomposition method (LOD) [HP13, MP14, KPY18, BGS21], Rough Polyharmonic Splines (RPS) [OZB14], and gamblets [Owh17].
The above-mentioned approaches can be distinguished into two classes.Methods like MS-GFEM and the AL-Basis first solve local spectral problems in the space of (locally) operator-harmonic functions.The respective ansatz spaces are then constructed by gluing together local eigenfunctions by means of a partition of unity [MB96,BM97].For such methods, the support of the basis functions is fixed by the choice of partition of unity.For convergence of optimal order, the number of local eigenfunctions taken into account needs to be increased logarithmically with the desired accuracy.In order to make these approaches computationally more efficient, one can use random sampling strategies as proposed, e.g., in [BS18].
The second class of methods includes the LOD, RPS, and gamblets.The idea is to construct problem-adapted ansatz spaces by applying the solution operator to specific classical finite element spaces with respect to some coarse mesh T H , which typically does not resolve the coefficient's oscillations.Due to their connection to isogeometric analysis in the case of constant coefficients, such methods are sometimes referred to as spline-type approaches.While optimal approximation orders of these methods are achieved by design, the true challenge is to construct a local basis of the problem-dependent ansatz space.An almost optimal solution is provided by the LOD, which constructs a fixed number of basis functions per mesh entity that decay exponentially fast with respect to the coarse mesh.This rapid decay enables a localization of the basis functions to -th order element patches with diameters of order H.For convergence of optimal order, the oversampling parameter needs to be increased logarithmically with the desired accuracy.
Recently, the Super-Localized Orthogonal Decomposition (SLOD) has been proposed in [HP22] (see also [FHP21,BFP22,HM22]).The key contribution of the SLOD is a novel localization strategy that enables a significantly improved localization compared to the LOD.The SLOD constructs rapidly decaying basis functions, yielding super-exponentially decaying localization errors.These improved localization leads to smaller local patch problems for the basis computation and a sparser coarse system matrix.Numerical experiments indicate that the SLOD outperforms the LOD, achieving similar magnitude errors for significantly smaller oversampling parameters.Until now, for the best practical realization of the SLOD in [HP22], the stability of the SLOD basis functions cannot be guaranteed a priori.For high-contrast channeled coefficients or convection-dominated regimes, these basis stability issues may deteriorate the method's approximation quality (see the numerical experiments in Section 7).
This paper proposes a novel multi-scale method that, on the one hand, preserves the unique localization properties of the SLOD and, on the other hand, resolves the aforementioned basis stability issues.This is achieved by combining the SLOD with a partition of unity approach.More precisely, locally on nodal patches, we apply the respective local solution operator to classical finite element source terms.Multiplying these spaces with the corresponding hat-functions yields local ansatz spaces with a low effective dimension.Consequently, low-dimensional optimally approximating spaces are constructed by solving local spectral problems.Compared to MS-GFEM methods, the proposed method has the major advantage that the local spectral problems are posed in a space spanned by a small number of deterministic snapshots.Hence, possible random sampling strategies can be avoided.Furthermore, due to their low dimension, the local spectral problems are easy to solve.The global problem-adapted ansatz space is then obtained by gluing together the local optimal approximation spaces using a partition of unity.We highlight that, from an application point of view, the proposed multi-scale method is conceptually simple and straightforward to implement.Further, by adapting the polynomial degree of the finite element source terms, one can easily construct higher-order versions of the method.Similarly as for the higher-order LOD [Mai21, DHM22], one obtains higher-order convergence rates using the regularity of the source only.
We prove that the proposed method possesses the advantageous localization and convergence properties of the SLOD which can be quantified a posteriori.Building on the well-understood theoretical foundation of the LOD, we additionally perform a pessimistic a priori error analysis, proving that the proposed method at least recovers the convergence and localization properties of the LOD.For the method's higher-order versions, we observe that solely increasing the polynomial degree significantly improves the localization properties.Numerical experiments even suggest that an (almost) local basis exists for sufficiently large polynomial degrees.Another noteworthy contribution of this work is the implementation of the proposed method as well as the SLOD in the Python-library gridlod [HK17].This serves the principles of open access and reproducibility while enabling computations on large parallel clusters.
The outline of this paper is as follows.In Section 2, we introduce the prototypical elliptic model problem.Section 3 recalls preliminary results, which are then used in Section 4 for the definition of the novel multi-scale method.An a posteriori error analysis is presented in Section 5 followed by a pessimistic a priori error analysis in Section 6.Finally, Section 7 presents a series of numerical experiments which confirm our theoretical findings.

Model problem
We consider the prototypical second-order elliptic PDE −divA∇u = f in weak form with homogeneous Dirichlet boundary conditions on a polygonal/polyhedral Lipschitz domain Ω ⊂ R d , d ∈ {2, 3}.Furthermore, without loss of generality, we assume that Ω is scaled such that its diameter is of order one.The coefficient function A ∈ L ∞ (Ω, R d×d ) may be matrix-valued and is assumed to be symmetric and positive definite almost everywhere.More specifically, we assume that there exist constants 0 < α ≤ β < ∞ such that (2.1) with |•| denoting the Euclidean norm.The weak formulation of the elliptic model problem uses the Sobolev space V := H 1 0 (Ω) and the bilinear form a : V × V → R, given by a(u, v) := ˆΩ(A∇u) • ∇v dx.
The symmetry and the condition (2.1) ensure that the above bilinear form is an inner product on V. Its induced norm is the energy norm • a,Ω := a(•, •), which is equivalent to the canonical Sobolev norm on V.The Lax-Milgram theorem ensures that, for all source terms f ∈ L 2 (Ω), there exists a unique weak solution u ∈ V to the boundary value problem, satisfying Note that the moderate restriction to source terms in L 2 (Ω) (rather than the dual space V = H −1 (Ω)) will be essential for the uniform convergence of the numerical homogenization method.We note that the possibly rough coefficient generally prevents H 2 (Ω)regularity of the solution, which would be required by classical finite elements.For a generalization to source terms with less regularity, we refer to [AHP21].Let us mention that the proposed method is not restricted to the class of elliptic PDEs with Dirichlet boundary conditions.Considering the extensions of the SLOD [FHP21,BFP22], also an extension of the proposed method to (non-symmetric) coercive operators and even Helmholtz-type problems seems possible.In particular, more general boundary conditions of Neumann and Robin type can be taken into account.Henceforth, we refer to A −1 : L 2 (Ω) → V as the solution operator that maps f ∈ L 2 (Ω) to the unique solution u ∈ V of (2.2).Moreover, for a subdomain ω ⊂ Ω, we denote by a ω (•, •) and • a,ω the restriction of the bilinear form a(•, •) to ω and the restricted energy norm, respectively.The restricted solution operator subject to homogeneous Dirichlet boundary conditions on ∂ω is denoted by A −1 ω : L 2 (ω) → H 1 0 (ω).

Preliminaries
Let T H denote a quasi-uniform coarse mesh of Ω consisting of closed, simplicial, or quadrilateral shape-regular elements.The subscript H denotes the maximal element diameter, i.e., H := max T ∈T H diam(T ) and by N H , we denote the set of all (interior and boundary) vertices of T H .For the ease of presentation, we henceforth only consider quadrilateral meshes; the extension to triangular meshes is straightforward.
The proposed multi-scale method utilizes the concept of patches.For any ∈ N 0 , we define the -th order (element) patch of a union of elements ω ⊂ Ω recursively by 3.1.Discontinuous finite element spaces.For a fixed (but arbitrary) polynomial degree p, we denote with the non-conforming space (with respect to V) consisting of element-wise defined polynomials.We define the restriction of P(T H ) to a subdomain ω ⊂ Ω by P(ω) := {v ∈ P(T H ) : supp(v) ⊂ ω}.
Let Π H : L 2 (Ω) → P(T H ) denote the L 2 -orthogonal projection, which for each v ∈ L 2 (Ω), is given by the element-wise equation for all w ∈ P(T ), T ∈ T H .
The projection satisfies the following stability and approximation estimates with constant C a > 0 depending only on the regularity of the mesh T H and the polynomial degree p, see, e.g., [HSS02].
In addition, we define the broken Sobolev space H k (T H ), k ∈ N, by 3.2.Conforming companion spaces.We next define local conforming companions of the functions Θ T,j , so-called bubble functions, which we denote by b T,j .For each element T ∈ T H , these functions fulfill, for 1 ≤ j ≤ J, We do not require an explicit characterization.However, it is important that such functions actually exist.This is guaranteed by [Mai21,Cor. 3.6] stating that, for all Θ T,j , there exist a corresponding bubble function b T,j such that with constant C b > 0 depending solely on the mesh regularity of T H and the polynomial degree p.By means of the bubble functions, we can define the operator B H mapping possibly nonconforming functions to T H -piecewise bubble functions with the same L 2 -projection.For any function in P(T H ), we uniquely define B H by setting B H Θ T,j := b T,j for all T ∈ T H , j = 1, . . ., J. We can extend the operator to L 2 (Ω) by setting Clearly, the kernels of the operators Π H and B H coincide and one can prove the local stability estimate with another constant C bo > 0 depending solely on the mesh regularity of T H and the polynomial degree p.By the definition of B H , we obtain, for all v ∈ L 2 (Ω), q ∈ P(T H ), i.e., B H is the L 2 -projection onto the space of bubble functions.
3.3.Partition of unity.The proposed multi-scale method is based on the framework of partition of unity methods, cf.[MB96,BM97].Although, there is great flexibility in the choice of such a partition, for simplicity, we restrict ourselves to the hat-functions {Λ z : z ∈ N H } corresponding to all (interior and boundary) nodes of T H . Recall that the hat-function Λ z associated with node z ∈ N H is a continuous T H -piecewise bilinear function uniquely defined by setting its nodal values for all y ∈ N H to Λ z (y) = δ yz with δ denoting the Kronecker symbol.By definition, the hat-functions have an L ∞ -norm of one and, due to the shape-regularity of T H , their gradients satisfy with constant C Λ > 0 depending solely on the mesh regularity of T H . Denoting ω z := supp(Λ z ), the shape-regularity of T H also implies that the supports {ω z : z ∈ N H } have a finite overlap, i.e., the maximal number of overlapping supports is uniformly bounded.Subsequently, we abbreviate the node patches around z ∈ N H and the element patches around T ∈ T H by ω z := N (ω z ) and ω T := N (T ), respectively.

Multi-scale method
This section introduces the proposed multi-scale method.The local ansatz spaces of the method are constructed by applying the local solution operator on an oversampling domain to piecewise polynomial source terms and by subsequent restriction to a subdomain.Due to the oversampling, the resulting local spaces have a low effective dimension, and thus, low-dimensional optimally approximating spaces are utilized.The ansatz space of the method is obtained by gluing together the low-dimensional local approximation spaces.Note that we consider a fixed polynomial degree p and do not track the dependence of constants on p. Explicitly tracking this dependence would make the analysis less clear and also add no value, as it relies on estimates that are pessimistic in p.

4.1.
Local approximation spaces.For any z ∈ N H , we aim to approximate the restriction of the solution space V| ωz using local approximation spaces.This is accomplished by choosing the local approximation space V H,z | ωz with being defined on the oversampling domain ω z .
Due to the oversampling, the restricted space contains many redundant functions.This holds, in particular, after the multiplication with the hat-function Λ z when gluing the local approximation spaces together.Hence, we investigate the optimal approximation of Λ z V H,z by n-dimensional subspaces Q(n) ⊂ H 1 0 (ω z ).Given the subspace Q(n), its worst-case best approximation error is defined as Typically, the minimal worst-case best approximation error is referred to as Kolmogorov n-width, cf.[Pin85], and is defined as Indeed, there exists a corresponding optimal local approximation space of dimension n, which we explicitly compute.For this, we solve the low-dimensional eigenvalue problem, which seeks eigenpairs We denote the eigenfunctions by {v i : i = 1, . . ., N := dim(V H,z )} assuming an ordering such that the corresponding eigenvalues satisfy H,z := span{v i : i = 1, . . ., n}, the optimal local approximation space of dimension n is given by Λ z V ,n H,z .4.2.Global approximation space.A global approximation space is obtained by gluing together the above local approximation spaces using the partition of unity, i.e., We measure the overall error when approximating Λ z V H,z by spaces of dimension n by (4.4) Finally, the proposed method seeks for all v ∈ V ,n H .Note that the mesh size H determines the accuracy of the approximation.The oversampling parameter specifies the size of the local patch problems and determines the method's localization error, whereas the number of local functions is given by n and needs to be chosen sufficiently large.For a precise choice of the parameters H, , and n, we refer to Remarks 5.6 and 6.3.The following two sections are devoted to the theoretical analysis of the proposed multi-scale method.In Section 5, we present an a posteriori error analysis, while in Section 6, we present a priori error bounds.

A posteriori error analysis
Subsequently, we derive an a posteriori error analysis of the proposed method by establishing a connection to the SLOD introduced in [HP22].The SLOD is conceptually related to the proposed method, as it also constructs its basis functions by applying the local solution operator to T H -piecewise polynomial source terms.5.1.Higher-order SLOD.For the a posteriori error analysis, we first briefly introduce a higher-order variant of the SLOD.Note that this variant only serves theoretical purposes and is not investigated numerically.Let us fix an arbitrary element T ∈ T H and oversampling parameter ∈ N. Henceforth, we drop all fixed indices and denote the -th order patch around T just by ω := ω T .Furthermore, we make the meaningful assumption that no patch ω coincides with the entire domain Ω.
An example of a continuous right-inverse of tr is the A-harmonic extension, henceforth denoted by tr −1 .Given w ∈ X, it satisfies tr tr −1 w = w and (5.2) a ω (tr −1 w, v) = 0, for all v ∈ H 1 0 (ω).Using definitions (5.1) and (5.2), and that (v This result yields, together with the definition of ϕ j and the local support of g j , the following key observation . Hence, we can rephrase the smallness of the localization error as the (almost) L 2 (ω)orthogonality of g j to the space (5.4) of A-harmonic functions on ω (which satisfy the homogeneous Dirichlet boundary condition on ∂Ω ∩ ∂ω).Since the restricted L 2 -projection Π H | Y has a finite rank of dimension less or equal to K := J • #(T H ∩ ω), there exists a singular value decomposition (SVD) such that where σ 1 ≥ • • • ≥ σ K ≥ 0 denote the singular values, {q 1 , . . ., q K } the L 2 (ω)-orthonormal left singular vectors, and {w 1 , . . ., w K } the H 1 (ω)-orthonormal right singular vectors.We choose the source terms g j as the left singular vectors corresponding to the J smallest singular values, i.e., (5.6)

This yields sup
which follows directly from the properties of the SVD.We define the quantity σ measuring the (quasi-)orthogonality between the g j and Y as (5.7) σ(H, ) := max The quantity σ is crucial for the error analysis of the SLOD as it determines the localization error.Note that the dependence of σ on the (fixed) polynomial degree p is not made explicit in (5.7).The following remark deals with the decay of σ with respect to the oversampling parameter and, for the sake of curiosity, also the polynomial degree p.
Remark 5.1 (Decay of σ).In [HP22], it has been numerically observed and conjectured that, for p = 0, the quantity σ decays super-exponentially as is increased, cf. Figure 5.1 (left).A similar decay in can also be observed for p > 0. In accordance with [HP22, FHP21], we state the following conjecture: there exists C σ > 0 depending algebraically on H, and C > 0 independent of H, such that Conversely, for fixed , a rapid decay of σ in p can be observed, cf. Figure 5.1 (right).The low level of magnitude of the last singular values may suggest that the respective source terms correspond to fully local basis functions.However, due to the low levels or singular values, this is difficult to verify numerically.
Choosing for each patch the basis functions corresponding to the J source terms (5.6), we obtain the SLOD ansatz space (5.9) V , SLOD H := span{ψ T,j : T ∈ T H , j = 1, . . ., J}.The Galerkin SLOD solution u , SLOD Note that a reasonable SLOD approximation requires a stable choice of the basis functions in (5.9).However, for large oversampling parameters and patches intersecting the boundary, the choice (5.6) may be insufficient, and a special treatment is required.In [HP22, App.B], such an algorithm is proposed, curing possible stability and uniqueness issues in practice.Since we still cannot guarantee stability in an a priori manner, we assume that the source terms corresponding to the basis function of (5.9) form a Riesz basis of P(T H ), which is formulated in the following.

5.2.
A posteriori error bound using SLOD.We provide an a posteriori error analysis of the proposed method based on SLOD techniques.Conceptually, it is similar to the one for the SLOD, cf.[HP22, Thm.6.1], but it additionally includes the local optimal approximation error d n defined in (4.4).
Theorem 5.3 (A posteriori error bound).Let Assumption 5.2 be satisfied and let u and u ,n H denote the solutions to (2.2) and (4.5), respectively.Then, there exists a constant C > 0 independent of H, , and n, such that, for any with s := min{k, p + 1} and the notation Proof.The application of Céa's Lemma yields, for arbitrary v ∈ V ,n H , (5.12) u − u ,n H a,Ω ≤ u − v a,Ω .Let the solution to the (higher-order) collocation variant of the SLOD, cf.[HP22, Rem.5.1], with oversampling parameter m := /2 be denoted by u m, SLOD H .Given a source term f ∈ L 2 (Ω), its solution is obtained by the following linear combination of localized basis functions ψ m T,j : with coefficients c T,j that are uniquely defined by (5.13) Adding and subtracting u m, SLOD H in (5.12) and employing the triangle inequality yields The first term is the error of the collocation variant of the SLOD.It can be bounded using a higher-order version of [HP22, Thm.6.1] stating the existence of a constant C slod > 0 independent of H and m such that . For the second term in (5.14), we choose v ∈ V ,n H as sum of functions v n z ∈ V ,n H,z to be specified later, i.e., v Using the partition of unity property of the hat-functions z∈N H Λ z ≡ 1, we obtain , we use , and performing the above-mentioned local replacement, we obtain a,ωz , (5.15) where we add and subtract v z and employ the triangle inequality.Using the product rule and the bound (3.7), we get for the first term Noting that by e 1 z ∈ H 1 0 (ω z ), we obtain for the first term in (5.16) (5.17) using Friedrichs' inequality on ω z with diam(ω z ) ≤ C p H, C p > 0. For the second term in (5.16), we infer the trivial estimate , which implies that, in order to bound (5.16), it suffices to estimate ∇e 1 z L 2 (ω z ) .Using the continuity estimate with a constant C tr > 0 independent of H, from the proof of [HP22, Thm.6.1], one can show that (5.18) (g m T,j , tr −1 . By (5.3) and (5.18), as well as the discrete Cauchy-Schwarz inequality, the finite overlap of the patches ω m T and Friedrichs' inequality, we get with constant C > 0 reflecting the overlap of the patches ω m T .Using (5.17), this yields an estimate for (5.16) and consequently bounds the first term in (5.15).
For the second expression in (5.15), using the definition of the Kolmogorov n-width (4.2), Friedrichs' inequality on the patch ω z , the discrete Cauchy-Schwarz inequality, and that the g m T,j are L 2 -normalized, we obtain with constant C > 0 appearing in the bound C 2 m d J of the number of terms in the above sum.Consequently, we conclude the estimate for the second term in (5.15).The assertion can be finalized by combining all estimates utilizing with constant C > 0 reflecting the overlap of the patches ω m z .Here, we used Assumption 5.2 and (3.1) and (5.13).Finally, for the sake of readability, we substitute m = 2 by , which may introduce additional constants that change the decay rate of σ by some constant factor.

5.3.
Local approximation error bound using SLOD.The error estimate from Theorem 5.3 incorporates the spectral approximation error d n defined in (4.4).Subsequently, we derive a bound for d n utilizing the SLOD for constructing bases of the local spaces V H,z defined in (4.1).For this purpose, we fix a node z ∈ N H and treat ω z as the whole domain.For the oversampling parameters m = 1, . . ., − 1, we denote the basis functions with a tilde to emphasize that, in general, they do not coincide with their global counterparts ψ m T,j , i.e., (5.19) { φm T,j := A −1 ω z gm T,j : T ⊂ ω z , j = 1, . . ., J} with source terms gm T,j ∈ P(ω m T ), where ωm T := ω m T ∩ ω z .We denote the corresponding localized basis by (5.20) { ψm T,j := A −1 ωm T gm T,j : T ⊂ ω z , j = 1, . . ., J}.
Similarly, as in the full domain setting, we need to measure, for all T ⊂ ω z , the (quasi-) orthogonality of the source terms {g m T,j : j = 1, . . ., J} on the corresponding space of A-harmonic functions The quantity σ is strongly related to its counterpart σ from (5.7) and, in numerical experiments, exhibits the same qualitative behavior as described in Remark 5.1.A local variant of Assumption 5.2 is required to ensure the stability of the local basis.
Theorem 5.5 (Bound on d n ).Let Assumption 5.4 be satisfied.Then, there exists a constant C > 0 independent of H, , and m such that, for m = 1, . . ., − 1 Proof.Let us consider a fixed node z ∈ N H and oversampling parameter .As approximation space Q(n) of dimension n ≈ m d , we choose . ., J} with basis functions defined in (5.20).For the approximation of v z ∈ V H,z , we choose the element where the c T,j are the coefficients of the expansion of v z in terms of the basis functions φm T,j defined in (5.19).Note that, by Assumption 5.4, the coefficients c T,j are uniquely determined.Thus, we can estimate the spectral approximation error (4.2) using v z and w z as we obtain for the numerator, using the product rule, the triangle inequality, and (3.7): We apply Friedrichs' inequality on the patch ω z using that diam(ω z ) ≤ C p H, with a constant C p > 0. Hence, we can bound the first term against the second term, i.e., We adapt estimate (5.18) to the local setting introducing a constant Ctr > 0. Using the finite overlap of the patches ωm T , the discrete Cauchy-Schwarz inequality and Friedrichs' inequality on ω z , we obtain where C > 0 reflects the overlap of the patches ωm T .Adding the remaining coefficients c T,j from the expansion of v z in terms of the φm T,j and using Assumption 5.4, we get C−1 r (H, , m) Here, we also employed that, by (3.5) and (3.6), we have, for all q ∈ P(ω z ), Combining the estimate yields the assertion.
Remark 5.6 (Choice of parameters).This remark specifies how to choose the oversampling parameter and the number of local functions n in order to preserve the optimal order of convergence of s + 1 in Theorem 5.3.For , the super-exponential decay (5.8) implies that it needs to be chosen of order |log H| (d−1)/d .Using Theorem 5.5 and that σ has similar decay properties as σ, we obtain that n needs to be chosen of order |log H| d−1 .Note that these choices require the validity of (5.8) and Assumptions 5.2 and 5.4.For a (pessimistic) estimate which is valid without additional assumptions, we refer to Remark 6.3 below.

(Pessimistic) a priori error analysis
This section presents an a priori error analysis of the proposed method, which is based on the LOD framework, cf.[MP14,HP13,AHP21].Note that the exponential localization properties of the LOD cannot match the practically observed super-exponential localization properties of the SLOD, cf.Remark 5.1.Nevertheless, the LOD construction has the decisive advantage that the basis stability is guaranteed by construction and that the exponential localization can be rigorously proved.This enables an a priori analysis without assumptions on the stability of the SLOD basis and without conjectures on the decay of singular values, cf.Assumptions 5.2 and 5.4 and Remark 5.1.6.1.Higher-order LOD.We briefly introduce a higher-order version of the LOD similar to the constructions in [Mai21,DHM22].The LOD constructs its problem-adapted basis functions by adding fine-scale information to coarse-scale finite element functions.We define the space of fine-scale functions as (6.1) The step of adding fine-scale information is called correction and utilizes the so-called correction operator C : V → W defined as the a-orthogonal projection onto W, i.e., a(Cv, w) = a(v, w), for all w ∈ W.
We split up the correction operator into a sum of element correction operators, i.e., C = T ∈T H C T with element correction operators C T : V → W defined by a(C T v, w) = a T (v, w), for all w ∈ W.
For any v ∈ V, the correction C T v decays exponentially fast away from its associated element T ∈ T H , cf. [DHM22, Lemma 5.1], which motivates a localization.For this purpose, we substitute the global space W by localized counterparts W T := W(ω T ), where, for a subdomain ω ⊂ Ω, we use the definition (6.2) W(ω) := {w ∈ W : supp(w) ⊂ ω}.
The localized element correction operator C T : V → W T is then defined by Similarly as for the correction operator C which can be decomposed into a sum of element correction operators, we define the localized correction operator by The ansatz space of the LOD is then obtained by adding (localized) corrections to the bubble functions {b T,j : T ∈ T H , j = 1, . . ., J} defined in (3.3), i.e.,

6.2.
A priori error bound using LOD.Using LOD techniques, we can prove the following a priori error estimate for the proposed multi-scale method.Numerical experiments show that this estimate is tentatively pessimistic.Nevertheless, compared to Theorem 5.3, it has the crucial advantage that it does not rely on additional assumptions or conjectures.
Theorem 6.1 (A priori error bound).Let u and u ,n H denote the solutions of (2.2) and (4.5), respectively.There exist constants C, C d > 0 independent of H, , and n such that, for any . with s := min{k, p + 1}.
Proof.We apply Céa's Lemma, which yields, for arbitrary v ∈ V ,n H , (6.4) For the oversampling parameter m := /2 , we define the approximation which is not the Galerkin LOD solution (6.3) but has the same approximation properties, cf.proof of [Mai21,Thm. 4.4].Adding and subtracting u m, LOD H in (6.4) and using the triangle inequality yields Using the above-mentioned approximation properties of u m, LOD

H
, we obtain the following estimate for the first term in (6.5) with constants C lod , C d > 0 independent of H, , and m.
For the second term in (6.5), we choose v ∈ V ,n H as sum of functions v n z ∈ V ,n H,z to be specified later, i.e., v Using the partition of unity property of the hat-functions z∈N H Λ z ≡ 1, we obtain for the second term in (6.5), with C ol defined in (3.8).For any z ∈ N H , we can, locally on ω z , substitute u m, LOD . Hence, we define v n z ∈ V ,n H,z as the (not necessarily unique) elements minimizing the expression Λ z v z − Λ z v n z a,ω z , where . We denote W z := W(ω z ) and define the above used correction operator C z : Note that it holds v z ∈ V H,z which is a non-trivial observation, cf.[AHP21, Rem 3.
, we obtain after performing the above-mentioned local substitution a,ωz .(6.7) Following (5.16), in order to bound the first term in (6.7), it suffices to bound e 1 z L 2 (ωz) and ∇e 1 z L 2 (ωz) .It holds z has vanishing element averages.Thus, by Poincare's inequality it is sufficient to estimate ∇e 2 z L 2 (ωz) in order to obtain a bound for the first term in (6.7).Given a function supported in ω m z (e.g., B H u| ω m z ) the correction operator C m coincides with the localization of C z to m-th order patches.Hence, we can apply the localization error estimate from the proof of [Mai21,Thm. 4.4], here in the oversampling parameter m and (3.5) to obtain . with constants C lo , C d > 0 independent of H, , and m.
For the second term in (6.7), we obtain by the definition of the Kolmogorov n-width, the stability of (1 . Combining the previous estimates, using Friedrichs' inequality on Ω (recall that Ω is scaled to unit size), we get Finally, we substitute m = 2 by , which introduces additional constants and changes the exponential decay rate C d by a factor of two.6.3.Local approximation error bound using LOD.Subsequently, we derive an a priori bound for d n , which is fully explicit in H and .This is the analog to Theorem 5.5, which does not rely on additional assumptions or conjectures.Numerical experiments show that this estimate is tentatively pessimistic.T : z , j = 1, . . ., J} and as approximation w z ∈ Q(n) of an element v z ∈ V H,z , we use Using the approximation space Q(n) and the above defined choice of w z , we can bound the Kolmogorov n-width as follows Abbreviating e z := ( Cm − C z )B H v z , we can estimate the numerator using (3.7) and It holds that e z ∈ W z , which implies that e 1 z has vanishing element averages.Thus, by Poincare's inequality e z L 2 (ωz) ≤ 2π −1 H ∇e z L 2 (ωz) .Applying the localization error estimate from the proof of [Mai21,Tehorem 4.4] to show that Cm approximates C z exponentially and using (3.5), and Friedrichs' inequality on the patch ω z with diam(ω z ) ≤ C p H, C p > 0, we obtain with constants C lo , C d > 0 independent of H, , and m.The assertion follows immediately.Remark 6.3 (Choice of parameters).This remark specifies how to choose the oversampling parameter and the number of local functions n in order to guarantee the optimal order of convergence of s+1 in Theorem 6.1.By Theorem 6.1, needs to be chosen of order |log H|.Using Theorem 6.2, we obtain that n needs to be chosen of order |log H| d .According to the experiments, these choices are pessimistic, cf.Remark 5.6.

Implementation and numerical experiments
In this section, we numerically investigate the proposed multi-scale method regarding the localization error, optimal convergence properties, high-contrast channeled coefficients, and higher-order polynomials using suitable benchmark problems.As a comparison, we use the SLOD from [HP22], which we consider as state-of-the-art.We refer to [HP22, Sec.8] for a comparison of the SLOD to other multi-scale methods such as the LOD.For underlining the origin of the proposed method and its super-localization properties, cf.Theorems 5.3 and 5.5, we subsequently refer to it as Super-Localized Generalized Finite Element Method (SL-GFEM).7.1.Implementation.For the practical implementation of the SL-GFEM, we need to perform a fine-scale discretization, i.e., we substitute the infinite-dimensional function space V by the finite element space V h := V ∩ P 1 (T h ).Here, T h denotes a fine-scale mesh obtained by uniform refinement of T H , where the number of refinements should be chosen such that the resulting mesh resolves all oscillations of A and f .For solving the patch problems (4.1) and (4.3), one considers patch-local subspaces of V h .
The SL-GFEM is straightforward to implement, as only very few technical details need to be addressed.The local spaces V H,z,h (discrete counterparts of V H,z ) can be computed in parallel.Their computation only requires the local stiffness and mass matrices on the respective patches.
In contrast to, e.g., MS-GFEM methods [BL11, EGH13, MSD22], the SL-GFEM solves local eigenvalue problems which are posed in the space spanned by a small number of deterministic snapshots.This results in a lower dimension of the eigenvalue problems (4.3) and, hence, makes them easier to solve numerically.After a multiplication with the respective hat-functions (partition of unity functions), we store the eigenfunctions corresponding to the n largest eigenvalues of the eigenvalue problems (4.3).These functions are then used as ansatz functions for computing the global approximation (4.5).Compared to the SLOD, by construction, no stability issues in the choice of basis can occur for the SL-GFEM and thus, no special treatment of boundary patches is required, cf.[HP22, App.B].
For the implementation, we use gridlod [HK17], which is a Python-library initially designed for the implementation of LOD-related methods.Although we do not require particular LOD-functionality from gridlod, it is convenient to use its flexible data structures for patches and its local discretization tools.Similarly, as in [KR21], our implementation can solve all local patch problems in parallel on an HPC cluster.As a comparison, we also implemented the SLOD from [HP22] in gridlod.All experiments are fully reproducible, and the corresponding source code can be found in [FHKP22]1 .7.2.Numerical experiments.We consider the domain Ω = (0, 1) 2 equipped with coarse Cartesian meshes T H and a fine Cartesian mesh T h obtained by uniform refinement of T H .Note that, for ease of presentation, H and h henceforth denote the elements side lengths instead of their diameters.For all numerical experiments, we use h = 2 −10 , which results in about one million degrees of freedom for the fine mesh.Note that our implementation also works for higher spatial dimensions d.For the numerical experiments, we consider two scalar diffusion coefficients A (realization of random field with short correlation length and high contrast channeled coefficient) and two source terms f (constant and non-polynomial).The precise definitions can be found in the respective experiments.Each configuration serves its own purpose for numerically investigating the SL-GFEM.For measuring the approximation quality, we use the relative energy error, i.e., e rel a (ũ) := where u h ∈ V h denotes the first order finite element approximation of (2.2) which we use as reference solution.Further, ũ is a placeholder for the SL-GFEM approximation (4.5) or the SLOD approximation (5.10).
7.2.1.Super-exponential localization.First, we investigate the localization properties of the SL-GFEM given several choices of the local approximation space size n.For the choice f ≡ 1, the optimal order term in Theorems 5.3 and 6.1 disappears and only the localization error and the approximation error d n are present.As coefficient A, we consider a realization of the random field taking piecewise constant values on T 2 −8 , which are independent and identically distributed in the interval [1, 100].This results in a maximum contrast of 100.We consider the fixed coarse mesh T 2 −5 and the polynomial degree p = 0.For several sizes of local approximation spaces n, Figure 7.1 depicts the relative energy errors of the SL-GFEM and the SLOD as a function of the oversampling parameter .
Clearly, the parameter n strongly impacts the approximation error of the SL-GFEM.One observes a large difference in the approximation quality, for instance, for n = 10 and n = 15.Conversely, choosing n = 20 does not yield a significantly better approximation than n = 15.This effect is related to the large jumps between the plateaus in Figure 5.1 and is also visible Theorems 5.5 and 6.2.In conclusion, for a sufficiently large n, the SL-GFEM shows a rapid decay of the localization error confirming Theorems 5.3 and 6.1 numerically.For the SLOD, the super-localization property [HP22, Sec.7] is visible in Figure 7.1.7.2.2.Optimal convergence.For investigating the convergence with respect to the coarse mesh size H, we use the same coefficient as in Section 7.2.1 but consider the non-polynomial source term f 1 (x 1 , x 2 ) := (x 1 + cos(3πx 1 )) • x 3 2 ∈ H 1 (Ω).Figure 7.2 depicts the errors of the SL-GFEM and SLOD for multiple choices of n and as a function of H.As a reference, a line with slope two indicates the expected order of convergence.For n and sufficiently large, one observes that the SL-GFEM converges with an order of two which numerically confirms Theorems 5.3 and 6.1.Notably, the errors of the SL-GFEM are smaller by nearly one order of magnitude than the errors of the SLOD.This effect only appears for non-trivial coefficients A, i.e., the effect is most probably related to the contrast of the coefficient.The contrast dependence is investigated more closely in the following subsection.7.2.3.High-contrast channeled coefficient.One of the major challenges for multi-scale methods is their sensitivity to high-contrast channeled coefficients.In this numerical experiment, we consider the coefficient A κ constructed by adding four channels of conductivity κ in the coefficient from Sections 7.2.1 and 7.2.2.Some of the channels touch the boundary, while others stop before.The number κ is the maximum contrast of A κ .Figure 7.3 illustrates the coefficient for κ = 10 5 , 10 7 .For this numerical experiment, we choose the same setup as in Section 7.2.1.For the SLOD, in contrast, the best practical realization known until now [HP22] yields a basis with deteriorating stability as κ is increased (growing constants in Assumptions 5.2 and 5.4).This explains the worse performance of the SLOD for κ = 10 7 compared to κ = 10 4 .Notably, when compared to Section 7.2.1, the SL-GFEM does not need more local functions n to attain a good approximation quality, which suggests that the choice of n is not by the contrast.
7.2.4.Higher-order polynomials.One key benefit of the proposed method is its flexibility with regard to the choice of polynomial degree, i.e., the construction of higher-order methods is straightforward.While the previous numerical experiments have investigated the performance of the SL-GFEM for p = 0, this experiment also considers higher polynomial degrees.Using the setup from Section 7.2.2, Figure 7.5 depicts the errors of the SL-GFEM for p = 0, 1, 2 as a function of H together with lines indicating the respective expected orders of convergence.For n and sufficiently large, it can be observed that the method of degree p converges with an order of p + 2 (recall that f is sufficiently smooth).This numerically confirms Theorems 5.3 and 6.1.Note that the choice of n needs to be adapted to p, which is related to the larger plateaus in Figure 5.1.We observed that n needs to be increased linearly as p is increased.It is left to future research to find a (possibly adaptive) choice of n such that pessimistic choices (that may result in unnecessarily many basis functions) can be avoided.

Figure 5 . 1 .
Figure 5.1.Singular values σ k of operator Π H | Y defined in (5.5) for an interior patch, for different pairs of , p.The singular values σ K−J+1 relevant for (5.7) are marked by dashed horizontal lines.
Theorem 6.2 (Bound on d n ).There exists C, C d > 0 independent of H, , and m such that, for m = 1, . . ., − 1 d n (H, ) ≤ C m d/2 exp(−C d m), where n ≈ m d .Proof.Let us consider a fixed node z ∈ N H and oversampling parameter .By [AHP21, Rem.3.7 ii], we can write any v ∈ V H,z as v = (1 − C z )B H v with the correction operator C z defined in (6.6).We define the patch-local localized correction operator by Cm := T ⊂ω z Cm T , where, denoting Wm T := W(ω m T ∩ ω z ), the element correctors Cm Figure 7.1.Localization errors of the SL-GFEM for multiple choices of n and of the SLOD for a fixed coarse mesh.

Figure 7 . 2 .
Figure 7.2.Convergence plot of the SL-GFEM and SLOD for multiple choices of n and .

Figure 7 . 5 .
Figure 7.5.Convergence plot of the higher-order SL-GFEM for multiple choices of n, , and p.
The work of Philip Freese, Moritz Hauck, and Daniel Peterseim is part of a project that has received funding from the European Research Council ERC under the European Union's Horizon 2020 research and innovation program (Grant agreement No. 865751).Tim Keil acknowledges funding by the Deutsche Forschungsgemeinschaft under Germany's Excellence Strategy EXC 2044 390685587, Mathematics Münster: Dynamics -Geometry -Structure.