<figure><img src="image_1.png" /></figure>

HAL
open science

# Fast discrete convolution in R2 using Sparse Bessel Decomposition

Martin Averseng

<figure><img src="image_2.png" /></figure>

## To cite this version:

Martin Averseng. Fast discrete convolution in R2 using Sparse Bessel Decomposition. 2017. hal-01643914

HAL Id: hal-01643914

https://hal.science/hal-01643914

Preprint submitted on 21 Nov 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Fast discrete convolution in $\mathbb{R}^2$ using Sparse Bessel Decomposition

Martin Averseng

November 21, 2017

## Abstract

We describe an efficient algorithm for computing the matrix vector products that appear in the numerical resolution of boundary integral equations in 2 space dimension. This work is an extension of the so-called Sparse Cardinal Sine Decomposition algorithm by Alouges et al., which is restricted to three-dimensional setups. Although the approach is similar, significant differences appear throughout the analysis of the method. Bessel decomposition, in particular, yield longer series for the same accuracy. We propose a careful study of the method that leads to a precise estimation of the complexity in terms of the number of points and chosen accuracy. We also provide numerical tests to demonstrate the efficiency of this approach. We give the compression performance for a $N \times N$ linear system for several values $N$ up to $10^7$ and report the computation time for the off-line and on-line parts of our algorithm. We also include a toy application to sound canceling to further illustrate the efficiency of our method.

## Introduction

The Boundary Element Method requires the resolution of fully populated linear systems $Au = b$. As the number of unknowns gets large, the storage of the matrix $A$ and the computational cost for solving the system through direct methods (e.g. LU factorization) become prohibitive. Instead, iterative methods can be used, which require very fast evaluations of matrix vector products. In the context of boundary integral formulations, this takes the form of discrete convolutions:

$$q_k = \sum_{l=1}^{N_z} G(z_k - z_l) f_l, \quad k \in \{1, \dots, N_z\}. \qquad (1)$$

Here, $G$ is the Green's kernel of the partial differential equation under consideration, $z = (z_k)_{1 \le k \le N_z}$ is a set of points in $\mathbb{R}^2$, with diameter $\delta_{\max}$, $|\cdot|$ is the Euclidian norm and $f = (f_k)_{1 \le k \le N_z}$ is a vector (typically the values of a function at the points $z_k$). For example, the resolution of the Laplace equation with Dirichlet boundary conditions leads to (1) with $G(x) = -\frac{1}{2\pi} \log|x|$ (the kernel of the single layer potential).

In principle, the effective computation of the $(q_k)_{1 \le k \le N_z}$ using (1) requires $O(N_z^2)$ operations. However, several more efficient algorithms have emerged to
compute an approximation of (1) with only quasilinear complexity in $N_z$. Among those are the celebrated Fast Multipole Method (see for example [7, 8, 13, 19, 20] and references therein), the Hierarchical Matrix [6], and more recently, the Sparse Cardinal Sine Decomposition (SCSD) [2].

One of the key ingredients in all those methods consists in writing the following local variable separation:

$$G(z_k - z_l) \approx \sum_j \lambda_j G_1^j(z_k) G_2^j(z_l),$$

which needs to be valid for $z_k$ and $z_l$ arbitrarily distant from each other, and up to a controlled accuracy. This eventually results in compressed matrix representations and accelerated matrix-vector product. Notice that, to be fully effective, the former separation is usually made locally with the help of a geometrical splitting of the cloud of points $z$ using a hierarchical octree.

Here we present an alternative compression and acceleration technique, which we call the Sparse Bessel Decomposition (SBD). This is an extension of the SCSD adapted to 2-dimensional problems. The SBD and SCSD methods achieve performances comparable to the aforementioned algorithms, they are flexible with respect to the kernel $G$, and do not rely on the construction of octrees, which makes them easier to implement. In addition, they express in an elegant way the intuition according to which a discrete convolution is nothing but a product of Fourier spectra.

The method heavily relies on the Non Uniform Fast Fourier Transform (NUFFT) of type III (see the seminal paper [12] and also [14, 17] and references therein for numerical aspects and open source codes). The NUFFT is a fast algorithm, which we denote by $\text{NUFFT}_{\pm}[z, \xi](\alpha)$ for later use, that returns, for arbitrary sets of $N_z$ points $z$ and $N_\xi$ points $\xi$ in $\mathbb{R}^2$ and a complex vector $\alpha \in \mathbb{C}^{N_\xi}$, the vector $q \in \mathbb{C}^{N_\xi}$ defined by:

$$q_v = \sum_{k=1}^{N_\xi} e^{\pm iz_k \cdot \xi_v} \alpha_k, \quad v \in \{1, \dots, N_\xi\}.$$

This algorithm generalizes the classical Fast Fourier Transform (FFT) [10], to nonequispaced data, preserving the quasi-linear complexity in $N_{z,\xi} := \max(N_z, N_\xi)$. The SBD method first produces a discrete and sparse approximation of the spectrum of $G$,

$$G(x) \approx G_{\text{approx}}(x) := \sum_{v=1}^{N_\xi} e^{ix \cdot \xi_v} \hat{w}_v, \quad |x| \le \delta_{\max}. \qquad (2)$$

This approximation is replaced in (1) to yield

$$q_k \approx \left( \sum_{v=1}^{N_\xi} e^{+iz_k \cdot \xi_v} \left[ \hat{w}_v \sum_{l=1}^{N_z} e^{-iz_l \cdot \xi_v} f_l \right] \right)_{1 \le k \le N_z} \qquad (3) \\ = \text{NUFFT}_+[z, \xi] (\hat{\omega} \odot \text{NUFFT}_-[z, \xi](f)).$$

where $\odot$ denotes the elementwise product between vectors. The decomposition (2) is obtained "off-line" and depends on the value of $\delta_{\max}$, but is independent of the choice of the vector $f$ and can thus be used for any evaluation of the matrix vector product (1). The approximation (3) reduces the complexity from $O(N_z^2)$ to $O(N_z \log(N_z, \xi))$, stemming from the NUFFT complexity.
The NUFFT has already been used in the literature for the fast evaluation of quantities of the form (1). In particular, our algorithm shares many ideas with the approach in [18]. The method presented therein also relies on an approximation of the form (2). However, we choose the set of frequencies $\xi$ in a different way, leading to a sparser representation of G (see Remark 2).

The kernel G is usually singular near the origin, in which case (2) can only be accurate for $|z|$ above some threshold $\delta_{min}$. Part of the SBD algorithm is thus dedicated to computing a local correction using a sparse matrix (which we call the close correction matrix, denoted by D in the sequel) to account for the closer interactions. This threshold must be chosen so as to balance the time spent for computing the far (3) and those close contributions. As a matter of fact, we shall prove the following:

**Theorem 1.** Assume the points z are uniformly distributed on a regular curve, and $G(x) = \log|x|$. Let $\varepsilon > 0$ the desired accuracy of the method, and assume $N_z > |\log \varepsilon|$. Fix

$$a = \frac{|\log \varepsilon|^{2/3}}{N_z^{2/3-\alpha}}$$

for some $\alpha \in [0, \frac{1}{6}]$, and choose

$$\delta_{\min} = a\delta_{\max}.$$

Then there exists a constant C > 0 independent of $N_z$, $\varepsilon$ and $\alpha$ such that:

(i) The number of operations required for the computation of the representation (2) valid for $|x| \in [\delta_{\min}, \delta_{\max}]$ is bounded by

$$C_{\text{off}}(N_z, \varepsilon, \alpha) \leq C |\log \varepsilon| N_z^{2-3\alpha};$$

(ii) The number of operations required for the assembling of the close correction matrix D is bounded by

$$C_{\text{assemble}}(N_z, \varepsilon, \alpha) \leq C |\log \varepsilon|^{2/3} C_{\text{NUFFT}}(\varepsilon) N_z^{4/3+\alpha} \log(N_z);$$

(iii) Once these two steps have been completed, (1) can be evaluated for any choice of vector f at a precision at least $\varepsilon \sum_l |f_l|$ in a number of operations bounded by

$$C_{\text{on}}(N_z, \varepsilon, \alpha) \leq C |\log \varepsilon|^{2/3} \left( N_z^{4/3+\alpha} + C_{\text{NUFFT}}(\varepsilon) N_z^{4/3-2\alpha} \log N_z \right).$$

We prove Theorem 1 using several steps. We first introduce the Fourier-Bessel series (section 2). In section 3, we present the Sparse Bessel Decomposition, and analyze its numerical stability. In section 4, we apply the SBD method to the Laplace kernel, and give estimates for the number of terms to reach a fixed accuracy. We also explain how to adapt this method for other kernels. In section 5, we show how to convert a SBD decomposition into an approximation of the form (2) through what we call "circular quadratures" and provide error estimates. In section 6, we summarize the complexities of each step and prove Theorem 1. We conclude with some numerical examples.

Before this, let us start by a brief overview of our algorithm.
# 1 Summary of the algorithm

The SBD algorithm can be summarized as follows:

## Off-line part

**Inputs:** A radially symmetric kernel $G$, a set of $N_z$ points $z$ in $\mathbb{R}^2$ of diameter $\delta_{\max}$, a value for the parameter $\delta_{\min}$, a tolerance $\varepsilon > 0$.

**Sparse spectrum sampling:** Compute a set of $N_\xi$ complex weights $\hat{\omega}$ and $N_\xi$ frequencies $\xi$ so that (2) is valid for $\delta_{\min} \le |x| \le \delta_{\max}$ up to tolerance $\varepsilon$.

**Correction Matrix:** Determine the set $\mathcal{P}$ of all pairs $(k, l)$ such that $|z_k - z_l| \le \delta_{\min}$ (fixed-radius neighbor search). Assemble the close correction sparse matrix:

$$D_{kl} = \delta_{(k,l)} \in \mathcal{P} \left( G(z_k - z_l) - \sum_{v=1}^{N_\xi} e^{i(z_k - z_l) \cdot \xi_v} \hat{\omega}_v \right). \quad (1.1)$$

Notice that the second term is a non-uniform discrete Fourier transform. Indeed, if we introduce $(y_p)_{p \in \mathcal{P}}$ given by $y_{(k,l)} = z_k - z_l$, and $Y = \text{NUFFT}_+[y, \xi](\hat{\omega})$, the non-zero entries of $D$ are given by

$$D_{k,l} = G(y_{(k,l)}) - Y_{(k,l)}.$$

**Outputs:** The set of weights $\hat{\omega}$, the frequencies $\xi$ and the sparse matrix $D$.

## On-line part

**Input:** All outputs of the off-line part, and a complex vector $f$ of size $N_z$.

**Far approximation:** Compute, for all $k$,

$$q^{\text{far}} = \sum_{l=1}^{N_r} G_{\text{approx}}(z_k - z_l) f_l. \quad (1.2)$$

For this, follow three steps

(i) **Space → Fourier:** Compute $\hat{f} = \text{NUFFT}_-[z, \xi](f)$.

(ii) **Fourier multiply:** Perform elementwise multiplication by $\hat{\omega}$:

$$\hat{g}_v = \hat{\omega}_v \hat{f}_v.$$

(iii) **Fourier → Space:** Compute $q^{\text{far}} = \text{NUFFT}_+[z, \xi](\hat{g})$.

**Close correction:** Compute the sparse matrix product:

$$q^{\text{close}} = Df.$$

**Output:** The vector $q = q^{\text{far}} + q^{\text{close}}$, with, for any $k \in \{1, \dots, N_z\}$,

$$\left| q_k - \sum_{l=1}^{N_r} G(z_k - z_l) f_l \right| \le \varepsilon \sum_{l=1}^{N_r} |f_l|.$$

## The sparse spectrum sampling step:

The main novelty in our algorithm is the method for producing an approximation of the form (2). We proceed in two steps, uncoupling the radial and circular coordinates.
**Sparse Dessel Decomposition:** In a first part, we compute a Bessel series approximating $G(r)$ on $[\delta_{\min}, \delta_{\max}]$, which we call a Sparse Bessel Decomposition. It takes the form

$$G(r) \approx G(\delta_{\max}) + \sum_{p=1}^{P} \alpha_p J_0(\delta_{\max} \rho_p r), \quad r \in [\delta_{\min}, \delta_{\max}],$$

where $J_0$ is the Bessel function of first kind and order zero and $(\rho_p)_{p \in \mathbb{N}^*}$ is the sequence of its roots (see Definition 1 for more details). In order to compute the coefficients $(\alpha_p)_{1 \le p \le P}$, we first choose a starting value for $P$ and compute the weights $\alpha_1, \dots, \alpha_p$ that minimize the least square error

$$\int_{\delta_{\min} \le |x| \le \delta_{\max}} \left| \nabla \left( G(x) - \sum_{p=1}^{P} \alpha_p J_0(\rho_p \delta_{\max} |x|) \right) \right|^2 dx.$$

This amounts to solving an explicit linear system. We keep increasing $P$ until the residual error goes below the required tolerance. To choose the stopping criterion, we suggest monitoring the $L^\infty$ error near $\delta_{\min}$ where it is usually the highest. The successive choices of $P$ are made using a dichotomy search.

**Circular quadrature:** In a second step, we use approximations of the form:

$$J_0(\rho_p |x|) \approx \frac{1}{M_p} \sum_{m=0}^{M_p-1} e^{i \rho_p \xi_p^m \cdot x}, \quad p \in \{1, \dots, P\},$$

which are discrete versions of the identity:

$$J_0(\rho_p |x|) = \frac{1}{2\pi} \int_{|\xi|=1} e^{i\rho_p x \cdot \xi} d\sigma(\xi).$$

We sum them to eventually obtain the formula (2).

**Remark 1.** In the SCSD method [2], the Bessel functions are replaced by cardinal sine functions, since, for $x \in \mathbb{R}^3$

$$\frac{1}{4\pi} \int_{|\xi|=1} e^{ix \cdot \xi} d\sigma(\xi) = \text{sinc}(|x|),$$

where the integral is now taken over $\mathbb{S}^2 \subset \mathbb{R}^3$.

## 2 Series of Bessel functions and error estimates

In this section, we give a short introduction to Fourier-Bessel series. A possible reference on this topic is [24], chapter XVIII. The main result needed for our purpose is Proposition 1, an equivalent statement of which can be found in to Theorem 1 in [22] chapter 8, section 20.

In subsection 2.1, we quickly recall some classical facts on the Laplacian eigenfunctions with Dirichlet boundary conditions on the unit ball and formulate a conjecture (Conjecture 1) for the normalization constant supported by strong numerical evidence (Figure 1). It is worth noting that the use Laplace eigenfunctions as the decomposition basis and the results we obtain hereafter do not rely on the
space dimension. For example, in $\mathbb{R}^3$ the radial eigenvalues of the Laplacian are proportional to

$$x \mapsto \frac{J_{1/2}(2\pi p |x|)}{|x|^{1/2}}, \quad p \in \mathbb{N}^*.$$

Therefore, our approach generalizes [2] to any dimension.

## 2.1 Radial eigenvalues of the Laplace operator with Dirichlet conditions

In the following, $B$ denotes the unit ball in $\mathbb{R}^2$, and $\mathscr{C}$ its boundary. We say that a function $u: \mathbb{R}^2 \to \mathbb{R}$ is radial if there exists $\tilde{u}: \mathbb{R}^+ \to \mathbb{R}$ such that for any $x \in \mathbb{R}^2$,

$$u(x) = \tilde{u}(|x|).$$

In this case, we use the notation $u(x) = u(r)$ where $|x| = r$. We note $L_{rad}^2(B)$ the closed subspace of $L^2$ that are radial functions and $C_{c,rad}^\infty(B) = \{\varphi \in C_c^\infty(B) | \varphi \text{ is radial}\}$. Similarly,

$$H_{rad}^1(B) = \{ u \in L_{rad}^2(B) \mid \forall 1 \le i \le 2, \frac{\partial u}{\partial x_i} \in L^2(B) \},$$

which is a Hilbert space for the norm

$$\|u\|_{H_{rad}^1(B)}^2 = \int_B |u|^2 + |\nabla u|^2 = 2\pi \int_0^1 r (u(r)^2 + u'(r)^2) dr.$$

Finally, we note $H_{0,rad}^1(B)$ the closure of $C_{c,rad}^\infty(B)$ in $H_{rad}^1(B)$, with the norm

$$\|u\|_{H_{0,rad}^1(B)}^2 = 2\pi \int_0^1 r u'(r)^2 dr.$$

We now briefly recall some facts on Bessel functions. All the results that we use on this topic can be found in the comprehensive book [1], most of which can be consulted handily on the digital library [16].

**Definition 1.** The Bessel function of the first kind and order $\alpha$, $J_\alpha$ is defined by the following series:

$$J_{\alpha}(r) := \sum_{m=0}^{\infty} \frac{(-1)^m}{m! \Gamma(m+1+\alpha)} \left(\frac{r}{2}\right)^{2m+\alpha} \quad (2.1)$$

$J_0$ is a $C^\infty$ solution of Bessel's differential equation

$$r^2 f''(r) + r f'(r) + r^2 f(r) = 0. \quad (2.2)$$

The roots $(\rho_p)_{p \in \mathbb{N}^*}$ of $J_0$, behave, for large $p$, as

$$\rho_p \sim \pi p.$$

More precisely, $\rho_p = \pi p - \frac{\pi}{4} + O(\frac{1}{p})$, and

$$\pi(p - 1/4) \le \rho_p \le \pi(p - 1/8). \quad (2.3)$$
For any $p \in \mathbb{N}^*$, we introduce:

$$e_p(x) = C_p J_0(\rho_p |x|),$$

where the normalization constant $C_p$ is chosen such that $\|e_p\|_{H_{\text{rad}}^1(B)} = 1$, that is

$$C_p = \frac{1}{(2\pi \int_B r \rho_p^2 J_1(\rho_p r)^2)^{1/2}} = \frac{1}{\sqrt{\pi \rho_p |J_1(\rho_p)|}}.$$

For any $p \in \mathbb{N}^*$, $e_p$ satisfies:

$$-\Delta e_p = \rho_p^2 e_p. \qquad (2.4)$$

The following result is well-known.

**Theorem 2.** The family $\{e_p, p \in \mathbb{N}^*\}$ is a Hilbert basis of $H_{0,\text{rad}}^1(B)$.

One can check, using asymptotic expansions of Bessel functions, that

$$C_p = \frac{1}{\sqrt{2\pi p}} + O\left(\frac{1}{p^{3/2}}\right), \qquad (2.5)$$

We will actually need a more precise knowledge on the constant $C_p$. We formulate the next conjecture, which seems hard to establish and for which we didn't find any element of proof in the literature. Numerical evidence exposed in Figure 1 strongly suggests it is true.

**Conjecture 1.** For all $p \in \mathbb{N}^*$:

$$\frac{1}{\sqrt{2\pi p}} \le C_p \le \frac{1}{\sqrt{2\pi(p-1/4)}}. \qquad (2.6)$$

<figure><img src="image_3.png" /><figcaption>Figure 1: Illustration of Conjecture 1 with numerical values of the first 1000 terms of the sequence $v_p = \sqrt{2\pi p}C_p - 1$ (blue circles) and $w_p = 1 - \sqrt{2\pi(p-1/4)}C_p$ (red crosses) in log-log scale</figcaption></figure>
## 2.2 Truncature error for Fourier-Bessel series of smooth functions

We now introduce the Fourier-Bessel series and prove a bound for the norm of the remainder. In Theorem 2, we have shown that any function $f \in H_{0,rad}^1(B)$ can be expanded through its so-called Fourier-Bessel series as

$$f = \sum_{p \in \mathbb{N}^*} c_p(f) e_p.$$

The generalized Fourier coefficients are obtained by the orthonormal projection:

$$c_p(f) = \int_B \nabla f(x) \cdot \nabla e_p(x) dx = \rho_k^2 \int_B f(x) e_p(x) dx.$$

Most references on this topic focus on proving pointwise convergence of the series even for not very regular functions $f$ (e.g. piecewise continuous, square integrable, etc.) [3, 9, 15, 21]. In such cases, the Fourier-Bessel series may exhibit a Gibb's phenomenon [26]. On the contrary here, we need to establish that the Fourier-Bessel series of very smooth functions converges exponentially fast. To this aim, we first introduce the following terminology:

**Definition 2.** We say that a radial function $f$ satisfies the multi-Dirichlet condition of order $n \in \mathbb{N}^*$ if $f$ is $H^{2n}$ in a neighborhood of $\mathcal{C}$ and if for all $s \le n-1$, the $s$-th iterate of $-\Delta$ on $f$, denoted by $(-\Delta)^s f$, vanishes on $\mathcal{C}$ (with the convention $(-\Delta)^0 f = f$).

**Proposition 1.** If $f \in H^{2n}(B)$ satisfies the multi-Dirichlet condition of order $n$, then for any $p \in \mathbb{N}^*$:

$$c_p(f) = \frac{1}{\rho_p^{2n-2}} \int_B (-\Delta)^n f(x) e_p(x) dx.$$

*Proof.* If $f$ satisfies the multi-Dirichlet condition of order $n=1$, then by integration by parts:

$$c_p(f) = \int_B (-\Delta) f(x) e_p(x) dx,$$

since $e_p$ vanishes on $\mathcal{C}$. Assume the result is true for some $n \ge 1$ and let $f$ satisfy the multi-Dirichlet condition of order $n+1$. Then, using the fact that $e_p$ is an eigenvector of $-\Delta$ associated to the eigenvalue $\rho_p^2$ we obtain:

$$c_p(f) = \frac{1}{\rho_p^{2n}} \int_B (-\Delta)^n f(x) (-\Delta) e_p(x) dx.$$

The result follows from integration by parts where we successively use that $(-\Delta)^n f$ and $e_p$ vanish on $\mathcal{C}$.

□

**Corollary 1.** If $f \in H^{2n}(B)$ satisfies the multi-Dirichlet condition of order $n$, there exists a constant $C$ independent of the function $f$ such that for all $p \in \mathbb{N}^*$,

$$|c_p(f)| \le C \frac{\|(-\Delta)^n f\|_{L_{\text{rad}}^2(B)}}{(\pi p)^{2n-1}}.$$

Notice that this is similar to the fact that the Fourier coefficients of smooth functions decay fast.
*Proof.* We apply the result of the previous proposition and remark that since $e_p$ is an eigenfunction of the Laplace operator with unit norm in $H_{0,rad}^1(B)$, $\|e_p\|_{L_{rad}^2(B)} = \frac{1}{\rho_p} \|e_p\|_{H_{0,rad}^1(B)} = \frac{1}{\rho_p}$. To conclude, recall $\rho_p \sim p\pi$ for large $p$. □

**Corollary 2.** Let the remainder be defined as

$$R_P(f) = \sum_{p=P+1}^{+\infty} c_p(f)e_p.$$

If $f \in H^{2n}(B)$ satisfies the multi-Dirichlet condition of order $n$, there exists a constant $C$ independent of $n$ and $P$ such that:

$$\|R_P(f)\|_{H_{0,rad}^1(B)} \le C \frac{\|(-\Delta)^n f\|_{L_{rad}^2(B)}}{(\pi P)^{2n}} \sqrt{\frac{P^3}{n}}.$$

*Proof.* Parseval's identity implies

$$\|R_P(f)\|_{H_{0,rad}^1(B)}^2 = \sum_{p=P+1}^{+\infty} |c_p(f)|^2.$$

According to the previous results, we find that:

$$\|R_P(f)\|_{H_{0,rad}^1(B)} \le C \|(-\Delta)^n f\|_{L_{rad}^2(B)} \sqrt{\sum_{p=P+1}^{+\infty} \frac{1}{(\pi p)^{4n-2}}}. $$

The announced result follows from $\sum_{p>P} \frac{1}{p^{\alpha}} \propto \frac{1}{(\alpha-1)P^{\alpha-1}}$ for $\alpha > 1$. □

## 2.3 Other boundary conditions

When we replace the Dirichlet boundary condition by the following Robin boundary conditions

$$\frac{\partial u}{\partial n} + Hu = 0 \qquad (2.7)$$

for some constant $H \ge 0$, the same analysis can be conducted, leading to Dini series (also covered in [24]). This time, we construct a Hilbert basis of $H_{rad}^1(B)$ with respect to the bilinear form

$$a_H(u,v) := \int_B \nabla u(x) \cdot \nabla v(x) dx + H \int_C u(x)v(x) d\sigma(x).$$

The following result holds.

**Theorem 3.** Let $(\rho_p^H)_{p \in \mathbb{N}_*}$ the sequence of positive solutions of

$$rJ_0'(r) + HJ_0(r) = 0.$$
(i) If $H > 0$, the functions

$$e_p^H(r) = C_p J_0(\rho_p^H r),$$

with $C_p$ such that $a_H(e_p^H, e_p^H) = 1$, form a Hilbert basis of $H_{\text{rad}}^1(B)$.

(ii) If $H = 0$, a constant function must be added to the previous family to form a complete set.

It can be checked that the truncature error estimates in Corollary 2 extend to this case, for functions satisfying multi-Robin conditions of order $n \ge 1$, that is for all $s \le n - 1, (-\Delta)^s u$ satisfies (2.7).

# 3 Sparse Bessel Decomposition

## 3.1 Definition of the SBD

Consider the kernel $G$ involved in (1). We can assume up to rescaling $G$ that the diameter $\delta_{\max}$ of $z$ is bounded by 1, and therefore, we need to approximate $G$ only on the unit ball $B$. If we wish to approximate $G$ in series of Bessel functions, two kinds of complications are encountered:

(i) $G$ is usually singular near the origin, therefore not in $H^{2n}(B)$ (even for $n = 1$).

(ii) The multi-Dirichlet conditions may not be fulfilled up to a sufficient order.

The point (ii) is crucial in order to apply the error estimates of the previous section. The first two kernels that we will study (Laplace and Helmholtz) satisfy the favorable property:

$$\Delta G = \lambda G$$

for some $\lambda \in \mathbb{C}$, which will be helpful to ensure (ii) at any order. For more general kernel, we propose in subsection 4.3 a simple trick to enforce multi-Dirichlet conditions up to a given order.

As for point (i), we will use the following method: for the approximation

$$G \approx \sum_{p=1}^{P} \alpha_p e_p,$$

we know by Theorem 2 that the minimal $H_0^1$ error on $B$ is reached for $\alpha_p = c_p(G)$. However, if the closest interaction in (1) are computed explicitly, it can be sufficient to approximate $G$ in a domain of the form $a \le r \le 1$ for some $a$. For this reason, we propose to define the coefficients $(\alpha_1, \dots, \alpha_P)$ as the minimizers of the quadratic form

$$Q^P(t_1, t_2, \dots, t_P) = \int_{\mathcal{S}(a)} \left| \nabla \left( G(x) - \sum_{p=1}^{P} t_p J_0(\rho_p |x|) \right) \right|^2 dx,$$

where $\mathcal{S}(a)$ is the annulus $\{x \in \mathbb{R}^2 \mid a < |x| < 1\}$. In the sequel, those coefficients will be called the SBD coefficients of $G$ of order $P$. Obviously, for any radial function $\tilde{G}$ defined on $B$ that coincides with $G$ on $\mathcal{S}(a)$, one has

$$Q^P(\alpha_1, \dots, \alpha_P) \le \int_B \left| \nabla \left( \tilde{G}(x) - \sum_{p=1}^P c_p(\tilde{G}) J_0(\rho_p |x|) \right) \right|^2 dx.$$
In particular, when $\tilde{G}$ is smooth up to the origin, this gives an error estimate via Corollary 2. If we choose a sufficiently high value for $a$, we ensure that smooth enough extensions exist, and ensure fast decay of the coefficients.

**Remark 2.** The SBD weights do not depend on any specific extension of $G$ outside the annulus. Therefore, they provide the sparsest approximation one can expect, contrary to the usual approach where an explicit regularization $\tilde{G}$ of the kernel is constructed and the coefficients $c_p(\tilde{G})$ are used (see, for example [18]).

The next result shows that the $H_0^1$ norm on $\mathcal{A}(a)$ controls the $L^\infty$ norm, thus ruling out any risk of Gibb's phenomenon.

**Lemma 1.** Let $a \in (0, 1)$ and $e \in H_{\text{rad}}^1(\mathcal{A}(a))$ that vanishes on $\mathcal{C}$. Then $e$ coincides almost everywhere with a continuous function with

$$|e(x_0)| \le \sqrt{\frac{-\log|x_0|}{2\pi}} \sqrt{\int_{\mathcal{A}(a)} |\nabla e(x)|^2 dx}, \quad \text{almost for all } x_0 \in \mathcal{A}(a).$$

*Proof.* It is sufficient to show the inequality for smooth $e$, the general result following by density. Let $x_0 \in \mathcal{A}(a)$, we have, since $e(1) = 0$:

$$|e(x_0)| \le \int_{|x_0|}^1 |e'(r)| dr, \quad (3.1)$$

$$\le \sqrt{2\pi \int_{|x_0|}^1 r |e'(r)|^2 dr} \sqrt{\int_{|x_0|}^1 \frac{1}{2\pi r} dr}. \quad (3.2)$$

□

## 3.2 Numerical computation of the SBD

For a given kernel $G$, the SBD coefficients $\alpha_p$ are obtained numerically by solving the following linear system:

$$\begin{align} & \sum_{q=1}^{P} \left( \int_{\mathcal{A}(a)} \rho_q J_1(\rho_p|x|) J_1(\rho_q|x|) dx \right) \alpha_q \tag{3.3} \\ & = - \int_{\mathcal{A}(a)} G'(x) J_1(\rho_p|x|) dx, \quad 1 \le p \le P, \nonumber \end{align}$$

Where $J_1$ is the Bessel function of first kind and order 1 (in fact, $J'_0 = -J_1$). We solve this system for increasing values of $P$ until a required tolerance is reached. It turns out that the matrix $A^P$ whose entries are given by

$$A_{k,l}^P = \int_{\mathcal{A}(a)} \nabla e_k \cdot \nabla e_l, \quad k, l \in \{1, \dots, P\},$$

is explicit: for $(i, j) \in \{1, \dots, P\}^2$, the non-diagonal entries of $A^P$ are

$$A_{i,j} = \frac{2\pi C_i C_j \rho_i \rho_j}{\rho_j^2 - \rho_i^2} \left[ F_{i,j}(1) - F_{j,i}(1) - F_{i,j}(a) + F_{j,i}(a) \right],$$

where

$$F_{i,j}(r) = \rho_i r J_0(\rho_i r) J_0'(\rho_j r),$$
while the diagonal entries are

$$A_{i,i} = 2\pi C_i^2 (F_i(1) - F_i(a)),$$

where

$$F_i(r) = \rho_i^2 r^2 \left[ \frac{1}{2} J_0(\rho_i r)^2 + \frac{1}{2} J_0'(\rho_i r)^2 \right] + \rho_i r J_0(\rho_i r) J_0'(\rho_i r).$$

Those formulas are obtained using Green's formulas together with the fact that $e_k$ are eigenfunctions of the Laplace operator. They are valid for any value of $\rho_k$ (not just the roots of $J_0$).

### 3.3 Conditioning of the linear system

The conditioning of A seems to depend almost exclusively on the parameter $\gamma := Pa$. We were only able to derive an accurate estimate of the conditioning of A when $\gamma$ is small enough. For large $\gamma$, we will show some numerical evidence for a conjectured bound on the conditioning of A.

#### Conditioning of A for small $\gamma$

**Theorem 4.** If Conjecture 1 is true, then, for $b = 1$, the eigenvalues of A lie in the interval $[F(\gamma) - \frac{\pi^4}{144} \frac{\gamma^4}{P}, 1]$ where

$$F(\gamma) = 1 - \int_{0}^{\pi\gamma} \frac{t}{2} (J_{1}(t)^{2} - J_{0}(t)J_{2}(t)) dt.$$

This estimate is only useful when $F(\gamma) > 0$, that is $\gamma < \gamma^*$ where $\gamma^*$ is the first positive root of $F(z)$. One has

$$\gamma^* \approx 1.471.$$

In particular, for $\gamma = 1$, The matrix A is well conditioned, the ratio of its largest to its smallest eigenvalues being of the order $F(1)^{-1} < 2$. A plot of F is provided below, Figure 2, and some numerical approximations of the minimal eigenvalue of A are shown in function of $\gamma$ for several values of P.

*Proof.* Let $u \in \operatorname{span}\{e_1, e_2, \dots, e_P\}$ and let $\alpha$ its coordinates on this basis. Then

$$\alpha^T A \alpha = \int_{\mathcal{S}(a)} \nabla u^2 < \int_B \nabla u^2 = \|\alpha\|_2^2,$$

showing that the eigenvalues of A are bounded by 1. Thus $I - A$ is a positive symmetric matrix and the smallest eigenvalue of A is bounded by

$$\lambda_{\min}(A) \geq 1 - \operatorname{tr}(I - A),$$

which yields:

$$\lambda_{\min}(A) \geq 1 - 2\pi \sum_{p=1}^{P} C_p^2 \int_{0}^{\rho_p a} u J_1(u)^2 du.$$

We now use Conjecture 1, which implies

$$2\pi C_p^2 \leq \frac{1}{p} + \frac{1}{3p^2},$$
<figure><img src="image_4.png" /><figcaption>Figure 2: Graph of F and numerical values of the minimal eigenvalue of A, $\lambda_{\min}(A)$ in function of $\gamma$, in the case $b = 1$, for $P = 50$ (blue circles) and $P = 500$ (red crosses)</figcaption></figure>

combined with (2.3), to get

$$ \lambda_{\min}(A) \geq 1 - \int_{0}^{P} \frac{dt}{t} \int_{0}^{\pi t a} u J_{1}(u)^{2} du - \frac{1}{3} \sum_{p=1}^{P} \frac{1}{p^{2}} \int_{0}^{\rho_{p} a} u J_{1}(u)^{2} du. $$

For the first term, we write

$$ \int_{0}^{t} u J_{1}(u)^{2} = \frac{t^{2}}{2} (J_{1}^{2}(t) - J_{0}(t)J_{2}(t)), $$

while for the second, we use the classical inequality $J_1(u) \le \frac{u}{2}$ to deduce

$$ \lambda_{\min}(A) \geq 1 - \int_{0}^{\pi\gamma} \frac{t}{2} (J_{1}(t)^{2} - J_{0}(t)J_{2}(t)) dt - \frac{1}{48} \sum_{p=1}^{P} \frac{(\rho_{p}a)^{4}}{p^{2}}. $$

We use again (2.3) to obtain:

$$ \lambda_{\min}(A) \geq 1 - F(\gamma) - \gamma^4 \frac{\pi^4}{48} \frac{1}{P^4} \sum_{p=1}^{P} p^2, $$

which obviously implies the claimed result. □

## Conditioning of A for large $\gamma$

The behavior of $\lambda_{\min}(A)$ is more difficult to study for large $\gamma$. Nevertheless, we observed the following exponential decay: for any $P \ge 10$, and $\gamma \ge 1.4$, the minimal eigenvalue of A is bounded below by

$$ \lambda_{\min}(A) \geq 180 \exp(-5.8\gamma). \qquad (3.4) $$
<figure><img src="image_5.png" /><figcaption>Figure 3: Lower bounds from Theorem 4 and Equation 3.4, and estimated minimal eigenvalue of A in function of $\gamma$ for $P = 50$ (blue circles), $P = 150$ (red crosses), $P = 500$ (yellow squares) and $P = 1500$ (purple diamonds). For $P = 50$, the computed minimal eigenvalue is negative due to numerical errors from $\gamma \approx 7.5$ so corresponding data cannot be displayed. The horizontal dashed line shows machine precision</figcaption></figure>

# 4 Application to Laplace and Helmholtz kernels

## 4.1 Laplace kernel

When solving PDE's involving the Laplace operator (for example in heat conduction or electrostatic problems), one is led to (1) with the Laplace kernel $G(r) = \log(r)$ (we have dropped the $-\frac{1}{2\pi}$ constant for simplicity). Here we show that its SBD converges exponentially fast:

**Theorem 5.** There exist two positive constants $L_1$ and $l_2$ such that

$$ \forall a \in (0, 1), \forall P \in \mathbb{N}^*, \forall r \in (a, 1), \quad \left| G(r) - \sum_{p=1}^{P} \alpha_p e_p(r) \right| \le L_1 e^{-l_2 a P} $$

where $\alpha_1, \dots, \alpha_P$ are the SBD coefficients of $G$ of order $P$.

We show this by exhibiting an extension $\tilde{G}$ for which we are able to estimate the remainder of the Fourier-Bessel series. For any $n \in \mathbb{N}^*$, let us define extensions $\tilde{G}_n$ of $G$ as

$$ \tilde{G}_n = \begin{cases} r^{2n} \sum_{k=0}^{2n} \frac{a_{k,n}}{k!} (r-a)^k & \text{if } r \le a, \\ G(r) & \text{otherwise.} \end{cases} \quad (4.1) $$

where the coefficients $a_{k,n}$ are chosen so that $\tilde{G}_n$ has continuous derivatives up to the order $2n$:

$$ a_{k,n} = \frac{d^k}{dr^k} \left( \frac{\log(r)}{r^{2n}} \right) \Big|_{r=a} $$
Notice that the $r^{2n}$ term ensures the boundedness of $(-\Delta)^n \tilde{G}_n$ near the origin. Also observe that for all $s \in \mathbb{N}$:

$$(-\Delta)^s \tilde{G}_n \text{ vanishes on } \mathcal{C},$$

since $\tilde{G}_n \equiv G$ in a vicinity of $\mathcal{C}$. We now go into some rather tedious computations to provide a crude bound for $\|(-\Delta)^n \tilde{G}_n\|_{L^2(B)}$ in terms of the coefficients $a_{k,n}$.

**Lemma 2.** There exists a constant $C$ independent of $n$ and a such that for $r < a$

$$|\Delta^n \tilde{G}_n(r)| \le C \left(\frac{16n}{e}\right)^{2n} \max_{k \in \{1, \dots, 2n\}} \left(\frac{|a_{k,n}|}{k!} a^k\right). \quad (4.2)$$

*Proof.* For $r \le a$, we have

$$\Delta^n \tilde{G}_n(r) = \sum_{k=0}^{2n} \sum_{l=0}^{k} \binom{k}{l} \frac{a_{k,n}}{k!} (-a)^{k-l} (2n+l)^2 (2(n-1)+l)^2 \times \dots \times (2+l)^2 r^l.$$

This result is obtained by expanding the sum in the definition of $\tilde{G}_n$ and using the fact that $\Delta r^k = k^2 r^{k-2}$. Hence, using triangular inequality

$$|(-\Delta)^n \tilde{G}_n(r)| \le \sum_{k=0}^{2n} \sum_{l=0}^{k} \binom{k}{l} \frac{|a_{k,n}|}{k!} a^{k-l} (2n+l)^2 (2(n-1)+l)^2 \times \dots \times (2+l)^2 r^l.$$

For $l \in \{1, \dots, 2n\}$, we apply the following (crude) inequality:

$$(2n+l)^2(2(n-1)+l)^2 \times \dots \times (2+l)^2 \le (4n)^2(4n-2)^2 \times \dots \times (2n+2)^2 \quad (4.3)$$

to obtain:

$$|(-\Delta)^n \tilde{G}_n(r)| \le (4n)^2 (4n-2)^2 \times \dots \times (2n+2)^2 \times \\ \times \max_{k \in \llbracket 0, 2n \rrbracket} \left( \frac{|a_{k,n}|}{k!} a^k \right) \sum_{k=0}^{2n} \sum_{l=0}^{k} \binom{k}{l} a^{-l} r^l$$

Since $r < a$, the last sum is bounded by $\sum_{k=0}^{2n} 2^k = 2^{2n+1} - 1 < 2^{2n+1}$, while

$$(4n)^2 (4n-2)^2 \times \dots \times (2n+2)^2 \sim 2 \left(\frac{8n}{e}\right)^{2n}$$

follows from Stirling formula.

<figure><img src="image_6.png" /></figure>

We are now able to prove the following, which implies Theorem 5.

**Theorem 6.** There exists a constant $C$ such that, for any $P \in \mathbb{N}^*$ and $a \in (0, 1)$, there exists a radial function $\tilde{G}$ which coincides with $G$ on $\mathcal{A}(a)$ satisfying:

$$\left\| \tilde{G} - \sum_{p=1}^{P} c_p(\tilde{G}) e_p \right\|_{H_{0,rad}^{1}(B)} \le C \sqrt{P} \exp \left( -\frac{aP\pi}{32} \right).$$
Proof. Let $n \in \mathbb{N}^*$, by Leibniz formula,

$$ \frac{d^k}{dr^k} \left( r^{-2n} \log(r) \right) = \frac{(-1)^k k!}{r^{2n+k}} \left( -\sum_{j=0}^{k-1} \frac{\binom{2n+j-1}{j}}{k-j} + \binom{2n+k-1}{k} \log(r) \right). $$

This leads to

$$ \frac{|a_{k,n}|}{k!} a^k \le a^{-2n} \binom{2n+k-1}{k} \left( \frac{k}{2n} - \log(a) \right), $$

where we have used the identity

$$ \sum_{j=0}^{k-1} \binom{j+2n-1}{j} = \frac{k}{2n} \binom{k+2n-1}{k}. $$

Observe that

$$ \binom{2n+k-1}{k} \le \binom{4n-1}{2n} = \frac{1}{2} \binom{4n}{2n} \le \frac{4^{2n}}{2\sqrt{2\pi n}} \quad k \in \{1, \dots, 2n\}, $$

and thus,

$$ \max_{0 \le k \le 2n} \left( \frac{|a_{k,n}|}{k!} a^k \right) \le \left( \frac{4}{a} \right)^{2n} \frac{1}{2\sqrt{2\pi n}} \left( \log \left( \frac{e}{a} \right) \right). \quad (4.4) $$

Combining (4.4) with estimation (4.2), we find that there exists a constant $C$ such that, for $r < a$

$$ |(-\Delta)^n \tilde{G}_n(r)| \le \frac{C}{\sqrt{n}} \left(\frac{16n}{e}\right)^{2n} \left(\frac{4}{a}\right)^{2n} \log\left(\frac{e}{a}\right). $$

Therefore, integrating on $B(0, a)$, we get

$$ \|(-\Delta)^n \tilde{G}_n\|_{L^2(B(0,a))} \le \frac{Ca^2}{\sqrt{n}} \log\left(\frac{e}{a}\right) \left(\frac{64n}{ae}\right)^{2n}, $$

and since

$$ (-\Delta)^n \tilde{G}_n(x) = (-\Delta)^n G(x) = 0 $$

for $|x| > a$, the same bound applies to $\|(-\Delta)^n \tilde{G}_n(x)\|_{L^2(B)}$. We now plug this estimate into the inequality of corollary 2, to get

$$ \left\| \tilde{G}_n - \sum_{p=1}^{P} c_p(\tilde{G}_n) e_p \right\|_{H_{0,rad}^1(B)} \le C \frac{P^{\frac{3}{2}}}{n} a^2 \log\left(\frac{e}{a}\right) \left(\frac{64n}{aeP\pi}\right)^{2n}. $$

The previous inequality holds true for any integer $n$ such that $n > 1$ and any $P \in \mathbb{N}$. Without loss of generality, one can assume that $\frac{aP\pi}{64} > 1$. In this case, let $n_P = \lfloor \frac{aP\pi}{64} \rfloor$, and $\tilde{G} = \tilde{G}_{nP}$. Using the fact that $x \mapsto x \log\left(\frac{e}{x}\right)$ is bounded on $(0, 1]$, we get

$$ \left\| \tilde{G} - \sum_{p=1}^{P} c_p(\tilde{G}) e_p \right\|_{H_{0,rad}^1(B)} \le C \sqrt{P} e^{-\frac{aP\pi}{32}}. $$

□
**Remark 3.** The convergence rate is indeed bounded by a function of the parameter $\gamma$. Figure 4d shows the decay of the $L^\infty$ error in function of $\gamma$ for different values of $P$. It can be seen that the error consistently decreases at an exponential rate of about $l_2 \approx 3.7$, and stagnates at the minimal error $e_{\min} \approx 10^{-10}$. We believe this stagnation is due to rounding errors related to the increasing condition number of the matrix $A$. However, note that the situation is invariant for constant $\gamma$, thus the stability of the linear system only depends on the target tolerance and not on the size of the problem. For example, to reach the error level $\varepsilon = 10^{-3}$, it is sufficient to take $\gamma \approx 1.8$ and the conditioning of $A$ is about 200, independent on the specific value of $P$ and $a$.

<figure><img src="image_7.png" /><figcaption>Figure 4: Approximation of $G(r) = \log r$ by SBD. (a): $L^\infty$ error between $\log(x)$ (dashed black) and its SBD approximation (solid blue) with $a = 0.05$ and $P = 5$. (b): Values of $\log r - \sum_p \alpha_p e_p(r)$ where $\alpha_p$ are the SBD coefficients for $a = 0.05$ of order $P = 5$ (Blue circles), $P = 30$, (red crosses) and $P = 50$ (yellow squares). (c): Numerical values of the SBD coefficients $\alpha_p$ in function of $p$, for $P = 5$ (blue circles), $P = 30$ (red crosses) and $P = 50$ (yellow squares). The black stars show the values of the exact (slowly decaying) Fourier-Bessel coefficients of $G$, that is $\alpha_p = c_p(G) = \frac{1}{\sqrt{\pi \rho_p |J_1(\rho_p)|}}$. (d): Evolution of the $L^\infty$ error over $[a, 1]$ associated to the SBD for different values of $P$ in function of $\gamma$. An exponential decay is indeed observed, at the roughly estimated rate of $\propto \exp(-3.7\gamma)$. The stops decreasing at $e_{\min} \approx 10^{-10}$, for a value of $\gamma \approx 6.7$ because of the ill conditioning of the linear system (3.3)</figcaption></figure>
## 4.2 Helmholtz kernel

Let $Y_0$ the classical Bessel function of second kind and of order 0. For any $k > 0$, the Helmholtz kernel, $r \mapsto -\frac{i}{4}H^{(1)}(kr)$, where $H^{(1)}(r) = J_0(r) + iY_0(r)$, is the fundamental Green's kernel associated to the harmonic wave operator $-\Delta - k^2$, that satisfies a Sommerfeld radiation condition at infinity, (see for example [25]). This kernel arises in various physical problems, such as sound waves scattering. To approximate $H^{(1)}(kr)$ as a sum of dilated $J_0$ functions, it is sufficient to produce a SBD decomposition of $r \mapsto Y_0(kr)$. We can obtain good approximations of $Y_0$ in series of dilated functions in the following way:

- **When $k$ is a root of $Y_0$**: In this case the multi-Dirichlet condition is satisfied at any order. Indeed, for any $n$,

$$(-\Delta)^n Y_0(kr)|_{r=1} = k^{2n} Y_0(k) = 0.$$

We thus produce a SBD decomposition of $Y_0$ on an interval $(a, 1)$. Just like for the Laplace kernel, it was observed that the approximation error converges exponentially fast to zero, as soon as $P$ is greater than $k$.

- **When $k$ is close to, or greater than the first root of $Y_0$**: we find a Sparse Bessel Decomposition for $r \mapsto Y_0(k'r)$ on $(a, 1)$, where $k'$ is the first root of $Y_0$ larger than $k$. This provides a decomposition for $r \mapsto Y_0(kr)$ valid on $(\frac{k'}{k}a, \frac{k'}{k})$.

- **When $k$ is much smaller than the first root of $Y_0$**: the previous idea might lead to unnecessary efforts. Indeed, to ensure that $\frac{k'}{k}a$ is small enough, one would have to choose a very small value of $a$ leading to a very long Bessel series. Instead, one can use the Bessel-Fourier series associated to the Robin condition (see subsection 2.3):

$$\frac{\partial u}{\partial n} + Hu = 0,$$

noticing that $H = -\frac{kY_0'(k)}{Y_0(k)} > 0$ in this region.

<figure><img src="image_8.png" /><figcaption>Figure 5: SBD of $Y_0(kr)$ with $k$ a root of $Y_0$ approximately equal to 7.086, with $P = 10$ terms and $a = 0.05$</figcaption></figure>
## 4.3 General kernel : enforcing the multi-Dirichlet condition

For general kernels $G$, the multi-Dirichlet conditions may not be fulfilled, even after rescaling. When applying the SBD method without any changes, this leads to the situation observed in the top panel of Figure 6. In this figure, the SBD is applied to the kernel $G_1(r) = \log(r) + \sin(r)$ (note that $\Delta G_1(1) \neq 0$), and we plot the error in function of $r$ for several values of $\gamma$ and $P = 10$. The error curve stagnates near $r = 1$. In the bottom panel, we apply the SBD to $G_2(r) = \log(r) + \sin(r) - \frac{1}{4}\Delta G_1(1)r^2$ where the last term enforces the Dirichlet condition, and the error decay is improved.

<figure><img src="image_9.png" /><figcaption>Figure 6: SBD of order 10 applied to the kernels $G_1(r) = \log(r) + \sin(r)$ (top panel) and $G_2(r) = \log(r) + \sin(kr) - \frac{1}{4}\Delta G_1(1)r^2$ (bottom panel), for $\gamma = 1$ (blue circles), $\gamma = 2$ (red crosses) and $\gamma = 3$ (purple squares). The vertical dashed lines show the position of $a$ for each $\gamma$.</figcaption></figure>

More generally, for any radial $G$, we can apply the SBD to a modified function $H = G - K$ where $K$ is chosen to enforce the multi-Dirichlet condition. Since we wish to obtain for $G$ a decomposition in sum of Bessel functions, we propose to
choose K in the form

$$K(r) = \sum_{t=1}^{n} \mu_t J_0(\omega_t r)$$

for some $(\omega_t)_{1 \le t \le n}$ that are not roots of $J_0$. It is not the aim of this paragraph to describe a systematic way of choosing $(\omega_t)_{1 \le t \le n}$, as this work is still in progress. However, when the first few iterates of the Laplace operator on $G$ are known, we suggest the following choice: let $\omega'$ the square root of the average ratio between too successive iterates, choose $\omega_1$ as the root of $J_1$ that is closest from $\omega'$. Then, assign $\omega_2, \dots, \omega_n$, successively to the closest roots of $J_1$. The coefficients $(\mu_t)_{1 \le t \le n}$ are finally found by inverting a small linear system

$$M\mu = \lambda,$$

where $\lambda$ is the vector given by

$$\lambda_t = (-\Delta)^t G|_{r=1}, \quad t \in \{1, \dots, n\},$$

and with

$$M = \begin{bmatrix} -\omega_1^2 & -\omega_2^2 & \cdots & -\omega_n^2 \\ \omega_1^4 & \omega_2^4 & \cdots & \omega_n^4 \\ \vdots & \vdots & \ddots & \vdots \\ (-1)^n \omega_1^{2n} & (-1)^n \omega_2^{2n} & \cdots & (-1)^n \omega_n^{2n} \end{bmatrix}$$

In Figure 7, we show the efficiency of this method by applying the SBD with 100 terms to some highly oscillating function $(x \mapsto \log(x) + \sin(250x))$, for $n \le 3$, and computing the maximal error of the decomposition in function of $\gamma$. The frequencies $(\omega_t)_{1 \le t \le n}$ are the roots of $J_1$ that are closest to 250.

<figure><img src="image_10.png" /><figcaption>Figure 7: Maximal error between the kernel $G_1(r) = \log(x) + \sin(250x)$ and its 100-terms SBD in function of $\gamma$ using the method described in this paragraph for several values of $n$.</figcaption></figure>
# 5 Circular quadrature

In this section, we study the approximation of the form

$$J_0(\rho_p|x|) \approx \frac{1}{M_p} \sum_{m=0}^{M_p-1} e^{i\rho_p \xi_m^p \cdot x},$$

for some integer $M_p$ and some quadrature points $(\xi_m^p)_{1 \le m \le M_p}$.

## 5.1 Theoretical bound

**Theorem 7.** There exists a constant $K$ such that for any $r > 0$, $M \in \mathbb{N}^*$, such that $M \ge \frac{\xi}{2}r$, and for any $\varphi \in \mathbb{R}$

$$\left| J_0(r) - \frac{1}{M} \sum_{m=0}^{M-1} e^{i r \sin\left(\frac{2m\pi}{M} - \varphi\right)} \right| \le K \left( \frac{er}{2M} \right)^M.$$

In order to prove this proposition, we first prove a result on Fourier series

**Lemma 3.** For any $2\pi$-periodic $\mathcal{C}^2$ complex-valued function $f$ one has

$$\frac{1}{2\pi} \int_{0}^{2\pi} f - \frac{1}{M} \sum_{m=0}^{M-1} f\left(\frac{2m\pi}{M}\right) = - \sum_{k \in \mathbb{Z}^*} c_{kM}(f),$$

where $c_n(f)$ denotes the Fourier coefficient of $f$ defined as

$$c_n(f) = \frac{1}{2\pi} \int_{0}^{2\pi} f(x)e^{-inx}dx.$$

*Proof.* Since $f$ is $\mathcal{C}^2$, it is equal to its Fourier Series, which converges normally:

$$\forall x \in \mathbb{R}, f(x) = \sum_{k \in \mathbb{Z}} c_k(f) e^{ikx}.$$

Using this expression, we obtain

$$\frac{1}{M} \sum_{m=0}^{M-1} f\left(\frac{2m\pi}{M}\right) = \sum_{k \in \mathbb{Z}^*} c_k(f) \left( \frac{1}{M} \sum_{m=0}^{M-1} e^{ik \frac{2m\pi}{M}} \right).$$

Now observe that

$$\frac{1}{M} \sum_{m=0}^{M-1} e^{ik \frac{2m\pi}{M}} = \begin{cases} 1 & \text{if } k \notin M\mathbb{Z}, \\ 0 & \text{otherwise.} \end{cases}$$

Therefore

$$\int_{0}^{2\pi} f(x)dx - \frac{1}{M} \sum_{m=0}^{M-1} f\left(\frac{2m\pi}{M}\right) = c_0(f) - \sum_{k \in M\mathbb{Z}} c_k(f) = - \sum_{k \in \mathbb{Z}^*} c_{kM}(f).$$

We now turn to the proof of Theorem 7:
*Proof.* The result is based on the following fact:

$$J_0(r) = \int_0^{2\pi} e^{ir\sin(x)} dx = \int_0^{2\pi} e^{ir\sin(x-\varphi)} dx.$$

Let $f: x \mapsto e^{ir\sin(x-\varphi)}$. Let us recall the integral representation of the Bessel function of the first kind and of order $k$ where $k$ is a relative integer:

$$J_k(r) = \int_0^{2\pi} e^{ir\sin(x)} e^{-ikx} dx = e^{-ik\varphi} \int_0^{2\pi} e^{ir\sin(x-\varphi)} e^{-ikx} dx.$$

Thus, one has $c_k(f) = e^{ik\varphi}J_k(r)$. Consequently, the former Lemma yields

$$J_0(r) - \frac{1}{M} \sum_{j=0}^{M-1} e^{ir\sin(\frac{2j\pi}{M}-\varphi)} = - \sum_{k \in \mathbb{Z}^*} e^{iNk\varphi} J_{Nk}(r).$$

For large $|k|,$

$$J_k(r) \sim \left( \frac{er}{2|k|} \right)^{|k|}.$$

Therefore, there exists a constant $C'$ such that:

$$\begin{aligned} \left| J_0(r) - \frac{1}{M} \sum_{m=0}^{M-1} e^{ir\sin(\frac{2m\pi}{M}-\varphi)} \right| &\le C' \sum_{k \in \mathbb{Z}^*} \left( \frac{er}{2M|k|} \right)^{M|k|} \\ &\le K \left( \frac{er}{2M} \right)^M \end{aligned}$$

for an appropriate choice of $K$.

<figure><img src="image_11.png" /></figure>

We conclude with the following result

**Proposition 2.** Let $\varepsilon > 0$, $r > 0$, and assume $M > \frac{e}{2}r + \log(\frac{K}{\varepsilon})$. Then

$$\left| J_0(r) - \frac{1}{M} \sum_{m=0}^{M-1} e^{ir\sin(\frac{2m\pi}{M}-\varphi)} \right| \le \varepsilon.$$

*Proof.* This result is a direct consequence of the previous proposition together with the following inequality: for any $(A,B) \in (\mathbb{R}_+^*)^2$ one has

$$\left( \frac{A}{A+B} \right)^{A+B} \le e^{-B}.$$

To prove it, we take the logarithm of this quantity,

$$f(A, B) = -B \left(1 + \frac{A}{B}\right) \log \left(1 + \frac{B}{A}\right)$$

and observe that for any positive $x$,

$$(1 + \frac{1}{x}) \log(1 + x) \ge 1.$$

<figure><img src="image_12.png" /></figure>
Hence, we can approximate very efficiently the functions $e_p$ of the previous paragraph as a finite sum as follows. We define the quadrature points $\xi_0^p, \xi_1^p, \dots, \xi_{M_p-1}^p$ by

$$ \xi_m^p := e^{i \frac{2\pi m}{M_p}} \quad 1 \le p \le P, \quad 0 \le m \le M_p - 1. \qquad (5.1) $$

With this definition, for any $x \in \mathbb{R}^2$

$$ e_p(|x|) = C_p J_0(\rho_p |x|) \approx \frac{C_p}{M_p} \sum_{m=0}^{M_p-1} e^{i \rho_p x \cdot \xi_m^p}, $$

and the approximation is valid at a precision $\epsilon$ as soon as $M_p > \frac{\epsilon}{2} \rho_p |x| + \log\left(\frac{KC_p}{\epsilon}\right)$.

# 6 Estimations of complexities

We now turn to the complexity estimate of the complete algorithm given in section 1. We fix $\alpha \in [0, 1/6]$ and let

$$ a = \frac{|\log \epsilon|^{2/3}}{N_z^{2/3} - \alpha} \qquad (6.1) $$

as in Theorem 1. We will give a bound for the number of operations of each step of the algorithm in function of $N_z$, $\epsilon$ and $\alpha$. We note $C_{SBD}$, $C_{circ}$, $C_{assemble}$, $C_{far}$ and $C_{close}$ respectively, the number of operations required to produce the SBD, the circular quadrature, to assemble the close correction matrix $D$ defined in (1.1), to compute the far approximation defined in (1.2), and to apply $D$ on a vector. We will denote by $C$ any positive constant that is independent of $N_z$, $\epsilon$ and $\alpha$.

## 6.1 Offline computations

The first part of the algorithm consists in combining the SBD with the circular quadratures detailed in the previous two sections to derive an approximation scheme for the log function in the following form:

$$ \log|x| = \sum_{v=1}^{N_z} \hat{w}_v e^{ix \cdot \xi_v} \quad a < |x| < 1 $$

valid to the accuracy $\epsilon$.

**Sparse Bessel Decomposition:** We first compute a SBD of log on the ring $\{a < r < 1\}$ to reach the accuracy $\epsilon/2$, as developed in Sections 3 and 4. We write this approximation

$$ \log r \approx \sum_{p=1}^{I} \alpha_p e_p(r), \quad a < r < 1. $$

Theorem 5 shows that the accuracy is reached for

$$ P = O\left(\frac{|\log(\epsilon)|}{a}\right). \qquad (6.2) $$
Since the coefficients $\alpha_1, \dots, \alpha_p$ are obtained through the inversion of a $P \times P$ matrix, the computation of the SBD requires $O(P^3)$ computations. Therefore, there exists a constant $C > 0$ independent of $N_z, \varepsilon$ and $\alpha$ such that

$$C_{\text{SBD}}(N_z, \varepsilon, \alpha) \leq C |\log \varepsilon| N_z^{2-3\alpha}. \quad (6.3)$$

**Circular quadrature:** We approximate each function $e_p$ using a circular quadrature as detailed in section 5. For each $p$, we choose the number $M_p$ of terms in the quadrature so that

$$\left| J_0(\rho_p |x|) - \frac{1}{M_p} \sum_{m=1}^{M_p} e^{i \rho_p x \cdot \xi_m^p} \right| \leq \frac{\varepsilon}{2P |\alpha_p| C_p} \quad a < |x| < 1, \quad (6.4)$$

where the quadrature points $\xi_m^p$ are defined in (5.1). Proposition 2 implies that taking

$$M_p > \frac{e}{2} \rho_p + \log \left( \frac{2KP |\alpha_p|}{\varepsilon} \right) \quad (6.5)$$

is sufficient to ensure (6.4). In this case, triangular inequality implies that for $a \leq |x| \leq 1$,

$$\left| \log|x| - \sum_{p=1}^{p} \sum_{m=1}^{M_p} C_p \frac{\alpha_p}{M_p} e^{ix \cdot \xi_m^p} \right| \leq \varepsilon.$$

Bessel's inequality

$$\sum_{p=1}^{p} |\alpha_p|^2 \leq \int_{B} |\log|x||^2$$

ensures the boundedness of $\alpha_p$. Moreover, $\rho_p = O(p)$, implying $M_p = O(P)$, and hence,

$$N_{\xi} = \sum_{p=1}^{p} M_p = O(P^2). \quad (6.6)$$

Since for any $p$, the computation of $(\xi_m^p)_{1 \leq m \leq M_p}$ has a linear complexity in $M_p$, we get:

$$C_{\text{circ}}(N_z, \varepsilon, \alpha) \leq C |\log \varepsilon|^{2/3} N_z^{4/3-2\alpha}. \quad (6.7)$$

Equations (6.3) and (6.7) yield the first part of Theorem 1.

**Close correction matrix:** Recall that $\delta_{\max}$ is defined as

$$\delta_{\max} = \max_{1 \leq k,l \leq N_z} |z_k - z_l|$$

where $z$ are the data points in (1). Let $\delta_{\min} = a\delta_{\max}$. We first determine the set $\mathcal{P}$ of all pairs $(k,l)$ such that $|z_k - z_l| \leq \delta_{\min}$. This is the classical "fixed-radius near neighbors search", and can be solved in $O(N_z \log N_z + \#P)$ operations (see for example [4, 5, 11, 23]). In order to compute an approximation of the close correction sparse matrix:

$$D_{kl} = \delta_{(k,l) \in \mathcal{P}} \left( \log |z_k - z_l| - \sum_{v=1}^{N_z} e^{i(z_k - z_l) \cdot \xi_v} \hat{\omega}_v \right),$$
we need $\#\mathcal{P}$ evaluations of log, and the computation of

$$ \text{NUFFT}_{-}[(z_k - z_l)_{(k,l) \in \mathcal{P}}, \xi](\hat{v}) $$

at precision $\varepsilon$. For data uniformly distributed on a curve, the number of close pairs scales as

$$ \#\mathcal{P} = O\left(\frac{\delta_{\min}}{\delta_{\max}} N_z\right) = O(N_z^2 a). \qquad (6.8) $$

One can check that $N_\xi \le N_z^2 a$ using (6.1), (6.2) and (6.6), so that

$$ C_{\text{close}} \le C C_{\text{NUFFT}}(\varepsilon) |\log \varepsilon|^{2/3} N_z^{4/3+\alpha} \log(N_z). \qquad (6.9) $$

This is the second part of Theorem 1.

## 6.2 On-line Computations

**Far approximation:** Recall that for all $k \in \{1, \dots, N_z\}$, the far approximation is defined by the following equation:

$$ q_k^{\text{far}} = \sum_{l=1}^{N_z} G_{\text{approx}}(z_k - z_l) f_l, $$

where, according to the previous subsection,

$$ G_{\text{approx}}(x) = \sum_{p=1}^{P} \sum_{m=1}^{M_p} \frac{\alpha_p C_p}{M_p} e^{ix \cdot \xi_m^p}. $$

Define $\hat{\omega} = (\hat{\omega}_v)_{1 \le v \le N_\xi}$ and $\xi = (\xi_v)_{1 \le v \le N_\xi}$ such that

$$ G_{\text{approx}}(x) = \sum_{v=1}^{N_\xi} \hat{\omega}_v e^{ix \cdot \xi_v}. $$

To compute $q^{\text{far}}$, recall the following three steps:

(i) **Space** → **Fourier:** Compute $\hat{f} = \text{NUFFT}_{-}[z, \xi](f)$,

(ii) **Fourier multiply** Perform elementwise multiplication by $\hat{\omega}$:

$$ \hat{g}_v = \hat{\omega}_v \hat{f}_v, $$

(iii) **Fourier** → **Space:** Compute $q^{\text{far}} = \text{NUFFT}_{+}[z, \xi](\hat{g})$.

One can check that (6.1) and (6.6) imply $N_\xi \ge N_z$, thus

$$ C_{\text{far}}(N_z, \varepsilon, \alpha) \le C |\log \varepsilon|^{2/3} C_{\text{NUFFT}}(\varepsilon) N_z^{4/3-2\alpha} \log(N_z). \qquad (6.10) $$
**Close correction:** D has #P non-zero entries so

$$C_{\text{close}}(N_z, \varepsilon, \alpha) \leq C |\log \varepsilon|^{2/3} N_z^{4/3+\alpha}. \qquad (6.11)$$

Summing (6.10) and (6.11) yields the second part of Theorem 1.

**Remark 4.** The extreme cases $\alpha = 0$ and $\alpha = 1/6$ correspond respectively to the situations where one wish to either minimize the total (off-line + on-line) computation time or just the on-line time. The complexities "off-line" and "on-line" then become (omitting the dependence in $\varepsilon$):

<table><caption>Table 1: Complexity of the algorithm (omitting dependence in $\varepsilon$) in the two extreme cases $\alpha = 0$ and $\alpha = 1/6$.</caption><thead><tr><th></th><th>Off-line</th><th>On-line</th></tr></thead><tbody><tr><td>$\alpha = 0$</td><td>$O(N_z^2)$</td><td>$O(N_z^{4/3})$</td></tr><tr><td>$\alpha = 1/6$</td><td>$O(N_z^{3/2} \log N_z)$</td><td>$O(N_z^{3/2})$</td></tr></tbody></table>

# 7 Numerical examples

## 7.1 System of $N \times N$ particles

To assess for the numerical performance of our method, we first generate two sets of $N$ points $(x_k)$ and $(y_k)$ uniformly distributed in a square, for $N$ ranging from 10 to $10^7$, and compute the discrete convolution

$$q_k = \sum_{l=1}^{N} \log |x_k - y_l| f_l$$

where $f_l$ is a random vector in $\mathbb{R}^N$. We measure both the time needed for off-line and on-line parts. We also measure the amount of memory occupied by the assembled operator (Memory usage goes a little bit above this value during the computation). The results are displayed in Table 2. The computer used for this test is a laptop cadenced to 1.6 GHz and possessing 32 GB of memory.

<table><caption>Table 2: Performances of the algorithm</caption><thead><tr><th></th><th>Off-line (s)</th><th>On-line(s)</th><th>Memory</th><th>Proportion of full matrix size</th></tr></thead><tbody><tr><td>$N = 10$</td><td>0.034</td><td>0.001</td><td>7 kb</td><td>800% (no compression)</td></tr><tr><td>$N = 10^2$</td><td>0.027</td><td>0.002</td><td>83 kb</td><td>105% (no compression)</td></tr><tr><td>$N = 10^3$</td><td>0.049</td><td>0.003</td><td>1.02 Mb</td><td>14%</td></tr><tr><td>$N = 10^4$</td><td>0.37</td><td>0.02</td><td>10.7 Mb</td><td>1.4%</td></tr><tr><td>$N = 10^5$</td><td>3.4</td><td>0.16</td><td>109 Mb</td><td>0.15%</td></tr><tr><td>$N = 10^6$</td><td>47.2</td><td>2.6</td><td>1.07 Gb</td><td>0.015%</td></tr><tr><td>$N = 10^7$</td><td>441</td><td>41</td><td>7.13 Gb</td><td>9.810<sup>-4</sup>%</td></tr></tbody></table>
## 7.2 Sound canceling

Consider $N_z$ punctual 2-dimensional sound sources located at $(z_l)_{1 \le l \le N_z}$ and emitting at a single frequency $f$, with unit amplitude and phases $(\varphi_l)_{1 \le l \le N_z}$. That is, for each $l \in \{1, \dots, N_z\}$, the source number $l$ generated an acoustic pressure at each point $x$ and time $t$ in $\mathbb{R}^2$ equal to

$$p_l(x,t) = \Re \left[ -\frac{i}{4} H^{(1)}(k|z_l - x|) e^{-i(\omega t - \varphi_l)} \right],$$

where $\omega = 2\pi f$, $k = \frac{2\pi f}{c}$, with $c$ the celerity of the sound waves, and $H^{(1)}$ is the Hankel function of first kind already defined in subsection 4.2. By superposition, the resulting pressure at $x$ is

$$p(x,t) = \Re \left[ -\frac{i}{4} e^{-i\omega t} \sum_{l=1}^{N_z} H^{(1)}(k|z_l - x|) e^{i\varphi_l} \right],$$

and the sound intensity is proportional to

$$\Pi(x, \varphi) = \left| \sum_{l=1}^{N_z} H^{(1)}(k|z_l - x|) e^{i\varphi_l} \right|^2.$$

Suppose one wishes to choose the phases that minimize the sound intensity in a prescribed zone $\Omega$ (called the silence zone), that is,

$$\varphi^* = \underset{\varphi \in [0, 2\pi]^{N_z}}{\operatorname{argmin}} \int_{\Omega} \Pi(x, \varphi).$$

If we approximate the integral over $\Omega$ by a uniform quadrature, this leads to solving

$$\varphi^* = \underset{\varphi \in [0, 2\pi]^{N_z}}{\operatorname{argmin}} \sum_{q=1}^{Q} \Pi(x_q, \varphi),$$

where $x_q$ are the coordinated of the quadrature points. If we let $A_{l,q} = H^{(1)}(k|x_q - z_l|)$ and $q(\varphi) = (e^{i\varphi_l})_{1 \le l \le N_z}$, this rewrites

$$\varphi^* = \underset{\varphi \in [0, 2\pi]^{N_z}}{\operatorname{argmin}} q(\varphi)^T A^T A q(\varphi).$$

Using our method, the cost function associated to this minimization problem can be evaluated rapidly, as well as its gradient. Thus, using black-box optimization procedures, we can find rapidly good candidates for $\varphi^*$. In Figure 8, we show the result of one such optimization, with $N_z = 100$ sound sources randomly located on a half circle and where the zone of silence is represented by the red circle. The silence zone is discretized using a mesh of $Q = 2.5 \times 10^5$ points. We stopped the optimization after 500 evaluations of the objective function and its gradient, which required a total computation time of 3 minutes on our computer. Note that since the silence zone is located far away from the sources, we don't need to use a close correction matrix to compute the objective function. To produce the image in Figure 8, the resulting field is evaluated at $2 \times 10^7$ grid points.
<figure><img src="image_13.png" /><figcaption>Figure 8: Sound level (in dB) created by 100 punctual 2D sound sources with wavenumber *k* ≈ 90 located at random spots on a half circle with radius *r* = 0.5. The sources locations are indicated by the small black circles. On the left, we show the sound level before the optimization of the phases. On the right, we show the result once the phases of the sources have been optimized to minimize the sound level in the region delimited by the red circle.</figcaption></figure>

References

[1] M. Abramowitz and I. A. Stegun. *Handbook of mathematical functions: with formulas, graphs, and mathematical tables*, volume 55. Courier Corporation, 1964.

[2] F. Alouges and M. Aussal. The sparse cardinal sine decomposition and its application for fast numerical convolution. *Numerical Algorithms*, 70(2):427–448, 2015.

[3] P. Balodis and A. Córdoba. The convergence of multidimensional fourier-bessel series. *Journal d'Analyse Mathématique*, 77(1):269–286, 1999.

[4] J. L. Bentley. Multidimensional binary search trees used for associative searching. *Communications of the ACM*, 18(9):509–517, 1975.

[5] J. L. Bentley, D. F. Stanat, and E. H. Williams. The complexity of finding fixed-radius near neighbors. *Information processing letters*, 6(6):209–212, 1977.

[6] S. Börm, L. Grasedyck, and W. Hackbusch. Introduction to hierarchical matrices with applications. *Engineering analysis with boundary elements*, 27(5):405–422, 2003.

[7] H. Cheng, L. Greengard, and V. Rokhlin. A fast adaptive multipole algorithm in three dimensions. *Journal of computational physics*, 155(2):468–498, 1999.

[8] R. Coifman, V. Rokhlin, and S. Wandzura. The fast multipole method for the wave equation: A pedestrian prescription. *IEEE Antennas and Propagation Magazine*, 35(3):7–12, 1993.

[9] L. Colzani, A. Crespi, G. Travaglini, and M. Vignati. Equiconvergence theorems for fourier-bessel expansions with applications to the harmonic analysis of radial functions in euclidean and non-euclidean spaces. *Transactions of the American Mathematical Society*, 338(1):43–55, 1993.
[10] J. W. Cooley and J. W. Tukey. An algorithm for the machine calculation of complex fourier series. *Mathematics of computation*, 19(90):297–301, 1965.

[11] M. T. Dickerson and R. S. Drysdale. Fixed-radius near neighbors search algorithms for points and segments. *Information Processing Letters*, 35(5):269–273, 1990.

[12] A. Dutt and V. Rokhlin. Fast fourier transforms for nonequispaced data. *SIAM J. Sci. Comput.*, 14(6):1368–1393, November 1993.

[13] L. Greengard. *The rapid evaluation of potential fields in particle systems*. MIT press, 1988.

[14] L. Greengard and J. Y. Lee. Accelerating the nonuniform fast fourier transform. *SIAM review*, 46(3):443–454, 2004.

[15] J. J. Guadalupe, M. Pérez, F. J. Ruiz, and J. L. Varona. Mean and weak convergence of fourier-bessel series. *Journal of mathematical analysis and applications*, 173(2):370–389, 1993.

[16] F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, and B. V. Saunders. *NIST Digital Library of Mathematical Functions*. http://dlmf.nist.gov/, Release 1.0.16 of 2017-09-18.

[17] D. Potts, G. Pöplau, and U. van Rienen. Calculation of 3d space-charge fields of bunches of charged particles by fast summation. In *Scientific Computing in Electrical Engineering*, volume 11, pages 241–246. Springer, 2006.

[18] D. Potts, G. Steidl, and A. Nieslony. Fast convolution with radial kernels at nonequispaced knots. *Numerische Mathematik*, 98(2):329–351, 2004.

[19] V. Rokhlin. Rapid solution of integral equations of scattering theory in two dimensions. *Journal of Computational Physics*, 86(2):414–439, 1990.

[20] V. Rokhlin. Diagonal forms of translation operators for the helmholtz equation in three dimensions. *Applied and Computational Harmonic Analysis*, 1(1):82–93, 1993.

[21] K. Stempak. On convergence and divergence of fourier-bessel series. *Electronic Transactions on Numerical Analysis*, 14:223–235, 2002.

[22] G. P. Tolstov. *Fourier series*. Courier Corporation, 2012.

[23] V. Turau. Fixed-radius near neighbors search. *Information processing letters*, 39(4):201–203, 1991.

[24] G. N. Watson. *A treatise on the theory of Bessel functions*. Cambridge university press, 1995.

[25] C. H. Wilcox. *Scattering theory for the d'Alembert equation in exterior domains*, volume 4. Springer Berlin, 1975.

[26] J. R. Wilton. The gibbs phenomenon in fourier-bessel series. *Journal f"ur die reine angewandte Mathematik*, 159:144–153, 1928.