Operator preconditioning: the simplest case
aa r X i v : . [ m a t h . NA ] F e b OPERATOR PRECONDITIONING: THE SIMPLEST CASE
ROB STEVENSON, RAYMOND VAN VENETI ¨EAbstract.
Using the framework of operator or Calder´on preconditioning, uniformpreconditioners are constructed for elliptic operators discretized with continuousfinite (or boundary) elements. The preconditioners are constructed as the compo-sition of an opposite order operator, discretized on the same ansatz space, and twodiagonal scaling operators. Introduction
This paper deals with the construction of uniform preconditioners for negativeand positive order operators, discretized by continuous piecewise polynomial trialspaces, using the framework of ‘operator preconditioning’ [Hip06], see also [SW98,Ste02, BC07, HJHUT20].For some d -dimensional closed domain (or manifold) Ω and an s ∈ [0 , , weconsider the (fractional) Sobolev space H s (Ω) and its dual that we denote by H − s (Ω) .Let ( S T ) T ∈ T be a family of continuous piecewise polynomials of some fixed degree ℓ w.r.t. uniformly shape regular, possibly locally refined, partitions.Given some families of uniformly boundedly invertible operators A T : (cid:0) S T , k · k H − s (Ω) (cid:1) → (cid:0) S T , k · k H − s (Ω) (cid:1) ′ ,B T : (cid:0) S T , k · k H s (Ω) (cid:1) → (cid:0) S T , k · k H s (Ω) (cid:1) ′ , we are interested in constructing a preconditioner for A T using operator precondi-tioning with B T , and vice versa. To this end, we introduce a uniformly bound-edly invertible operator D T : (cid:0) S T , k · k H − s (Ω) (cid:1) → (cid:0) S T , k · k H s (Ω) (cid:1) ′ , yielding pre-conditioned systems D − T B T ( D ′T ) − A T and ( D ′T ) − A T D − T B T that are uniformlyboundedly invertible.In earlier research, [SvV19, SvV20], we already constructed such precondition-ers in a more general setting where different ansatz spaces were used to define A T and B T . The setting studied in the current work, however, allows for precondition-ers with a remarkably simple implementation.A typical setting is that for some A : H s (Ω) → H − s (Ω) and B : H − s (Ω) → H s (Ω) , both boundedly invertible and coercive, it holds that ( A T u )( v ) := ( Au )( v ) and ( B T u )( v ) := ( Bu )( v ) with u, v ∈ S T . An example for s = is that A is the Sin-gle Layer Integral operator and B is the Hypersingular Integral operator. For thiscase, continuity of piecewise polynomial trial functions is required for discretizing Date : February 25, 2021.2010
Mathematics Subject Classification.
Key words and phrases.
Operator preconditioning, uniform preconditioners, finite- and boundaryelements.The second author has been supported by the Netherlands Organization for Scientific Research(NWO) under contract. no. 613.001.652. B , but not for A , for which often discontinuous piecewise polynomials are em-ployed. Nevertheless, when the solution of the Single Layer Integral equation isexpected to be smooth, e.g., when Ω is a smooth manifold, then it is advantageousto take an ansatz space of continuous (or even smoother) functions also for A .An obvious choice for D T would be to consider ( D T u )( v ) := h u, v i L (Ω) . How-ever, a problem becomes apparent when one considers the matrix representation D T of D T in the standard basis: the inverse matrix D − T , that appears in the pre-conditioned system, is densely populated. In view of application cost, this inversematrix has to be approximated, where it generally can be expected that, in orderto obtain a uniform preconditioner, approximation errors have to decrease with adecreasing (minimal) mesh size, which will be confirmed in a numerical experi-ment. To circumvent this issue, we will introduce a D T that has a diagonal matrixrepresentation, so that its inverse can be exactly evaluated.1.1. Notation.
In this work, by λ . µ we mean that λ can be bounded by a multipleof µ , independently of parameters which λ and µ may depend on, with the soleexception of the space dimension d , or in the manifold case, on the parametrizationof the manifold that is used to define the finite element spaces on it. Obviously, λ & µ is defined as µ . λ , and λ h µ as λ . µ and λ & µ .For normed linear spaces Y and Z , in this paper for convenience over R , L ( Y , Z ) will denote the space of bounded linear mappings Y → Z endowed with the op-erator norm k·k L ( Y , Z ) . The subset of invertible operators in L ( Y , Z ) with inversesin L ( Z , Y ) will be denoted as L is( Y , Z ) .For Y a reflexive Banach space and C ∈ L ( Y , Y ′ ) being coercive , i.e., inf = y ∈ Y ( Cy )( y ) k y k Y > , both C and ℜ ( C ) := ( C + C ′ ) are in L is( Y , Y ′ ) with kℜ ( C ) k L ( Y , Y ′ ) ≤ k C k L ( Y , Y ′ ) , k C − k L ( Y ′ , Y ) ≤ kℜ ( C ) − k L ( Y ′ , Y ) = (cid:16) inf = y ∈ Y ( Cy )( y ) k y k Y (cid:17) − . The set of coercive C ∈ L is( Y , Y ′ ) is denoted as L is c ( Y , Y ′ ) . If C ∈ L is c ( Y , Y ′ ) ,then C − ∈ L is c ( Y ′ , Y ) and kℜ ( C − ) − k L ( Y , Y ′ ) ≤ k C k L ( Y , Y ′ ) kℜ ( C ) − k L ( Y ′ , Y ) .Given a family of operators C i ∈ L is( Y i , Z i ) ( L is c ( Y i , Z i ) ), we will write C i ∈L is( Y i , Z i ) ( L is c ( Y i , Z i ) ) uniformly in i , or simply ‘uniform’, when sup i max( k C i k L ( Y i , Z i ) , k C − i k L ( Z i , Y i ) ) < ∞ , or sup i max( k C i k L ( Y i , Z i ) , kℜ ( C i ) − k L ( Z i , Y i ) ) < ∞ . Construction of D T in the domain case For some d -dimensional domain Ω and an s ∈ [0 , , we consider the Sobolevspaces H s (Ω) := [ L (Ω) , H (Ω)] s, , H − s (Ω) := H s (Ω) ′ , which form the Gelfand triple H s (Ω) ֒ → L (Ω) ≃ L (Ω) ′ ֒ → H − s (Ω) . PERATOR PRECONDITIONING: THE SIMPLEST CASE 3
Remark . In this work, for convenience we restrict ourselves to Sobolev spaceswith positive smoothness index which do not incorporate homogeneous Dirich-let boundary conditions and their duals. The proofs given below can howeverbe extended to the setting with boundary conditions, see the arguments foundin [SvV19, SvV20].Let ( T ) T ∈ T be a family of conforming partitions of Ω into (open) uniformly shaperegular d -simplices. Thanks to the conformity and the uniform shape regularity,for d > we know that neighbouring T, T ′ ∈ T , i.e. T ∩ T ′ = ∅ , have uniformlycomparable sizes. For d = 1 , we impose this uniform ‘ K -mesh property ’ explicitly.Fix ℓ > . For T ∈ T , let S T denote the space of continuous piecewise polyno-mials of degree ℓ w.r.t. T , i.e., S T := { u ∈ H (Ω) : u | T ∈ P ℓ ( T ∈ T ) } . Additionally, for r ∈ [ − , , we will write S T ,r as shorthand notation for thenormed linear space (cid:0) S T , k · k H r (Ω) (cid:1) .Denote N T for the set of the usual Lagrange evaluation points of S T , and equipthe latter space with Φ T = { φ T ,ν : ν ∈ N T } , being the canonical nodal basis definedby φ T ,ν ( ν ′ ) := δ νν ′ ( ν, ν ′ ∈ N T ). For T ∈ T , set h T := | T | /d and let N T := T ∩ N T be the set of evaluation points in T . We will omit notational dependence on T if itis clear from the context, e.g., we will simply write φ ν .2.1. Operator preconditioning.
Given some family of opposite order operators A T ∈ L is c ( S T , − s , ( S T , − s ) ′ ) and B T ∈ L is c ( S T ,s , ( S T ,s ) ′ ) , both uniformly in T ∈ T ,we are interested in constructing optimal preconditioners for both A T and B T , us-ing the idea of opposite order preconditioning ([Hip06]).That is, if one has an additional family of operators D T ∈ L is( S T , − s , ( S T ,s ) ′ ) uniformly in T ∈ T , then uniformly preconditioned systems for A T and B T aregiven by(2.1) D − T B T ( D ′T ) − A T ∈L is( S T , − s , S T , − s ) , ( D ′T ) − A T D − T B T ∈L is( S T ,s , S T ,s ) , see the following diagram: S T , − s ( S T , − s ) ′ ( S T ,s ) ′ S T ,sA T ( D ′T ) − D − T B T . In the following we shall be concerned with constructing a suitable family D T .2.1.1. An obvious but unsatisfactory choice for D T . An option would be to consider ( D T u )( v ) := h u, v i L (Ω) ( u, v ∈ S T ) , being uniformly in L ( S T , − s , ( S T ,s ) ′ ) . Forshowing boundedness of its inverse, let Q T be the L (Ω) -orthogonal projector onto S T then k D − T k − L (( S T ,s ) ′ , S T , − s ) = inf = u ∈ S T , − s sup = v ∈ H s (Ω) h u, v i L (Ω) k u k H − s (Ω) k Q T v k H s (Ω) ≥ k Q T k − L ( H s (Ω) ,H s (Ω)) , ROB STEVENSON, RAYMOND VAN VENETI ¨E
As follows from [SvV19, Prop. 2.3], the converse is also true, i.e., uniform bound-edness of k D − T k L (( S T ,s ) ′ , S T , − s ) is actually equivalent to uniform boundedness of k Q T k L ( H s (Ω) ,H s (Ω)) .This uniform boundedness of k Q T k L ( H s (Ω) ,H s (Ω)) is well-known for families of quasi-uniform , uniformly shape regular conforming partitions of Ω into say d -simplices.It has also been demonstrated for families of locally refined partitions, for d = 2 in-cluding those that are generated by the newest vertex bisection (NVB) algorithm,see [Car02, GHS16, DST20]. On the other hand, in [BY14] a one-dimensional coun-terexample was presented in which the L (Ω) -orthogonal projector on a familyof sufficiently strongly graded, although uniform K meshes, is not H (Ω) -stable.Thus, in any case uniform H (Ω) -stability cannot hold without any restrictions onthe grading of the meshes.Aside from this latter theoretical shortcoming, more importantly, there is a com-putational problem with the current choice of D T . The matrix representation of D T w.r.t. Φ T is the ‘mass matrix’ D T := h Φ T , Φ T i L (Ω) . Its inverse D − T , appearing inthe preconditioner, is densely populated, and therefore has to be approximated,where generally the error in such approximations has to decrease with a decreas-ing (minimal) mesh-size in order to arrive at a uniform preconditioner.2.2. Constructing a practical D T . To avoid the aforementioned problems, we shallconstruct D T ∈ L is( S T , − s , ( S T ,s ) ′ ) with a diagonal matrix representation. To thisend, we require some auxiliary space f S T ⊂ H (Ω) equipped with a local basis e Φ T that is L (Ω) -biorthogonal to Φ T and that has ‘approximation properties’. To beprecise, let e Φ T := { e φ ν ∈ H (Ω) : ν ∈ N T } be some collection that satisfies:(2.2) h e φ ν , φ ν ′ i L (Ω) = δ νν ′ h , φ ν i L (Ω) , X ν ∈ N T e φ ν = Ω , k e φ ν k H k (Ω) . k φ ν k H k (Ω) (cid:0) k ∈ { , } (cid:1) , supp e φ ν ⊆ supp φ ν . We will take D T := I ′T e D T with e D T and I T being defined and analyzed in the nexttwo theorems. Theorem 2.2.
The operator e D T : S T , − s → ( f S T ,s ) ′ , defined by ( e D T u )( v ) := h u, v i L (Ω) ,satisfies e D T ∈ L is( S T , − s , ( f S T ,s ) ′ ) uniformly in T ∈ T .Proof. This proof largely follows [SvV19, Sect. 3.1], but because here we considera Sobolev space H s (Ω) that does not incorporate homogeneous boundary condi-tions, it allows for an easier proof.From the assumptions (2.2), it follows that the biorthogonal ‘Fortin’ projector P T : L (Ω) → H (Ω) onto f S T with ran(Id − P T ) = S ⊥ L T exists, and is given by P T u = X ν ∈ N T h u, φ ν i L (Ω) h e φ ν , φ ν i L (Ω) e φ ν . Let T ∈ T , by (2.2) and the fact that h , φ ν i L (Ω) h k φ ν k L (Ω) , we find for k ∈ { , } (2.3) k P T u k H k ( T ) . X ν ∈ N T k e φ ν k H k ( T ) k φ ν k L (Ω) k u k L (supp φ ν ) . h − kT k u k L ( ω T ( T )) , This last condition can be replaced by e φ ν having (uniformly) local support. PERATOR PRECONDITIONING: THE SIMPLEST CASE 5 with ω T ( T ) := S { ν ∈ N T } supp φ ν . This shows sup T ∈ T k P T k L ( L (Ω) ,L (Ω)) < ∞ .From the above inequality, and P ν ∈ N T e φ ν = , we deduce that k (Id − P T ) u k H ( T ) = inf p ∈P k (Id − P T )( u − p ) k H ( T ) . inf p ∈P k u − p k H ( T ) + h − T k u − p k L ( ω T ( T )) . inf p ∈P h − T k u − p k L ( ω T ( T )) + | u | H ( T ) . | u | H ( ω T ( T )) , with the last step following from the Bramble-Hilbert lemma. We conclude that sup T ∈ T k P T k L ( H (Ω) ,H (Ω)) < ∞ , and consequently by the Riesz-Thorin interpola-tion theorem, that sup T ∈ T k P T k L ( H s (Ω) ,H s (Ω)) < ∞ . This latter property guarantees that e D T is uniformly boundedly invertible: k e D T k L ( S T , − s , ( f S T ,s ) ′ ) = sup = u ∈ S T , − s sup = v ∈ f S T ,s h u, v i L (Ω) k u k H − s (Ω) k v k H s (Ω) ≤ , k e D − T k − L (( f S T ,s ) ′ , S T , − s ) = inf = u ∈ S T , − s sup = v ∈ f S T ,s h u, v i L (Ω) k u k H − s (Ω) k v k H s (Ω) = inf = u ∈ S T , − s sup = v ∈ H s (Ω) h u, v i L (Ω) k u k H − s (Ω) k P T v k H s (Ω) ≥ k P T k − L ( H s (Ω) ,H s (Ω)) . (cid:3) Theorem 2.3.
For I T : S T ,s → f S T ,s being the bijection given by I T φ ν = e φ ν ( ν ∈ N T ) ,it holds that I T ∈ L is( S T ,s , f S T ,s ) uniformly in T ∈ T .Proof. Note that we may write I T u = X ν ∈ N T h u, e φ ν i L (Ω) h φ ν , e φ ν i L (Ω) e φ ν and I − T u = X ν ∈ N T h u, φ ν i L (Ω) h e φ ν , φ ν i L (Ω) φ ν . Equivalently to (2.3), we see for k ∈ { , } that k I T u k H k ( T ) . X ν ∈ N T k e φ ν k H k ( T ) k e φ ν k L (Ω) k φ ν k L (Ω) k u k L (supp φ ν ) . h − kT k u k L ( ω T ( T )) . Following the same arguments as in the proof of Theorem 2.2, using that I T = ,then reveals that I T is uniformly bounded. Uniformly boundedness of I − T followssimilarly. (cid:3) As announced earlier, we define D T ∈ L ( S T , − s , ( S T ,s ) ′ ) by D T := I ′T e D T , so ( D T u )( v ) := h u, I T v i L (Ω) ( u, v ∈ S T ). Combining the previous theorems givesthe following corollary. Corollary 2.4.
The operator D T is in L is( S T , − s , ( S T ,s ) ′ ) uniformly in T ∈ T .Remark . The matrix representation of D T w.r.t. Φ T given by D T = h Φ T , I T Φ T i L (Ω) = diag {h , φ ν i L (Ω) : ν ∈ N T } , ROB STEVENSON, RAYMOND VAN VENETI ¨E which is diagonal and therefore easily invertible. The matrix D T is known as thelumped mass matrix. Remark . The operator D T depends merely on the existence of a biorthogonal ba-sis e Φ T that satisfies (2.2). Indeed, this basis does not appear in the implementationof D T .A possible construction of e Φ T can be given using techniques from [SvV19]. Con-sider some collection of local ‘bubble’ functions Θ T = { θ ν ∈ H (Ω) : ν ∈ N T } thatsatisfy: (cid:12)(cid:12) h θ ν , φ ν ′ i L (Ω) (cid:12)(cid:12) h δ νν ′ k φ ν k L (Ω) , k θ ν k H k (Ω) . k φ ν k H k (Ω) ( k ∈ { , } ) , and supp θ ν ⊆ supp φ ν . Existence of such a collection can be shown by a constructionon a reference d -simplex, and then using an affine bijection to transfer it to generalelements, see [SvV19, Sect. 4.1]. A suitable e Φ T that satisfies (2.2) is then given by e φ ν := φ ν + h , φ ν i L (Ω) h θ ν , φ ν i L (Ω) θ ν − X ν ′ ∈ N T h φ ν , φ ν ′ i L (Ω) h θ ν ′ , φ ν ′ i L (Ω) θ ν ′ . We emphasize that the construction of a uniform preconditioner outlined in thesubsection does require any assumptions on the mesh grading.2.2.1.
Implementation.
Taking Φ T as basis for both S T , − s and S T ,s , the matrix rep-resentation of the preconditioned systems from (2.1) read as D − T B T D −⊤T A T and D −⊤T A T D − T B T , where A T := ( A T Φ T )(Φ T ) , B T := ( B T Φ T )(Φ T ) , D T := ( D T Φ T )(Φ T ) = diag {h , φ ν i L (Ω) : ν ∈ N T } . Alternatively, we could equip the spaces with the scaled nodal basis ˘Φ T := D − T Φ T ,so that the L (Ω) -norm of any basis function is proportional to , yielding ˘ A T := ( A T ˘Φ T )( ˘Φ T ) = ( D − T ) ⊤ A T D − T , ˘ B T := ( B T ˘Φ T )( ˘Φ T ) = ( D − T ) ⊤ B T D − T , ˘ D T := ( D T ˘Φ T )( ˘Φ T ) = ( D − T ) ⊤ D T D − T = Id , showing that ˘ B T is a uniform preconditioner for ˘ A T (and vice versa). To the bestof our knowledge, so far this most easy form of operator preconditioning, wherethe stiffness matrix of some operator w.r.t. some basis is preconditioned by stiffnessmatrix of an opposite order operator w.r.t. the same basis, has not been shown tobe optimal. 3. Manifold case
Let Γ be a compact d -dimensional Lipschitz, piecewise smooth manifold in R d ′ for some d ′ ≥ d without boundary ∂ Γ . For s ∈ [0 , , we consider the Sobolevspaces H s (Γ) := [ L (Γ) , H (Γ)] s, , H − s (Γ) := H s (Γ) ′ . PERATOR PRECONDITIONING: THE SIMPLEST CASE 7
We assume that Γ is given as the closure of the disjoint union of ∪ pi =1 χ i (Ω i ) , with,for ≤ i ≤ p , χ i : R d → R d ′ being some smooth regular parametrization, and Ω i ⊂ R d an open polytope. W.l.o.g. assuming that for i = j , Ω i ∩ Ω j = ∅ , we define χ : Ω := ∪ pi =1 Ω i → ∪ pi =1 χ i (Ω i ) by χ | Ω i = χ i . Let T be a family of conforming partitions T of Γ into ‘panels’ such that, for ≤ i ≤ p , χ − ( T ) ∩ Ω i is a uniformly shape regular conforming partition of Ω i into d -simplices (that for d = 1 satisfies a uniform K -mesh property).Fix ℓ > , we set S T := { u ∈ H (Γ) : u ◦ χ | χ − ( T ) ∈ P ℓ ( T ∈ T ) } , equipped with the canonical nodal basis Φ T = { φ ν : ν ∈ N T } .For construction of an operator D T ∈ L is( S T , − s , ( S T ,s ) ′ ) one can proceed as inthe domain case. A suitable collection e Φ T that is L (Γ) -biorthogonal to Φ T exists.Moreover, the analysis from the domain case applies verbatim by only changing h· , ·i L (Ω) into h· , · , i L (Γ) . A hidden problem , however, is that the computation of D T = diag {h , φ ν i L (Γ) : ν ∈ N T } involves integrals over Γ that generally have tobe approximated using numerical quadrature.In [SvV19] we solved this issue by defining an additional ‘mesh-dependent’scalar product h u, v i T := X T ∈T | T || χ − ( T ) | ˆ χ − ( T ) u ( χ ( x )) v ( χ ( x )) dx. This is constructed by replacing on each χ − ( T ) , the Jacobian | ∂χ | by its average | T || χ − ( T ) | over χ − ( T ) .By considering e Φ T that is biorthogonal to Φ T with respect to h· , ·i T , and thelinear bijection I T given by I T φ ν = e φ ν , one is able show that the operator D T defined as ( D T u )( v ) := h u, I T v i T satisfies the necessary requirements. For detailswe refer to [SvV19]. The resulting matrix representation of D T w.r.t. Φ T is thengiven by D T = diag {h , φ ν i T : ν ∈ N T } .4. Numerical results
Let
Γ = ∂ [0 , ⊂ R be the two-dimensional manifold without boundary givenas the boundary of the unit cube, s = , and S T the space of continuous piecewisepolynomials of degree ℓ w.r.t. a partition T . We will evaluate preconditioning ofthe discretized Single Layer Integral operator A T ∈ L is c ( S T , − s , ( S T , − s ) ′ ) and an(essentially) discretized Hypersingular Integral operator B T ∈ L is c ( S T ,s , ( S T ,s ) ′ ) .The Hypersingular Integral operator ˜ B ∈ L ( H (Γ) , H − (Γ)) , is only-semi co-ercive, but solving ˜ Bu = f for f with f ( ) = 0 is equivalent to solving Bu = f with B given by ( Bu )( v ) = ( ˜ Bu )( v ) + α h u, i L (Γ) h v, i L (Γ) , for some fixed α > .This operator B is in L is c ( H (Γ) , H − (Γ)) , and we shall consider discretizations B T ∈ L is c ( S T ,s , ( S T ,s ) ′ ) of B . We found α = 0 . to give good results in ourexamples.Equipping both S T ,s and S T , − s with the standard nodal basis Φ T = { φ ν : ν ∈ N T } , the matrix representations of the preconditioned systems from Sect. 2.2 readas D − T B T D −⊤T A T and D −⊤T A T D − T B T , ROB STEVENSON, RAYMOND VAN VENETI ¨E for D T = diag {h , φ ν i L (Γ) : ν ∈ N T } , A T = ( A T Φ T )(Φ T ) and B T := ( B T Φ T )(Φ T ) .We calculated (spectral) condition numbers of these preconditioned systems,where this condition number is given by κ S ( X ) := ρ ( X ) ρ ( X − ) with ρ ( · ) denotingthe spectral radius. Note that the condition numbers of the preconditioned systemscoincide, i.e., κ S ( D − T B T D −⊤T A T ) = κ S ( D −⊤T A T D − T B T ) , so we may restrict ourselves to results for preconditioning of A T .We used the BEM++ software package [ ´SBA +
15] to approximate the matrixrepresentation of A T and B T by hierarchical matrices based on adaptive cross ap-proximation [Hac99, Beb00].As initial partition T ⊥ = T of Γ we take a conforming partition consisting of triangles per side, so triangles in total, with an assignment of the newest verticesthat satisfies the so-called matching condition. We let T be the sequence {T k } k ≥ where the (conforming) partition T k is found by applying both uniform and localrefinements. To be precise, T k is constructed by first applying k uniform bisectionsto T ⊥ , and then k local refinements by repeatedly applying NVB to all trianglesthat touch a corner of the cube.4.1. Comparison preconditioners.
Write G D T := D − T B T D −⊤T for the precondi-tioner constructed in Sect. 2.2. We will compare this with the preconditioner de-scribed in Sect. 2.1.1, for which the matrix representation is given by G M T := M − T B T M −⊤T with mass matrix M T = h Φ T , Φ T i L (Γ) . Because our partitions of the two-dimensionalsurface are created with NVB, we know that also the latter preconditioner providesuniformly bounded condition numbers. In contrast to D − T , the inverse M − T can-not be evaluated in linear complexity. We implemented the application of M − T bycomputing an LU-factorization of M T .Table 4.1 compares the spectral condition numbers for the preconditioned SingleLayer systems with trial spaces given by continuous piecewise linears and those bycontinuous piecewise cubics. The condition numbers κ S ( G D T A T ) are uniformlybounded, but quantitatively the condition numbers κ S ( G M T A T ) are better.4.2. Improving the preconditioner quality.
As observed in Table 4.1, the precon-ditioner G M T appears to be of superior quality, but it has unfavourable computa-tional complexity. It does suggest a way for improving G D T : by replacing D − T with a better approximation of M − T , one may hope to improve the quality. Tothis end, we introduce damped (preconditioned) Richardson. Let < λ − ≤ λ min ( D − T M T ) , λ max ( D − T M T ) ≤ λ + , R (0) T := 0 and for k ≥ define R ( k +1) T := R ( k ) T + ω D − T (Id − M T R ( k ) T ) , ω = 2 λ − + λ + , being the result of k Richardson iterations. Correspondingly define(4.1) G ( k ) T := R ( k ) T B T R ( k ) T . It follows that G (1) T = G D T and lim k →∞ G ( k ) T = G M T . Although we have no proof,we suspect that G ( k ) T provides a uniform preconditioner for A T due to the factthat R ( k ) T approximates M − T , while preserving constant functions, being a keyingredient in the proofs of Theorems 2.2 and 2.3. PERATOR PRECONDITIONING: THE SIMPLEST CASE 9
Table 1.
Spectral condition numbers, κ S ( G ◦T A T ) for ◦ ∈ { D, M } ,of the preconditioned Single Layer system discretized on {T k } k ≥ ,by continuous piecewise linears ( ℓ = 1) in the middle columns anddiscretized by continuous piecewise cubics ( ℓ = 3) in the rightcolumns. Here G D T is the preconditioner introduced in Sect. 2.2,whereas G M T is the preconditioner described in Sect. 2.1.1 whoseapplication requires an application of M − T , which we imple-mented using an LU-factorization.Partition T Linears ( ℓ = 1 ) Cubics ( ℓ = 3 ) h min h max dofs G D T A T G M T A T dofs G D T A T G M T A T . · . · . .
20 56 90 . . . · − . · −
218 14 . .
91 1946 87 . . . · − . · −
482 14 . .
04 4322 86 . . . · − . · −
962 14 . .
10 8642 85 . . . · − . · − . .
14 20738 84 . . . · − . · − . .
16 63938 84 . . . · − . · − . .
17 231554 84 . . . · − . · − . .
17 896834 84 . . Values for λ − and λ + can be found by calculating the extremal eigenvalues of thecorresponding preconditioned mass matrix on a reference simplex, see e.g. [Wat87].For ℓ = 1 this gives ω = d +2) d +3 , whereas for ℓ = 3 and d = 2 we computed ω = 0 . .Table 4.2 compares the condition numbers κ S ( G ( k ) T A T ) for k ∈ { , , } . We seethat a few Richardson iterations drastically improves our preconditioner, makingits quality on par with that of G M T while having a favourable linear application cost.Finally, to show that one cannot simply use any (iterative) method for approx-imating M − T , we consider the case where one approximates this inverse using aJacobi preconditioner. The resulting preconditioner is then given by(4.2) G J T := (diag M T ) − B T (diag M T ) −⊤ . Table 4.2 clearly displays that this is not a uniformly bounded preconditioner, whichwe assume is due to the fact that (diag M T ) − does not preserve constant functionsfor ℓ > . 5. Conclusion
Considering discretized opposite order operators A T and B T using the sameansatz space of continuous piecewise polynomial w.r.t. a possibly locally refinedpartition T , we consider matrices D T such that D − T B T D −⊤T is a uniform precon-ditioner for A T , and D −⊤T A T D − T for B T . The obvious choice for D T would bethe mass matrix, however, it yields uniformly bounded condition numbers onlyunder a (mild) grading assumption on the mesh, and more importantly, it has thedisadvantage that its inverse is dense. We proved that when taking D T as thelumped mass matrix the condition numbers are uniformly bounded, remarkablywithout any gradedness assumption on the mesh, while obviously its inverse canbe applied in linear cost. Table 2.
Spectral condition numbers κ S ( G ( k ) T A T ) with G ( k ) T thepreconditioner from (4.1) that incorporates k Richardson itera-tions. The systems are discretized by continuous piecewise linearsin the left columns and discretized by continuous piecewise cubicsin the right columns.Linears ( ℓ = 1 ) Cubics ( ℓ = 3 )dofs k = 2 k = 4 k = 6 dofs k = 2 k = 4 k = 68 2 .
26 1 .
29 1 .
22 56 10 . .
99 2 . .
05 2 .
07 1 .
94 1946 8 .
96 3 .
57 2 . .
53 2 .
28 2 .
08 4322 8 .
80 3 .
59 2 . .
79 2 .
44 2 .
19 8642 8 .
63 3 .
59 2 . .
98 2 .
52 2 .
24 20738 8 .
54 3 .
59 2 . .
18 2 .
57 2 .
27 63938 8 .
54 3 .
59 2 . .
35 2 .
61 2 .
28 231554 8 .
54 3 .
59 2 . .
47 2 .
65 2 .
29 896834 8 .
54 3 .
59 2 . Table 3.
Spectral condition numbers κ S ( G J T A T ) with G J T from (4.2), and systems discretized by continuous piecewise cu-bics ( ℓ = 3) . dofs G J T A T
56 62 . . . . . . In our experiments with locally refined meshes generated by Newest Vertex Bi-section, the condition numbers with D T as the mass matrix are quantitatively bet-ter than those found with D T as the lumped mass matrix though. Constructing D − T as an approximation for the inverse mass matrix by a few preconditioneddamped Richardson steps with the lumped mass matrix as a preconditioner, boththe resulting matrix can be applied at linear cost and the observed condition num-bers are essentially as good as with the inverse mass matrix. References [BC07] A. Buffa and S.H. Christiansen. A dual finite element complex on the barycentric refine-ment.
Math. Comp. , 76(260):1743–1769, 2007.[Beb00] M. Bebendorf. Approximation of boundary element matrices.
Numerische Mathematik ,86(4):565–589, 2000.[BY14] R.E. Bank and H. Yserentant. On the H -stability of the L -projection onto finite elementspaces. Numer. Math. , 126(2):361–381, 2014.[Car02] C. Carstensen. Merging the Bramble-Pasciak-Steinbach and the Crouzeix-Thom´ee cri-terion for H -stability of the L -projection onto finite element spaces. Math. Comp. ,71(237):157–163, 2002.
PERATOR PRECONDITIONING: THE SIMPLEST CASE 11 [DST20] L. Diening, J. Storn, and T. Tscherpel. On the Sobolev and L p -stability of the L -projection,2020.[GHS16] F. D. Gaspoz, C.-J. Heine, and K. G. Siebert. Optimal grading of the newest vertex bisectionand H -stability of the L -projection. IMA J. Numer. Anal. , 36(3):1217–1241, 2016.[Hac99] W. Hackbusch. A sparse matrix arithmetic based on H -matrices. Part i: Introduction to H -matrices. Computing , 62(2):89–108, 1999.[Hip06] R. Hiptmair. Operator preconditioning.
Comput. Math. Appl. , 52(5):699–706, 2006.[HJHUT20] R. Hiptmair, C. Jerez-Hanckes, and C. Urz ´ua-Torres. Optimal operator preconditioningfor Galerkin boundary element methods on 3-dimensional screens.
SIAM J. Numer. Anal. ,58(1):834–857, 2020.[´SBA +
15] W. ´Smigaj, T. Betcke, S. Arridge, J. Phillips, and M. Schweiger. Solving boundary integralproblems with BEM++.
ACM Transactions on Mathematical Software (TOMS) , 41(2):1–40,2015.[Ste02] O. Steinbach. On a generalized L projection and some related stability estimates inSobolev spaces. Numer. Math. , 90(4):775–786, 2002.[SvV19] R.P. Stevenson and R. van Veneti¨e. Uniform preconditioners for problems of negative order.
Math. Comp. , 2019.[SvV20] R.P. Stevenson and R. van Veneti¨e. Uniform preconditioners for problems of positive order.
Comput. Math. Appl. , 79(12):3516–3530, 2020.[SW98] O. Steinbach and W. L. Wendland. The construction of some efficient preconditioners in theboundary element method.
Adv. Comput. Math. , 9(1-2):191–216, 1998. Numerical treatmentof boundary integral equations.[Wat87] A. J. Wathen. Realistic eigenvalue bounds for the Galerkin mass matrix.
IMA Journal ofNumerical Analysis , 7(4):449–457, 1987.
Korteweg-de Vries Institute for Mathematics, University of Amsterdam, P.O. Box 94248, 1090 GEAmsterdam, The Netherlands
Email address ::