⬇ Download Markdown  ·  ⬇ Download LaTeX  ·  ⬇ original  ·  view original  ·  CSV: table 1, table 2, table 3

H 2 -matrices for translation-invariant kernel functions

Steffen B¨orm and Janne Henningsen

arXiv:2210.16609v2 [math.NA]

7

Jun 2024

June 10, 2024 Boundary element methods for elliptic partial differential equations typ- ically lead to boundary integral operators with translation-invariant kernel functions. Taking advantage of this property is not straightforward if general unstructured meshes and general basis functions are used, since we need the supports of these basis functions to be contained in a hierarchy of subdomains with translational symmetry. In this article, we present a modified construction for H2 -matrices on un- structured quasi-uniform meshes that uses translation-invariance to signifi- cantly reduce the storage requirements for the farfield representation. We construct a nested hierarchy of axis-parallel boxes so that translational symmetry is preserved and prove optimal-order complexity estimates under moderate assumptions. In particular, we need only one weak assumption for proving that the entire farfield requires only O(log(n)) coefficients. It should be mentioned that, since we are working with an unstructured mesh and general basis functions, the nearfield of the matrix still requires O(n) units of storage. This work was funded by the DFG in project BO 3289/7-1. Keywords: integral equations, data-sparse approximation, hierarchical matrices MSC codes: 65N38, 65R20, 45B05, 65D05, 65D15, 41A10, 41A63

1 Introduction

Boundary element methods are an attractive technique for handling homogeneous linear partial differential equations, e.g., the Laplace, Lam´e, Helmholtz, or Maxwell equations [24, 27, 30], on a domain Ω ⊆ Rd with d ∈ {2, 3}. A Galerkin discretization of related boundary integral operators based on test functions (ϕi ) i ∈I and trial functions (ψj )j ∈J leads to matrices G ∈ RI ×J given by Gi j := Γ ϕi ( x) Γ g(x, y)ψj (y) dy dx for all i ∈ I, j ∈ J, (1) 1

where Γ := ∂Ω denotes the boundary of Ω and g : Rd × Rd → R is a kernel function. In typical applications, the kernel function g is non-zero almost everywhere, and the matrix G is therefore densely populated. In order to reduce the storage requirements and computational complexity, compression techniques are employed, e.g., the panel- clustering technique [26], the fast multipole method [28, 21, 1, 22, 33], interpolation [18, 10, 12, 13, 17], algebraic approximations [31, 32, 3, 4, 14], or hybrid methods [19, 9, 8]. All of these techniques split the matrix G into submatrices that can belong either to the nearfield or the farfield. Nearfield submatrices are small and can be stored directly, while farfield submatrices can be large and have to be approximated, e.g., by low-rank matrices that can be stored efficiently in factorized form. Unfortunately, even the factorized form still requires a large amount of storage, and while recompression techniques [11, 20, 7] can help, they require additional computational work. In this article, we pursue an alternative approach frequently used in the context of fast multipole methods for particle systems: if the kernel function g is invariant under translation, i.e., if g(x, y) = g(x + c, y + c) for all x, y, c ∈ Rd , (2) we can modify interpolation and certain hybrid methods in a way that reduces the storage requirements for the farfield matrices to O(k2 log(n)), where k denotes the rank of the approximation and n := max{|I|, |J|} the maximal dimension of the matrix G. This task is straightforward for particle systems with approximately uniformly dis- tributed particles, since each particle is represented by just one point in space, but it is significantly more challenging for Galerkin discretizations, since every matrix entry depends on the supports of the basis functions ϕi and ψj , and the supports of different basis functions may overlap or have a non-trivial shape. In some implementations, this problem is circumvented by replacing the integrals in eq. (1) by quadrature, thus reduc- ing the computation to a sum of kernel evaluations closely related to particle methods [5]. While this is certainly an elegant approach and allows BEM codes to take advan- tage of highly sophisticated implementations of the fast multipole method, the number of quadrature points may grow excessively large if higher quadrature orders are used to keep up with the discretization error as the grid is refined. In this article, we present a different approach: Initially, we construct a hierarchy of axis-parallel boxes that cover representative points of the supports of the basis functions, and we ensure that all boxes on a given level are translations of the same reference box with a fixed displacement step size. By suitably enlarging the reference boxes and applying the aforementioned translations to these enlarged reference boxes we then construct a corresponding hierarchy of axis-parallel boxes that cover the entire supports. If we apply standard techniques like interpolation [18, 10] or Green quadrature [8] to this structure, matrices corresponding to translation-equivalent pairs of boxes are identical and therefore have to be stored only once. As long as the enlargement of the reference boxes is sufficiently bounded, this allows us to reduce the corresponding storage requirements to O(k2 ) for each level. No further assumption is needed for proving this result. 2

This approach also reduces the assembly time, since the matrices have to be set up only once. The run-time for the matrix-vector multiplication may benefit, too, since fewer data has to be moved between main memory and the processor, although the number of arithmetic operations usually grows since guaranteeing translational symmetry limits our flexibility when constructing the boxes. That is why a major part of this article is devoted to proving that our proposed construction does not change the asymptotically optimal complexity of the H2 -matrix method. In section 2 we introduce H2 -matrices, a suitable representation of our matrix approx- imation, and outline the modifications required to take full advantage of translation- invariance. Section 3 contains the key result of this article: we prove that the entire farfield of the modified H2 -matrices requires only O(k2 log(n)) coefficients if the enlarge- ment of the reference boxes is sufficiently bounded. The nearfield matrices and the leaf matrices are directly connected to the unstructured surface mesh and therefore cannot take advantage of the translation-invariance property, so we store these matrices explic- itly using O(nk) units of storage. Section 4 illustrates the advantage of the modified representation in a series of numerical experiments.

2 H 2 -matrices for translation-invariant kernel functions

In order for our approximation technique to achieve the desired accuracy, the kernel function g needs to be asymptotically smooth, i.e., there have to be constants C as ∈ R ≥0, c 0 ∈ R >0 and a singularity degree σ ∈ N such that for all (x, y) ∈ Rd × Rd with x ̸= y the following condition is fulfilled: |∂νι g(x, y)| ≤ C as (σ − 1 + ν)! cν 0 ∥x − y∥ν +σ 2 for all ν ∈ N 0, ι ∈ {1, ..., 2d}. (3) Kernel functions occurring in typical applications are known to be asymptotically smooth, see for example [25, Appendix E] and [2]. Additionally, we assume that g satisfies the translation-invariance property (2).

Definition 1 (Tree notations) Let T be a tree. We use the notation t ∈ T for ”t is a node in T ”. For each node t ∈ T we denote the set of its sons by sons T (t). The set of

a node in T ”. For each node t ∈ T we denote the set of its sons by sonsT (t). The set of leaves of T is denoted by ℒ T := {t ∈ T : sonsT (t) = ∅}. We denote the root of T by ℛ T

and define the level of a node by

levelT (t) := 0 if t = ℛ T , l evelT (ˇt ) + 1 if t has a father ˇt ∈ T for all t ∈ T . For all ℓ ∈ N 0 we define T (ℓ) := {t ∈ T : levelT (t) = ℓ}. The first phase of our approximation procedure relies on characteristic points (xi ) i ∈I and (yj )j ∈J satisfying xi ∈ supp(ϕi ) for all i ∈ I and yj ∈ supp(ψj ) for all j ∈ J. 3

These points are split hierarchically into a hierarchy of boxes that will then give rise to a decomposition of the matrix into submatrices. We choose a maximal level ℓm ax ∈ N ≥d and construct trees TI and TJ of closed axis- parallel boxes in Rd satisfying t = r( ℓ) + δ( ℓ) ⊙ pt =: m t = r( ℓ) + mt for all t ∈ T (ℓ) I ∪ T (ℓ) J (4) on every level ℓ ∈ {0, ..., ℓm ax }, where (r( ℓ) )ℓ m ax ℓ =0 is a family of reference boxes, δ( ℓ) ∈ Rd consists of the interval lengths of r( ℓ) for each ℓ ∈ {0, ..., ℓm ax } and (pt ) t ∈TI ∪TJ are suitable integer vectors. We start with computing an axis-parallel box r( 0) ⊆ Rd containing Γ. In order to create the isotropic boxes required by our complexity analysis, we then compute (r( ℓ) )ℓ m ax ℓ =1 recursively by consecutive splitting in each coordinate direction. The corresponding cyclic sequence of splitting directions is given by ι( 0) := 1, ι( ℓ) := ι ( ℓ−1) + 1 if ι( ℓ−1) < d, 1 if ι( ℓ−1) = d for all ℓ ∈ N. (5) Assuming that r( ℓ) = [a( ℓ) 1 , b( ℓ) 1 ] × ... × [a( ℓ) d , b( ℓ) d ] on a level ℓ ∈ {0, ..., ℓm ax − 1} is already given, we compute r( ℓ+1) := {z ∈ r( ℓ) : zι ( ℓ) ≤ c( ℓ) ι( ℓ) } with c( ℓ) ι( ℓ) := a( ℓ) ι( ℓ) + b( ℓ) ι( ℓ) 2 . (6) If ℓ + 1 < ℓm ax holds, we proceed with r( ℓ+1) by recursion. Otherwise, we stop. Next, we create the tree TI by shifting the reference boxes r( ℓ) in a regular pattern. In this stage, we also construct subsets (It ) t ∈TI of I such that xi ∈ t for all i ∈ It . (7) Beginning with r( 0) as the root ℛ TI of TI , pr ( 0) := 0 and Ir ( 0) := I, we construct TI also recursively: Assuming that a box t = [at ,1, bt ,1] × ... × [at ,d, bt ,d] ∈ T (ℓ) I on a level ℓ ∈ {0, ..., ℓm ax − 1} with t = r( ℓ) + δ( ℓ) ⊙ pt for a vector pt ∈ Nd 0 and a corresponding non-empty subset It of I with (7) are already given, we compute the midpoint ct ,ι( ℓ) := (at ,ι( ℓ) + bt ,ι( ℓ) )/2 and t1 := {z ∈ Rd : zι ( ℓ) ≤ ct ,ι( ℓ) } = r( ℓ+1) + δ( ℓ+1) ⊙ pt 1 , t2 := {z ∈ Rd : zι ( ℓ) ≥ ct ,ι( ℓ) } = r( ℓ+1) + δ( ℓ+1) ⊙ pt 2 with pt 1 , pt 2 ∈ Nd 0 given by pt 1 , ι := 2 pt ,ι if ι = ι( ℓ) p t ,ι if ι ̸= ι( ℓ) , pt 2 , ι := 2 pt ,ι + 1 if ι = ι( ℓ) p t ,ι if ι ̸= ι( ℓ) for all ι ∈ {1, ..., d}. (8) 4

We let It 1 := It \ It 2 , It 2 := ∅ if {xi : i ∈ It } ⊆ t1 { i ∈ It : xi ∈ t2 } otherwise (9) and define sonsT I (t) := {t′ ∈ {t1 , t2 } : It ′ ̸= ∅}. If ℓ + 1 < ℓm ax holds, we proceed with all sons t′ ∈ sonsT I (t) by recursion. Otherwise, we stop the recursion. By replacing I with J and (xi ) i ∈I with (yj )j ∈J we analogously compute the tree TJ along with a corresponding hierarchical family (Js )s ∈TJ of non-empty subsets of J. In conjunction with the families (It ) t ∈TI and (Js )s ∈TJ , our trees TI and TJ form cluster trees for I and J, respectively, which implies (cf. [6, Corollary 3.9]) I = ˙ t∈ℒT I It , J = ˙ s∈ℒT J Js . (10) Furthermore, by construction, we have ℒ TI = T (ℓm ax) I , ℒ TJ = T (ℓm ax) J . (11) In the second phase we create corresponding support bounding boxes (Bt ) t ∈TI and (Cs ) s ∈TJ with (see fig. 1) a property based on eq. (4): For every level ℓ ∈ {0, ..., ℓm ax } we construct the smallest closed axis-parallel box B (ℓ) containing r( ℓ) and the shifted supports (supp(ϕi ) − mt ) i ∈It for all t ∈ T (ℓ) I . Likewise, for every level ℓ ∈ {0, ..., ℓm ax } we construct the smallest closed axis-parallel box C (ℓ) containing r( ℓ) and the shifted supports (supp(ψj ) − ms ) j ∈Js for all s ∈ T (ℓ) J . According to (4), the properties r( ℓ) ∪ i∈It supp(ϕi ) − mt ⊆ B (ℓ) for all t ∈ T (ℓ) I , (12a) r( ℓ) ∪ j∈Js supp(ψj ) − ms ⊆ C (ℓ) for all s ∈ T (ℓ) J (12b) then lead to t ∪ i∈I t supp(ϕi ) ⊆ B (ℓ) + mt =: Bt for all t ∈ T (ℓ) I , (13a) s ∪ j∈Js supp(ψj ) ⊆ C (ℓ) + ms =: Cs for all s ∈ T (ℓ) J . (13b) In the third phase we construct an approximation of G by combining the trees, boxes and index sets created before with tensor interpolation. For this purpose we choose a degree θ ∈ N and define Θ := {0, ..., θ}d . For every axis-parallel box Q ⊂ Rd let I Q : C(Q) → Πθ , f → ν∈Θ f (ξQ ,ν )LQ ,ν , 5

r( 6) mt t Γ B (6) mt Bt Γ Figure 1: Relation between the boxes of TI and their support bounding boxes. be the corresponding tensor Chebyshev interpolation operator mapping to the set Πθ of tensor polynomials of degree θ with associated interpolation points (ξQ ,ν )ν ∈Θ and Lagrange polynomials (LQ ,ν )ν ∈Θ .

We use the admissibility condition

Aη : TI × TJ → {true, false}, (t, s) → t rue if max{diam(Bt ) , diam(Cs ) } ≤ η dist(Bt , Cs ) f alse if max{diam(Bt ) , diam(Cs ) } > η dist(Bt , Cs ) (14) with a given parameter η ∈ R >0 to decide whether the kernel function g can be approx- imated in a domain Bt × Cs . For every pair (t, s) ∈ TI × TJ with Aη (t, s) = true the corresponding tensor interpolation error satisfies ∥g − (IB t ⊗ I Cs )[g]∥∞ ,Bt ×Cs ≤ C in dist(Bt , Cs ) σ qθ +1 η for a constant C in ∈ R >0 and qη := min c 0η c 0η +2 , c0 η 4 (cf. [6, Remark 4.23]) and therefore converges exponentially to zero with respect to the degree θ. The actual convergence rates that appear in practice are often considerably better than qη . For every t ∈ TI let now Vt ∈ RI t × Θ be given by (Vt ) i ν := Γ ϕi ( x)LB t , ν (x) dx for all i ∈ It , ν ∈ Θ, and for every s ∈ TJ let Ws ∈ RJ s × Θ be given by (Ws ) j µ := Γ ψj (y)LC s , µ( y) dy for all j ∈ Js , µ ∈ Θ. Then for every pair (t, s) ∈ TI × TJ with Aη (t, s) = true we get, by using (1) and (13), 6

for all i ∈ It , j ∈ Js Gi j ≈ Γ ϕi ( x) Γ (IB t ⊗ I Cs )[g](x, y)ψj (y) dy dx = ν∈Θ µ∈Θ Γ ϕi ( x)LB t ,ν (x) dx =(Vt ) i ν g(ξB t ,ν , ξC s , µ) Γ L Cs , µ( y)ψj (y) dy =(Ws ) j µ , (15) which leads to G|I t × Js ≈ Vt S ( t,s)W ∗ s (16) with S( t,s) ∈ RΘ ×Θ given by (S( t,s)) ν µ := g(ξB t ,ν , ξC s , µ) for all ν ∈ Θ, µ ∈ Θ, i.e., we have a factorized low-rank approximation of admissible submatrices G|I t ×Js . We compute a tree TI ×J ⊆ TI × TJ defined by the following properties: • ℛ TI ×J = (ℛT I , ℛ TJ ) = (r( 0) , r( 0) ); • For all (t, s) ∈ TI ×J the sons are given by sonsT I ×J (t, s) = ∅ if Aη (t, s) = true, s onsT I (t) × sonsT J (s) otherwise. Note that this implies sonsT I ×J (t, s) = ∅ if sonsT I (t) = ∅ or sonsT J (s) = ∅. Its admissible leaves and inadmissible leaves are denoted by ℒ + TI ×J := {(t, s) ∈ ℒ TI ×J : Aη (t, s) = true}, ℒ − TI ×J := {(t, s) ∈ ℒ TI ×J : Aη (t, s) = false}. Drawing on the properties of the tree TI and TJ , we find (cf. [6, Corollary 3.15]) I × J = ˙ (t,s)∈ℒT I×J (It × Js ) = ˙ (t,s)∈ℒ + T I ×J (It × Js ) ˙∪ ˙ (t,s)∈ℒ − T I ×J (It × Js ). We can therefore use (16) to construct an approximation G ∈ RI ×J of G by G|I t × Js := Vt S ( t,s)W ∗ s for all (t, s) ∈ ℒ + TI ×J , G|I t × Js := G|I t × Js for all (t, s) ∈ ℒ − TI ×J . This means that for admissible blocks (t, s) ∈ ℒ + TI ×J , only a small coupling matrix S( t,s) ∈ RΘ ×Θ has to be stored. 7

The row cluster basis (Vt ) t ∈TI and the column cluster basis (Ws ) s ∈TJ can be stored efficiently using the identity theorem for polynomials: we have L Bt , ν = µ∈Θ L Bt , ν (ξB t ′ ,µ) LB t ′ ,µ for all t ∈ TI \ ℒ TI , t′ ∈ sonsT I (t), ν ∈ Θ. Consequently, for every t ∈ TI \ ℒ TI , t′ ∈ sonsT I (t), i ∈ It ′ , ν ∈ Θ we get (Vt ) i ν = Γ ϕi ( x)LB t ,ν (x) dx = Γ ϕi ( x) µ∈Θ L Bt , ν (ξB t ′ ,µ) LB t ′ ,µ( x) dx = µ∈Θ Γ ϕi ( x)LB t ′ ,µ( x) dx L Bt , ν (ξB t ′ ,µ) = µ∈Θ (Vt ′ )i µ L Bt , ν (ξB t ′ ,µ) and therefore Vt | I t ′ ×Θ = Vt ′ Et ′ (17) with Et ′ ∈ RΘ ×Θ given by (Et ′ )µ ν := L Bt , ν (ξB t ′ ,µ) for all µ ∈ Θ, ν ∈ Θ. Since we also have It = ˙t ′ ∈sonsT I (t) It ′ for all t ∈ TI \ ℒ TI by construction, we do not have to store all the matrices (Vt ) t ∈TI explicitly. It is sufficient to store only the leaf matrices (Vt ) t ∈ℒT I and the transfer matrices (Et ) t ∈TI \{ℛT I } instead. So far, we have not taken advantage of the translation-invariance eq. (2) of the kernel function. It allows us to avoid storing all the coupling matrices (Sb ) b ∈ℒ + TI×J individually: eq. (13) implies ξB t , ν = ξB ( ℓ), ν + mt for all t ∈ T (ℓ) I , ν ∈ Θ, (18a) ξC s , µ = ξC ( ℓ), µ + ms for all s ∈ T (ℓ) J , µ ∈ Θ (18b) on every level ℓ ∈ {0, ..., ℓm ax } while, due to eq. (11), the tree TI ×J satisfies levelT I ×J (t, s) = levelT I (t) = levelT J (s) for all (t, s) ∈ TI ×J . (19) By combining eq. (19) with eq. (18) and eq. (2) we get (S( t,s)) ν µ = g(ξB t ,ν , ξC s , µ) = g(ξB ( ℓ) ,ν + mt , ξC ( ℓ), µ + ms ) = g(ξB ( ℓ) ,ν + mt − mt , ξC ( ℓ), µ + ms − mt ) = g(ξB ( ℓ) ,ν , ξC ( ℓ), µ − (mt − ms ) ) for all (t, s) ∈ ℒ + TI ×J ∩ T (ℓ) I ×J , ν, µ ∈ Θ on every level ℓ ∈ {0, ..., ℓm ax }. Hence, on every level ℓ ∈ {0, ..., ℓm ax } the implication mt − ms = m˜t − m˜s =⇒ S( t,s) = S( ˜t ,˜s) (20) holds for all (t, s), (˜t , ˜s) ∈ ℒ + TI ×J ∩ T (ℓ) I ×J . 8

Moreover, due to our regular construction of the cluster tree, we also do not have to store all the transfer matrices individually: On every level ℓ ∈ {0, ..., ℓm ax − 1} we have, according to eq. (18a), L Bt , ν (x) = L B( ℓ), ν (x − mt ) for all t ∈ T (ℓ) I , ν ∈ Θ, x ∈ Rd , which implies (Et ′ )µ ν = L Bt , ν (ξB t ′ ,µ) = L B( ℓ), ν (ξB t ′ ,µ − mt ) = L B( ℓ), ν (ξB ( ℓ+1) ,µ + mt ′ − mt ) = L B( ℓ), ν (ξB ( ℓ+1) ,µ − (mt − mt ′ )) for all t ∈ T (ℓ) I , t′ ∈ sonsT I (t), ν, µ ∈ Θ. Hence, on every level ℓ ∈ {0, ..., ℓm ax − 1}, we have mt − mt ′ = m˜t − m˜t ′ =⇒ Et ′ = E˜t ′ (21) for all t, ˜t ∈ T (ℓ) I , t′ ∈ sonsT I (t), ˜t′ ∈ sonsT I (˜t ). The column cluster basis (Ws ) s ∈TJ and the corresponding transfer matrices (Fs )s ∈TJ have similar properties. The implications eq. (20) and eq. (21) allow us to significantly reduce the storage requirements of the H2 -matrix representation by taking advantage of the translation- invariance of the kernel function.

3 Complexity estimates

In the following, let n := max{|I|, |J|}, k := |Θ| = (θ + 1)d . (22) For every ι ∈ {1, ..., d} let eι ∈ Rd denote the ι-th canonical unit vector.

Lemma 1 (Storage requirements of the leaf matrices) Storing the leaf matrices

Lemma 1 (Storage requirements of the leaf matrices) Storing the leaf matrices (Vt ) t ∈ℒT I and (Ws ) s ∈ℒT J requires not more than 2kn units of storage. Proof. Combining eq. (22) with eq. (10) implies that storing the matrices (Vt ) t ∈ℒT I and (Ws ) s ∈ℒT J directly requires not more than t∈ℒT I |It | |Θ| + s∈ℒT J |Js | |Θ| = k t∈ℒT I |It | + s∈ℒT J |Js | = k(|I| + |J|) ≤ 2kn units of storage. □ Lemma 2 On every level ℓ ∈ {0, ..., ℓm ax − 1} we have {mt − mt ′ : t ∈ T (ℓ) I , t′ ∈ sonsT I (t)} ⊆ {0, −δ( ℓ+1) ι( ℓ) eι ( ℓ) }, {ms − ms ′ : s ∈ T (ℓ) J , s′ ∈ sonsT J (s)} ⊆ {0, −δ( ℓ+1) ι( ℓ) eι ( ℓ) }. 9

Proof. Let ℓ ∈ {0, ..., ℓm ax − 1}, let t ∈ T (ℓ) I and let t′ ∈ sonsT I (t). Then pt ′ ∈ Nd 0 satisfies pt ′ ,ι = pt ,ι for all ι ∈ {1, ..., d}\{ι( ℓ) } and pt ′ ,ι( ℓ) ∈ {2pt ,ι( ℓ) , 2pt ,ι( ℓ) +1}, according to eq. (8). Due to eq. (4) and eq. (6) we therefore get, as illustrated in Figure 2, (mt − mt ′ )ι = 0 for all ι ∈ {1, ..., d} \ {ι( ℓ) } and (mt − mt ′ )ι ( ℓ) = pt ,ι( ℓ) δ( ℓ) ι( ℓ) − pt ′ ,ι( ℓ) δ( ℓ+1) ι( ℓ) ∈ {pt ,ι( ℓ) δ( ℓ) ι( ℓ) − 2pt ,ι( ℓ) δ( ℓ+1) ι( ℓ) , pt ,ι( ℓ) δ( ℓ) ι( ℓ) − 2pt ,ι( ℓ) δ( ℓ+1) ι( ℓ) − δ( ℓ+1) ι( ℓ) } = {pt ,ι( ℓ) δ( ℓ) ι( ℓ) − pt ,ι( ℓ) δ( ℓ) ι( ℓ) , pt ,ι( ℓ) δ( ℓ) ι( ℓ) − pt ,ι( ℓ) δ( ℓ) ι( ℓ) − δ( ℓ+1) ι( ℓ) } = {0, −δ( ℓ+1) ι( ℓ) }, which proves the first statement. The second statement can be proven similarly. □ r( 5) r( 4) Γ t mt mt 1 mt 2 r( 5) r( 4) Γ m˜t m˜t 1 m˜t 2 ˜t Figure 2: The translations of the boxes of TI and the translations of their sons.

Theorem 1 (Storage requirements of the transfer matrices) The transfer matri-

Theorem 1 (Storage requirements of the transfer matrices) The transfer matri- ces (Et ) t ∈TI \{ℛT I } and (Fs ) s ∈TJ \{ℛT J } require not more than 4k2 ℓm ax units of storage. Proof. Combining the implication eq. (21) with lemma 2 and eq. (22) directly implies that for representing (Et ) t ∈TI \{ℛT I } and (Fs ) s ∈TJ \{ℛT J } not more than |Θ|2 ℓm ax− 1 ℓ=0 {mt − mt ′ : t ∈ T (ℓ) I , t′ ∈ sonsT I (t)} + |Θ|2 ℓm ax− 1 ℓ=0 {ms − ms ′ : s ∈ T (ℓ) J , s′ ∈ sonsT J (s)} ≤ |Θ|2 ℓm ax− 1 ℓ=0 4 = 4k2 ℓm ax units of storage are needed. □ 10

In order to obtain further results, we now consider not only the boxes appearing as nodes in our trees TI and TJ , but the infinitely many boxes belonging to T := ℓ m ax ℓ =0 T ( ℓ) given by (cf. eq. (4)) T ( ℓ) := {r( ℓ) + δ( ℓ) ⊙ p : p ∈ Zd } for all ℓ ∈ {0, ..., ℓm ax }. (23) For all ℓ ∈ {0, ..., ℓm ax } let (pt ) t∈T ( ℓ) be defined by pt ∈ Zd , t = r( ℓ) + δ( ℓ) ⊙ pt =:mt = r( ℓ) + mt for all t ∈ T ( ℓ) (24) and let (cf. eq. (13)) Bt := B (ℓ) + δ( ℓ) ⊙ pt = B (ℓ) + mt for all t ∈ T ( ℓ) , (25a) Cs := C (ℓ) + δ( ℓ) ⊙ ps = C (ℓ) + ms for all s ∈ T ( ℓ) . (25b) Based on (14), let Aη : T × T → {true, false}, ( t, s) → t rue if max{diam(Bt ) , diam(Cs ) } ≤ η dist(Bt , Cs ) f alse if max{diam(Bt ) , diam(Cs ) } > η dist(Bt , Cs ) . (26) We base the following analysis on one key assumption.

Assumption 1 (Diameters of supports) We assume that there is a constant C bb ∈

Assumption 1 (Diameters of supports) We assume that there is a constant C bb ∈ R ≥1 satisfying diam(B (ℓ) ) ≤ C bb diam(r( ℓ) ) and diam(C (ℓ) ) ≤ C bb diam(r( ℓ) ) (27) for all ℓ ∈ {0, ..., ℓm ax }. We will see later on that under normal circumstances, due to our choice of ℓm ax, this condition is actually fulfilled since we do not allow the boxes to become too small (see Remark 1). Let λd denote the d-dimensional Lebesgue measure.

Lemma 3 There is a constant C dv ∈ R >1 such that

diam(r( ℓ) )d ≤ C dv λd ( r( ℓ) ) for all ℓ ∈ {0, ..., ℓm ax }. (28) Proof. According to eq. (5) and eq. (6), we have diam(r( ℓ) ) = 12 diam(r( ℓ−d) ) and λd ( r( ℓ) ) = 1 2d λd ( r( ℓ−d) ) for all ℓ ∈ {d, ..., ℓm ax } and therefore diam(r( ℓ) )d λd ( r( ℓ) ) = diam(r( ℓ−d) )d λd ( r( ℓ−d) ) for all ℓ ∈ {d, ..., ℓm ax }. 11

This implies diam(r( ℓ) )d λd ( r( ℓ) ) : ℓ ∈ {0, ..., ℓm ax } = diam(r( ℓ) )d λd ( r( ℓ) ) : ℓ ∈ {0, ..., d − 1} =: R. Hence, C dv := max R satisfies eq. (28). □

Lemma 4 (Sparsity) We define

C sp := 2− d 3 + 2η− 1 d Cd bbω d C dv , where ωd := λd ( {z ∈ Rd : ∥z∥2 ≤ 1}) is the volume of the d-dimensional unit ball. For every ℓ ∈ {0, ..., ℓm ax } we have {s ∈ T ( ℓ) : Aη (t, s) = false} ≤ C sp for all t ∈ T ( ℓ) , {t ∈ T ( ℓ) : Aη (t, s) = false} ≤ C sp for all s ∈ T ( ℓ) . Proof. Let ℓ ∈ {0, ..., ℓm ax } , let t ∈ T ( ℓ) , let Ht := {s ∈ T ( ℓ) : Aη (t, s) = fa lse} and let βt denote the midpoint of Bt . For every s ∈ Ht we have according to eq. (26) dist(Bt , Cs ) < η− 1 max{diam(Bt ) , diam(Cs ) }. Thus for every s ∈ Ht there is at least one pair (xs , ys ) ∈ Bt × Cs satisfying ∥xs − ys ∥ 2 < η− 1 max{diam(Bt ) , diam(Cs ) }, which, according to (25), implies ∥z − βt ∥2 ≤ ∥z − ys ∥ 2 + ∥ys − xs ∥ 2 + ∥xs − βt ∥ 2 < diam(Cs ) + η− 1 max{diam(Bt ) , diam(Cs ) } + 12 diam(Bt ) ≤ 1 2 (3 + 2η− 1 ) max{diam(Bt ) , diam(Cs ) } = 1 2 (3 + 2η− 1 ) max{diam(B (ℓ) ), diam(C (ℓ) )} =: ρ( ℓ) for all z ∈ Cs , which, combined with eq. (12b), eq. (24) and eq. (25b), implies (cf. Figure 3) s ⊆ Cs ⊂ B 2[ βt , ρ( ℓ) ] := {z ∈ Rd : ∥z − βt ∥ 2 ≤ ρ( ℓ) }. (29) Moreover, according to eq. (23), we have (cf. Figure 3) λd z ∈ B 2[ βt , ρ( ℓ) ] : |{s ∈ Ht : z ∈ s}| > 1 = 0. (30) By using eq. (28), eq. (23), eq. (29), eq. (30) and eq. (27) we now get diam(r( ℓ) )d |Ht | ≤ C dv λd ( r( ℓ) )|Ht | = C dv s∈Ht λd ( s) = C dv s∈Ht B 2 [ βt , ρ( ℓ)] 1 s ( z) dz = C dv B 2 [ βt , ρ( ℓ)] s∈Ht 1 s ( z) dz ≤ C dv B 2 [ βt , ρ( ℓ)] 1 B 2[ βt , ρ( ℓ)] ( z) dz 12

t Bt ·βt ρ( 6) B 2[ βt , ρ( 6) ] Figure 3: For t ∈ T ( 6) the ball B 2[ βt , ρ( 6) ] contains all inadmissible s ∈ T ( 6) . = C dv λd ( B2 [ βt , ρ( ℓ) ]) = C dv ωd ( ρ( ℓ) )d = C dv ωd 2 − d (3 + 2η− 1 )d max{diam(B (ℓ) ), diam(C (ℓ) )}d ≤ C dv ωd 2 − d (3 + 2η− 1 )d Cd bb diam(r( ℓ) )d , which proves the first statement. The second statement can be proven similarly. □

Lemma 5 (Difference bound for T I×J ) On every level ℓ ∈ {0, ..., ℓ max } we have

{mt − ms : (t, s) ∈ T (ℓ) I ×J } ≤ 3Cs p. Proof. By definition, we have T (0) I ×J = {(ℛT I , ℛ TJ )} = {(r( 0) , r( 0) )} and therefore {mt − ms : (t, s) ∈ T (0) I ×J } = {mr ( 0) − mr ( 0) } = {0}. Let ℓ ∈ {1, ..., ℓm ax } and let (t, s) ∈ T (ℓ) I ×J . Then, by the definition of TI ×J and due to eq. (19), there is a pair (ˇt , ˇs) ∈ T (ℓ−1) I ×J ⊆ T (ℓ−1) I × T (ℓ−1) J satisfying t ∈ sonsT I (ˇt ), s ∈ sonsT J (ˇs) and Aη (ˇt , ˇs) = false. According to eq. (13) and eq. (14), we therefore get max{diam(B (ℓ−1) ), diam(C (ℓ−1) )} = max{diam(Bˇt ) , diam(Cˇs ) } > η dist(Bˇt , Cˇs ) = η inf{∥x + mˇt − (y + mˇs ) ∥2 : x ∈ B (ℓ−1) , y ∈ C (ℓ−1) } = η inf x − y − (mˇt − mˇs ) 2 : x ∈ B (ℓ−1) , y ∈ C (ℓ−1) = η dist B (ℓ−1) , C (ℓ−1) − (mˇt − mˇs ) , which, in view of eq. (23), eq. (24), eq. (25) and eq. (26), implies r( ℓ−1) −(mˇt − mˇs ) = r( ℓ−1) +δ( ℓ−1) ⊙(pˇs − pˇt ) ∈ u ∈ T ( ℓ−1) : Aη (r( ℓ−1) , u) = false =:H (ℓ− 1) . (31) Moreover, lemma 2 implies mt − ms ∈ {mˇt − mˇs , mˇt − mˇs + δ( ℓ) ι( ℓ−1) eι ( ℓ−1) , mˇt − mˇs − δ( ℓ) ι( ℓ−1) eι ( ℓ−1) }, 13

which, combined with eq. (31) and eq. (24), leads to mt − ms ∈ u∈H (ℓ− 1) {−mu , −mu + δ( ℓ) ι( ℓ−1) eι ( ℓ−1) , −mu − δ( ℓ) ι( ℓ−1) eι ( ℓ−1) }. Since we have chosen (t, s) ∈ T (ℓ) I ×J arbitrarily, this leads to |{mt − ms : (t, s) ∈ T (ℓ) I ×J }| ≤ 3|H (ℓ−1) |. Furthermore, due to lemma 4, we have |H (ℓ−1) | ≤ C sp. □

Theorem 2 (Storage requirements of the coupling matrices) The coupling ma-

Theorem 2 (Storage requirements of the coupling matrices) The coupling ma- trices (Sb ) b ∈ℒ + TI×J require not more than 3Cs pk 2 (ℓm ax − d) units of storage. Proof. On every level ℓ ∈ {0, ..., d} all boxes t ∈ T (ℓ) I and s ∈ T (ℓ) J contain, by con- struction, the midpoint of r( 0) , which implies dist(Bt , Cs ) = 0 for all (t, s) ∈ T (ℓ) I ×J and therefore ℒ + TI ×J ∩ T (ℓ) I ×J = {(t, s) ∈ T (ℓ) I ×J : Aη (t, s) = true} = ∅. Combining this with the implication eq. (20), lemma 5 and eq. (22) directly implies that for representing all coupling matrices (Sb ) b ∈ℒ + TI×J not more than |Θ|2 ℓm ax ℓ=d+1 {mt − ms : (t, s) ∈ ℒ + TI ×J ∩ T (ℓ) I ×J } ≤ k2 ℓm ax ℓ=d+1 3Cs p = 3Cs pk 2 (ℓm ax − d) units of storage are needed. □ In order for arithmetic operations like the matrix-vector multiplication to have the desired complexity, TI ×J has to be sparse. This property is guaranteed by our construc- tion.

Lemma 6 (Sparsity of T I×J ) We have

|{s ∈ TJ : (t, s) ∈ TI ×J }| ≤ 2Cs p for all t ∈ TI , |{t ∈ TI : (t, s) ∈ TI ×J }| ≤ 2Cs p for all s ∈ TJ . Proof. The statements are implied by the definition of TI ×J , eq. (19), lemma 4 and the fact that every node of TI or TJ has at most two sons. □ For the sake of completeness, we also have to consider the storage requirements of the

nearfield matrices

N( t,s) := G|I t × Js for all (t, s) ∈ ℒ − TI ×J that capture the part of the matrix that cannot be compressed. Since we want to maintain the ability to work with unstructured surface meshes, we cannot avoid storing all of these matrices explicitly. 14

Lemma 7 (Storage requirements of the nearfield matrices) We define the reso-

Lemma 7 (Storage requirements of the nearfield matrices) We define t lutions of the cluster trees TI and TJ by γI := max{|It | : t ∈ ℒ TI }, γJ := max{|Js | : s ∈ ℒ TJ }. Storing (Nb ) b ∈ℒ − TI×J requires not more than C sp min{γI , γJ }n units of storage.

Proof. By construction and due to (11), (4) and (23) we have ℒ − TI ×J ⊆ ℒ TI × ℒ TJ = T (ℓm ax) I × T (ℓm ax) J ⊆ T ( ℓm ax) × T ( ℓm ax) . Lemma 4 and (10) therefore imply that storing all nearfield matrices (Nb ) b ∈ℒ − TI×J requires not more than (t,s)∈ℒ − T I ×J |It | |Js | = t∈ℒT I s∈ℒT J (t,s)∈ℒ − T I ×J |It | |Js | ≤ t∈ℒT I s∈T ( ℓma x) Aη (t,s)=false γJ |It | ≤ t∈ℒT I C spγ J |It | = C spγ J t∈ℒT I |It | = C spγ J |I| units of storage and similarly not more than C spγ I |J| units of storage, which, in view of (22), implies the statement. □ Theorem 1, lemma 7 and theorem 2 suggest that we still have to discuss how to properly choose the maximal level ℓm ax. On the one hand, we have to choose ℓm ax large enough to keep at least one of the resolutions γI and γJ sufficiently small. On the other hand, since the complexity of the matrix-vector multiplication mainly depends on |TI | and |TJ | (cf. [6, Theorem 3.42]), we have to choose ℓm ax small enough to prevent |TI | and |TJ | from getting too large. In order to be able to analyze γI and γJ , we assume that the relative number of characteristic points assigned to a box (cf. (7) and (9)) can essentially be estimated by its size, i.e., we assume that there is a constant C nd ∈ R >0 satisfying |It | ≤ C nd diam(t)d −1 |I|, |Js | ≤ C nd diam(s)d −1 |J| for all t ∈ TI , s ∈ TJ . (32) This condition is, for example, usually fulfilled on shape-regular and quasi-uniform meshes (cf. [30, Remark 4.1.14]). Due to (5) and (6), we have diam(r( ℓ+d) ) = 12 diam(r( ℓ) ) for all ℓ ∈ {0, ..., ℓm ax − d}. (33) For the sake of simplicity, we therefore choose ℓm ax ∈ d · N. According to (32), (4) and (33), we thus get |It | ≤ C nd diam(t)d −1 |I| = C nd diam(r( ℓm ax ) )d −1 |I| 15

= C nd2 − d−1 d ℓm ax diam(r( 0) )d −1 |I| for all t ∈ T (ℓm ax) I = ℒ TI and likewise |Js | ≤ C nd2 − d−1 d ℓm ax diam(r( 0) )d −1 |J| for all s ∈ T (ℓm ax) J = ℒ TJ . Hence, with C rs := C nd diam(r( 0) )d −1 we have γI ≤ C rs2 − d−1 d ℓm ax |I|, γJ ≤ C rs2 − d−1 d ℓm ax |J|. (34) We choose a constant C rk ∈ R >0 satisfying C rk k ≤ min{|I|, |J|} and observe 2− d−1 d ℓm ax |I| ≤ C rk k ⇐⇒ ℓm ax ≥ d d−1 log2 ( |I|) − log2 ( Cr k k) , 2− d−1 d ℓm ax |J| ≤ C rk k ⇐⇒ ℓm ax ≥ d d−1 log2 ( |J|) − log2 ( Cr k k) . Based on this, we choose ℓm ax := min ℓ ∈ d · N : ℓ ≥ d d−1 min{log2 ( |I|), log2 ( |J|)} − log2 ( Cr k k) (35) to gain, according t o eq. (34), min{γI , γJ } ≤ C rsC rk k. (36) In order to be able to analyze |TI | and |TJ |, we assume that the number of boxes per level which contain only a relatively small part of the boundary Γ does not grow too fast, i.e., we assume that there are constants C da, C sa ∈ R >0 satisfying |{t ∈ T (ℓ) I : diam(t)d −1 > C da| t ∩ Γ|}| ≤ C sa diam(r( ℓ) )− (d−1) (37a) |{s ∈ T (ℓ) J : diam(s)d −1 > C da| s ∩ Γ|}| ≤ C sa diam(r( ℓ) )− (d−1) (37b) on every level ℓ ∈ {0, ..., ℓm ax }, where |X| := Γ 1 X (y) dy for all measurable subsets X ⊆ Γ. Lemma 8 On every level ℓ ∈ {0, ..., ℓm ax } we have |{t ∈ T (ℓ) I : diam(t)d −1 ≤ C da| t ∩ Γ|}| ≤ 2Cd a| Γ| diam(r( ℓ) )− (d−1) , |{s ∈ T (ℓ) J : diam(s)d −1 ≤ C da| s ∩ Γ|}| ≤ 2Cd a| Γ| diam(r( ℓ) )− (d−1) . Proof. Let ℓ ∈ {0, ..., ℓm ax } and let T (ℓ) I ,Cd a := {t ∈ T (ℓ) I : diam(t)d −1 ≤ C da| t ∩ Γ|}. Due to (4), the surface area of the set x ∈ Γ : |{t ∈ T (ℓ) I ,Cd a : x ∈ t}| > 2 equals 0. Thus diam(r( ℓ) )d −1 T (ℓ) I ,Cd a = t∈T (ℓ)I,Cd a diam(t)d −1 ≤ t∈T (ℓ)I,Cd a C da| t ∩ Γ| 16

= C da t∈T (ℓ)I,Cd a Γ 1 t ∩Γ( y) dy = C da Γ t∈T (ℓ)I,Cd a 1 t ∩Γ (y) dy ≤ C da Γ 2 · 1 Γ ( y) dy = 2Cd a |Γ| holds, which proves the first statement. The second statement can be proven similarly. □ Defining C nb := 2Cd a |Γ| + C sa and combining assumption (37) with Lemma 8 yield |T (ℓ) I | ≤ C nb diam(r( ℓ) )− (d−1) for all ℓ ∈ {0, ..., ℓm ax }, (38a) |T (ℓ) J | ≤ C nb diam(r( ℓ) )− (d−1) for all ℓ ∈ {0, ..., ℓm ax }. (38b) Moreover, according to (35), we also have ℓm ax ≤ d d−1 min{log2 ( |I|), log2 ( |J|)} − log2 ( Cr k k) + d. (39) By using (38), (33) and (39) we get |TI | = ℓm ax ℓ=0 |T (ℓ) I | ≤ ℓm ax ℓ=0 C nb diam(r( ℓ) )− (d−1) = C nb ℓm ax ℓ=0 diam(r( ℓ) )− (d−1) = C nb diam(r( 0) )− (d−1) + ℓmax d− 1 ℓ=0 d u=1 diam(r( dℓ+u) )− (d−1) < C nb diam(r( 0) )− (d−1) + ℓmax d− 1 ℓ=0 d diam(r( d(ℓ+1)) )− (d−1) = C nb diam(r( 0) )− (d−1) + d ℓmax d− 1 ℓ=0 2( ℓ+1)(d−1) diam(r( 0) )− (d−1) < C nb d diam(r( 0) )− (d−1) + d diam(r( 0) )− (d−1) ℓmax d− 1 ℓ=0 2( ℓ+1)(d−1) = C nbd diam(r( 0) )− (d−1) 1 + ℓmax d ℓ=1 2d −1 ℓ = C nbd diam(r( 0) )− (d−1) ℓmax d ℓ=0 2d −1 ℓ = C nbd diam(r( 0) )− (d−1) 1 − (2d −1 ) ℓmax d +1 1 − 2d −1 < C nbd diam(r( 0) )− (d−1) (2d −1 − 1)− 1 2 d−1 d ℓm ax+ d−1 ≤ C nbd diam(r( 0) )− (d−1) (2d −1 − 1)− 1 2 d−1 d ( d d−1 (log2 ( |I|)−log2 ( Cr k k))+d)+d−1 = C nbd diam(r( 0) )− (d−1) (2d −1 − 1)− 1 2l og2 ( |I|)−log2 ( Cr k k)+2(d−1) 17

= C nbd diam(r( 0) )− (d−1) (2d −1 − 1)− 1 22 (d−1) C− 1 rk k− 1 |I|. Due to x2 x−1 = x2 (x+2) (x−1)(x+2) = x2 (x+2) x2 + x−2 ≤ x + 2 for x ∈ R ≥2 we therefore have |TI | < C nbd diam(r( 0) )− (d−1) (2d −1 + 2)C− 1 r k k− 1 |I| and analogously |TJ | < C nbd diam(r( 0) )− (d−1) (2d −1 + 2)C− 1 r k k− 1 |J|. Finally, based on our choice of the maximal level ℓm ax , we can now specify the storage requirements for our H2 -matrix approximation: Theorem 3 (Total storage requirements) G requires not more than • 2kn units of storage for the leaf matrices, • 4k2 d d−1 log2 ( n) − log2 ( Cr k k) + d units of storage for the transfer matrices, • C spC rsC rk kn units of storage for the nearfield matrices, • 3Cs pk 2 d d−1 log2 ( n) − log2 ( Cr k k) units of storage for the coupling matrices. Proof. The statement is a direct consequence of Lemma 1, Theorem 1, Lemma 7, The- orem 2, (36) and (39). □ Compared to the standard H2 -approach, these storage complexity bounds lead to a significant reduction of the total storage requirements, as we will see and discuss in chapter 4.

Remark 1 (Bounded enlargement) Due to (33) and (39), we have

diam(r( ℓ) ) ≥ diam(r( ℓm ax) ) ∼ 2− ℓmax d ≳ 2− log2( n) d−1 = n− 1 d−1 for all ℓ ∈ {0, ..., ℓm ax }. n many situations, the diameters of the supports behave like n− 1 d−1 , for example on

shape-regular and quasi-uniform meshes. In these cases, the condition (27) is therefore fulfilled.

4 Numerical experiments

In the following, let g : Rd × Rd → R, (x, y) → 1 4π 1 ∥x−y∥2 if x ̸= y 0 if x = y and for all y ∈ Γ = ∂Ω let ny denote the corresponding outward-pointing unit normal vector. Under suitable conditions, a solution u of the interior Dirichlet boundary value

problem

−∆u = 0 in Ω, u = gD on Γ (40) 18

is given by Green’s representation formula (cf. [23, Theorem 2.2.2], [30, Theorem 3.1.6]), which reduces obtaining u to obtaining the corresponding Neumann boundary values gN : Γ → R, y → ∂u ∂ny (y). Depending on the situation, gN as well as the given Dirichlet boundary values gD might only exist in a generalized sense. The Neumann values gN can be obtained by solving a boundary integral equation resulting from combining Green’s representation formula with suitable trace operators (cf. [30, Section 3.4.2.1]). Applying a Galerkin discretization to this integral equation on a triangulation of Γ with discontinuous piecewise constant basis functions (ϕi ) i ∈I and replacing gD with the corresponding L2 -orthogonal projection ˜gD into a space spanned by continuous piecewise linear nodal basis functions (ψj )j ∈J leads to a linear system Gx = K + 1 2 M b , (41) where b ∈ RJ contains the coefficients of ˜gD = j ∈J bj ψj , x ∈ RI contains the coefficients of the approximate Neumann values ˜gN := i ∈ I xi ϕ i , and where G ∈ RI ×I , K ∈ RI ×J and M ∈ RI ×J are given by Gi j := Γ ϕi ( x) Γ g(x, y)ϕj (y) dy dx for all i ∈ I, j ∈ I, Ki j := Γ ϕi ( x) Γ ∂g ∂ny (x, y) ψj (y) dy dx for all i ∈ I, j ∈ J, Mi j := Γ ϕi ( x)ψj (x) dx for all i ∈ I, j ∈ J. The derivative ∂ ∂ny is applied with respect to the y-variable. We construct an H2 -matrix approximation G of G as described in section 2. In order to approximate K, we replace ∂ ∂ny g(x, y) with ∂ ∂ny ( IB t ⊗IC s )[g](x, y) for all (x, y) ∈ Bt ×Cs , (t, s) ∈ ℒ + TI ×J (cf. eq. (15)). This approach leads to an H2 -matrix approximation K of K which differs from the H2 -matrix described in section 2 only in terms of the leaf matrices for the column cluster basis, which are given by Ws ∈ RJ s × Θ , (Ws ) j µ := Γ ψj (y) ∂ ∂ny L Cs , µ( y) dy for all s ∈ ℒ TJ , j ∈ Js , µ ∈ Θ. Hence, all results from section 3 also hold for K. Moreover, we can expect a similar convergence behavior as for G (cf. [6, Chapter 4]). The sparse matrix 1 2 M can simply be added to the nearfield entries o f K and therefore does not require any additional storage. In this way, we approximate the Neumann values of the function u0 := 4πg(·, y0 ) with y0 := (1.2, 1.2, 1.2)T on an aproximation ΓS of the unit sphere S 2 := {x ∈ R3 : ∥x∥2 = 1} = ∂{x ∈ R3 : ∥x∥2 < 1} consisting of plane triangles. S 2 and its approximation are both contained in the box [−1, 1]3 , which we use as our starting box r( 0) . We choose the centers of gravity of the triangles as characteristic points for (ϕi ) i ∈I , the vertices of the triangles as characteristic points for (ψj )j ∈J , C rk = 2 and η = 2. We denote the 19

maximal level of the trees used for the construction of G by ℓI ×I m ax and the maximal level of the trees used for the construction of K by ℓI ×J m ax . W e compute the matrix entries by tensor Gauss quadrature as described in [16, 29],[30, Chapter 5] with 3 quadrature points per dimension for the regular integrals and 5 quadrature points per dimension for the singular integrals. In order to preserve the convergence behavior of the standard Galerkin solution, we increase the interpolation degree θ by 1 whenever the mesh width h is halved, i.e., whenever the number n = max{|I|, |J|} = |I| of triangles is quadrupled (cf. [6, Section 10.1 and Section 10.2]). We solve the linear system using the cg method with a relative residual accuracy of 10− 6 (cf. [30, Section 6.1]). n θ k ℓI ×I m ax ℓI ×J m ax ǫL 2 8192 4 125 9 9 1.194− 2 32768 5 216 12 9 5.683− 3 131072 6 343 12 12 2.870− 3 524288 7 512 15 15 1.412− 3 2097152 8 729 18 15 7.050− 4 Table 1: Parameters and resulting L2 -error ǫL 2 := ∥gN − ˜gN ∥L 2 ( ΓS ). n leaf matrices transfer matrices nearfield matrices coupling matrices 8192 15.6 4.3 74.4 51.8 32768 108.1 17.1 303.8 218.6 131072 686.1 43.1 3913.3 509.9 524288 4096.3 120.0 16090.5 1496.1 2097152 23329.2 291.9 65051.2 3762.8 Table 2: Storage requirements of G on ΓS in MB. n leaf matrices transfer matrices nearfield matrices coupling matrices 8192 11.7 4.3 38.8 54.2 32768 81.0 12.8 476.1 143.9 131072 514.6 43.1 1956.3 524.3 524288 3072.3 120.0 8047.3 1528.1 2097152 17496.3 243.3 128696.6 2765.3 Table 3: Storage requirements of K on ΓS in MB. It is clearly visible in the Tables 1, 2 and 3 that the storage requirements of our H2 -matrices G and K behave as predicted in Chapter 3: The storage for the leaf ma- trices grows like O(nk) (cf. Lemma 1), the storage for the transfer matrices grows like O(ℓm ax k2 ) (cf. Theorem 1), the storage for the nearfield matrices grows like O(2− 2 3 ℓm ax n2 ) (cf. Lemma 7, (34)) and the storage for the coupling matrices grows like O((ℓm ax − 3)k2 ) (cf. Theorem 2). 20

21 3 21 5 21 7 21 9 22 1 22 3 23 25 27 29 21 1 21 3 21 5 21 7 21 9 n S t o r a g e i n M B O(nk) nearfield matrices leaf matrices 21 3 21 5 21 7 21 9 22 1 22 3 22 24 26 28 21 0 21 2 21 4 n S t o r a g e i n M B O(log2 ( n)k2 ) coupling matrices transfer matrices Figure 4: Storage requirements of G on ΓS in MB. 21 3 21 5 21 7 21 9 22 1 22 3 23 25 27 29 21 1 21 3 21 5 21 7 21 9 n S t o r a g e i n M B O(nk) nearfield matrices leaf matrices 21 3 21 5 21 7 21 9 22 1 22 3 22 24 26 28 21 0 21 2 n S t o r a g e i n M B O(log2 ( n)k2 ) coupling matrices transfer matrices Figure 5: Storage requirements of K on ΓS in MB. Figure 4 and Figure 5 show that the storage requirements of the matrix components can be bounded as predicted by Theorem 3. The expected theoretical convergence rate of O(h) for the Neumann values in the L2 -norm (cf. [30, Theorem 4.1.33] and [15, Theorem 4.6]) is preserved (see Table 1). We want to compare these results with the corresponding results obtained by conven- tional H2 -matrix approximations Gc onv and Kc onv , which are based on different cluster trees and support bounding boxes ( cf. [6, Section 3.3]): The construction is focused on reducing the diameters of the boxes as fast as p ossible instead of keeping them uniform on every level. Hence, one has to store an individual transfer matrix for almost every 21

box and an individual coupling matrix for every admissible leaf of the block cluster tree. Since approximating a block with less than k2 entries would therefore be inefficient, the recursion of the algorithm for the cluster tree construction stops when it reaches a box with not more than 2k characteristic points. On the finest mesh (n = 2097152), Gc onv requires 23328.1 MB for the leaf matrices, 33718.1 MB for the transfer matrices, 241218.8 MB for the nearfield matrices and 206819.6 MB for the coupling matrices. Kc onv requires 17496.1 MB for the leaf matrices, 25187.2 MB for the transfer matrices, 2 00100.7 MB for the nearfield matrices and 126874.0 MB for the coupling matrices. The resulting L2 -error is 7.022− 4 . In order to test our translation-invariant compression method on a more challenging boundary containing edges and corners, we approximate the Neumann values of the function u1 : R3 → R, x → x21 −x23 on the boundary ΓC := ∂([−1, 1]×[− 34 , 3 4 ]×[− 12 , 12 ]) of a rectangular cuboid by the same method as before, i.e., we replace the matrices G and K appearing in (41) with their respective translation-invariant H2 -matrix approximations G and K and then solve the resulting perturbed linear system with the cg method. This time, we use r( 0) = [−1, 1] × [− 3 4 , 3 4 ] × [− 12 , 12 ] as the starting box, C rk = 1 for the determination of ℓI ×J m ax and increasing quadrature orders. The number of quadrature points per dimension for the regular integrals is denoted by qr eg and the number of quadrature points per dimension for the singular integrals is denoted by qs ing . The characteristic points, η, ℓI ×I m ax and the accuracy of the cg method are chosen as before. n qr eg qs ing θ k ℓI ×I m ax ℓI ×J m ax ǫL 2 8112 4 6 4 125 9 9 3.895− 2 32448 4 7 5 216 12 12 1.484− 2 129792 4 8 6 343 12 12 7.408− 3 519168 5 9 7 512 15 15 2.626− 3 2076672 5 10 8 729 18 18 1.332− 3 Table 4: Parameters and resulting L2 -error ǫL 2 := ∥gN − ˜gN ∥L 2 ( ΓC ). n leaf matrices transfer matrices nearfield matrices coupling matrices 8112 15.5 4.3 55.9 48.0 32448 107.0 17.1 192.4 215.8 129792 679.4 43.1 1906.6 484.8 519168 4056.4 120.0 6915.7 1488.1 2076672 23101.6 291.9 26143.3 3843.8 Table 5: Storage requirements of G on ΓC in MB. As in the previous experiment, the storage requirements of the translation-invariant H2 -matrix approximations and the resulting L2 -error of the Neumann values behave as predicted by the theory (cf. Tables 4, 5, 6 and Figures 6, 7). 22

n leaf matrices transfer matrices nearfield matrices coupling matrices 8112 11.6 4.3 31.0 59.0 32448 80.3 17.1 106.7 269.2 129792 509.6 43.1 1237.3 545.8 519168 3042.4 120.0 4480.0 1624.1 2076672 17326.6 291.9 16913.0 4119.6 Table 6: Storage requirements of K on ΓC in MB. 21 3 21 5 21 7 21 9 22 1 23 25 27 29 21 1 21 3 21 5 21 7 n S t o r a g e i n M B O(nk) nearfield matrices leaf matrices 21 3 21 5 21 7 21 9 22 1 22 24 26 28 21 0 21 2 n S t o r a g e i n M B O(log2 ( n)k2 ) coupling matrices transfer matrices Figure 6: Storage requirements of G on ΓC in MB. 21 3 21 5 21 7 21 9 22 1 23 25 27 29 21 1 21 3 21 5 21 7 n S t o r a g e i n M B O(nk) nearfield matrices leaf matrices 21 3 21 5 21 7 21 9 22 1 22 24 26 28 21 0 21 2 n S t o r a g e i n M B O(log2 ( n)k2 ) coupling matrices transfer matrices Figure 7: Storage requirements of K on ΓC in MB. 23

We again want to compare the results of the translation-invariant H2 -matrix compres- sion with the corresponding results of the conventional H2 -matrix compression on the finest mesh (n = 2076672): Gc onv requires 23100.3 MB for the leaf matrices, 26630.6 MB for the transfer matrices, 2 47629.2 MB for the nearfield matrices and 167992.2 MB for the coupling matrices. Kc onv requires 17325.2 MB for the leaf matrices, 20321.7 MB for the transfer matrices, 193189.7 MB for the nearfield matrices and 115508.9 MB for the coupling matrices. The r esulting L2 -error is 1.503− 3 . As the complexity bounds suggest and the results demonstrate, the translation-invari- ant approach reduces the storage requirements for the transfer and coupling matrices drastically. Furthermore, the translation-invariant H2 -matrix approximations require remarkably less storage for the nearfield matrices than their respective conventional counterparts. This is due to the fact that the leaves of the cluster trees used in the translation-invariant approach contain significantly less characteristic points than the leaves of the conventional cluster trees. The resulting increase of the farfield can be handled very efficiently, since no individual but only very few coupling matrices have to be stored. We conclude that by properly exploiting the translation-invariance property (2) of the kernel function the storage requirements of H2 -matrix approximations can be greatly reduced without affecting the desired accuracy.

References

[1] C. R. Anderson. An implementation of the fast multipole method without multi- poles. SIAM J. Sci. Stat. Comp., 13:923–947, 1992. [2] L. Banjai and W. Hackbusch. H- and H2 -matrices for low and high frequency Helmholtz equations. IMA J. Numer. Anal., 28:46–79, 2008. [3] M. Bebendorf. Approximation of boundary element matrices. Numer. Math., 86(4):565–589, 2000.

[4] M. Bebendorf and R. Venn. Constructing nested bases approximations from the
entries of non-local operators. Numer. Math., 121(4):609–635, 2012.
[5] T. Betcke and M. W. Scroggs. Bempp-cl: A fast Python based just-in-time com-
piling boundary element library. J. Open Source Software, 6(59), 2021. available at
https://doi.org/10.21105/joss.02879.

[6] S. Börm. Efficient Numerical Methods for Non-local Operators: H 2 -Matrix Com-

pression, Algorithms and Analysis, volume 14 of EMS Tracts in Mathematics. EMS,

pression, Algorithms and Analysis, volume 14 of EMS Tracts in Mathematics. EMS,
2010.
[7] S. B¨orm and C. B¨orst. Hybrid matrix compression for high-frequency problems.
SIAM J. Matrix Anal. Appl., 41(4):1704–1725, 2020.
24

[8] S. B¨orm and S. Christophersen. Approximation of integral operators by Green quadrature and nested cross approximation. Numer. Math., 133(3):409–442, 2016. [9] S. B¨orm and L. Grasedyck. Hybrid cross approximation of integral operators. Nu- mer. Math., 101:221–249, 2005. [10] S. B¨orm and W. Hackbusch. H2 -matrix approximation of integral operators by interpolation. Appl. Numer. Math., 43:129–143, 2002. [11] S. B¨orm and W. Hackbusch. Approximation of boundary element operators by adaptive H2 -matrices. Found. Comp. Math., 312:58–75, 2004. [12] S. B¨orm, M. L¨ohndorf, and J. M. Melenk. Approximation of integral operators by variable-order interpolation. Numer. Math., 99(4):605–643, 2005. [13] S. B¨orm and J. M. Melenk. Approximation of the high-frequency Helmholtz kernel by nested directional interpolation: error analysis. Numer. Math., 137(1):1–34, 2017. [14] E. Corona, A. Rahimian, and D. Zorin. A tensor-train accelerated solver for integral equations in complex geometries. J. Comp. Phys, 334:145–169, 2015. [15] W. Dahmen, B. Faermann, I. G. Graham, W. Hackbusch, and S. A. Sauter. Inverse inequalities on non-quasiuniform meshes and applications to the mortar element method. Math. Comp., 73:1107–1138, 2004. [16] S. Erichsen and S. A. Sauter. Efficient automatic quadrature in 3-d Galerkin BEM. Comput. Meth. Appl. Mech. Eng., 157:215–224, 1998. [17] W. Fong and E. Darve. The black-box fast multipole method. J. Comp. Phys., 228:8712–8725, 2009. [18] K. Giebermann. Multilevel approximation of boundary integral operators. Com- puting, 67:183–207, 2001. [19] Z. Gimbutas and V. Rokhlin. A generalized fast multipole method for nonoscillatory kernels. SIAM J. Sci. Comput., 24(3):796–817, 2002. [20] L. Grasedyck. Adaptive recompression of H-matrices for BEM. Computing, 74(3):205–223, 2004.

[21] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. J. Comp.
Phys., 73:325–348, 1987.
[22] L. Greengard and V. Rokhlin. A new version of the fast multipole method for
the Laplace equation in three dimensions. In Acta Numerica 1997, pages 229–269.
Cambridge University Press, 1997.

[23] W. Hackbusch. Elliptic Differential Equations. Theory and Numerical Treatment.

W. Hackbusch. Elliptic DifferentialEqua
Springer-Verlag Berlin, 1992.
25
[24] W. Hackbusch. Integral equations. Birkh¨auser Basel, 1995.
[25] W. Hackbusch. Hierarchical Matrices: Algorithms and Analysis. Springer, 2015.
[26] W. Hackbusch and Z. P. Nowak. On the fast matrix multiplication in the boundary
element method by panel clustering. Numer. Math., 54(4):463–491, 1989.
[27] G. C. Hsiao and W. L. Wendland. Boundary Integral Equations. Number 164 in
Appl. Math. Sci. Springer, 2008.
[28] V. Rokhlin. Rapid solution of integral equations of classical potential theory. J.
Comp. Phys., 60:187–207, 1985.
[29] S. A. Sauter. Cubature techniques for 3-d Galerkin BEM. In W. Hackbusch and

G. Wittum, editors, Boundary Elements: Implementation and Analysis of Advanced

G. Wittum, editors, Boundary Elements: Implementation and Analysis of Advanced
Algorithms, pages 29–44. Vieweg-Verlag, 1996.
[30] S. A. Sauter and C. Schwab. Boundary Element Methods. Springer, 2011.
[31] E. E. Tyrtyshnikov. Mosaic-skeleton approximation. Calcolo, 33:47–57, 1996.
[32] E. E. Tyrtyshnikov. Incomplete cross approximation in the mosaic-skeleton method.
Computing, 64:367–380, 2000.
[33] L. Ying, G. Biros, and D. Zorin. A kernel-independent adaptive fast multipole
algorithm in two and three dimensions. J. Comp. Phys., 196(2):591–626, 2004.
26