An analysis of a butterfly algorithm

TL;DR

This paper refines the stability analysis of a butterfly algorithm that uses Chebyshev interpolation to efficiently compress oscillatory integral operators, enhancing understanding of its error bounds.

Contribution

It provides an improved stability analysis of the butterfly algorithm with Chebyshev interpolation, leading to better error estimates.

Findings

Enhanced stability estimates for the butterfly algorithm

Improved error bounds for oscillatory kernel compression

Refined analysis supports more reliable application in practice

Abstract

Butterfly algorithms are an effective multilevel technique to compress discretizations of integral operators with highly oscillatory kernel functions. The particular version of the butterfly algorithm considered here realizes the transfer between levels by Chebyshev interpolation. We present a refinement of the analysis that improves the stability estimates underlying the error bounds.

Tables

Table 1. Table 1: Estimated spectral errors for the butterfly approximation

Table 2. Table 2: Frobenius errors for the butterfly approximation

Equations

$$k(x,y)=\exp (\mathbf{i}\kappa \Phi (x,y))A(x,y)$$

$$\underset{i}{\mathrm{diam}}F:=\sup \{x_i|x\in F\}-\inf \{x_i|x\in F\}.$$

$${\mathcal{E}}_{\rho }^{[a,b]}:=\frac{a+b}{2}+\frac{b-a}{2}{\mathcal{E}}_{\rho }.$$

$${\mathcal{E}}_{\mathbf{\rho }}$$

$$[\mathbf{a},\mathbf{b}]:=B^{[\mathbf{a},\mathbf{b}]}:=\prod _{i=1}^d[a_i,b_i].$$

$$M:=(m+1{)}^d=\dim {\mathcal{Q}}_m.$$

$$\Vert u-I_m^{[a,b]}[u]{\Vert }_{L^{\infty }([a,b])}$$

$$\Vert u-I_m^{[\mathbf{a},\mathbf{b}]}[u]{\Vert }_{L^{\infty }([\mathbf{a},\mathbf{b}])}$$

$$t\mapsto u_{x,i}(t)=u(x_1,\dots ,x_{i-1},t,x_{i+1},\dots ,x_d).$$

$$I_m^{[\mathbf{a},\mathbf{b}]}[f]=\sum _{i=1}^{M}f({\xi }_i^{[\mathbf{a},\mathbf{b}]})L_{i,[\mathbf{a},\mathbf{b}]},$$

$$R_{x_0,y_0}(x,y):=\hat{\Phi }(x,y)-\hat{\Phi }(x,y_0)-\hat{\Phi }(x_0,y)+\hat{\Phi }(x_0,y_0)$$

$$R_{x_0,y_0}(x,y)=(x-x_0{)}^{\top }\hat{G}(x,y)(y-y_0),$$

$$|{\hat{G}}_{ij}(x,y)|\le \frac{1}{d_K^2}\Vert \hat{\Phi }{\Vert }_{L^{\infty }(\Omega )}\text{ for all }(x,y)\in K.$$

$$R_{x_0,y_0}(x,y)$$

$$={\int }_{s=0}^1{\int }_{t=0}^1\sum _{j=1}^d\sum _{i=1}^d{\partial }_{x_j}{\partial }_{y_i}\hat{\Phi }(x_s,y_t)(x-x_0{)}_j(y-y_0{)}_{i}dtds;$$

$$R_{x_0,y_0}(x,y)=\Phi (x,y)-\Phi (x,y_0)-\Phi (x_0,y)+\Phi (x_0,y_0)$$

$$\exp (\mathbf{i}\kappa \Phi (x,y))$$

$$=\exp (\mathbf{i}\kappa \Phi (x,y_0))\exp (\mathbf{i}\kappa \Phi (x_0,y))\exp (\mathbf{i}\kappa (R_{x_0,y_0}(x,y)-\Phi (x_0,y_0))),$$

$$k(x,y)=\exp (\mathbf{i}\kappa \Phi (x,y_0))\exp (\mathbf{i}\kappa \Phi (x_0,y))\underset{=:k_{x_0,y_0}(x,y)}{\underbrace{k(x,y)\exp (-\mathbf{i}\kappa (\Phi (x,y_0)+\Phi (x_0,y)))}}.$$

$$k(x,y)\approx \exp (\mathbf{i}\kappa \Phi (x,y_0))\exp (\mathbf{i}\kappa \Phi (x_0,y))I_m^{B_0^X\times B_0^Y}[k_{x_0,y_0}](x,y).$$

$$k(x,y)$$

$$=\sum _{p,q=1}^M{L}_{p,B_0^X,y_0}^x(x)L_{q,B_0^Y,x_0}^y(y)k_{x_0,y_0}({\xi }_p^{B_0^X},{\xi }_q^{B_0^Y}),$$

$$L_{p,B_0^X,y_0}^x$$

$$\left({\mathfrak{I}}_{y_0}^{B_0^X,x}\otimes {\mathfrak{I}}_{x_0}^{B_0^Y,y}\right)[k],$$

$${\mathfrak{I}}_z^{B,x}[f]$$

$$L_{p,B_0^X,y_0}^x{|}_{B_1^X}$$

$$L_{p,B_0^X,y_0}^x{|}_{B_L^X}$$

$$k^{BF}:=\left({\mathfrak{I}}_{y_{-L}}^{B_L^X,x}\circ \cdots \circ {\mathfrak{I}}_{y_0}^{B_0^X,x}\right)\otimes \left({\mathfrak{I}}_{x_{-L}}^{B_L^Y,y}\circ \cdots \circ {\mathfrak{I}}_{x_0}^{B_0^Y,y}\right)[k].$$

$${\mathbf{K}}_{i,j}={\int }_{x\in {\Gamma }^{\prime }}{\int }_{y\in \Gamma }k(x,y){\phi }_i(x){\psi }_j(y)dydx,i\in \mathcal{I},j\in \mathcal{J}.$$

$${\sigma }_{\ell +1}\in \mathrm{sons}({\sigma }_{\ell })\text{ for all }\ell \in \{0,\dots ,n-1\}.$$

$$\mathrm{supp}{\phi }_i\subset B_{\sigma }\text{ for all }i\in \sigma ,$$

$${\mathbf{K}}_{i,j}$$

Full text

add1]sb@informatik.uni-kiel.de

add1]chb@informatik.uni-kiel.de add2]melenk@tuwien.ac.at

1 Introduction

Nonlocal operators with highly oscillatory kernel functions arise in many applications. Prominent examples of such operators include the Fourier transform, special function transforms, and integral operators whose kernels are connected to the high-frequency Helmholtz or Maxwell’s equations. Many algorithms with log-linear complexity have been proposed in the past. A very incomplete list includes: for the non-uniform Fourier transform [16] (and references therein), [1, 6, 17]; for special function transforms [21]; for Helmholtz and Maxwell integral operators [5, 23, 20, 13, 8, 7, 11, 19, 2, 4].

Underlying the log-linear complexity of the above algorithms for the treatment of high-frequency kernels is the use of a multilevel structure and in particular efficient transfer operators between the levels. In the language of linear algebra and with $L$ denoting the number of levels, the matrix realizing the operator is (approximately) factored into $O(L)$ matrices, each of which can be treated with nearly optimal complexity $O(N\log^{\alpha}N)$, where $N$ denotes the problem size. This observation is formulated explicitly in [18] in connection with a class of butterfly algorithms.

A central issue for a complete error analysis of such algorithms is that of stability. That is, since the factorization into $O(L)$ factors results from an approximation process, the application of each of the $O(L)$ steps incurs an error whose propagation in the course of a realization of the matrix-vector multiplications needs to be controlled. For some algorithms, a full error analysis is possible. Here, we mention the fast multipole method [23] (and its stable realizations [8, 13]) for the Helmholtz kernel, which exploits that a specific kernel is analyzed. One tool to construct suitable factorizations for more general high-frequency kernels is polynomial interpolation (typically tensor product Chebyshev interpolation) as proposed, e.g., in [17, 6] for Fourier transforms and generalizations of the Fourier transform and in [19, 4] for Helmholtz kernels. A full analysis of the case of the (classical) Fourier transform is developed in [17] and for the algorithm proposed in [6, Sec. 4] in [9]. For the procedure proposed in [19] for the Helmholtz kernel, a detailed analysis is given in [4]. Based on generalizations of the tools developed in the latter work [4] the novel contribution of the present work is a sharper analysis of the butterfly algorithm proposed in [6, Sec. 4] for general integral operators with oscillatory kernel functions. Indeed, the present analysis of the butterfly algorithm of [6, Sec. 4] improves over [9] and [17] in that the necessary relation between the underlying polynomial degree $m$ and the depth $L$ of the cluster tree is improved from $m\geq CL\log(\log L)$ to $m\geq C\log L$ (cf. Theorem 2.6 and Theorem 3.1 in comparision to [9, Sec. 4.3]); we should, however, put this improvement in perspective by mentioning that a requirement $m\geq CL$ typically results from accuracy considerations.

It is worth stressing that, although techniques that base the transfer between levels on polynomial interpolation are amenable to a rigorous stability analysis and may even lead to the desirable log-linear complexity, other techniques are available that are observed to perform better in practice. We refer the reader to [18] for an up-to-date discussion of techniques associated with the name “butterfly algorithms” and to [2] for techniques related to adaptive cross approximation (ACA).

In the present work, we consider kernel functions $k$ of the form

[TABLE]

on a product $B_{X}\times B_{Y}$, where $B_{X}$, $B_{Y}\subset{\mathbb{R}}^{d}$ are axis-parallel boxes. The phase function $\Phi$ is assumed to be real analytic on $B_{X}\times B_{Y}$, in particular, it is assumed to be real-valued in this box. The amplitude function $A$ is likewise assumed to be analytic, although it could be complex-valued or even vector-valued. The setting of interest for the parameter $\kappa\in{\mathbb{R}}$ is that of $\kappa\gg 1$.

An outline of the paper is as follows. We continue the introduction in Section 1.1 with notation that will be used throughout the paper. In Section 1.2, we discuss how a butterfly representation for kernels of the form (1.1) can be obtained with the aid of an iterated Chebyshev interpolation procedure. The point of view taken in Section 1.2 is one of functions and their approximation in separable form. The following Section 1.3, therefore, focuses on the realization of the butterfly structure on the matrix level. As alluded to in the introduction, butterfly structures are one of several techniques for highly oscillatory kernel. In Section 1.3.3, we discuss the relation of butterfly techniques with directional ${\mathcal{H}}^{2}$-matrices, [2, 3, 4]. The discussion in Section 1.2 concentrated on the case of analytic kernel functions. Often, however, kernel functions are only asymptotically smooth, e.g., if they are functions of the Euclidean distance $\|x-y\|$. We propose in Section 1.3.4 to address this issue by combining the butterfly structure with a block clustering based on the “standard” admissibility condition that takes the distance of two clusters into account. We mention in passing that alternative approaches are possible to deal with certain types of singularities (see, e.g., the transform technique advocated in [6, Sec. 1.1] to account for the Euclidean distance $\|x-y\|$).

Section 2 is concerned with an analysis of the errors incurred by the approximations done to obtain the butterfly structure, which is enforced by an iterated interpolation process. The stability of one step of this process is analyzed in the univariate case in Lemma 2.2; the multivariate case is inferred from that by tensor product arguments in Lemma 2.4. The final stability analysis of the iterated process is given in Theorem 2.5.

Section 3 specializes the analysis of Theorem 2.5 to the single layer operator for the Helmholtz equation. This operator is also used in the numerical examples in Section 4.

1.1 Notation and preliminaries

We start with some notation: $B_{\varepsilon}(z)$ denotes the (Euclidean) ball of radius $\varepsilon>0$ centered at $z\in{\mathbb{C}}^{d}$. For a bounded set $F\subset{\mathbb{R}}^{d}$, we denote by $\operatorname*{diam}\nolimits_{i}F$, $i\in\{1,\ldots,d\}$, the extent of $F$ in the $i$-th direction:

[TABLE]

For $\rho>1$ the Bernstein elliptic discs are given by ${\mathcal{E}}_{\rho}:=\{z\in{\mathbb{C}}\,|\,|z-1|+|z+1|<\rho+1/\rho\}.$ More generally, for $[a,b]\subset{\mathbb{R}}$, we also use the scaled and shifted elliptic discs

[TABLE]

In a multivariate setting, we collect $d$ values $\rho_{i}>1$, $i=1,\dots,d$, in the vector $\boldsymbol{\rho}$ and define, for intervals $[a_{i},b_{i}]\subset{\mathbb{R}}$, $i=1,\ldots,d$, the elliptic polydiscs (henceforth simply called “polydiscs”)

[TABLE]

We will write ${\mathcal{E}}_{{\rho}}$ and ${\mathcal{E}}^{[\boldsymbol{a},\boldsymbol{b}]}_{{\rho}}$ if $\boldsymbol{\rho}_{i}=\rho$ for $i=1,\ldots,d$. Axis-parallel boxes are denoted

[TABLE]

In fact, throughout the text, a box $B$ is always understood to be a set of the form (1.3). For vector-valued objects such as $\boldsymbol{\rho}$ we will use the notation $\boldsymbol{\rho}>1$ in a componentwise sense.

We employ a univariate polynomial interpolation operator $I_{m}:C([-1,1])\rightarrow{\mathcal{P}}_{m}$ with Lebesgue constant $\Lambda_{m}$ (e.g., the Chebyshev interpolation operator). Here, ${\mathcal{P}}_{m}=\operatorname*{span}\{x^{i}\,|\,0\leq i\leq m\}$ is the space of (univariate) polynomials of degree $m$. We write ${\mathcal{Q}}_{m}:=\operatorname*{span}\{x_{1}^{i_{1}}x_{2}^{i_{2}}\cdots x_{d}^{i_{d}}\,|\,0\leq i_{1},\ldots,i_{d}\leq m\}$ for the space of $d$-variate polynomials of degree $m$ (in each variable). Throughout the text, we will use

[TABLE]

The interpolation operator $I_{m}$ may be scaled and shifted to an interval $[a,b]\subset{\mathbb{R}}$ and is then denoted by $I^{[a,b]}_{m}$. Tensor product interpolation on the box $[\boldsymbol{a},\boldsymbol{b}]$ is correspondingly denoted $I^{[\boldsymbol{a},\boldsymbol{b}]}_{m}:C([\boldsymbol{a},\boldsymbol{b}])\rightarrow{\mathcal{Q}}_{m}$. We recall the following error estimates:

Lemma 1.1**.**

[TABLE]

where, for given $x\in[\boldsymbol{a},\boldsymbol{b}]$, the univariate functions $u_{x,i}$ are defined by

[TABLE]

It will sometimes be convenient to write the interpolation operator $I^{[\boldsymbol{a},\boldsymbol{b}]}_{m}$ explicitly as

[TABLE]

where $\xi^{[\boldsymbol{a},\boldsymbol{b}]}_{i}$, $i=1,\ldots,M$, are the interpolation points and $L_{i,[\boldsymbol{a},\boldsymbol{b}]}$, $i=1,\ldots,M$, are the associated Lagrange interpolation polynomials. The following lemma is a variant of a result proved in [6]:

Lemma 1.2**.**

Let $\Omega\subset{\mathbb{C}}^{2d}$ be open. Let $(x_{0},y_{0})\in\Omega$. Let $(x,y)\mapsto\widehat{\Phi}(x,y)$ be analytic on $\Omega$. Then the function

[TABLE]

can be written in the form

[TABLE]

where the entries $\widehat{G}_{ij}$, $i$, $j=1,\dots,d$, of the matrix $\widehat{G}$ are analytic on $\Omega$. Furthermore, for any convex $K\subset\Omega$ with $(x_{0},y_{0})\in K$ and $d_{K}:=\sup\{\varepsilon>0\,|\,B_{\varepsilon}(x)\times B_{\varepsilon}(y)\subset\Omega\quad\text{ for all }(x,y)\in K\}>0$ there holds

[TABLE]

Proof.

Let us first show the estimate for the convex $K$. Let $(x,y)\in K$. We denote the parametrizations of straight lines from $x_{0}$ to $x$ and from $y_{0}$ to $y$ by $x_{s}:=x_{0}+s(x-x_{0})$ and $y_{t}:=y_{0}+t(y-y_{0})$ for $s,t\in[0,1]$. By integrating along the second line for fixed $x$ and fixed $x_{0}$ and along the first line for fixed $y_{t}$, we arrive at

[TABLE]

here, we abbreviated, for example, $(y-y_{0})_{i}$ for the $i$-th component of the vector $(y-y_{0})$. Using [15, Cor. 2.2.5] one can show that the double integral indeed represents an analytic function $K\ni(x,y)\mapsto\widehat{G}(x,y)$. In order to bound $\widehat{G}$, one observes that the integrand involves the partial derivatives of $\widehat{\Phi}$ with respect to the variables $x_{j}$, $y_{i}$, $i$, $j=1,\ldots,d$, which can be estimated in view of the Cauchy integral formula (cf., e.g., [15, Thm. 2.2.1]).

To see that the coefficients $G_{ij}$ are analytic on $\Omega$ (and not just on $K$), we note that we have to define $G_{ij}$ on $\Omega\setminus\operatorname*{closure}K$ as $G_{ij}(x,y)=R_{x_{0},y_{0}}(x,y)/((x-x_{0})_{i}(y-y_{0})_{i})$. Now that $G_{ij}$ is defined on $\Omega$, we observe that by the analyticity of $R_{x_{0},y_{0}}$, is it analytic in each variable separately and thus, by Hartogs’ theorem (see, e.g., [15, Thm. 2.2.8]) analytic on $\Omega$. ∎

1.2 Butterfly representation: the heart of the matter

Typically, the key ingredient of fast summation schemes is the approximation of the kernel function by short sums of products of functions depending only on $x$ or $y$. Following [6, Sec. 4] we achieve this by applying a suitable modification to the kernel function $k$ and then interpolating the remaining smooth term on a domain $B^{X}_{0}\times B^{Y}_{0}$, where $B^{X}_{0}$ and $B^{Y}_{0}$ are two boxes and $x_{0}\in B^{X}_{0}$, $y_{0}\in B^{Y}_{0}$.

According to Lemma 1.2, we can expect

[TABLE]

to be “small” if the product $\|x-x_{0}\|\,\|y-y_{0}\|$ is small. For fixed $x_{0}$ and $y_{0}$ we obtain

[TABLE]

and we observe that the first term on the right-hand side depends only on $x$, the second one only on $y$, while the third one is smooth if the product $\kappa\|x-x_{0}\|\,\|y-y_{0}\|$ is small, since the modified phase function $R_{x_{0},y_{0}}(x,y)$ takes only small values and $\Phi(x_{0},y_{0})$ is constant. This observation allows us to split the kernel function $k$ into oscillatory factors depending only on $x$ and $y$, respectively, and a smooth factor $k_{x_{0},y_{0}}$ that can be interpolated:

[TABLE]

Applying the polynomial interpolation operator to $k_{x_{0},y_{0}}$ yields

[TABLE]

Written with the $M=(m+1)^{d}$ interpolation points $(\xi_{p}^{B^{X}_{0}})_{p=1}^{M}\subset B^{X}_{0}$ and $(\xi_{q}^{B^{Y}_{0}})_{q=1}^{M}\subset B^{Y}_{0}$ and the corresponding Lagrange polynomials $L_{p,B^{X}_{0}}$, $L_{q,B^{Y}_{0}}$, we have

[TABLE]

where the expansion functions $L^{x}_{p,B^{X}_{0},y_{0}}$ and $L^{y}_{q,B^{Y}_{0},x_{0}}$ are defined by

[TABLE]

A short form of the approximation is given by

[TABLE]

where, for a box $B\subset{\mathbb{R}}^{d}$, a point $z\in{\mathbb{R}}^{d}$, and a polynomial degree $m$, we have introduced the operators

[TABLE]

We have seen that the product $\kappa\|x-x_{0}\|\,\|y-y_{0}\|$ controls the smoothness of $k_{x_{0},y_{0}}$, so we can move $y$ away from $y_{0}$ as long as we move $x$ closer to $x_{0}$ without changing the quality of the approximation. This observation gives rise to a multilevel approximation of (1.8) that applies an additional re-interpolation step to the expansion functions $L^{x}_{p,B^{X}_{0},y_{0}}$ and $L^{y}_{q,B^{Y}_{0},x_{0}}$. We only describe the process for $L^{x}_{p,B^{X}_{0},y_{0}}$, since $L^{y}_{q,B^{Y}_{0},x_{0}}$ is treated analogously. Let $(B^{X}_{\ell})_{\ell=0}^{L}$ be a nested sequence of boxes and $(y_{-\ell})_{\ell=0}^{L}$ be a sequence of points. The first step of the iterated re-interpolation process is given by

[TABLE]

The actual approximation of $L^{x}_{p,B^{X}_{0},y_{0}}$ is then given by iteratively applying the operators ${\mathfrak{I}}^{B^{X}_{\ell},x}_{y_{-\ell}}$, $\ell=0,\ldots,L$, i.e.,

[TABLE]

This process is formalized in the following algorithm:

Algorithm 1.3** (Butterfly representation by interpolation).**

Let two sequences $(B_{\ell}^{X})_{\ell=0}^{L}$, $(B_{\ell}^{Y})_{\ell=0}^{L}$, of nested boxes and sequences of points $(y_{-\ell})_{\ell=0}^{L}$, $(x_{-\ell})_{\ell=0}^{L}$ be given. The butterfly representation of $k$ on $B^{X}_{L}\times B^{Y}_{L}$ is defined by

[TABLE]

In Theorem 2.6 below we will quantify the error $k|_{B^{X}_{L}\times B^{Y}_{L}}-k^{BF}$.

Remark 1.4**.**

*In Algorithm 1.3, the number of levels $L$ is chosen to be the same for the first argument “$x$” and the second argument “$y$”. The above developments show that this is not essential, and Algorithm 1.3 naturally generalizes to a setting with boxes $(B_{\ell}^{X})_{\ell=0}^{L^{X}}$, $(B_{\ell}^{Y})_{\ell=0}^{L^{Y}}$ and corresponding point sequences $(y_{-\ell})_{\ell=0}^{L^{X}}$, $(x_{-\ell})_{\ell=0}^{L^{Y}}$ with $L^{X}\neq L^{Y}$.  *

1.3 Butterfly structures on the matrix level

It is instructive to formulate how the approximation described in Algorithm 1.3 is realized on the matrix level. To that end, we consider the Galerkin discretization of an integral operator $K:L^{2}(\Gamma)\rightarrow L^{2}(\Gamma^{\prime})$ defined by $(K\varphi)(x):=\int_{y\in\Gamma}k(x,y)\varphi(y)$. Let $(\varphi_{i})_{i\in\mathcal{I}}\subset L^{2}(\Gamma^{\prime})$, $(\psi_{j})_{j\in\mathcal{J}}\subset L^{2}(\Gamma)$ be bases of trial and test spaces. We have to represent the Galerkin matrix ${\mathbf{K}}$ with entries

[TABLE]

Following the standard approach for fast multipole [22] and panel-clustering methods [14], the index sets $\mathcal{I}$ and $\mathcal{J}$ are assumed to be organized in cluster trees ${\mathcal{T}}_{\mathcal{I}}$ and ${\mathcal{T}}_{\mathcal{J}}$, where the nodes of the tree are called clusters. We assume that the maximal number of sons of a cluster is fixed. A cluster tree ${\mathcal{T}}_{\mathcal{I}}$ can be split into levels; the root, which is the complete index set $\mathcal{I}$, is at level [math]. We employ the notation ${\mathcal{T}}^{\ell}_{\mathcal{I}}$ and ${\mathcal{T}}^{\ell}_{\mathcal{J}}$ for the clusters on level $\ell$. We use the notation $\operatorname*{sons}(\sigma)$ for the collection of sons of a cluster $\sigma$ (for leaves $\sigma$, we set $\operatorname*{sons}(\sigma)=\emptyset$). We let $\operatorname*{father}(\sigma)$ denote the (unique) father of a cluster $\sigma$ that is not the root. We use $\operatorname*{descendants}(\sigma)$ for the set of descendants of the cluster $\sigma$ (including $\sigma$ itself).

A tuple $(\sigma_{0},\sigma_{1},\ldots,\sigma_{n})$ of clusters is called a cluster sequence if

[TABLE]

We introduce for each cluster $\sigma\in{\mathcal{T}}_{\mathcal{I}}$ a bounding box $B_{\sigma}$, which is an axis-parallel box such that

[TABLE]

and similarly for clusters $\tau\in{\mathcal{T}}_{\mathcal{J}}$ and basis functions $\psi_{j}$.

1.3.1 Butterfly structure in a model situation

We illustrate the butterfly structure (based on interpolation as proposed in [6, Sec. 4]) for a model situation, where the leaves of the cluster trees ${\mathcal{T}}_{\mathcal{I}}$ and ${\mathcal{T}}_{\mathcal{J}}$ are all on the same level. In particular, this assumption implies $0pt({\mathcal{T}}_{\mathcal{I}})=0pt({\mathcal{T}}_{\mathcal{J}})$.

We fix points $x_{\sigma}\in B_{\sigma}$ and $y_{\tau}\in B_{\tau}$ for all $\sigma\in\mathcal{T}_{\mathcal{I}}$ and $\tau\in\mathcal{T}_{\mathcal{J}}$.

For a given pair $(\sigma,\tau)\in\mathcal{T}_{\mathcal{I}}\times\mathcal{T}_{\mathcal{J}}$, combining the intermediate decomposition (1.8) with (1.13) yields

[TABLE]

If $\sigma$ is not a leaf of $\mathcal{T}_{\mathcal{I}}$, we employ the additional approximation (1.12). We choose a middle level $L\in{\mathbb{N}}_{0}$ and cluster sequences $(\sigma_{-L},\ldots,\sigma_{0},\ldots,\sigma_{L})$ and $(\tau_{-L},\ldots,\tau_{0},\ldots,\tau_{L})$ of clusters such that $\sigma_{0}=\sigma$, $\tau_{0}=\tau$. Since each cluster has at most one father, the clusters $\sigma_{-L},\ldots,\sigma_{0}$ are uniquely determined by $\sigma=\sigma_{0}$ and $\tau_{-L},\ldots,\tau_{0}$ are uniquely determined by $\tau=\tau_{0}$. The approximation (1.12) implies for $i\in\sigma_{L}$

[TABLE]

The first interpolation step is given by

[TABLE]

where the transfer matrix $\mathbf{E}^{\sigma_{1},\sigma_{0},\tau_{0},\tau_{-1}}$ is given by

[TABLE]

The re-expansion (1.15) implies

[TABLE]

so we can avoid storing $\mathbf{V}^{\sigma_{\ell},\tau_{-\ell}}$ by storing the $M\times M$-matrix $\mathbf{E}^{\sigma_{\ell+1},\sigma_{\ell},\tau_{-\ell},\tau_{-\ell-1}}$ and $\mathbf{V}^{\sigma_{\ell+1},\tau_{-\ell-1}}$ instead. A straightforward induction yields that we only have to store $\mathbf{V}^{\sigma_{L},\tau_{-L}}$ and the transfer matrices.

Algorithm 1.5** (Butterfly representation of matrices).**

Let $L=\lfloor 0pt(\mathcal{T}_{\mathcal{I}})/2\rfloor$ and $L^{\text{middle}}=0pt(\mathcal{T}_{\mathcal{I}})-L=\lceil 0pt(\mathcal{T}_{\mathcal{I}})/2\rceil\geq L$.

Compute, for all $\sigma\in\mathcal{T}_{\mathcal{I}}^{L^{\text{middle}}}$, $\tau\in\mathcal{T}_{\mathcal{J}}^{L^{\text{middle}}}$, the coupling matrices

[TABLE] 2. 2.

Compute the transfer matrices

[TABLE]

for all $\ell\in\{0,\ldots,L-1\}$, $\sigma_{\ell}\in\mathcal{T}_{\mathcal{I}}^{L^{\text{middle}}+\ell}$, $\sigma_{\ell+1}\in\operatorname*{sons}(\sigma_{\ell})$, $\tau_{-\ell-1}\in\mathcal{T}_{\mathcal{J}}^{L^{\text{middle}}-\ell-1}$, $\tau_{-\ell}\in\operatorname*{sons}(\tau_{-\ell-1})$, and

[TABLE]

for all $\ell\in\{0,\ldots,L-1\}$, $\sigma_{-\ell-1}\in\mathcal{T}_{\mathcal{I}}^{L^{\text{middle}}-\ell-1}$, $\sigma_{-\ell}\in\operatorname*{sons}(\sigma_{-\ell-1})$, $\tau_{\ell}\in\mathcal{T}_{\mathcal{J}}^{L^{\text{middle}}+\ell}$, $\tau_{\ell+1}\in\operatorname*{sons}(\tau_{\ell})$. 3. 3.

Compute the leaf matrices

[TABLE]

for all $\sigma_{L}\in\mathcal{T}^{L^{\text{middle}}+L}_{\mathcal{I}}$, $\tau_{-L}\in\mathcal{T}^{L^{\text{middle}}-L}_{\mathcal{J}}$, $i\in\sigma_{L}$, $p\in\{1,\ldots,M\}$, and

[TABLE]

for all $\tau_{L}\in\mathcal{T}^{L^{\text{middle}}+L}_{\mathcal{J}}$, $\sigma_{-L}\in\mathcal{T}^{L^{\text{middle}}-L}_{\mathcal{I}}$, $j\in\tau_{L}$, $q\in\{1,\ldots,M\}$. 4. 4.

For leaf clusters $\sigma_{L}\in\mathcal{T}_{\mathcal{I}}$ and $\tau_{L}\in\mathcal{T}_{\mathcal{J}}$ there are uniquely determined cluster sequences $(\tau_{-L},\ldots,\tau_{0},\ldots,\tau_{L})$ and $(\sigma_{-L},\ldots,\sigma_{0},\ldots,\sigma_{L})$. The matrix $\mathbf{K}$ is approximated by

[TABLE]

Remark 1.6**.**

The costs of representing a butterfly matrix are as follows (for even $\operatorname*{depth}{\mathcal{T}}_{\mathcal{I}}$ and $L^{\text{middle}}=L=\operatorname*{depth}({\mathcal{T}}_{\mathcal{I}}/2$):

For the coupling matrices ${\mathbf{S}}^{\sigma\times\tau}$ on the middle level $L^{\text{middle}}$: $M^{2}|{\mathcal{T}}^{L^{\text{middle}}}_{\mathcal{I}}||{\mathcal{T}}^{L^{\text{middle}}}_{\mathcal{J}}|$ 2. 2.

For the transfer matrices $\mathbf{E}^{\sigma_{\ell+1},\sigma_{\ell},\tau_{-\ell-1},\tau_{\ell}}$: $\sum_{\ell=0}^{L-1}M^{2}|{\mathcal{T}}^{L^{\text{middle}}+\ell}_{\mathcal{I}}||{\mathcal{T}}^{L^{\text{middle}}-\ell}_{\mathcal{J}}|$ 3. 3.

For the transfer matrices $\mathbf{E}^{\tau_{\ell+1},\tau_{\ell},\sigma_{-\ell-1},\sigma_{-\ell}}$: $\sum_{\ell=0}^{L-1}M^{2}|{\mathcal{T}}^{L^{\text{middle}}-\ell}_{\mathcal{I}}||{\mathcal{T}}^{L^{\text{middle}}+\ell}_{\mathcal{J}}|$ 4. 4.

For the leaf matrices ${\mathbf{V}}^{\sigma_{L},\tau_{-L}}$ and ${\mathbf{W}}^{\tau_{L},\sigma_{-L}}$ with leaves $\sigma_{L}\in\mathcal{T}_{\mathcal{I}}$ and $\tau_{L}\in\mathcal{T}_{\mathcal{J}}$: $\sum_{\sigma\in{\mathcal{T}}_{\mathcal{I}}}M|\sigma|$ and $\sum_{\tau\in{\mathcal{T}}_{\mathcal{J}}}M|\tau|$.

In a model situation with $|\mathcal{I}|=|\mathcal{J}|=N$ and balanced binary trees ${\mathcal{T}}_{\mathcal{I}}$, ${\mathcal{T}}_{\mathcal{J}}$ of depth $2L=O(\log N)$ and leaves of ${\mathcal{T}}_{\mathcal{I}}$, ${\mathcal{T}}_{\mathcal{J}}$ that have at most $n_{\mathrm{leaf}}$ elements we get

[TABLE]

*We expect for approximation-theoretical reasons that the polynomial degree satisfies $m=O(\log N)$. Since $M=(m+1)^{d}$, the total complexity is then $O(N\log^{2d+1}N)$.  *

Remark 1.7**.**

*The butterfly structure presented here is suitable for kernel functions with analytic phase function $\Phi$ and amplitude function $A$. When these functions are “asymptotically smooth”, for example, when they are functions of the Euclidean distance $(x,y)\mapsto\|x-y\|$, a modification is necessary to take care of the singularity at $x=y$. For example, one could create a block partition that applies the approximation scheme only to pairs $(\sigma,\tau)$ of clusters that satisfy the standard admissibility condition $\max\{\operatorname*{diam}\nolimits(B_{\sigma}),\operatorname*{diam}\nolimits(B_{\tau})\}\leq\operatorname*{dist}\nolimits(B_{\sigma},B_{\tau})$. Each block $\mathbf{K}|_{\sigma\times\tau}$ that satisfies this condition is treated as a butterfly matrix in the above sense. We illustrate this procedure in Section 1.3.4 below for general asymptotically smooth kernel functions $k$ and specialize to the 3D Helmholtz kernel in Section 3.  *

1.3.2 ${\mathcal{D}}{\mathcal{H}}^{2}$-matrices

It is worth noting that the above butterfly structure can be interpreted as a special case of directional ${\mathcal{H}}^{2}$-matrices (short: ${\mathcal{D}}{\mathcal{H}}^{2}$-matrices) as introduced in [2, 3, 4] in the context of discretizations of Helmholtz integral operators.

Let us recall the definition of a ${\mathcal{D}}{\mathcal{H}}^{2}$-matrix ${\mathbf{K}}\in{\mathbb{C}}^{\mathcal{I}\times\mathcal{J}}$ with cluster trees ${\mathcal{T}}_{\mathcal{I}}$, ${\mathcal{T}}_{\mathcal{J}}$.

Definition 1.8** (Directional cluster basis for ${\mathcal{T}}_{\mathcal{I}}$).**

For each cluster $\sigma\in{\mathcal{T}}_{\mathcal{I}}$, let ${\mathcal{D}}_{\sigma}$ be a given index set. Let ${\mathcal{V}}=({\mathbf{V}}^{\sigma,c})_{\sigma\in{\mathcal{T}}_{\mathcal{I}},c\in{\mathcal{D}}_{\sigma}}$ be a two-parameter family of matrices. This family is called a directional cluster basis with rank $M$ if

${\mathbf{V}}^{\sigma,c}\in{\mathbb{C}}^{\sigma\times M}$* for all $\sigma\in{\mathcal{T}}_{\mathcal{I}}$ and $c\in\mathcal{D}_{\sigma}$, and*

• 2.

there is, for every $\sigma$ that is not a leaf of ${\mathcal{T}}_{\mathcal{I}}$ and every $\sigma^{\prime}\in\operatorname*{sons}(\sigma)$ and every $c\in{\mathcal{D}}_{\sigma}$, an element $c^{\prime}\in{\mathcal{D}}_{\sigma^{\prime}}$ and a matrix ${\mathbf{E}}^{\sigma^{\prime},\sigma,c^{\prime},c}\in{\mathbb{C}}^{M\times M}$ such that

[TABLE]

The matrices ${\mathbf{E}}^{\sigma^{\prime},\sigma,c^{\prime},c}$ are called transfer matrices for the directional cluster basis.

${\mathcal{D}}{\mathcal{H}}^{2}$-matrices are blockwise low-rank matrices. To describe the details of this structure, let ${\mathcal{T}}_{\mathcal{I}\times\mathcal{J}}$ be a block tree based on the cluster trees ${\mathcal{T}}_{\mathcal{I}}$ and ${\mathcal{T}}_{\mathcal{J}}$. Specifically, we assume that a) the root of ${\mathcal{T}}_{\mathcal{I}\times\mathcal{J}}$ is $\mathcal{I}\times\mathcal{J}$, b) each node of ${\mathcal{T}}_{\mathcal{I}\times\mathcal{J}}$ is of the form $(\sigma,\tau)\in{\mathcal{T}}_{\mathcal{I}}\times{\mathcal{T}}_{\mathcal{J}}$, and c) for every node $(\sigma,\tau)\in{\mathcal{T}}_{\mathcal{I}\times\mathcal{J}}$ we have

[TABLE]

We denote the leaves of the block tree ${\mathcal{T}}_{\mathcal{I}\times\mathcal{J}}$ by

[TABLE]

The leaves form a disjoint partition of $\mathcal{I}\times\mathcal{J}$, so a matrix ${\mathbf{G}}$ is uniquely determined by the submatrices ${\mathbf{G}}|_{\sigma\times\tau}$ for $b=(\sigma,\tau)\in\mathcal{L}_{\mathcal{I}\times\mathcal{J}}$. The set of leaves $\mathcal{L}_{\mathcal{I}\times\mathcal{J}}$ is written as the disjoint union $\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{+}\dot{\cup}\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{-}$ of two sets, which are called the admissible leaves $\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{+}$, corresponding to submatrices that can be approximated, and the inadmissible leaves $\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{-}$, corresponding to small submatrices that have to be stored explicitly. We are now in a position to define $\mathcal{DH}^{2}$-matrices as in [3, 4]:

Definition 1.9** (Directional $\mathcal{H}^{2}$-matrix).**

Let ${\mathcal{V}}$ and ${\mathcal{W}}$ be directional cluster bases of rank $M$ for ${\mathcal{T}}_{\mathcal{I}}$ and ${\mathcal{T}}_{\mathcal{J}}$, respectively. A matrix ${\mathbf{G}}\in{\mathbb{C}}^{\mathcal{I}\times\mathcal{J}}$ is called a directional $\mathcal{H}^{2}$-matrix (or simply: a $\mathcal{DH}^{2}$-matrix) if there are families ${\mathcal{S}}=({\mathbf{S}}_{b})_{b\in\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{+}}$ and $(c_{b}^{\sigma})_{b\in\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{+}}$, $(c_{b}^{\tau})_{b\in\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{+}}$ such that

${\mathbf{S}}_{b}\in{\mathbb{C}}^{M\times M}$* for all $b=(\sigma,\tau)\in\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{+}$, and*

• 2.

${\mathbf{G}}|_{\sigma\times\tau}={\mathbf{V}}^{\sigma,c^{\sigma}_{b}}{\mathbf{S}}_{b}({\mathbf{W}}^{\tau,c^{\tau}_{b}})^{\top}$* with $c^{\sigma}_{b}\in\mathcal{D}_{\sigma}$, $c^{\tau}_{b}\in\mathcal{D}_{\tau}$ for all $b=(\sigma,\tau)\in\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{+}$.*

The elements of the family ${\mathcal{S}}$ are called coupling matrices. The cluster bases ${\mathcal{V}}$ and ${\mathcal{W}}$ are called the row cluster basis and column cluster basis, respectively. A $\mathcal{DH}^{2}$-matrix representation of a $\mathcal{DH}^{2}$-matrix ${\mathbf{G}}$ consists of ${\mathcal{V}}$, ${\mathcal{W}}$, ${\mathcal{S}}$ and the family $({\mathbf{G}}|_{\sigma\times\tau})_{b=(\sigma,\tau)\in\mathcal{L}_{\mathcal{I}\times\mathcal{I}}^{-}}$ of nearfield matrices corresponding to the inadmissible leaves of ${\mathcal{T}}_{\mathcal{I}\times\mathcal{J}}$.

1.3.3 The butterfly structure as a special $\mathcal{DH}^{2}$-matrix

We now show that the butterfly structure discussed in Section 1.3.1 can be understood as a ${\mathcal{D}}{\mathcal{H}}^{2}$-matrix: for $L=\lfloor 0pt(\mathcal{T}_{\mathcal{I}})/2\rfloor$ and the middle level $L^{\text{middle}}:=0pt(\mathcal{T}_{\mathcal{I}})-L=\lceil 0pt(\mathcal{T}_{\mathcal{I}})/2\rceil$, we let

[TABLE]

The key is to observe that the sets ${\mathcal{D}}_{\sigma}$ associated with a cluster $\sigma\in{\mathcal{T}}^{L^{\text{middle}}+\ell}_{\mathcal{I}}$ on level $L^{\text{middle}}+\ell$ are taken to be points $y_{\tau}$ with $\tau\in{\mathcal{T}}^{L^{\text{middle}}-\ell}_{\mathcal{J}}$:

[TABLE]

Analogously, we define the sets ${\mathcal{D}}_{\tau}$ for $\tau\in{\mathcal{T}}_{\mathcal{J}}$. The transfer matrices ${\mathbf{E}}^{\sigma^{\prime},\sigma,c^{\prime},c}$ appearing in Definition 1.8 are those of Section 1.3.1, and the same holds for the leaf matrices $\mathbf{V}^{\sigma_{L},c}$ and $\mathbf{W}^{\tau_{L},c}$.

1.3.4 Butterfly structures for asymptotically smooth kernels in a model situation

The model situation of Section 1.3.1 is appropriate when the kernel function $k$ is analytic. Often, however, the kernel function $k$ is only “asymptotically smooth”, i.e., it satisfies estimates of the form

[TABLE]

for some $C$, $\gamma>0$, $\delta\in{\mathbb{R}}$. Prominent examples include kernel function such as the 3D Helmholtz kernel (3.1) where the dependence on $(x,y)$ is through the Euclidean distance $\|x-y\|$.

We describe the data structure for an approximation of the stiffness matrix ${\mathbf{K}}$ given by (1.13) for asymptotically smooth kernels. We will study a restricted setting that focuses on the essential points and is geared towards kernels such as the Helmholtz kernel (3.1).

Again, let the index sets $\mathcal{I}$ and $\mathcal{J}$ be organized in trees ${\mathcal{T}}_{\mathcal{I}}$ and ${\mathcal{T}}_{\mathcal{J}}$ with a bounded number of sons. We assume that the trees are balanced and that all leaves are on the same level $0pt({\mathcal{T}}_{\mathcal{I}})=0pt({\mathcal{T}}_{\mathcal{J}})$. Recall the notion of bounding box in (1.14). It will also be convenient to introduce for $\sigma\in{\mathcal{T}}_{\mathcal{I}}$ and $\tau\in{\mathcal{T}}_{\mathcal{J}}$ the subtrees

[TABLE]

with roots $\sigma$ and $\tau$, respectively, and the clusters on level $\ell$:

[TABLE]

Concerning the block cluster tree ${\mathcal{T}}_{\mathcal{I}\times\mathcal{J}}$, we assume that its leaves $\mathcal{L}_{\mathcal{I}\times\mathcal{J}}=\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{+}\dot{\cup}\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{-}$ are created as follows:

Apply a clustering algorithm to create a block tree ${\mathcal{T}}^{standard}_{\mathcal{I}\times\mathcal{J}}$ based on the trees ${\mathcal{T}}_{\mathcal{I}}$ and ${\mathcal{T}}_{\mathcal{J}}$ according to the standard admissibility condition

[TABLE]

for a fixed admissibility parameter $\eta_{1}>0$. 2. 2.

The leaves of ${\mathcal{T}}^{standard}_{\mathcal{I}\times\mathcal{J}}$ are split as $\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{standard,+}\dot{\cup}\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{standard,-}$ into admissible leaves $\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{standard,+}$ and inadmissible leaves $\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{standard,-}$. 3. 3.

Set $\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{-}:=\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{standard,-}$. 4. 4.

For all $(\widehat{\sigma},\widehat{\tau})\in\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{standard,+}$, we define a set $\mathcal{L}^{parabolic,+}(\sigma,\tau)$ of sub-blocks satisfying a stronger admissibility condition: Let $(\widehat{\sigma},\widehat{\tau})\in\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{standard,+}$ and $\ell:=\operatorname*{level}(\widehat{\tau})$. Our assumptions imply $\ell=\operatorname*{level}(\widehat{\sigma})$.

Set $L_{\ell}:=\lfloor(0pt(\mathcal{T}_{\mathcal{I}})-\ell)/2\rfloor$ and $L_{\ell}^{\text{middle}}:=0pt(\mathcal{T}_{\mathcal{I}})-L_{\ell}$.

Define the parabolically admissible leaves corresponding to $(\widehat{\sigma},\widehat{\tau})$ by

[TABLE] 5. 5.

Define the set of admissible leaves by

[TABLE]

In order to approximate $\mathbf{K}$, we consider each block $(\widehat{\sigma},\widehat{\tau})\in\mathcal{L}^{standard}_{\mathcal{I}\times\mathcal{J}}$ individually: if it is an inadmissible block, we store $\mathbf{K}|_{\widehat{\sigma}\times\widehat{\tau}}$ directly. If it is an admissible block, we apply the butterfly representation described in the previous section to the sub-clustertrees $\mathcal{T}_{\mathcal{I}}(\sigma)$ and $\mathcal{T}_{\mathcal{J}}(\tau)$. This is equivalent to approximating $\mathbf{K}|_{\widehat{\sigma}\times\widehat{\tau}}$ by a local $\mathcal{DH}^{2}$-matrix.

Remark 1.10**.**

In Section 1.3.3 we argued that a matrix with a butterfly structure can be understood as a ${\mathcal{D}}{\mathcal{H}}^{2}$-matrix. The situation is different here, where only submatrices are endowed with a butterfly structure. While the submatrices are ${\mathcal{D}}{\mathcal{H}}^{2}$-matrices, the global matrix is not a ${\mathcal{D}}{\mathcal{H}}^{2}$-matrix. To see this, let $p:=0pt(\mathcal{T}_{\mathcal{I}})$. If we start with an admissible block $\sigma\times\tau$ (with respect to the standard admissibility condition) on level $\ell$, choose a middle level

[TABLE]

and re-interpolate $p-\ell^{\text{middle}}$ times until we reach leaf clusters $\sigma_{p},\tau_{p}$, we will use points $x_{\sigma_{\hat{\ell}}}$ and $y_{\tau_{\hat{\ell}}}$ on level

[TABLE]

*for the approximation, i.e., the point sets $\mathcal{D}_{\sigma}$ and $\mathcal{D}_{\tau}$ depend on the level $\ell$ of the admissible block. For a $\mathcal{DH}^{2}$-matrix, these sets are only allowed to depend on $\sigma$ and $\tau$, but not on the level $\ell$ of the admissible block.  *

The error analysis of the resulting matrix approximation will require some assumptions. The following assumptions will be useful in Section 3 for the analysis of the Helmholtz kernel (3.1).

Assumption 1.11**.**

Blocks $(\widehat{\sigma},\widehat{\tau})\in\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{standard,+}$ satisfy the admissibility condition (1.19). 2. 2.

For blocks $(\widehat{\sigma},\widehat{\tau})\in\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{standard,+}$ there holds for all $\widehat{\ell}\in\{L_{\ell}^{\text{middle}}-L_{\ell},\ldots,L_{\ell}^{\text{middle}}+L_{\ell}\}$

[TABLE]

for some fixed parameter $\eta_{2}$. 3. 3.

(shrinking condition) There is a constant $\overline{q}\in(0,1)$ such that for all blocks $(\widehat{\sigma},\widehat{\tau})\in\mathcal{L}_{\mathcal{I}\times\mathcal{J}}^{standard,+}$ there holds for all $\sigma\in{\mathcal{T}}_{\mathcal{I}}(\widehat{\sigma})$, $\sigma^{\prime}\in\operatorname*{sons}(\sigma)$ as well as all $\tau\in{\mathcal{T}}_{\mathcal{J}}(\widehat{\tau})$, $\tau^{\prime}\in\operatorname*{sons}(\tau)$

[TABLE]

2 Analysis

A key point of the error analysis is the understanding of the re-interpolation step, which hinges on the following question: Given, on an interval $[-1,1]$, an analytic function that is the product of an analytic function and a polynomial, how well can we approximate it by polynomials on a subinterval $[a,b]\subset[-1,1]$? In turn, sharp estimates for polynomial approximation of analytic functions rely on bounds of the function to be approximated on Bernstein’s elliptic discs. The following Lemma 2.1, which is a refinement of [4, Lemma 5.4], shows that, given $\rho_{1}>1$, it is possible to ensure $\mathcal{E}^{[a,b]}_{\rho_{0}}\subset\mathcal{E}^{[-1,1]}_{\rho_{1}}$ in conjunction with $\rho_{0}>\rho_{1}$:

Lemma 2.1** (Inclusion).**

Fix $\overline{q}\in(0,1)$ and $\underline{\rho}>1$. Then there exists $\widehat{q}\in(0,1)$ such that for any $\rho_{0}\geq\underline{\rho}$ there exists $\rho_{1}\in(1,\widehat{q}\rho_{0}]$ such that for any interval $[a,b]\subset[-1,1]$ with $(b-a)/2\leq\overline{q}$ there holds

[TABLE]

In fact, the smallest $\rho_{1}$ satisfying (2.1) is given by the solution $\rho_{1}>1$ of the quadratic equation (A.6).

Proof.

We remark that [4, Lemma 5.4] represents a simplified version of Lemma 2.1 that is suitable for large $\rho_{0}$; in particular, for $\rho_{0}\rightarrow\infty$, the ratio $\rho_{1}/\rho_{0}$ tends to $(b-a)/2$. The proof of the present general case is relegated to A. ∎

In view of our assumption (1.21), we can fix a “shrinking factor” $\overline{q}\in(0,1)$ in the following. We also assume $\underline{\rho}>1$; the parameter $\widehat{q}\in(0,1)$ appearing in the following results will be as in Lemma 2.1.

We study one step of re-interpolation in the univariate case:

Lemma 2.2**.**

Let $J_{0}=[a_{0},b_{0}]$ and $J_{1}=[a_{1},b_{1}]$ be two intervals with $J_{1}\subset J_{0}$. Let $x_{1}\in J_{1}$. Set $h_{0}:=(b_{0}-a_{0})/2$ and $h_{1}=(b_{1}-a_{1})/2$ and assume

[TABLE]

Let $G$ be holomorphic on ${\mathcal{E}}^{[a_{1},b_{1}]}_{\rho_{1}}$ for some $\rho_{1}\geq\underline{\rho}>1$. Assume that $G|_{J_{1}}$ is real-valued. Let $\kappa\geq 0$. Then there exists $\widehat{q}\in(0,1)$ depending solely on $\overline{q}$ and $\underline{\rho}$ such that

[TABLE]

Hence, for the interpolation error we get

[TABLE]

Proof.

Let $x_{m}:=(b_{1}+a_{1})/2$ denote the midpoint of $J_{1}$. Since $x_{1}\in J_{1}$, we have $|x-x_{1}|\leq|x-x_{m}|+|x_{m}-x_{1}|\leq((\rho_{1}+1/\rho_{1})/2+1)h_{1}$ for any $x\in{\mathcal{E}}^{[a_{1},b_{1}]}_{\rho_{1}}\subseteq B_{(\rho_{1}+1/\rho_{1})/2}(x_{m})$. We estimate (generously) with the abbreviations $M:=\|G\|_{L^{\infty}({\mathcal{E}}^{[a_{1},b_{1}]}_{\rho_{1}})}$ and $R:=(\rho_{1}+1/\rho_{1})/2+1$

[TABLE]

We conclude with a polynomial approximation result (cf. [10, eqn. (8.7) in proof of Thm. 8.1, Chap. 7])

[TABLE]

Let $\rho_{0}$ be the smallest value such that ${\mathcal{E}}^{[a_{1},b_{1}]}_{\rho_{1}}\subset{\mathcal{E}}^{[a_{0},b_{0}]}_{\rho_{0}}$ as given by Lemma 2.1. In particular, we have $\rho_{1}/\rho_{0}\leq\widehat{q}$ with $\widehat{q}$ given by Lemma 2.1. By Bernstein’s estimate ([10, Chap. 4, Thm. 2.2]) we can estimate

[TABLE]

and arrive at

[TABLE]

Recalling $\rho_{1}/\rho_{0}\leq\widehat{q}$ allows us to finish the proof of (2.2). The estimate (2.3) then follows from Lemma 1.1. ∎

Remark 2.3**.**

*The limiting case $\rho_{1}\rightarrow\infty$ corresponds to an entire bounded $G$, which is, by Liouville’s theorem a constant. This particular case is covered in [4, Lemma 5.4].  *

It is convenient to introduce, for $z\in{\mathbb{R}}^{d}$, the function $E_{z}$ and the operator $\widehat{\mathfrak{I}}^{B}_{z}$ by

[TABLE]

The multivariate version of Lemma 2.2 is as follows:

Lemma 2.4**.**

Let $\overline{q}\in(0,1)$ and $\underline{\rho}>1$. Let $B_{1}^{x}:=[\boldsymbol{a}_{1}^{x},\boldsymbol{b}_{1}^{x}]\subset B_{0}^{x}:=[\boldsymbol{a}_{0}^{x},\boldsymbol{b}_{0}^{x}]\subset{\mathbb{R}}^{d}$ and $B_{0}^{y}:=[\boldsymbol{a}_{0}^{y},\boldsymbol{b}_{0}^{y}]\subset B_{-1}^{y}:=[\boldsymbol{a}_{-1}^{y},\boldsymbol{b}_{-1}^{y}]\subset{\mathbb{R}}^{d}$. Let $y_{0}\in B_{0}^{y}$, $y_{-1}\in B_{-1}^{y}$ be given. Assume

[TABLE]

Assume furthermore for a $\rho\geq\underline{\rho}$ and an open set $\Omega\subset{\mathbb{C}}^{2d}$ with ${\mathcal{E}}^{[\boldsymbol{a}^{x}_{0},\boldsymbol{b}^{x}_{0}]}_{{\rho}}\times{\mathcal{E}}^{[\boldsymbol{a}^{y}_{-1},\boldsymbol{b}^{y}_{-1}]}_{{\rho}}\subset\Omega$ that the function $\Phi\in L^{\infty}(\Omega)$ is analytic on $\Omega$. Set

[TABLE]

Assume, for some $\gamma>0$,

[TABLE]

Then there holds for a $\widehat{q}\in(0,1)$ that depends solely on $\overline{q}$ and $\underline{\rho}$

[TABLE]

Proof.

Fix $x_{1}\in B_{1}^{x}$ and $x_{0}\in B_{0}^{x}$. Notice that with the function $R_{x_{1},y_{-1}}$ of Lemma 1.2

[TABLE]

By Lemma 1.2, we have $R_{x_{1},y_{-1}}(x,y_{0})=(x-x_{1})^{\top}G(x,y_{0})(y_{0}-y_{-1})$ for a function $G$ that is analytic on $\Omega$ and satisfies

[TABLE]

Noting that $\exp({\mathbf{i}}\kappa(\Phi(x_{1},y_{0})-\Phi(x_{1},y_{-1})))$ is a constant of modulus $1$, we obtain by combining the univariate interpolation estimate (2.3) of Lemma 2.2 with the multivariate interpolation error estimate (1.6) the bound

[TABLE]

∎

Lemma 2.4 handles one step of a re-interpolation process. The following theorem studies the question of iterated re-interpolation:

Theorem 2.5**.**

Let $\overline{q}\in(0,1)$ and $\underline{\rho}>1$. Let $B_{\ell}^{x}=[\boldsymbol{a}_{\ell}^{x},\boldsymbol{b}_{\ell}^{x}]$, $\ell=0,\ldots,L$, be a nested sequence with $B_{L}^{x}\subset B_{L-1}^{x}\subset\cdots\subset B_{0}^{x}$ satisfying the shrinking condition

[TABLE]

Let $B^{y}=[\boldsymbol{a}^{y},\boldsymbol{b}^{y}]\subset{\mathbb{R}}^{d}$. Assume that $\Phi$ is analytic on the open set $\Omega\subset{\mathbb{C}}^{2d}$ with $\mathcal{E}^{[\boldsymbol{a}_{0}^{x},\boldsymbol{b}_{0}^{x}]}_{\rho}\times\mathcal{E}^{[\boldsymbol{a}^{y},\boldsymbol{b}^{y}]}_{\rho}\subset\Omega$ for some $\rho\geq\underline{\rho}>1$. Define

[TABLE]

Let $(y_{-i})_{i=0}^{L}\subset B^{y}$ be a sequence of points. Assume furthermore

[TABLE]

Abbreviate the operators

[TABLE]

Then, for a $\widehat{q}\in(0,1)$ depending solely on $\overline{q}$ and $\underline{\rho}$ and with the constant

[TABLE]

there holds for $\ell=1,\ldots,L$:

[TABLE]

Proof.

Step 1: By Lemma 2.4, we have the following approximation property for the operators ${\mathfrak{I}}_{\ell}$: with $C_{1}$ given by (2.13) (cf. the definition of $C_{T}$ in (2.9)) we find

[TABLE]

where the last equality follows from the fact that $\Phi(x,z)$ is real for real arguments $x$ and $z$, so we have $|E_{y_{-\ell+1}}|=1$.

Step 2: Observe the telescoping sum

[TABLE]

We claim the following estimates:

[TABLE]

This is proven by induction on $\ell$. For $\ell=1$, the estimate (2.19) expresses (2.17), and (2.20) follows with an additional application of the triangle inequality since ${\mathfrak{I}}_{1}=\operatorname*{I}-\widetilde{E}_{1}$. The case $\ell=0$ is trivial as ${\mathfrak{I}}_{\ell}\circ\cdots\circ{\mathfrak{I}}_{1}$ is understood as the identity. To complete the induction argument, assume that there is an $n\in\mathbb{N}$ such that (2.19), (2.20) hold for all $\ell\in\{0,\ldots,\min\{n,L-1\}\}$. Let $\ell\in\{0,\ldots,\min\{n,L-1\}\}$ and $\pi\in\mathcal{Q}_{m}$. We observe that there is a $\widetilde{\pi}\in\mathcal{Q}_{m}$ such that ${\mathfrak{I}}_{\ell}\circ\cdots\circ{\mathfrak{I}}_{1}(E_{y_{0}}\pi)=E_{y_{-\ell}}\widetilde{\pi}$. The induction hypothesis and (2.20) imply

[TABLE]

Now let $\ell=\min\{n,L-1\}$. We get from (2.18), (2.20), (2.21), and the geometric series

[TABLE]

which is the desired induction step for (2.19). The induction step for (2.20) now follows with the triangle inequality.

Step 3: The estimate (2.19) is the desired estimate (2.14). The bound (2.16) follows from (2.20) and the stability properties of $I^{B^{X}_{0}}_{m}$: For $u\in C(B^{X}_{0})$ we compute

[TABLE]

We are now in a position to prove our main result, namely, an error estimate for the butterfly approximation of the kernel function $k$ given in (1.1):

Theorem 2.6** (Butterfly approximation by interpolation).**

Let $B_{L}^{X}\subset B_{L-1}^{X}\subset\cdots\subset B_{-L}^{X}$ and $B_{L}^{Y}\subset B_{L-1}^{Y}\subset\cdots\subset B_{-L}^{Y}$ be two sequences of the form $B_{\ell}^{X}=[\boldsymbol{a}_{\ell}^{x},\boldsymbol{b}_{\ell}^{x}]$, $B_{\ell}^{Y}=[\boldsymbol{a}_{\ell}^{y},\boldsymbol{b}_{\ell}^{y}]$. Let $(x_{\ell})_{\ell=-L}^{L}$, $(y_{\ell})_{\ell=-L}^{L}$ be two sequences with $x_{\ell}\in B^{X}_{\ell}$ and $y_{\ell}\in B^{Y}_{\ell}$, $\ell=-L,\ldots,L$. Assume:

(analyticity of $\Phi$ and $A$) Let ${\rho}_{\Phi}>1$ and ${\rho}_{A}>1$ be such that $A$ is holomorphic on ${\mathcal{E}}^{[\boldsymbol{a}_{0}^{x},\boldsymbol{b}_{0}^{x}]}_{\rho_{A}}\times{\mathcal{E}}^{[\boldsymbol{a}_{0}^{y},\boldsymbol{b}_{0}^{y}]}_{\rho_{A}}$ and the phase function $\Phi$ is holomorphic on $\Omega\supset{\mathcal{E}}^{[\boldsymbol{a}_{-L}^{x},\boldsymbol{b}_{-L}^{x}]}_{\rho_{\Phi}}\times{\mathcal{E}}^{[\boldsymbol{a}_{-L}^{y},\boldsymbol{b}_{-L}^{y}]}_{\rho_{\Phi}}\supset{\mathcal{E}}^{[\boldsymbol{a}_{0}^{x},\boldsymbol{b}_{0}^{x}]}_{\rho_{A}}\times{\mathcal{E}}^{[\boldsymbol{a}_{0}^{y},\boldsymbol{b}_{0}^{y}]}_{\rho_{A}}$. Define

[TABLE]

• 2.

(shrinking condition) Let $\overline{q}\in(0,1)$ be such that

[TABLE]

• 3.

Let $\gamma>0$ be such that

[TABLE]

• 4.

(polynomial growth of Lebesgue constant) Let $C_{\Lambda}>0$ and $\lambda>0$ be such that the Lebesgue constant of the underlying interpolation process satisfies

[TABLE]

Then:* There exist constants $C$, $b$, $b^{\prime}>0$ depending only on $\gamma$, $\rho_{A}$, $\rho_{\Phi}$, $d$, $C_{\Lambda}$, $\lambda$ such that under the constraint*

[TABLE]

the following approximation result holds for $k$ given by (1.1):

[TABLE]

Proof.

Step 1: We may assume $\rho_{\Phi}\geq\rho_{A}$. It is convenient to abbreviate

[TABLE]

We note

[TABLE]

where the function $R_{x_{0},y_{0}}$ is defined in Lemma 1.2. Using the representation of $R_{x_{0},y_{0}}$ given there, we write

[TABLE]

where the function $G$ is holomorphic on the domain $\Omega$. We estimate

[TABLE]

and get with analogous arguments

[TABLE]

Lemma 1.1 in conjunction with Lemma 1.2 implies together with the univariate polynomial approximation result that led to (2.4)

[TABLE]

Recall the definition of $H$ in (2.29). Since $\Phi$ is real for real arguments, the estimate (2.30) yields

[TABLE]

Step 2: We quantify the effect of $(\widehat{\mathfrak{I}}^{x}_{L}\circ\cdots\circ\widehat{\mathfrak{I}}^{x}_{1})\otimes(\widehat{\mathfrak{I}}^{y}_{L}\circ\cdots\circ\widehat{\mathfrak{I}}^{y}_{1})$. The key is to observe that Theorem 2.5 is applicable since this operator is applied to the function $H$, which is a tensor product of functions of the form suitable for an application of Theorem 2.5. We note with the constants $C_{1}$, $\widehat{q}\in(0,1)$ of Theorem 2.5

[TABLE]

We get, noting the trivial bound $\|k\|_{L^{\infty}(B^{X}_{0}\times B^{Y}_{0})}\leq M_{A}$,

[TABLE]

Step 3: (2.33) is valid for arbitrary $m$ and $L$. We simplify (2.33) by making a further assumption on the relation between $m$ and $L$: The assumption (2.27) on $\Lambda_{m}$ implies that for any chosen $\widetilde{q}\in(\widehat{q},1)$ we have for sufficiently large $m$

[TABLE]

Hence, we obtain for a suitable constant $C>0$ that is independent of $m$

[TABLE]

Using the estimate $\exp(x)-1\leq ex$, which is valid for $x\in[0,1]$, and assuming that $C\widetilde{q}^{m}L\leq 1$ (note that this holds for $m\geq K\log(L+2)$ for sufficiently large $K$), we obtain

[TABLE]

where $b>0$ if we assume that $K$ is selected sufficiently large. Inserting this estimate in (2.33) and noting that $C_{m,L}=2+\widehat{\varepsilon}_{m,L}$ allows us to conclude the proof. ∎

3 Application: the 3D Helmholtz kernel

The case of the 3D Helmholtz kernel

[TABLE]

corresponds to the phase function $\Phi(x,y)=\|x-y\|$ and the amplitude function $A(x,y)=1/(4\pi\|x-y\|)$. We illustrate the butterfly representation for a Galerkin discretization of the single layer operator, i.e.,

[TABLE]

where $\Gamma$ is a bounded surface in ${\mathbb{R}}^{3}$. Given a family of shape functions $(\varphi_{i})_{i=1}^{N}$, the stiffness matrix ${\mathbf{K}}$ is given by

[TABLE]

We place ourselves in the setting of Section 1.3.4 with $\mathcal{I}=\mathcal{J}=\{1,\ldots,N\}$.

Theorem 3.1**.**

Assume the setting of Section 1.3.4 and let Assumption 1.11 be valid. Then there are constants $C$, $b$, $b^{\prime}>0$ that depend solely on the admissibility parameters $\eta_{1}$, $\eta_{2}$, and the parameter $q\in(0,1)$ of Assumption 1.11 such that for the stiffness matrix ${\mathbf{K}}\in{\mathbb{C}}^{\mathcal{I}\times\mathcal{I}}$ given by (3.2) and its approximation $\widetilde{\mathbf{K}}\in{\mathbb{C}}^{\mathcal{I}\times\mathcal{I}}$ that is obtained by the butterfly representation as described in Section 1.3.4 the following holds: If $m\geq b^{\prime}\log(2+0pt{\mathcal{T}}_{\mathcal{I}})$ then

[TABLE]

Proof.

We apply Theorem 2.6 for blocks $(\widehat{\sigma},\widehat{\tau})\in\mathcal{L}_{\mathcal{I}\times\mathcal{I}}^{standard,+}$. To that end, we note that Lemma 3.3 gives us the existence of $\varepsilon>0$ and $\rho>1$ (depending only on the admissibility parameter $\eta_{1}$) such that phase function $(x,y)\mapsto\Phi(x,y)=\|x-y\|$ is holomorphic on

[TABLE]

and satisfies

[TABLE]

Hence, the constants $M_{\Phi}$, $M_{A}$, and $d_{\Omega}$, $\rho_{A}$, $\rho_{\Phi}$ appearing in Theorem 2.6 can be bounded by

[TABLE]

We observe

[TABLE]

so that the conditions (2.26) of Theorem 2.6 are satisfied in view of our Assumption in (1.20). The result now follows from Theorem 2.6. ∎

We conclude this section with a proof of the fact that the Euclidean norm admits a holomorphic extension.

Lemma 3.2**.**

Let $\omega\subset{\mathbb{R}}^{d}$ be open. Define the set

[TABLE]

Then the function

[TABLE]

has an analytic extension to ${\mathcal{C}}_{\omega}$. Furthermore,

[TABLE]

Proof.

The assertion of analyticity will follow from Hartogs’ theorem (cf., e.g., [15, Thm. 2.2.8]), which states that a function that is analytic in each variable separately is in fact analytic. In order to apply Hartogs’ theorem, we ascertain that ${\mathcal{C}}_{\omega}$ is chosen in such a way that

[TABLE]

To see this, abbreviate $D:=\sqrt{2}-1$ and write $(z_{1},\ldots,z_{d})\in{\mathcal{C}}_{\omega}$ in the form $z_{i}=x_{i}+\zeta_{i}$ with $x\in\omega$ and $\zeta_{i}\in{\mathbb{C}}$ with $\sum_{i=1}^{d}|\zeta_{i}|^{2}<D^{2}|x|^{2}$. Then, with Young’s inequality with $\delta:=D=\sqrt{2}-1$:

[TABLE]

Since the square root function is well-defined on the right half plane $\{z\in{\mathbb{C}}\,|\,\operatorname*{Re}z>0\}$, the function $\mathfrak{n}$ is naturally defined on ${\mathcal{C}}_{\omega}$ and analytic in each variable separately. The equalities $|\mathfrak{n}(z_{1},\ldots,z_{d})|=\sqrt{|\sum_{i=1}^{d}z_{i}^{2}|}$ in (3.4) follow from the equation $|\sqrt{z}|=\sqrt{|z|}$ for $z\in{\mathbb{C}}$ with $\operatorname*{Re}z>0$, and the two inequalities in (3.4) are straightforward. ∎

Lemma 3.3**.**

Let $\eta>0$. Then there exist $\varepsilon>0$ and $\rho>1$ depending solely on $\eta$ and the spatial dimension $d$ such that the following is true for any ${[\boldsymbol{a}^{x},\boldsymbol{b}^{x}]}$ and ${[\boldsymbol{a}^{y},\boldsymbol{b}^{y}]}$ satisfying the admissibility condition

[TABLE]

Set $\delta_{B}:=\operatorname*{dist}({[\boldsymbol{a}^{x},\boldsymbol{b}^{x}]},{[\boldsymbol{a}^{y},\boldsymbol{b}^{y}]})$ and define

[TABLE]

Then the function $(x,y)\mapsto\|x-y\|$ has an analytic extension $(x,y)\mapsto\mathfrak{n}(x-y)$ on $\Omega$, and this extension satisfies, for a constant $C>0$ that also depends solely on $\eta$ and $d$,

[TABLE]

Proof.

It is convenient to introduce the abbreviations

[TABLE]

We identify $\operatorname*{Re}\Omega_{x}$ and $\operatorname*{Re}\Omega_{y}$. We start by observing that $\operatorname*{Re}{\mathcal{E}}_{\boldsymbol{\rho}}={\mathcal{E}}_{\boldsymbol{\rho}}\cap{\mathbb{R}}^{d}=B^{[\boldsymbol{a},\boldsymbol{b}]}$ with

[TABLE]

More generally, $\operatorname*{Re}{\mathcal{E}}^{[\boldsymbol{a},\boldsymbol{b}]}_{\boldsymbol{\rho}}={\mathcal{E}}^{[\boldsymbol{a},\boldsymbol{b}]}_{\boldsymbol{\rho}}\cap{\mathbb{R}}^{d}$ is again a box obtained from the box ${[\boldsymbol{a},\boldsymbol{b}]}$ by stretching the $i$-th direction by a factor $1/2(\boldsymbol{\rho}_{i}+1/\boldsymbol{\rho}_{i})$. We now restrict to the case that $\boldsymbol{\rho}_{i}=\rho$ for all $i=1,\ldots,d$. We note that

[TABLE]

Using (3.8) and a triangle inequality, we obtain from (3.5) for $\rho>1$ sufficiently small

[TABLE]

Consider now the set

[TABLE]

and ${\mathcal{C}}_{\omega}$ as defined by (3.3). Note that for $z\in\Omega_{x}$ we have

[TABLE]

with an analogous statement about $\zeta\in\Omega_{y}$. We conclude for $z\in\Omega_{x}$ and $\zeta\in\Omega_{y}$ that the difference

[TABLE]

satisfies

[TABLE]

Hence, $z-\zeta\in{\mathcal{C}}_{\omega}$ provided

[TABLE]

Selecting first $\varepsilon$ sufficient small and then $\rho$ sufficiently close to $1$, this last condition can be ensured. By Lemma 3.2, we conclude the desired analyticity assertion as well as the upper bound (3.6) on $|\mathfrak{n}(x-y)|$. For the lower bound (3.7), we use

[TABLE]

where $C>0$ depends only on $\eta$ and $d$. ∎

Remark 3.4**.**

The condition (1.20) can be met. To illustrate this, assume that the basis functions $(\varphi_{i})_{i=1}^{N}$ all have support of size $O(h)$. For each index $i$, fix a “proxy” $x_{i}\in\operatorname*{supp}\varphi_{i}$, e.g., the barycenter of $\operatorname*{supp}\varphi_{i}$. Consider a tree $\mathcal{T}_{\mathcal{I}}$ for the index set $\mathcal{I}=\{1,\ldots,N\}$ that results from organizing the proxies $(x_{i})_{i=1}^{N}$ in a tree based on a standard octree. In particular, to each cluster $\sigma\in\mathcal{T}_{\mathcal{I}}$ we can associate a box $B^{oct}_{\sigma}$ of the octree such that

[TABLE]

The tree $\mathcal{T}_{\mathcal{I}}$, which was created using the proxies, is also a cluster tree for the shape functions $\varphi_{i}$. The bounding box $B_{\sigma}\supset B^{oct}_{\sigma}$ for $\sigma$ can be chosen close to $B^{oct}_{\sigma}$ in the sense that

[TABLE]

This allows us to show that (1.20) and (1.21) can be met if the leaf size is sufficiently large: For $\sigma$ and $\sigma^{\prime}\in\operatorname*{sons}(\sigma)$ we compute

[TABLE]

This last expression can be made $<1$ if the leaves are not too small, i.e., if the smallest boxes of the octree are large compared to $h$. Let us consider (1.20). We assume that for $\sigma\in{\mathcal{T}}^{L^{\text{middle}}_{\ell}+i}_{\mathcal{I}}(\widehat{\sigma})$ and $\tau\in{\mathcal{T}}^{L^{\text{middle}}_{\ell}-i}_{\mathcal{I}}(\widehat{\tau})$ we have

[TABLE]

Furthermore, we observe $h\lesssim\operatorname*{diam}B_{\widehat{\sigma}}^{oct}$ and $h\lesssim\operatorname*{diam}B_{\widehat{\tau}}^{oct}$ so that we can estimate

[TABLE]

•  *

4 Numerical experiment

In view of the main result of Theorem 2.6, we expect the butterfly approximation $\widetilde{\mathbf{K}}$ to converge exponentially to $\mathbf{K}$ as the degree $m$ of the interpolation polynomials is increased.

In order to get an impression of the convergence properties of the approximation scheme, we apply the butterfly approximation to the discretization (3.2) of the Helmholtz single layer operator. The surface $\Gamma$ is taken to be the polyhdral approximation of the unit sphere $\{x\in{\mathbb{R}}^{3}\ :\ \|x\|_{2}=1\}$ that is obtained by applying regular refinement to the sides of the double pyramid $\{x\in{\mathbb{R}}^{3}\ :\ \|x\|_{1}=1\}$ and projecting the resulting vertices to the sphere. These polyhedra form quasi-uniform meshes. The test and trial spaces consist of piecewise constant functions on these meshes, taking the characteristic functions of the elements are the basis functions $\varphi_{i}$. When forming the stiffness matrix the singular integrals are evaluated by the quadratures described in [12, 24], while regular integrals are evaluated by tensor quadrature in combination with the Duffy transformation.

The cluster tree $\mathcal{T}_{\mathcal{I}}$ is constructed by finding an axis-parallel bounding box containing the entire surface $\Gamma$ and bisecting it simultaneously in all coordinate directions. Due to this construction, clusters can have at most eight sons (empty boxes are discarded) and the diameters of the son boxes are approximately half the diameter of their father. The subdivision algorithm stops on the first level containing a cluster with not more than $32$ indices. The block tree is constructed by the standard admissibility condition using the parameter $\eta_{1}=1$.

The butterfly approximation is constructed by tensor product Chebyshev interpolation. Table 1 lists the spectral errors for $n\in\{32\,768,73\,728,131\,072\}$ triangles with wave numbers $\kappa\in\{16,24,32\}$, corresponding to $\kappa h\approx 0.6$, i.e., approximately ten mesh elements per wavelength. The spectral errors are estimated by applying the power iteration to approximate the spectral radius of the Gramian of the error. We can see that the error reduction factors are quite stable and close to $10$. Table 2 lists the error in the Frobenius norm. The Frobenius error is computed by direct comparison with the exact matrix $\mathbf{K}$. We observe convergence at a rate close to $10$.

Appendix A Proofs of auxiliary results

Proof of Lemma 2.1:

We consider $h\in(0,1]$ and the ellipse ${\mathcal{E}}^{[-1,-1+2h]}_{\rho_{0}}$. We note that once we find $\rho_{1}$ such that ${\mathcal{E}}_{\rho_{1}}\supset{\mathcal{E}}^{[-1,-1+2h]}_{\rho_{0}}$, then by symmetry we also have ${\mathcal{E}}_{\rho_{1}}\supset{\mathcal{E}}^{[1,1-2h]}_{\rho_{0}}$ and then, by convexity of ${\mathcal{E}}_{\rho_{1}}$ also ${\mathcal{E}}_{\rho_{1}}\supset{\mathcal{E}}^{[x-h,x+h]}_{\rho_{0}}$ for any $x\in[-1+h,1-h]$. This justifies our restricting to ${\mathcal{E}}^{[-1,-1+2h]}_{\rho_{0}}$. In Cartesian coordinates, this ellipse is given by

[TABLE]

and $x\in[-1+h-a,-1+h+a]$. The value $\rho_{1}>1$ such that ${\mathcal{E}}^{[-1,-1+2h]}_{\rho_{0}}\subset{\mathcal{E}}_{\rho_{1}}$ has to satisfy

[TABLE]

where

[TABLE]

Claim: The maximum value $\widetilde{M}$ in (A.1) is attained at the left endpoint $x=-1+h-a$ and is given by

[TABLE]

To compute the supremum in (A.1) we introduce the function $f:=s_{1}+s_{2}$ with

[TABLE]

The special structure of the values $a$ and $b$ implies that $s_{2}$ is actually a polynomial:

[TABLE]

We compute $f$ at the endpoints:

[TABLE]

From $h\in(0,1]$ and (A.3a), (A.3b) we obtain

[TABLE]

We are now ready for a further analysis, for which we distinguish the cases that $s_{1}$ is convex or concave. Indeed, only these two cases can occur since the function $s_{1}$ is the square root of a polynomial of degree $2$, and a calculation shows that

[TABLE]

so that $s_{1}^{\prime\prime}$ has a sign. We write $s_{1}(x)=\sqrt{\alpha(x+1-h)^{2}+\beta(x+1-h)+\gamma}$ with

[TABLE]

The case of $s_{1}$ concave: This case is characterized by $4\alpha\gamma-\beta^{2}\leq 0$, i.e.,

[TABLE]

Since $s_{2}$ is affine, the function $f=s_{1}+s_{2}$ is concave. We compute $f^{\prime}(-1+h-a)$:

[TABLE]

In order to see that $f^{\prime}(-1+h-a)\leq 0$, we observe that this last difference is the sum of two terms of opposite sign. For $\xi$, $\eta\geq 0$ we have $-\xi+\eta=(\eta^{2}-\xi^{2})/(\eta+\xi)$ so that the sign of $-\xi+\eta$ is the same as the sign of $\eta^{2}-\xi^{2}$. Hence,

[TABLE]

We conclude that $f$ has its maximum at the left endpoint $x=-1+h-a$.

The case of $s_{1}$ convex: Since $s_{1}$ is convex and $s_{2}$ affine (and thus convex), the function $f$ is convex and therefore attains its maximum at one of the endpoints. We get

[TABLE]

The condition on $\rho_{1}$ is therefore $\widetilde{M}\leq\rho_{1}+1/\rho_{1}$, and hence the smallest possible $\rho_{1}$ is given by the condition

[TABLE]

Note that the right-hand side is $<1$ for every $\rho_{0}>1$ and every $h\in(0,1)$. One can solve for $\rho_{1}$ for given $(\rho_{0},h)$, i.e., solve the quadratic equation (A.6) for $\rho_{1}=\rho_{1}(\rho_{0},h)$. We first observe the asymptotic behavior of $\rho_{1}$: we have $\lim_{\rho_{0}\rightarrow\infty}\rho_{1}(\rho_{0},h)/\rho_{0}=h<1$. One can check the sightly stronger statement that for every $\widehat{q}\in(\overline{q},1)$ there is $\overline{\rho_{0}}>1$ such that $\rho_{1}(\rho_{0},h)/\rho_{0}\leq\widehat{q}$ for all $h\in[0,\overline{q}]\subset[0,1)$ and all $\rho_{0}>\overline{\rho_{0}}$. Hence, we are left with checking the finite range $(\rho_{0},h)\in[\underline{\rho},\overline{\rho_{0}}]\times[0,\overline{q}]$. For that, we note that function $g:x\mapsto x+1/x$ is strictly monotone increasing. Noting $g(\rho_{0})=2r_{x}$, the equation (A.6) takes the form

[TABLE]

from which we get in view of the strict monotonicity of $g$ that $\rho_{1}(\rho_{0},h)<\rho_{0}$ for every $\rho_{0}>1$ and every $h\in[0,1)$. The continuity of the mapping $(\rho_{0},h)\mapsto\rho_{1}$, then implies the desired bound

[TABLE]

$\Box$