# Localization of point scatterers via sparse optimization on measures

Giovanni S. Alberti, Romain Petit*, Matteo Santacesaria

MaLGa Center, Department of Mathematics, University of Genoa
Via Dodecaneso 35, 16146 Genoa, Italy

July 24, 2024

## Abstract

We consider the inverse scattering problem for time-harmonic acoustic waves in a medium with pointwise inhomogeneities. In the Foldy-Lax model, the estimation of the scatterers' locations and intensities from far field measurements can be recast as the recovery of a discrete measure from nonlinear observations. We propose a "linearize and locally optimize" approach to perform this reconstruction. We first solve a convex program in the space of measures (known as the Beurling LASSO), which involves a linearization of the forward operator (the far field pattern in the Born approximation). Then, we locally minimize a second functional involving the nonlinear forward map, using the output of the first step as initialization. We provide guarantees that the output of the first step is close to the sought-after measure when the scatterers have small intensities and are sufficiently separated. We also provide numerical evidence that the second step still allows for accurate recovery in settings that are more involved.

## 1 Introduction

The localization of inhomogeneities included in a medium has numerous applications in imaging. These inclusions might model defects in material science or tumors in medical imaging. A popular method for detecting them is to exploit scattering phenomena, which originate from the interaction of waves propagating in the medium with the inhomogeneities. By sending incident waves and measuring the scattered ones away from the medium, an observer can aim at estimating the localization and characteristics of the inclusions. This reconstruction problem is called the *inverse scattering problem* [37]. In this work, we focus on the case of time-harmonic acoustic incident waves.

In order to solve the inverse scattering problem with few measurements, one needs to make strong assumptions on the type of inhomogeneities included in the medium. The simplest case is arguably the one of very small inclusions, which can be approximately modeled by point-like inhomogeneities. This particular setting has been the subject of numerous works; see for example [40, 67, 2]. The classical model for the scattering of acoustic waves by point-like inclusions, known as *Foldy-Lax model*, was formally introduced in [50] (see also [64]). In this framework, the inverse scattering problem can be recast as the recovery of a discrete measure, which encodes the locations and the intensities of the scatterers, from nonlinear measurements.

In this article, we propose to connect this inverse problem to the *sparse spikes problem*, which has attracted a lot of attention in the past ten years (see e.g. [39, 20] and the review [63]). These works investigated the possibility of recovering a discrete measure from noisy linear observations. Although the measurements we consider in the context of inverse scattering are nonlinear, we propose to investigate how the guarantees and reconstruction methods developed for the sparse spikes problem can still be leveraged.

<footnote>*Corresponding author (email: romain.petit@edu.unige.it).</footnote>
## 1.1 Problem formulation

Before introducing the Foldy-Lax model for the scattering of acoustic waves by point-like inhomogeneities [2], we briefly describe the case of a piecewise continuous perturbation of the refractive index $q$. Let $d \in \{2, 3\}$ be the ambient dimension. Considering time-harmonic acoustic incident waves with wavenumber $\kappa > 0$, the incident field $u^{\text{in}}$ associated to an incident direction $\theta \in \mathbb{S}^{d-1}$ is given by

$$u^{\text{in}}(y) = e^{i\kappa\theta \cdot y}. \qquad (1)$$

The total field $u^{\text{tot}}$ measured by the observer, which is the sum of the incident field $u^{\text{in}}$ and the scattered field $u^s$, solves the Helmholtz equation (2) with Sommerfeld radiation condition (3)

$$\Delta u^{\text{tot}}(y) + \kappa^2(1+q(y))u^{\text{tot}}(y) = 0 \quad \text{in } \mathbb{R}^d, \qquad (2)$$

$$\lim_{|y| \to +\infty} |y|^{(d-1)/2} \left( \frac{\partial u^s}{\partial |y|}(y) - i\kappa u^s(y) \right) = 0. \qquad (3)$$

Now, let us denote the intensities and the locations of a set of $s$ point scatterers by $a \in (\mathbb{C}^*)^s$ and $x \in (\mathbb{R}^d)^s$, respectively, where $\mathbb{C}^* \stackrel{\text{def.}}{=} \mathbb{C} \setminus \{0\}$ denotes the set of nonzero complex numbers. The Foldy-Lax model formally corresponds to taking $q = \sum_{i=1}^s a_i \delta_{x_i}$ in (2). In this model, the total field $u^{\text{tot}}$ is given by

$$u^{\text{tot}}(y) = u^{\text{in}}(y) + \sum_{i=1}^{s} G(y, x_i) a_i u_i, \quad y \in \mathbb{R}^d,$$

where $G$ is the Green's function of the Helmholtz equation (see Section 1.4) and the coefficients $u = (u_i)_{1 \le i \le s}$ are the solutions of the Foldy-Lax system

$$u_i = u^{\text{in}}(x_i) + \kappa^2 \sum_{\substack{j=1 \\ j \ne i}}^{s} G(x_i, x_j) a_j u_j, \quad i \in \{1, \dots, s\}. \qquad (4)$$

We assume that $u^{\text{tot}}$ can be measured by the observer at infinity, i.e. that we have access to the far field pattern $u^{\infty}(\hat{x}, \theta)$ for several incident directions $\theta$ and observation directions $\hat{x}$, given by

$$u^{\infty}(\hat{x}, \theta) \stackrel{\text{def.}}{=} \frac{\kappa^2}{4\pi} \sum_{i=1}^{s} a_i u_i e^{-i\kappa \hat{x} \cdot x_i}, \quad (\hat{x}, \theta) \in \mathbb{S}^{d-1} \times \mathbb{S}^{d-1}. \qquad (5)$$

Our goal is to recover $a$ and $x$ from the knowledge of a possibly noisy version of $(u^{\infty}(\hat{x}, \theta))_{(\hat{x}, \theta) \in \Omega}$, with $\Omega$ a finite subset of $\mathbb{S}^{d-1} \times \mathbb{S}^{d-1}$. In this work, we consider the case where $\Omega = \{(\hat{x}_k, \theta_k)\}_{k=1}^m$, for a suitable choice of incident directions $\theta_k$ and corresponding observation directions $\hat{x}_k$. In other words, we consider $m$ incident waves, and for each incident wave we measure the corresponding far-field pattern at a single point at infinity. We stress that the difficulty of this task stems from the nonlinear dependence of $u^{\infty}$ on $a$ and $x$ and on the availability of only finitely many measurements.

**Derivation of the Foldy-Lax model as a vanishing size limit.** As a rule of thumb, the point approximation for small inhomogeneities is valid as soon as they have sub-wavelength size. The rigorous derivation of the Foldy-Lax model as a limit case of scattering by inhomogeneities of vanishing size has been investigated in several works. If the homogeneities have a moderate contrast<fnref n="1" />, the first term in the asymptotic expansion of the scattered field is indeed the one given by the Foldy-Lax model, but it goes to zero in the vanishing size limit (see e.g. [28, 33, 15]). On the other hand, there exist critical scalings between the size and the contrast

<footnote n="1">¹See also [26], in which the same behavior arises from a particular scaling between the (high) contrast and the size of the inhomogeneities.</footnote>
of the inhomogeneities under which resonating frequencies appear. When using incident frequencies close to those, the Foldy-Lax model is valid and the limit scattered field is non-trivial (see e.g. [7, 4, 5, 3, 66, 9]). In the acoustic setting, such resonance phenomena (called Minnaert resonances) have important applications in medical ultrasonic imaging [47], where the small scatterers model micro bubbles. We also mention their use in optics, where the scatterers model dielectric nanoparticles [6].

## 1.2 Related works

**Reconstruction methods in inverse acoustic scattering.** There is a wide literature on the numerical resolution of the inverse acoustic scattering problem. We mention here the three main approaches, namely iterative methods, decomposition methods, and sampling methods (see for example [37, Sections 4-6 of Chapter 5 and Sections 3-5 of Chapter 10]). Iterative methods (see for example [54, 12, 55, 56]) solve the equation associated to the inverse problem (in the case of noiseless observations) or a relevant optimization problem (in the case of noisy observations) by iteratively updating an initial guess of the inhomogeneity. Their main drawback is that, due to the nonlinearity of the forward map, convergence guarantees are only local. As a consequence, in order to provide a good initial guess, one needs extensive a priori information on the distribution of the scatterers. Decomposition methods, introduced in [59, 60, 61], split the inverse problem into an ill-posed linear problem and a well-posed nonlinear problem. Just as iterative methods, they only offer local recovery guarantees, but avoid the high cost of solving the full forward problem at each iteration. Sampling methods (see e.g. [36, 22, 58]) do not require extensive a priori knowledge of the inhomogeneity, but are of qualitative type: they mainly provide information on the location of the scatterers, but not on the perturbation of the refraction index. Finally, let us mention [13], which introduces a qualitative reconstruction method based on the notion of topological derivative.

**Reconstruction of small inhomogeneities.** The reconstruction of small inhomogeneities has also attracted a lot of attention outside the Foldy-Lax model and the specific case of acoustic waves. We refer the reader to [8] and to the long list of references therein concerning the localization of small conductivity, elastic and electromagnetic inclusions. In these works, the contrast of the inhomogeneities is moderate, and the measurements are hence close to those corresponding to the medium without the inclusions. Reconstruction methods are mainly based on asymptotic expansions of the perturbation resulting from the presence of the inhomogeneities.

**Born approximation.** When the scatterers have small intensities and are sufficiently separated, the solution of the Foldy-Lax system is close to $(u^{\text{in}}(x_i))_{1\le i\le s}$ (see Section 4.1 for more details). As a result, the far field pattern is close to the Fourier transform of $\mu = \sum_{i=1}^s a_i \delta_{x_i}$. Consequently, the recovery of $a$ and $x$ is a particular instance of the sparse spikes problem mentioned above. In the context of inverse scattering, most reconstruction methods are based on the Multiple Signal Classification (MUSIC) algorithm. This algorithm, first used in signal processing, was introduced for inverse scattering in [65] (see also [44, Chapter 10]). As pointed out in e.g. [34, 58], it has strong connections with sampling methods. Its noise robustness has been extensively studied in the literature on the sparse spikes problem. However, in the context of inverse scattering, we are not aware of any work leveraging these results to provide quantitative reconstruction guarantees accounting for the linearization error.

**Inverse scattering in the nonlinear case.** To our knowledge, in the Foldy-Lax setting, most works also rely on the MUSIC algorithm [53, 45, 32]. Exact recovery guarantees exist in the noiseless case, but we are not aware of any sensitivity analysis in the case of noisy measurements. In [52], the scatterers are constrained to belong to a fixed finite grid on the domain. The convergence of a nonlinear iterative hard thresholding algorithm is then investigated. The authors obtain global convergence guarantees, under the assumption that the iterates remain bounded and some coherence estimates are satisfied by the forward operator. Checking the validity of these assumptions is however involved, and is the subject of an intensive numerical investigation.
## 1.3 Contributions

We propose a “linearize and locally optimize” approach to solve the considered inverse problem. First, we solve a convex program on the space of measures, which involves the linearized forward operator (the far field pattern in the Born approximation). Then, we use this first estimate to initialize a local optimization procedure involving the nonlinear forward operator.

This approach allows us to leverage several results on the recovery of a discrete measure from noisy linear measurements (or sparse spikes problem). To our knowledge, their application to inverse scattering is new. In particular, they allow us to obtain quantitative recovery guarantees, as well as lower bounds on the number of required measurements. To this aim, we show that the assumptions of [72, Theorem 3] and [71, Theorem 1] are satisfied in the case of continuous Fourier sampling with frequencies uniformly distributed in a Euclidean ball of fixed radius. This setting was not covered by previous works, and is a result of independent interest.

We also extensively study the error committed by linearizing the forward operator. We treat the special case of two scatterers, and derive two different bounds in the general case. Numerical evidence suggests that their dependence on the wave number and on the minimal separation between the spikes is sharp.

We developed a JAX-based Python package implementing the proposed reconstruction method, which is publicly available at https://github.com/rpetit/pointscat. It can run seamlessly on CPUs, GPUs and TPUs, and benefits from just-in-time compilation.

**Structure of the paper.** In Section 2, we provide a brief introduction to the sparse spikes problem. Then, we give an overview of the approach we propose to solve the inverse problem in Section 3. Section 4 is dedicated to the derivation of theoretical guarantees for the linear part of our recovery procedure. Finally, we introduce our reconstruction algorithm in Section 5 and present numerical results on a few examples in Section 6.

## 1.4 Notations

**General notations.** Let $\mathcal{X}$ be a connected bounded open subset of $\mathbb{R}^d$ and denote the space of complex-valued Radon measures on $\mathcal{X}$ by $\mathcal{M}(\mathcal{X})$. The total variation of a measure $\nu \in \mathcal{M}(\mathcal{X})$ is denoted by $|\nu|(\mathcal{X})$. If $a \in \mathbb{C}^s$ and $1 \le i \le s$, we define

$$a_{-i} \stackrel{\text{def.}}{=} (a_1, \dots, a_{i-1}, a_{i+1}, \dots, a_s) \in \mathbb{C}^{s-1}.$$

When $x = (x_1, \dots, x_s) \in \mathcal{X}^s$, we sometimes use the notation $u^{\text{in}}(x) \stackrel{\text{def.}}{=} (u^{\text{in}}(x_i))_{1 \le i \le s}$. Finally, we use the notation $|\cdot|$ for the Euclidean norm in $\mathbb{R}^d$.

**Forward operators.** Given $s$ scatterers with intensities $a \in (\mathbb{C}^*)^s$ and locations $x \in \mathcal{X}^s$, we denote the far field pattern in the Born approximation by $u^{\infty, b}$. Recalling (1) and (5), the Born approximation $u \approx u^{\text{in}}(x)$ yields

$$u^{\infty, b}(\hat{x}, \theta) \stackrel{\text{def.}}{=} \frac{\kappa^2}{4\pi} \sum_{i=1}^{s} a_i e^{-i\kappa(\hat{x}_i - \theta) \cdot x_i}, \quad (\hat{x}, \theta) \in \mathbb{S}^{d-1} \times \mathbb{S}^{d-1}. \qquad (6)$$

Given $m$ incident directions $(\theta_k)_{1 \le k \le m}$ and corresponding observation directions $(\hat{x}_k)_{1 \le k \le m}$, we denote the associated forward operator on $\mathcal{M}(\mathcal{X})$ by $\Phi^b$, which is

$$\Phi^b: \mathcal{M}(\mathcal{X}) \to \mathbb{C}^m \\ \mu \mapsto \frac{1}{\sqrt{m}} \left[ \int_{\mathcal{X}} e^{-i\kappa(\hat{x}_k - \theta_k) \cdot y} d\mu(y) \right]_{k=1}^m.$$
Note that $\Phi^b$ is a linear map. Finally, we define the Born and Foldy forward operators associated to the arrangement of scatterers $x$:

$$ \Phi_x^b: \mathbb{C}^s \to \mathbb{C}^m \qquad \Phi_x^f: \mathbb{C}^s \to \mathbb{C}^m \\ a \mapsto \frac{1}{\sqrt{m}} [u^{\infty,b}(\hat{x}_k, \theta_k)]_{k=1}^m, \qquad a \mapsto \frac{1}{\sqrt{m}} [u^{\infty}(\hat{x}_k, \theta_k)]_{k=1}^m. $$

With this notation, taking $\mu = \sum_{i=1}^s a_i \delta_{x_i}$ yields $\Phi_x^b a = \Phi^b \mu$.

**Green functions.** We denote the Green function associated to the Helmholtz equation with wavenumber $\kappa$ by $G$, i.e.

$$ G(x, y) \stackrel{\text{def.}}{=} \begin{cases} \frac{i}{4} H_0^{(1)}(\kappa|x-y|) & \text{if } d=2, \\ \frac{1}{4\pi|x-y|} e^{i\kappa|x-y|} & \text{if } d=3, \end{cases} \quad x, y \in \mathbb{R}^d, $$

where $H_0^{(1)}$ denotes the Hankel function of the first kind and order zero. As we often need to upper bound the modulus of $G$, by [77, page 446] and [51, Theorem A] we have

$$ \forall x, y \in \mathbb{R}^d, |G(x, y)| \le \phi(|x - y|), \qquad (7) $$

where $\phi: \mathbb{R}_+^* \to \mathbb{R}$ is the non-increasing function given by

$$ \phi(t) \stackrel{\text{def.}}{=} \begin{cases} \frac{1}{4} \min \left( \sqrt{1 + \frac{4}{\pi^2} (\gamma + \log(\kappa t/2))^2}, \sqrt{\frac{2}{\pi \kappa t}} \right) & \text{if } d=2, \\ \frac{1}{4\pi t} & \text{if } d=3, \end{cases} \quad t \in \mathbb{R}_+^*, \qquad (8) $$

where $\gamma$ is Euler's constant. We stress that the function $\phi$ depends on $\kappa$ when $d=2$.

**Foldy matrix.** Given $a \in \mathbb{C}^s$ and $x \in (\mathbb{R}^d)^s$, we define the matrix $T_{a,x} \in \mathbb{C}^{s \times s}$ by

$$ (T_{a,x})_{ij} = \begin{cases} 0 & \text{if } i=j, \\ G(x_i, x_j) a_j & \text{otherwise,} \end{cases} \quad (i,j) \in \{1, \dots, s\}^2. $$

With this definition, the Foldy system (4) writes

$$ (Id - \kappa^2 T_{a,x}) u = u^{\text{in}}(x). \qquad (9) $$

## 2 Background on the sparse spikes problem

Since our approach heavily relies on existing results on the sparse spikes problem, in this section we provide a brief introduction to this topic. We also refer the reader to [41, 30, 63] and to the lecture notes [70].

**Problem formulation.** The sparse spikes problem consists in recovering a discrete measure $\mu = \sum_i a_i \delta_{x_i}$ from the knowledge of $y = \Phi\mu + w$, where $\Phi: \mathcal{M}(\mathcal{X}) \to \mathbb{C}^m$ is linear and $w \in \mathbb{C}^m$ is an additive noise. This task naturally appears in a wide range of applications. Let us mention the localization of point light sources in biological and astronomical imaging [43], the estimation of the parameters of mixture models in statistics [57, 29], and the training of two-layer neural networks in machine learning [11]. The most popular approaches for solving it are either of variational or of parametric type. In this work, we focus on the former. The latter is however briefly discussed below.

**Parametric approaches.** Broadly speaking, parametric approaches consist in encoding the locations of the atoms in the roots of some polynomial. The most popular are Prony's method, Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT), and the MUSIC algorithm. We refer the reader to [30, Section 0.2] for a detailed presentation and a comparison with variational approaches.
**The Beurling LASSO.** In this work, we focus on a variational approach. It consists in looking for an estimate of $\mu$ among the solutions of a convex optimization problem in the space of measures, called the Beurling LASSO or BLASSO. It was introduced in [39, 20], and writes

$$ \min_{\nu \in \mathcal{M}(\mathcal{X})} \frac{1}{2} \|\Phi_{\nu} - y\|_{2}^{2} + \lambda |\nu|(\mathcal{X}), \qquad (10) $$

where $\lambda$ is a regularizer parameter that should be chosen according to an estimate of the noise level $\|w\|_2$. This problem can be seen as an infinite-dimensional version of the celebrated LASSO, which is a popular method for the recovery of sparse finite-dimensional signals.

**Stable recovery.** Under suitable assumptions, the solutions of (10) can be guaranteed to be close to $\mu$. In [20, 23, 48, 10], some convergence guarantees were proved in the low noise regime. In [46, 42, 71], stronger exact support recovery results were derived. Outside the low noise regime, guarantees were also obtained in [74, 18, 21] under a Gaussian noise assumption.

**Numerical methods.** The algorithm we focus on for solving (10), introduced in [43] (see also the earlier works [20, 17]), is based on the Frank-Wolfe or conditional gradient algorithm. It is known as the *sliding Frank-Wolfe algorithm*. It consists in iteratively adding atoms to the support of the zero measure based on the residual observations, and comes with strong convergence guarantees. Let us mention that there exist other approaches for solving (10). These include conic particle gradient descent [35], projected gradient descent [14], proximal methods [76], exchange algorithms [49], as well as semidefinite relaxations [75, 24, 38, 31].

## 3 Overview of the proposed approach

The inverse problem we wish to solve aims to recover $a \in (\mathbb{C}^*)^s$ and $x \in \mathcal{X}^s$ from $y = \Phi_x^f a + w^{\text{noise}}$, where $w^{\text{noise}} \in \mathbb{C}^m$ is an additive noise. In this section, we describe our two-step “linearize and locally optimize” approach to this problem. We first focus on the first step, hereby called the *linear step*, and then turn to the second step, hereby called the *nonlinear step*.

**Linear step.** The key idea of our approach is to decompose $\Phi_x^f a$ by writing $\Phi_x^f a = \Phi_x^b a + w^{\text{lin}}$, where $w^{\text{lin}}$ is the linearization error. As a result, we can recast the recovery of $a$ and $x$ from $y$ as the recovery of the discrete measure $\mu = \sum_{i=1}^s a_i \delta_{x_i}$ from the noisy linear measurements $\Phi^b \mu + w$, where $w = w^{\text{lin}} + w^{\text{noise}}$:

$$ y = \Phi_x^f a + w^{\text{noise}} = \Phi_x^b a + w^{\text{lin}} + w^{\text{noise}} = \Phi_x^b \mu + w^{\text{lin}} + w^{\text{noise}}. $$

This is a particular instance of the sparse spikes problem. We hence propose to estimate $\mu$ by solving the BLASSO program (10) with linear forward map $\Phi = \Phi^b$:

$$ \min_{\nu \in \mathcal{M}(\mathcal{X})} \frac{1}{2} \|\Phi^{b}\nu - y\|_{2}^{2} + \lambda^{b}|\nu|(\mathcal{X}) \qquad (11) $$

for some regularization parameter $\lambda^b > 0$.

**Nonlinear step.** Approximating $\mu$ by solving (11) yields reconstructions whose quality deteriorates as the norm of the linearization error $\|w^{\text{lin}}\|_2$ grows. We hence propose to improve the estimate obtained in the linear step by performing a local minimization step involving the nonlinear forward operator $\Phi^f$. More precisely, we wish to locally optimize the non-convex functional $J^f$ defined by

$$ J^f: \mathbb{C}^s \times \mathcal{X}^s \to \mathbb{R} \\ (\tilde{a}, \tilde{x}) \mapsto \frac{1}{2} \|\Phi_x^f \tilde{a} - y\|_2^2 + \lambda^f \|\tilde{a}\|_1, \qquad (12) $$

by using the output of the linear step as initialization. Provided the linearization error is sufficiently small, we can hope that the linear estimate lies in a basin of attraction of the true intensities and locations.
# 4 Recovery guarantees

In this section, we first study the error $w^{\text{lin}} = \Phi_x^f a - \Phi_x^b a$ committed by linearizing the forward operator, i.e. the difference between the far field patterns associated to the Foldy model and to the Born approximation. We show that the error can be controlled provided the scatterers are sufficiently separated and their intensities are small. Then, we leverage existing results on the sparse spikes problem, and show that, in this regime, solving (10) allows us to stably recover $\mu$.

## 4.1 Bounding the linearization error

In this subsection, we study the linearization error $\|u^{\infty} - u^{\infty,b}\|_{\infty}$, see (5) and (6). Classical error bounds are derived in e.g. [68], which, however, does not treat the case of point-like scatterers. In short, we show that the linearization error is small as soon as $\kappa^2\phi(\Delta)\|a\|_1$ is small, where $\phi$ is given in (8) and

$$ \Delta = \min_{i \neq j} |x_i - x_j|, \qquad (13) $$

i.e. the scatterers are sufficiently separated and have small intensities.

We first investigate the special case of two scatterers.

**Proposition 4.1** Let $a \in (\mathbb{C}^*)^2$ and $x \in \mathcal{X}^2$. Define $\alpha = \kappa^2|G(x_1, x_2)|\sqrt{|a_1||a_2|}$. If $\alpha < 1$, then we have

$$ \|u^{\infty} - u^{\infty,b}\|_{\infty} \le \frac{\kappa^2}{4\pi} \frac{2\alpha}{1-\alpha^2} \left( \alpha \frac{|a_1| + |a_2|}{2} + \sqrt{|a_1||a_2|} \right). \qquad (14) $$

*Proof.* For $s = 2$ scatterers, the Foldy matrix $Id - \kappa^2 T_{a,x}$ can be explicitly inverted. Indeed, we have that $\det(\text{Id} - \kappa^2 T_{a,x}) = 1 - \beta^2$ with

$$ \beta \stackrel{\text{def.}}{=} \kappa^2 G(x_1, x_2) \sqrt{|a_1||a_2|} e^{i(\varphi_1 + \varphi_2)/2} \text{ and } a_i = |a_i| e^{i\varphi_i}. $$

Therefore, provided $\beta^2 \neq 1$, we obtain that $Id - \kappa^2 T_{a,x}$ is invertible and

$$ (\text{Id} - \kappa^2 T_{a,x})^{-1} = \frac{1}{1 - \beta^2} \begin{pmatrix} 1 & \kappa^2 G(x_1, x_2) a_2 \\ \kappa^2 G(x_1, x_2) a_1 & 1 \end{pmatrix}. $$

In view of (1), (5) and (9), this yields

$$ u^{\infty}(\hat{x}, \theta) = \frac{\kappa^2}{4\pi} \left( a_1 e^{-i\kappa\hat{x}\cdot x_1} \frac{e^{i\kappa\theta\cdot x_1} + \kappa^2 G(x_1, x_2) a_2 e^{i\kappa\theta\cdot x_2}}{1 - \beta^2} + a_2 e^{-i\kappa\hat{x}\cdot x_2} \frac{e^{i\kappa\theta\cdot x_2} + \kappa^2 G(x_1, x_2) a_1 e^{i\kappa\theta\cdot x_1}}{1 - \beta^2} \right), $$

which, by (6), finally allows us to get the following bound:

$$ \|u^{\infty} - u^{\infty,b}\|_{\infty} \le \frac{\kappa^2}{4\pi} \frac{2|\beta|}{|1 - \beta^2|} \left( |\beta| \frac{|a_1| + |a_2|}{2} + \sqrt{|a_1||a_2|} \right). \qquad (15) $$

When $\alpha = |\beta| < 1$, we finally obtain (14). □

Next, we turn to the general case, in which we provide two different bounds. The second one is finer and explicitly depends on the strength of the pairwise interactions between the scatterers.

**Proposition 4.2** Let $a \in (\mathbb{C}^*)^s$ and $x \in \mathcal{X}^s$. Define $\Delta \stackrel{\text{def.}}{=} \min_{i \neq j} |x_i - x_j|$ and $\alpha \stackrel{\text{def.}}{=} \kappa^2\phi(\Delta) \max_i \|a_{-i}\|_1$. If $\alpha < 1$, then we have

$$ \|u^{\infty} - u^{\infty,b}\|_{\infty} \le \frac{\kappa^2}{4\pi} \|a\|_1 \frac{\alpha}{1 - \alpha} \qquad (16) $$
and

$$ \|u^{\infty} - u^{\infty,b}\|_{\infty} \le \frac{\kappa^2}{4\pi} \left( \|a\|_1 \frac{\alpha^2}{1-\alpha} + 2 \sum_{i<j} |a_i||a_j| \kappa^2 \phi(|x_i - x_j|) \right). \quad (17) $$

*Proof.* The infinity operator norm of $T_{a,x}$, defined by

$$ \|T_{a,x}\|_{\infty} \stackrel{\text{def}}{=} \sup \{ \|T_{a,x}u\|_{\infty} : \|u\|_{\infty} = 1 \}, $$

can be bounded as follows (see (7) and (13)):

$$ \|T_{a,x}\|_{\infty} = \max_{1 \le i \le s} \sum_{\substack{j=1 \\ j \ne i}}^{s} |G(x_i, x_j)| |a_j| \le \phi(\Delta) \max_{1 \le i \le s} \|a_{-i}\|_1 \le \phi(\Delta) \|a\|_1. $$

Provided $\kappa^2 \|T_{a,x}\|_{\infty} < 1$, by (1), (5), (6), and (9) this yields

$$ \begin{aligned} |u^{\infty}(\hat{x}, \theta) - u^{\infty,b}(\hat{x}, \theta)| &= \frac{\kappa^2}{4\pi} \left| \sum_{i=1}^{s} a_i e^{-i\kappa\hat{x}\cdot x_i} \left[ \left( (Id - \kappa^2 T_{a,x})^{-1} - Id \right) u^{\text{in}}(x) \right]_i \right| \\ &\le \frac{\kappa^2}{4\pi} \|a\|_1 \|(Id - \kappa^2 T_{a,x})^{-1} - Id\|_{\infty} \\ &\le \frac{\kappa^2}{4\pi} \|a\|_1 \left\| \sum_{n=1}^{+\infty} \kappa^{2n} T_{a,x}^n \right\|_{\infty} \\ &\le \frac{\kappa^2}{4\pi} \|a\|_1 \sum_{n=1}^{+\infty} (\kappa^2 \|T_{a,x}\|_{\infty})^n \\ &= \frac{\kappa^2}{4\pi} \|a\|_1 \frac{\kappa^2 \|T_{a,x}\|_{\infty}}{1 - \kappa^2 \|T_{a,x}\|_{\infty}}, \end{aligned} $$

where we have used that $\|u^{\text{in}}(x)\|_{\infty} \le 1$. We therefore obtain (16) provided $\alpha < 1$.

Now, we turn to the proof of (17). When $s > 2$, obtaining a simple expression of the powers of $T_{a,x}$ appears difficult. However, the contribution of first order interactions between pairs of scatterers can be treated separately from higher order interactions. This is done as follows:

$$ \begin{aligned} |u^{\infty}(\hat{x}, \theta) - u^{\infty,b}(\hat{x}, \theta)| &= \frac{\kappa^2}{4\pi} \left| \sum_{i=1}^{s} a_i e^{-i\kappa\hat{x}\cdot x_i} \left[ \left( (Id - \kappa^2 T_{a,x})^{-1} - Id \right) u^{\text{in}}(x) \right]_i \right| \\ &= \frac{\kappa^2}{4\pi} \left| \sum_{i=1}^{s} a_i e^{-i\kappa\hat{x}\cdot x_i} \left[ \left( \kappa^2 T_{a,x} + \sum_{n=2}^{+\infty} \kappa^{2n} T_{a,x}^n \right) u^{\text{in}}(x) \right]_i \right| \\ &\le \frac{\kappa^2}{4\pi} \left| \sum_{i=1}^{s} a_i e^{-i\kappa\hat{x}\cdot x_i} \left[ \kappa^2 T_{a,x} u^{\text{in}}(x) \right]_i \right| + \left| \frac{\kappa^2}{4\pi} \left| \sum_{i=1}^{s} a_i e^{-i\kappa\hat{x}\cdot x_i} \left[ \left( \sum_{n=2}^{+\infty} \kappa^{2n} T_{a,x}^n \right) u^{\text{in}}(x) \right]_i \right| \right| \end{aligned} $$

We bound the first term as follows:

$$ \left| \sum_{i=1}^{s} a_i e^{-i\kappa\hat{x}\cdot x_i} \left[ \kappa^2 T_{a,x} u^{\text{in}}(x) \right]_i \right| = \kappa^2 \left| \sum_{i=1}^{s} \sum_{\substack{j=1 \\ j \ne i}}^{s} a_i a_j G(x_i, x_j) e^{-i\kappa(\hat{x}\cdot x_i - \theta \cdot x_j)} \right| \le 2 \sum_{i<j} |a_i||a_j| \kappa^2 \phi(|x_i - x_j|), $$

and treat the second term as above to obtain

$$ \left| \sum_{i=1}^{s} a_i e^{-i\kappa\hat{x}\cdot x_i} \left[ \left( \sum_{n=2}^{+\infty} \kappa^{2n} T_{a,x}^n \right) u^{\text{in}} \right]_i \right| \le \|a\|_1 \frac{\alpha^2}{1-\alpha}. $$

This finally yields (17). □
<figure><img src="image_1.png" /><figcaption>Figure 1: Dependence of the linearization error on $\kappa$ and $\Delta$. The dashed line corresponds to the bound (14), and the plain line to the empirical error.</figcaption></figure>

We notice that, in the case $s = 2$ and $a_1 = a_2$, bounding $|G(x_1, x_2)|$ by $\phi(\Delta)$ in (16) yields (14). However, (14) is significantly better when the two amplitudes are different. The second bound (17), which explicitly depends on the strength of the pairwise interactions between the scatterers, allows us to bridge this gap. It is finer than (16), since $\frac{\alpha}{1-\alpha} = \frac{\alpha^2}{1-\alpha} + \alpha$ and

$$ 
\begin{aligned}
2 \sum_{i<j} |a_i||a_j| \kappa^2 \phi(|x_i - x_j|) &= \sum_{i=1}^{s} |a_i| \sum_{\substack{j=1 \\ j \neq i}}^{s} |a_j| \kappa^2 \phi(|x_i - x_j|) \\
&\leq \kappa^2 \phi(\Delta) \sum_{i=1}^{s} |a_i| \sum_{\substack{j=1 \\ j \neq i}}^{s} |a_j| \\
&\leq \kappa^2 \phi(\Delta) \|a\|_1 \max_i \|a_{-i}\|_1 = \alpha \|a\|_1.
\end{aligned}
$$

**Numerical experiments.** We end this section by studying how tight our bounds are in practice. We are especially interested in knowing whether the true error has the same dependence on $\kappa^2\phi(\Delta)$. Since (16) and (17) reduce to (14) when $s = 2$ and $a_1 = a_2$, we focus on the study of (14). We take $d = 2$ and consider 100 uniformly sampled incident and observation directions. We compute the linearization error $\|\Phi_x^f a - \Phi_x^b a\|_2$ where $a = (1, 1)$ and $x = (x_1, x_2)$ with $|x_1 - x_2| = \Delta$ for multiple values of $\Delta$ and $\kappa \in \{0.1, 1, 10\}$. For each value of $\Delta$, we average the errors corresponding to 20 configurations obtained as follows: we draw $x_1$ uniformly in $[-1, 1]^2$, then draw $e$ uniformly in $\mathbb{S}^1$ and define $x_2 \stackrel{\text{def.}}{=} x_1 + \Delta e$. The results are presented in Figure 1. We choose not to display error bars for the empirical error since we observe its standard deviation to be very small. We see that (14) is almost tight when $\alpha$ is sufficiently bigger than 1, but deteriorates as $\alpha$ approaches 1. One could also wonder how (15) behaves outside the regime $|\beta| < 1$. However, Figure 2 suggests that, in this case, $\|\Phi_x^f a - \Phi_x^b a\|_2$ is almost equal to $\|\Phi_x^f a\|_2 + \|\Phi_x^b a\|_2$, so that even if (15) remains valid, this regime makes the linearization approach unsuitable for solving the inverse problem.

## 4.2 Stable reconstruction

We now provide guarantees that solutions of (10) are close to $\mu$, provided a minimal separation condition holds and the sum $w$ of the linearization error and the noise is small. Our results are direct consequences of existing results on the sparse spikes problem (specifically, [72, Theorem 3] and [71, Theorem 1]). We focus on the “off-the-grid compressed sensing” line of works [75, 72], which exhibit conditions under which a discrete
<figure><img src="image_2.png" /><figcaption>Figure 2: True linearization error $||\Phi_x^f a - \Phi_x^b a||_2$ and naive upper bound $||\Phi_x^f a||_2 + ||\Phi_x^b a||_2$ for $\kappa = 1$. On the left of the vertical black line we have $|\beta| = \alpha > 1$, and on the right $|\beta| = \alpha < 1$.</figcaption></figure>

measure can be stably recovered from randomly sampled measurements. They consider forward operators of the form<fnref n="2" />

$$ \mathcal{M}(\mathcal{X}) \to \mathbb{C}^m \\ \mu \mapsto \frac{1}{\sqrt{m}} \left[ \int_{\mathcal{X}} \varphi_{\omega_k} d\mu \right]_{k=1}^m \qquad (18) $$

where $\omega_1, ..., \omega_m$ are independent and identically distributed according to a certain probability distribution $\Lambda$ and $\varphi_\omega$ are continuous functions on $\mathcal{X}$.

Recall that the far field pattern in the Born approximation (6) gives access to the Fourier transform of $\mu = \sum_{i=1}^s a_i \delta_{x_i}$ in $\kappa \mathbb{S}^{d-1} - \kappa \mathbb{S}^{d-1} = \bar{B}(0, 2\kappa)$, with $B(0, 2\kappa)$ the open ball of radius $2\kappa$ centered at 0. Thus, we take $\varphi_\omega \stackrel{\text{def.}}{=} e^{-i\langle\omega, \cdot\rangle}$ and we let $\Lambda$ be the uniform law on $B(0, 2\kappa)$.

**Sampling of the incident and observation directions.** Now, we explain how to sample the incident and observation directions in order to obtain uniform Fourier samples in $B(0, 2\kappa)$. For simplicity, we consider the case $d = 2$. From (6), we have that, for an incident direction $\theta \in \mathbb{S}^1$ and an observation direction $\hat{x} \in \mathbb{S}^1$, the far field pattern in the Born approximation is proportional to the Fourier transform of the unknown measure at the frequency $\kappa(\hat{x} - \theta)$. We hence want to find $\pi_{(\hat{x}, \theta)}$ such that, if $(\hat{X}, \Theta) \sim \pi_{(\hat{x}, \theta)}$, then $\kappa(\hat{X} - \Theta) \sim \pi_\omega$ with $\pi_\omega = \mathcal{U}(B(0, 2\kappa))$. In order to do this, let us identify $\mathbb{R}^2$ with $\mathbb{C}$ and define

$$ \begin{gathered} f: B(0, 2\kappa) \to \mathbb{S}^1 \\ g: B(0, 2\kappa) \to \mathbb{S}^1 \\ \omega = re^{i\theta} \mapsto e^{i(\theta+\arccos(\frac{r}{2\kappa}))}, \qquad \omega = re^{i\theta} \mapsto -e^{i(\theta-\arccos(\frac{r}{2\kappa}))}. \end{gathered} $$

A simple computation shows that, for any $\omega \in B(0, 2\kappa)$, we have $\omega = \kappa(f(\omega) - g(\omega))$. As a consequence, we obtain that taking $\pi_{(\hat{x}, \theta)} = (f, g)_{\#} \pi_\omega$ yields $\kappa(\hat{X} - \Theta) \sim \pi_\omega$, where $(f, g)_{\#} \pi_\omega$ denotes the pushforward of $\pi_\omega$ by $(f, g)$. In other words, one can sample $\omega_1, ..., \omega_m$ uniformly in $B(0, 2\kappa)$ and take $\hat{x}_k = f(\omega_k)$ and $\theta_k = g(\omega_k)$ for all $k \in \{1, ..., m\}$.

**Stable recovery.** We first focus on the application of [72, Theorem 3]. This result shows the stability of the reconstruction in terms of a partial optimal transport distance $\mathcal{T}_\delta$, which allows us to compare measures with different masses. It is defined, for any pair of positive Radon measures $\nu_1, \nu_2 \in \mathcal{M}_+(\mathcal{X})$, by

$$ \mathcal{T}_{\delta}(\nu_1, \nu_2)^2 \stackrel{\text{def.}}{=} \inf_{\tilde{\nu}_1, \tilde{\nu}_2} \delta W_2^2(\tilde{\nu}_1, \tilde{\nu}_2) + |\nu_1 - \tilde{\nu}_1|(\mathcal{X}) + |\nu_2 - \tilde{\nu}_2|(\mathcal{X}), $$

<footnote n="2">²Recall that $\mathcal{X}$ is a connected bounded open subset of $\mathbb{R}^d$.</footnote>
where $\delta \stackrel{\text{def.}}{=} 2\kappa/\sqrt{d+2} \sim \kappa$, $W_2$ denotes the 2-Wasserstein distance and the infimum is taken over all measures $\tilde{\nu}_1, \tilde{\nu}_2$ satisfying $|\tilde{\nu}_1|(\mathcal{X}) = |\tilde{\nu}_2|(\mathcal{X})$. We refer the reader to e.g. [73, 69] for more details on optimal transport and its extensions.

**Theorem 4.3** Let $a \in (\mathbb{C}^*)^s$ and $x \in \mathcal{X}^s$. Assume that $(\hat{x}_k, \theta_k)_{1\le k\le m}$ are sampled so that $(\kappa(\hat{x}_k - \theta_k))_{1\le k\le m}$ are i.i.d. samples from $\mathcal{U}(B(0, 2\kappa))$, and that $y = (u^{\infty,b}(\hat{x}_k, \theta_k) + w_k)_{1\le k\le m}$ with $w \in \mathbb{C}^m$. Then there exists $r \sim 1/\kappa$ such that, assuming

$$m \gtrsim s \left( \log(s) \log\left(\frac{s}{\rho}\right) + \log\left(\frac{s\kappa}{\rho}\right) \right)$$

for some $\rho \in (0, 1)$ and

$$\min_{i \neq j} |x_i - x_j| \gtrsim \frac{s^{2/(d+1)}}{\kappa},$$

with probability at least $1 - \rho$, every solution $\hat{\mu}$ of (11) with $\lambda^b \sim \|w\|_2/\sqrt{s}$ is $\sqrt{s}\|w\|_2$-close to $\sum_{i=1}^s a_i \delta_{x_i}$, i.e.

$$\mathcal{T}_{\delta}^{2} \left( \sum_{i=1}^{s} \hat{A}_{i} \delta_{x_{i}}, |\hat{\mu}| \right) \lesssim \sqrt{s} \|w\|_{2} \quad \text{and} \quad \max_{1 \le i \le s} |\hat{a}_{i} - a_{i}| \lesssim \sqrt{s} \|w\|_{2}, \qquad (19)$$

where $\hat{A}_i = |\hat{\mu}|(B(x_i, r))$ and $\hat{a}_i = \hat{\mu}(B(x_i, r))$.

Let us stress that (19) tells that the support of $\hat{\mu}$ is essentially located on small neighborhoods of the true locations $\{x_i\}_{1\le i\le s}$, and that, in each of these neighborhoods, the average value $\hat{a}_i$ of $\hat{\mu}$ is close to the true amplitude $a_i$. The number of scalar measurements $m$ is, up to logarithmic factors, linearly proportional to the number of spikes $s$, which is the best one can hope for.

**Exact support recovery.** In their earlier work [71], the authors of [72] showed stronger exact support recovery results under more restrictive assumptions. Their application to our problem is the object of Theorem 4.4. Contrary to Theorem 4.3, this result does not hold for arbitrary $w$. It holds only if $\|w\|_2$ is sufficiently small, in terms of the sparsity $s$ and the intensities $(a_i)_{1\le i\le s}$, and with a larger number of measurements $m$.

**Theorem 4.4** Let $a \in (\mathbb{C}^*)^s$, $x \in \mathcal{X}^s$ and define $a \stackrel{\text{def.}}{=} \min_{1\le i\le s}\{|a_i|^2, |a_i|^{-2}\}$ and $D \stackrel{\text{def.}}{=} a \min(\sqrt{s}, \sqrt{s}/\|a\|_2, 1)$. Assume that $(\hat{x}_k, \theta_k)_{1\le k\le m}$ are sampled so that $(\kappa(\hat{x}_k - \theta_k))_{1\le k\le m}$ are i.i.d. samples from $\mathcal{U}(B(0, 2\kappa))$, and that $y = (u^{\infty,b}(\hat{x}_k, \theta_k) + w_k)_{1\le k\le m}$ with $w \in \mathbb{C}^m$. If

$$m \gtrsim s^{3/2} \log(\kappa/\rho), \quad \lambda^b \lesssim D/s, \quad \|w\|_2 \lesssim \lambda^b,$$

for some $\rho \in (0, 1)$ and

$$\min_{i \neq j} |x_i - x_j| \gtrsim \frac{s^{2/(d+1)}}{\kappa},$$

then, with probability at least $1 - \rho$, the solution $\hat{\mu}$ of (11) is unique and made of exactly $s$ spikes. Moreover, writing $\hat{\mu} = \sum_{i=1}^s \hat{a}_i \delta_{\hat{x}_i}$, we have:

$$\|\hat{a} - a\|_2 + \kappa \|\hat{x} - x\|_2 \lesssim \frac{\sqrt{s}(\lambda^b + \|w\|_2)}{\min_{1 \le i \le s} |a_i|}.$$

The proof of Theorems 4.3 and 4.4 is postponed to Appendix A. These results are a consequence of [72, Theorem 3] and [71, Theorem 1]. We prove their applicability in the case of continuous Fourier sampling with frequencies uniformly distributed in a Euclidean ball of fixed radius.

We end this section with two remarks on the above results.
Remark 4.5. Our bound on $||u^{\infty} - u^{\infty,b}||_{\infty}$ (see Section 4.1) yields a bound on $||w^{\text{lin}}||_2$, which can in turn be used in Theorems 4.3 and 4.4 by taking $w = w^{\text{lin}} + w^{\text{noise}}$. Indeed, we have

$$ \begin{aligned} \|w^{\text{lin}}\|_2 &= \|\Phi_x^b a - \Phi_x^f a\|_2 = \sqrt{\sum_{k=1}^m \left( \frac{1}{\sqrt{m}} u^{\infty}(\hat{x}_k, \theta_k) - \frac{1}{\sqrt{m}} u^{\infty,b}(\hat{x}_k, \theta_k) \right)^2} \\ &= \sqrt{\frac{1}{m} \sum_{k=1}^m |u^{\infty}(\hat{x}_k, \theta_k) - u^{\infty,b}(\hat{x}_k, \theta_k)|^2} \le \|u^{\infty} - u^{\infty,b}\|_{\infty}. \end{aligned} $$

Remark 4.6. One could wonder how changing the sampling scheme $\Lambda$ would impact the results of Theorems 4.3 and 4.4. Owing to the proofs of [72, Theorem 3] and [71, Theorem 1], we only expect this choice to impact the minimal separation distance condition (at least for reasonable choices of sampling schemes). In fact, we believe that sampling the lower frequencies more often than the higher frequencies (instead of sampling them uniformly as above) would yield a kernel with a faster decay (see Appendix A). Given that $d = 2$ or $3$ in our case, this could allow for a minimal separation of the form

$$ \min_{i \neq j} |x_i - x_j| \gtrsim \frac{1}{\kappa}, $$

without any dependence on the sparsity $s$ (see [72, Section 2.1 and Lemma 17]).

# 5 Reconstruction algorithm

**Description.** As discussed in Section 3, the algorithm we propose to recover $\mu$ from $y$ is made of two main steps. The first one consists in solving (11) using the sliding Frank-Wolfe algorithm introduced in [43]. The second one consists in locally minimizing (12) using the output of the first step as initialization. This procedure is precisely described in Algorithm 1.

**Stopping condition and convergence of the linear step.** The condition of Line 6 is the standard stopping criterion for the sliding Frank-Wolfe algorithm (see [43]). In particular, it guarantees that the estimated measure is at most $\epsilon$-suboptimal in terms of the objective value for (11). By [43, Lemma 5, Proposition 5], we have that the sequence of iterates produced by the linear step (i.e. without performing Lines 6-10) converges in objective value towards the value of (11) at a rate $O(1/n)$. Moreover, it has an accumulation point for the weak-* topology on $M(\mathcal{X})$, which is a solution of (11). We also stress that, under a non-degeneracy assumption, Theorem 3 in the above reference yields the weak-* convergence of the sequence towards the unique minimizer of (11) in a finite number of steps.

**Implementation.** The implementation of Algorithm 1 requires several oracles to carry out the operations of Lines 5, 7, 13 and 14. We refer the reader to [43, Remark 9] regarding the operations of Lines 5, 13 and 14. The local minimization of Line 7 is performed as in Line 14. To solve the various optimization problems appearing in Algorithm 1, we rely on the Python libraries JAX [19] and JAXopt [16].

**Hyperparameters.** Algorithm 1 involves several hyperparameters. The main ones are the regularization parameters for the linear and nonlinear steps $\lambda^b$ and $\lambda^f$, as well as the stopping criterion tolerance $\epsilon$. When carrying out Lines 9 and 16, we also need to define two thresholds: one below which the amplitudes are set to zero, and one below which neighboring spikes are merged$^3$. The most critical parameters are $\lambda^b$ and $\lambda^f$. In the literature on variational regularization, they are usually selected after performing a grid search, or using more advanced techniques under statistical assumptions on the noise. We leave the design of rules for automatically choosing them to future works, and provide in Section 6 the set of parameters used for each

<footnote>$^3$We stress that, due to the linearity of $\Phi^b$ and the nonlinearity of $\Phi^f$, we merge neighboring spikes only in the linear step.</footnote>
<caption>**Algorithm 1:** Reconstruction of a discrete measure from far field measurements</caption>

**Data:** measurements $y$, regularization parameters $(\lambda^b, \lambda^f)$, stopping criterion tolerance $\epsilon$
**Result:** estimated measure $\hat{\mu}$
1 & $k \leftarrow 0$
2 & $\mu^{[0]} \leftarrow 0$
3 & **while** true **do**
4 & & $\eta^{[n]} \leftarrow (\Phi^b)^*(y - \Phi^b \mu^{[n]})$
5 & & $x_* \leftarrow \text{Argmax}_{x \in \mathcal{X}} |\eta^{[n]}(x)|$ // new atom finding
6 & & **if** $|\eta^{[n]}(x_*)| \le \lambda^b(1 + \epsilon)$ and $n \ge 1$ **then**
7 & & & & // nonlinear sliding step
8 & & & & find a critical point $(a^*, x^*)$ of $(a, x) \mapsto \frac{1}{2}\|\Phi_x^f a - y\|_2^2 + \lambda^f \|a\|_1$ by performing a local descent
9 & & & & initialized with $(a^{[n]}, x^{[n]})$
10 & & & & $\hat{\mu} \leftarrow \sum_i a_i^* \delta_{x_i^*}$
11 & & & & remove atoms with zero amplitude
12 & & & & output $\hat{\mu}$
13 & & **else**
14 & & $x^{[n+1/2]} \leftarrow (x^{[n]}, x_*)$
15 & & $a^{[n+1/2]} \leftarrow \text{Argmin}_a \frac{1}{2}\|\Phi_x^b a - y\|_2^2 + \lambda^b \|a\|_1$ // weights update
16 & & & & // linear sliding step
17 & & & & find a critical point $(a^{[n+1]}, x^{[n+1]})$ of $(a, x) \mapsto \frac{1}{2}\|\Phi_x^b a - y\|_2^2 + \lambda^b \|a\|_1$ by performing a local descent initialized with $(a^{[n+1/2]}, x^{[n+1/2]})$
18 & & & & $\mu^{[n+1]} \leftarrow \sum_i a_i^{[n+1]} \delta_{x_i^{[n+1]}}$
19 & & & & remove atoms with zero amplitude and merge spikes with same locations
20 & & & & $n \leftarrow n + 1$
21 & & **end**
22 & & **end**
23 & &

experiment. Let us point out that $\lambda^b$ should be chosen according to an estimate of the linearization error, which could be difficult to know a priori in practice.

# 6 Numerical results

In this section, we investigate the performance of Algorithm 1. We consider the two-dimensional case. For visualization purposes, we restrict ourselves to the case of real amplitudes, i.e. $a \in \mathbb{R}^s$. All the experiments below can be reproduced using the code available at https://github.com/rpetit/pointscat. The domain $\mathcal{X}$ is taken to be $(-r/2, r/2)^2$ with $r = 5$ in the first setting and $r = 10$ in the second one.

**Two scatterers.** In Figure 3, we present the output of Algorithm 1 in a first simple setting. We refer to the output of Algorithm 1 without performing Lines 6-10 as the *linear estimate* and to the output of the full algorithm as the *nonlinear estimate*. The unknown measure is $\mu = \delta_{x_1} + \delta_{x_2}$ with $x_1 = (-\Delta/2, 0)$ and $x_2 = (\Delta/2, 0)$, the wave number is $\kappa = 1$, and the number of randomly sampled frequencies is $m = 20$. In Figure 3(a), there is no noise on the observations and we fix $\lambda^b = 0.5$ and $\lambda^f = 10^{-3}$. In Figure 3(b), the measurement noise is Gaussian and has standard deviation 0.1 (corresponding to a relative $\ell_2$ noise level of 11%, 11%, 10%, 8%, 6% respectively) and we fix $\lambda^b = 1.0$ and $\lambda^f = 0.1$. We observe that, in both cases, the quality of the linear and nonlinear estimates deteriorates as $\Delta$ decreases. This is expected from the results of Section 4: a smaller $\Delta$ calls for a larger $\kappa$ (Theorems 4.3 and 4.4), yielding in turn a larger
<figure><img src="image_3.png" /><figcaption>(b) Gaussian measurement noise with standard deviation 0.1.</figcaption></figure>

<caption>Figure 3: Linear and nonlinear estimates for two unknown scatterers with several separation distances $Δ$, obtained with $m = 20$ measurements.</caption>
linearization error (Section 4.1), also because $\phi$ is non-increasing. In some cases, the nonlinear step allows us to significantly improve the quality of the reconstruction. For very close scatterers, neither the linear nor the nonlinear estimation allow us to accurately recover the unknown measure.

**A more challenging example.** In Figure 4, we show the output of Algorithm 1 for a measure $\mu$ made of 9 spikes and $m = 100$ sampled frequencies. The measurement noise is Gaussian and has standard deviation 0.1 (corresponding to a relative $\ell_2$ noise level of 6%). This setting is already highly nonlinear, since the relative $\ell_2$ error between the noiseless measurements and those corresponding to the Born approximation is of 24%. Yet, the linear estimate is close to the unknown measure and the nonlinear step allows further improvement in the quality of the reconstruction.

<figure><img src="image_4.png" /><figcaption>Figure 4: Linear and nonlinear estimates for nine unknown scatterers with $m = 100$ measurements and Gaussian measurement noise with standard deviation 0.1.</figcaption></figure>

**Influence of the initialization of the nonlinear step.** As a final experiment, we investigate the outcome of the nonlinear step if, instead of being initialized with the linear estimate, we initialize it with a discrete measure whose support is a uniform discretization of the domain. We use the same experimental setup as in Figure 4. In Figure 5, we show that even a rather fine initialization might fail to recover one of the scatterers. Still, increasing the number of initial spikes allows for an accurate recovery. This experiment suggests to further analyze the nonlinear dynamic, and in particular the existence of spurious critical points.

# 7 Conclusion

In this paper, we proposed to connect the inverse scattering problem associated to pointwise inhomogeneities with the sparse spikes problem. We leveraged existing results on the sparse spikes problem to obtain guarantees for a reconstruction method involving the linearized forward operator. We showed that the error is controlled by the strength of the scatterers and the minimal separation distance between them. Our analysis suggests the existence of a tradeoff in the choice of the incident frequency: a low frequency allows us to reduce the linearization error but imposes a larger separation distance between the scatterers. Finally, we proposed a refinement step involving the nonlinear forward operator, and provided numerical evidence of improvements in the reconstructions with respect to the linear estimation procedure.
<figure><img src="image_5.png" /><figcaption>Figure 5: Output of the nonlinear step initialized with a measure supported on a 4 × 4 (top row) and a 5 × 5 (bottom row) regular grid. The experimental setup is the same as in Figure 4.</figcaption></figure>
A natural question that we left open is to understand when the nonlinear step initialized with the output of the linear step allows us to recover the unknown accurately. One could start by investigating a simpler setting, where the scatterers are constrained to lie on a uniform discretization of the domain.

Future works could also include taking into account higher order terms in the first estimation step, for example by using the lifting and relaxation strategies developed for quadratic inverse problems (see for instance [25]). This would yield an estimation error of the form $\alpha^2/(1-\alpha)$ instead of $\alpha/(1-\alpha)$ as in (16), but would not allow to lift the requirement $\alpha < 1$. Another possible direction could be to rely on consecutive approximations. By inserting the estimates of $a$ and $x$ obtained in the linear step in the system (4), its solution $u$ (instead of $u^{in}(x)$ in the Born approximation) could be in turn inserted in (5), preserving the linearity of the forward map. Repeating this procedure could lead to more accurate estimates.

Finally, one could also wish to consider a more realistic measurement acquisition scenario, i.e. observations of the form $u^{\infty}(\hat{x}_{k_1}, \theta_{k_2})$ with $1 \le k_1 \le m_1$ and $1 \le k_2 \le m_2$ instead of $u^{\infty}(\hat{x}_k, \theta_k)$ with $1 \le k \le m$ [52, 1]. This would correspond to fixing the observation directions (and hence the set of sensing devices used), and, for every incident wave sent in the medium, acquire a full set of measurements. Our reconstruction algorithm applies to this setting, but the results of Section 4.2 correspond instead to acquiring a single measurement for every incident wave. We stress that treating the former scenario would require the development of different tools, as the corresponding frequencies would no longer be independent and identically distributed. In a similar spirit, it would be interesting to investigate the problem when the Foldy-Lax system is replaced by more advanced models for the scattering by small inhomogeneities [27].

## Acknowledgments

The authors thank Mourad Sini for valuable discussions about this work. Co-funded by the European Union (ERC, SAMPDE, 101041040). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. GSA and MS are members of the “Gruppo Nazionale per l’Analisi Matematica, la Probabilità e le loro Applicazioni”, of the “Istituto Nazionale di Alta Matematica”. The research was supported in part by the MIUR Excellence Department Project awarded to Dipartimento di Matematica, Università di Genova, CUP D33C23001110001. Funded by European Union – Next Generation EU.

## References

[1] Giovanni S. Alberti and Matteo Santacesaria. Infinite-Dimensional Inverse Problems with Finite Measurements. *Archive for Rational Mechanics and Analysis*, 243(1):1–31, January 2022.

[2] Sergio Albeverio, Friedrich Gesztesy, Raphael Hoegh-Krohn, and Helge Holden. *Solvable Models in Quantum Mechanics*. Springer Science & Business Media, December 2012.

[3] Abdulaziz Alsenafi, Ahcene Ghandriche, and Mourad Sini. The Foldy-Lax approximation is valid for nearly resonating frequencies. *Zeitschrift für angewandte Mathematik und Physik*, 74(1):11, December 2022.

[4] Habib Ammari, Durga Prasad Challa, Anupam Pal Choudhury, and Mourad Sini. The point-interaction approximation for the fields generated by contrasted bubbles at arbitrary fixed frequencies. *Journal of Differential Equations*, 267(4):2104–2191, August 2019.

[5] Habib Ammari, Durga Prasad Challa, Anupam Pal Choudhury, and Mourad Sini. The Equivalent Media Generated by Bubbles of High Contrasts: Volumetric Metamaterials and Metasurfaces. *Multiscale Modeling & Simulation*, 18(1):240–293, January 2020.
[6] Habib Ammari, Alexander Dabrowski, Brian Fitzpatrick, Pierre Millien, and Mourad Sini. Subwave-
length resonant dielectric nanoparticles with high refractive indices. *Mathematical Methods in the Ap-
plied Sciences*, 42(18):6567–6579, 2019.

[7] Habib Ammari, Brian Fitzpatrick, David Gontier, Hyundae Lee, and Hai Zhang. Minnaert resonances for acoustic waves in bubbly media. *Annales de l'Institut Henri Poincaré, Analyse non linéaire*, 35(7):1975–1998, November 2018.

[8] Habib Ammari and Hyeonbae Kang. *Reconstruction of Small Inhomogeneities from Boundary Measurements*. Springer Science & Business Media, August 2004.

[9] Habib Ammari, Bowen Li, and Jun Zou. Mathematical analysis of electromagnetic scattering by di-
electric nanoparticles with high refractive indices. *Transactions of the American Mathematical Society*,
376(01):39–90, January 2023.

[10] Jean-Marc Azaïs, Yohann de Castro, and Fabrice Gamboa. Spike detection from inaccurate samplings.
*Applied and Computational Harmonic Analysis*, 38(2):177–195, March 2015.

[11] Francis Bach. Breaking the Curse of Dimensionality with Convex Neural Networks. *Journal of Machine Learning Research*, 18(19):1–53, 2017.

[12] A. B. Bakushinsky and M. Yu Kokurin. *Iterative Methods for Approximate Solution of Inverse Problems*. Springer Science & Business Media, February 2005.

[13] Cédric Bellis, Marc Bonnet, and Fioralba Cakoni. Acoustic inverse scattering using topological derivative of far-field measurements-based L2 cost functionals. *Inverse Problems*, 29(7):075012, June 2013.

[14] Pierre-Jean Bénard, Yann Traonmilin, Jean-François Aujol, and Emmanuel Soubies. Estimation of off-the-grid sparse spikes with over-parametrized projected gradient descent: Theory and application, September 2023.

[15] Abderrahmane Bendali, Pierre-Henri Cocquet, and Sébastien Tordeux. Approximation by Multipoles of the Multiple Acoustic Scattering by Small Obstacles in Three Dimensions and Application to the Foldy Theory of Isotropic Scattering. *Archive for Rational Mechanics and Analysis*, 219(3):1017–1059, March 2016.

[16] Mathieu Blondel, Quentin Berthet, Marco Cuturi, Roy Frostig, Stephan Hoyer, Felipe Llinares-Lopez, Fabian Pedregosa, and Jean-Philippe Vert. Efficient and Modular Implicit Differentiation. *Advances in Neural Information Processing Systems*, 35:5230–5242, December 2022.

[17] Nicholas Boyd, Geoffrey Schiebinger, and Benjamin Recht. The Alternating Descent Conditional Gradient Method for Sparse Inverse Problems. *SIAM Journal on Optimization*, 27(2):616–639, January 2017.

[18] Claire Boyer, Yohann De Castro, and Joseph Salmon. Adapting to unknown noise level in sparse deconvolution. *Information and Inference: A Journal of the IMA*, 6(3):310–348, September 2017.

[19] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, and Skye Wanderman-Milne. JAX: Composable transformations of Python+ NumPy programs. 2018.

[20] Kristian Bredies and Hanna Katriina Pikkarainen. Inverse problems in spaces of measures. *ESAIM: Control, Optimisation and Calculus of Variations*, 19(1):190–218, 2013.

[21] Cristina Butucea, Jean-François Delmas, Anne Dutfoy, and Clément Hardy. Off-the-grid learning of sparse mixtures from a continuous dictionary, June 2022.
[22] Fioralba Cakoni and David Colton. *Qualitative Methods in Inverse Scattering Theory: An Introduction*. Springer Science & Business Media, December 2005.

[23] Emmanuel J. Candès and Carlos Fernandez-Granda. Super-Resolution from Noisy Data. *Journal of Fourier Analysis and Applications*, 19(6):1229–1254, December 2013.

[24] Emmanuel J. Candès and Carlos Fernandez-Granda. Towards a Mathematical Theory of Super-resolution. *Communications on Pure and Applied Mathematics*, 67(6):906–956, 2014.

[25] Emmanuel J. Candès, Thomas Strohmer, and Vladislav Voroninski. PhaseLift: Exact and Stable Signal Recovery from Magnitude Measurements via Convex Programming. *Communications on Pure and Applied Mathematics*, 66(8):1241–1274, 2013.

[26] Yves Capdeboscq and Shaun Chen Yang Ong. Extending representation formulas for boundary voltage perturbations of low volume fraction to very contrasted conductivity inhomogeneities. *Comptes Rendus. Mathématique*, 360(G2):127–150, 2022.

[27] Rémi Carminati and John C. Schotland. *Principles of Scattering and Transport of Light*. Cambridge University Press, Cambridge, 2021.

[28] Maxence Cassier and Christophe Hazard. Multiple scattering of acoustic waves by small sound-soft obstacles in two dimensions: Mathematical justification of the Foldy–Lax model. *Wave Motion*, 50(1):18–28, January 2013.

[29] Y. De Castro, S. Gadat, C. Marteau, and C. Maugis-Rabusseau. SuperMix: Sparse regularization for mixtures. *The Annals of Statistics*, 49(3):1779–1809, June 2021.

[30] Paul Catala. *Relaxations semi-definies positives pour l'imagerie*. PhD thesis, Université Paris sciences et lettres, October 2020.

[31] Paul Catala, Vincent Duval, and Gabriel Peyré. A Low-Rank Approach to Off-the-Grid Sparse Super-resolution. *SIAM Journal on Imaging Sciences*, 12(3):1464–1500, January 2019.

[32] Durga Prasad Challa and Mourad Sini. Inverse scattering by point-like scatterers in the Foldy regime. *Inverse Problems*, 28(12):125006, November 2012.

[33] Durga Prasad Challa and Mourad Sini. On the Justification of the Foldy–Lax Approximation for the Acoustic Scattering by Small Rigid Bodies of Arbitrary Shapes. *Multiscale Modeling & Simulation*, 12(1):55–108, January 2014.

[34] Margaret Cheney. The linear sampling method and the MUSIC algorithm. *Inverse Problems*, 17(4):591, August 2001.

[35] Lénaïc Chizat. Sparse optimization on measures with over-parameterized gradient descent. *Mathematical Programming*, 194(1):487–532, July 2022.

[36] David Colton and Andreas Kirsch. A simple method for solving inverse scattering problems in the resonance region. *Inverse Problems*, 12(4):383, August 1996.

[37] David Colton and Rainer Kress. *Inverse Acoustic and Electromagnetic Scattering Theory*. Springer Science & Business Media, October 2012.

[38] Yohann De Castro, F. Gamboa, Didier Henrion, and J.-B. Lasserre. Exact Solutions to Super Resolution on Semi-Algebraic Domains in Higher Dimensions. *IEEE Transactions on Information Theory*, 63(1):621–630, January 2017.

[39] Yohann De Castro and Fabrice Gamboa. Exact reconstruction using Beurling minimal extrapolation. *Journal of Mathematical Analysis and Applications*, 395(1):336–354, November 2012.
[40] Pedro de Vries, David V. van Coevorden, and Ad Lagendijk. Point scatterers for classical waves. *Reviews of Modern Physics*, 70(2):447-466, April 1998.

[41] Quentin Denoyelle. *Theoretical and Numerical Analysis of Super-Resolution Without Grid*. PhD thesis, Université Paris sciences et lettres, July 2018.

[42] Quentin Denoyelle, Vincent Duval, and Gabriel Peyré. Support Recovery for Sparse Super-Resolution of Positive Measures. *Journal of Fourier Analysis and Applications*, 23(5):1153-1194, October 2017.

[43] Quentin Denoyelle, Vincent Duval, Gabriel Peyré, and Emmanuel Soubies. The sliding Frank-Wolfe algorithm and its application to super-resolution microscopy. *Inverse Problems*, 36(1):014001, December 2019.

[44] Anthony J. Devaney. *Mathematical Foundations of Imaging, Tomography and Wavefield Inversion*. Cambridge University Press, Cambridge, 2012.

[45] Anthony J. Devaney, Edwin A. Marengo, and Fred K. Gruber. Time-reversal-based imaging and inverse scattering of multiply scattering point targets. *The Journal of the Acoustical Society of America*, 118(5):3129-3138, November 2005.

[46] Vincent Duval and Gabriel Peyré. Exact Support Recovery for Sparse Spikes Deconvolution. *Foundations of Computational Mathematics*, 15(5):1315-1355, October 2015.

[47] Claudia Errico, Juliette Pierre, Sophie Pezet, Yann Desailly, Zsolt Lenkei, Olivier Couture, and Mickael Tanter. Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging. *Nature*, 527(7579):499-502, November 2015.

[48] Carlos Fernandez-Granda. Support detection in super-resolution. In *Proceedings of the 10th International Conference on Sampling Theory and Applications*, pages 145-148, 2013.

[49] Axel Flinth, Frédéric de Gournay, and Pierre Weiss. On the linear convergence rates of exchange and continuous methods for total variation minimization. *Mathematical Programming*, 190(1):221-257, November 2021.

[50] Leslie L. Foldy. The Multiple Scattering of Waves. I. General Theory of Isotropic Scattering by Randomly Distributed Scatterers. *Physical Review*, 67(3-4):107-119, February 1945.

[51] Pedro Freitas. Sharp bounds for the modulus and phase of Hankel functions with applications to Jaeger integrals. *Mathematics of Computation*, 87(309):289-308, January 2018.

[52] Anna C. Gilbert, Howard W. Levinson, and John C. Schotland. Nonlinear Iterative Hard Thresholding for Inverse Scattering. *SIAM Journal on Imaging Sciences*, 13(1):108-140, January 2020.

[53] Fred K. Gruber, Edwin A. Marengo, and Anthony J. Devaney. Time-reversal imaging with multiple signal classification considering multiple scattering between the targets. *The Journal of the Acoustical Society of America*, 115(6):3042-3047, June 2004.

[54] Eldad Haber, Uri M. Ascher, and Doug Oldenburg. On optimization techniques for solving nonlinear inverse problems. *Inverse Problems*, 16(5):1263, October 2000.

[55] Barbara Kaltenbacher, Andreas Neubauer, and Otmar Scherzer. Iterative Regularization Methods for Nonlinear Ill-Posed Problems. In *Iterative Regularization Methods for Nonlinear Ill-Posed Problems*. De Gruyter, September 2008.

[56] Barbara Kaltenbacher, Frank Schöpfer, and Thomas Schuster. Iterative methods for nonlinear ill-posed problems in Banach spaces: Convergence and applications to parameter identification problems. *Inverse Problems*, 25(6):065003, April 2009.
[57] Nicolas Keriven, Anthony Bourrier, Rémi Gribonval, and Patrick Pérez. Sketching for large-scale learning of mixture models. *Information and Inference: A Journal of the IMA*, 7(3):447–508, September 2018.

[58] Andreas Kirsch and Natalia Grinberg. *The Factorization Method for Inverse Problems*. OUP Oxford, December 2007.

[59] Andreas Kirsch and Rainer Kress. On an Integral Equation of the first Kind in Inverse Acoustic Scattering. In John Rozier Cannon and Ulrich Hornung, editors, *Inverse Problems: Proceedings of the Conference Held at the Mathematical Research Institute at Oberwolfach, Black Forest, May 18–24, International Series of Numerical Mathematics*, pages 93–102. Birkhäuser, Basel, 1986.

[60] Andreas Kirsch and Rainer Kress. A numerical method for an inverse scattering problem. In Heinz W. Engl and C. W. Groetsch, editors, *Inverse and Ill-Posed Problems*, pages 279–289. Academic Press, January 1987.

[61] Andreas Kirsch and Rainer Kress. An optimization method in inverse acoustic scattering. *Boundary elements IX*, 3:3–18, 1987.

[62] L. J. Landau. Bessel functions: monotonicity and bounds. *Journal of the London Mathematical Society*, 61(1):197–215, February 2000.

[63] Bastien Laville, Laure Blanc-Féraud, and Gilles Aubert. Off-The-Grid Variational Sparse Spike Recovery: Methods and Algorithms. *Journal of Imaging*, 7(12):266, December 2021.

[64] Melvin Lax. Multiple Scattering of Waves. *Reviews of Modern Physics*, 23(4):287–310, October 1951.

[65] H. Lev-Ari and A.J. Devancy. The time-reversal technique re-interpreted: Subspace-based signal processing for multi-static target location. In *Proceedings of the 2000 IEEE Sensor Array and Multichannel Signal Processing Workshop. SAM 2000 (Cat. No.00EX410)*, pages 509–513, March 2000.

[66] Andrea Mantile, Andrea Posilicano, and Mourad Sini. On the origin of Minnaert resonances. *Journal de Mathématiques Pures et Appliquées*, 165:106–147, September 2022.

[67] P. A. Martin. *Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles*. Encyclopedia of Mathematics and Its Applications. Cambridge University Press, Cambridge, 2006.

[68] F. Natterer. An error bound for the Born approximation. *Inverse Problems*, 20(2):447, January 2004.

[69] Gabriel Peyré and Marco Cuturi. Computational Optimal Transport: With Applications to Data Science. *Foundations and Trends® in Machine Learning*, 11(5-6):355–607, February 2019.

[70] Clarice Poon. An introduction to sparse spikes recovery via the BLASSO, 2019.

[71] Clarice Poon, Nicolas Keriven, and Gabriel Peyré. Support Localization and the Fisher Metric for off-the-grid Sparse Regularization. In *Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics*, pages 1341–1350. PMLR, April 2019.

[72] Clarice Poon, Nicolas Keriven, and Gabriel Peyré. The Geometry of Off-the-Grid Compressed Sensing. *Foundations of Computational Mathematics*, 23(1):241–327, February 2023.

[73] Filippo Santambrogio. *Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling*. Birkhäuser, October 2015.

[74] Gongguo Tang, Badri Narayan Bhaskar, and Benjamin Recht. Near Minimax Line Spectral Estimation. *IEEE Transactions on Information Theory*, 61(1):499–512, January 2015.
[75] Gongguo Tang, Badri Narayan Bhaskar, Parikshit Shah, and Benjamin Recht. Compressed Sensing Off the Grid. *IEEE Transactions on Information Theory*, 59(11):7465-7490, November 2013.

[76] Tuomo Valkonen. Proximal methods for point source localisation. *Journal of Nonsmooth Analysis and Optimization*, Volume 4(Original research articles), September 2023.

[77] G. N. Watson. *A Treatise on the Theory of Bessel Functions*. Cambridge University Press, August 1995.

## A Proof of Theorems 4.3 and 4.4

In this section, we prove Theorems 4.3 and 4.4 by showing that the assumptions of [72, Theorem 3] and [71, Theorem 1] are satisfied in our setting. This amounts to studying the properties of a kernel $K: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ associated to the measurement operator. For a forward operator of the form (18), $K$ is defined by

$$ K(x, x') \stackrel{\text{def.}}{=} \int_{\mathcal{X}} \overline{\varphi_{\omega}(x)\varphi_{\omega}(x')} d\Lambda(\omega), \quad (x, x') \in \mathcal{X} \times \mathcal{X}. \qquad (20) $$

We first begin by deriving the expression of $K$ in our setting, and then turn to the study of its properties.

**Proposition A.1** If $\varphi_{\omega} = e^{-i\langle\omega, \cdot\rangle}$ and $\Lambda$ is the uniform law on $B(0, 2\kappa)$, then for every $x, x' \in \mathcal{X}$ we have $K(x, x') = \rho(\sigma^{-1}|x - x'|)$ with $\sigma^{-1} \stackrel{\text{def.}}{=} 2\kappa$ and $\rho(s) \stackrel{\text{def.}}{=} (2/s)^{d/2}\Gamma(d/2 + 1)J_{d/2}(s)$, where $J_{\alpha}$ denotes the Bessel function of the first kind and order $\alpha$.

*Proof.* A simple computation shows that $K(x, x')$ is proportional to the Fourier transform of $1_{B(0,\sigma^{-1})}$ at $x - x'$, which allows us to conclude. $\square$

We notice that, since $J_{\alpha}(s) \sim (s/2)^{\alpha}/\Gamma(\alpha + 1)$ at $0$, we have $K(x, x) = \mathbb{E}_{\omega \sim \Lambda}[|\varphi_{\omega}(x)|^2] = \rho(0) = 1$ for every $x \in \mathcal{X}$. For convenience, we define $\rho_{\sigma} \stackrel{\text{def.}}{=} \rho(\sigma^{-1}.)$. Let us now list a few notations and definitions introduced in [71, 72], before stating the various results we need to prove.

**The Fisher metric.** In [71, 72], a special metric on $\mathcal{X}$ related to $K$ is introduced. In our case, the metric tensor $g_x \stackrel{\text{def.}}{=} \nabla_1\nabla_2 K(x, x)$ is constant and $g_x = -\rho_{\sigma}''(0) \text{Id} = -\sigma^{-2}\rho_{\sigma}''(0) \text{Id}$ with $\rho_{\sigma}''(0) = -1/(d+2)$. Hence, the associated distance is given by

$$ \mathfrak{D}_{g}(x, x') = \sqrt{|\rho_{\sigma}''(0)|} |x - x'| = \sqrt{-\rho_{\sigma}''(0)} \sigma^{-1} |x - x'|, \quad x, x' \in \mathcal{X}. $$

For every $u \in \mathbb{C}^d$, we define $|u|_{\sigma} \stackrel{\text{def.}}{=} \sqrt{u^* g_x u} = \sigma^{-1} \sqrt{-\rho_{\sigma}''(0)} |u|$.

**Covariant derivatives.** For $j \in \{0, 1, 2\}$, the “covariant derivative” $D_j[f](x): (\mathbb{C}^d)^j \to \mathbb{C}$ are defined as

$$ D_0[f](x) \stackrel{\text{def.}}{=} f(x), \quad D_1[f](x)[v] \stackrel{\text{def.}}{=} v^T \nabla f(x), \quad D_2[f](x)[v_1, v_2] \stackrel{\text{def.}}{=} v_1^T \nabla^2 f(x)v_2. $$

The associated operator norms are defined by

$$ \|D_1[f]\|_{x} \stackrel{\text{def.}}{=} \|g_x^{-1/2} \nabla f(x)\|_{2}, \quad \|D_2[f]\|_{x} \stackrel{\text{def.}}{=} \|g_x^{-1/2} \nabla^2 f(x) g_x^{-1/2}\|_{2}, $$

and we also define $L_0(\omega) \stackrel{\text{def.}}{=} \|\varphi_{\omega}\|_{\infty}$, $L_j(\omega) \stackrel{\text{def.}}{=} \sup_{x \in \mathcal{X}} \|D_j[\varphi_{\omega}]\|_{x}$, ($j = 1, 2$). Given $0 \le i, j \le 2$, let $K^{(ij)}(x, x')$ be defined, for $Q \in (\mathbb{C}^d)^i$ and $V \in (\mathbb{C}^d)^j$, by

$$ [Q] K^{(ij)}(x, x') [V] \stackrel{\text{def.}}{=} \mathbb{E} \left[ \overline{D_i[\varphi_{\omega}](x)[Q]} D_j[\varphi_{\omega}](x') [V] \right]. $$
The associated operator norms are defined by

$$ \left\| K^{(ij)}(x, x') \right\|_{x,x'} \stackrel{\text{def}}{=} \sup_{Q,V} [Q] K^{(ij)}(x, x') [V], \qquad (21) $$

where the supremum is taken over all $Q = [q_1, ..., q_i]$ and $V = [v_1, ..., v_j]$ such that $|q_k|_x \le 1$ for every $k \in \{1, ..., i\}$ and $|v_l|_{x'} \le 1$ for every $l \in \{1, ..., j\}$. In (21), we rather write $\| \cdot \|_x$ instead of $\| \cdot \|_{x,x'}$ when the dependence is only on $x$ (that is to say when $j = 0$). We also define $B_{ij} \stackrel{\text{def.}}{=} \sup_{x,x' \in \mathcal{X}} \|K^{(ij)}(x, x')\|_{x,x'}$.

**Kernel width.** Finally, we define the kernel width of $K$, defined in [72, Definition 3] as:

$$ W(h, s) \stackrel{\text{def.}}{=} \inf \left\{ \Delta \left| \sum_{k=2}^{s} \|K^{(ij)}(x_1, x_k)\|_{x_1, x_k} \le h, (i, j) \in \{0, 1\} \times \{0, 2\}, (x_k)_{k=1}^{s} \in \mathcal{S}_{\Delta} \right. \right\} $$

where $\mathcal{S}_{\Delta} \stackrel{\text{def.}}{=} \{(x_k)_{k=1}^s \in \mathcal{X}^s | \forall k \ne l, \mathfrak{D}_g(x_k, x_l) \ge \Delta\}$.

**Outline of the appendix.** First, we state useful properties of Bessel functions in Appendix A.1. In Appendix A.2, we give the expressions of the derivatives of the kernel and the operator norm of its covariant derivatives. The remaining of the appendix is dedicated to the proof that assumptions (i) - (v) below hold.

(i) $B_{ij} = \mathcal{O}(1)$ for every $(i, j) \in \{0, 1, 2\}$.

(ii) There exist $r_{\text{near}} > 0$ and $0 < \bar{\epsilon}_2 < r_{\text{near}}^{-2}$ such that $-K^{(02)}(x, x') [v, v] \ge \bar{\epsilon}_2 |v|_x^2$ for every $v \in \mathbb{C}^d$ as soon as $\mathfrak{D}_g(x, x') \le r_{\text{near}}$.

(iii) There exist $\bar{\epsilon}_0 < 1$ such that $|K(x, x')| \le 1 - \bar{\epsilon}_0$ as soon as $\mathfrak{D}_g(x, x') \ge r_{\text{near}}$.

(iv) $W(h, s) = \mathcal{O}(s^{2/(d+1)})$ as soon as $h = \mathcal{O}(1)$.

(v) $\|L_j\|_{\infty} = \mathcal{O}(1)$ for every $j \in \{0, ..., 3\}$.<fnref n="4" />

We stress that the constants appearing in the various bounds do not depend on $\sigma$. Provided (i) - (v) hold, [72, Theorem 3] and [71, Theorem 1] can be applied to our setting, which yields precisely Theorems 4.3 and 4.4.

## A.1 Bessel functions

In this subsection, we give two properties of Bessel functions that we use in the following. First, we have

$$ \forall \alpha > 0, \forall x > 0, J'_{\alpha}(x) - (\alpha/x)J_{\alpha}(x) = -J_{\alpha+1}(x). \qquad (22) $$

Then, we have the following bound, which implies the boundedness of $x \mapsto |J_{\alpha}(x)|/x^{\beta}$ for any $\beta \le \alpha$.

**Lemma A.2** For every $\alpha > 0$ and $x > 0$, we have:

$$ |J_{\alpha}(x)| \le \min \left( \frac{1}{\Gamma(\alpha + 1)} \left(\frac{x}{2}\right)^{\alpha} \left( \left|1 - \frac{1}{\alpha + 1} \frac{x^2}{4}\right| + \frac{1}{(\alpha + 1)(\alpha + 2)} \left(e^{x^2/4} - 1 - \frac{x^2}{4}\right) \right), \frac{0.8}{x^{1/3}} \right). $$

<footnote n="4">⁴See [72, Section 5.1] for the definition of $L_3$.</footnote>
Proof. The inequality $|J_\alpha(x)| \le 0.8x^{-1/3}$ is proved in [62]. For the other inequality, we use the series expansion of $J_\alpha$:

$$J_{\alpha}(x) = \frac{1}{\Gamma(\alpha + 1)} \left(\frac{x}{2}\right)^{\alpha} \left(1 - \frac{1}{\alpha + 1} \frac{x^2}{4} + \Gamma(\alpha + 1) \sum_{p=2}^{+\infty} \frac{(-1)^p}{p! \Gamma(p + \alpha + 1)} \left(\frac{x}{2}\right)^{2p}\right).$$

Then, we notice that

$$\begin{aligned} \left| \Gamma(\alpha + 1) \sum_{p=2}^{+\infty} \frac{(-1)^p}{p! \Gamma(p + \alpha + 1)} \left(\frac{x}{2}\right)^{2p} \right| & \le \Gamma(\alpha + 1) \sum_{p=2}^{+\infty} \frac{1}{p! \Gamma(p + \alpha + 1)} \left(\frac{x}{2}\right)^{2p} \\ & \le \frac{\Gamma(\alpha + 1)}{\Gamma(\alpha + 3)} \sum_{p=2}^{+\infty} \frac{1}{p!} \left(\frac{x^2}{4}\right)^p \\ & = \frac{1}{(\alpha + 1)(\alpha + 2)} \left(e^{x^2/4} - 1 - \frac{x^2}{4}\right), \end{aligned}$$

which yields the result.

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

## A.2 Derivatives of the kernel

Let $x, x' \in \mathcal{X}$. We define $t = x - x'$ and, abusing notation, denote $K(t) = \rho_\sigma(|t|)$. We have the following:

$$\begin{aligned} \nabla K(t) &= \frac{\rho'_{\sigma}(|t|)}{|t|} t, \\ \nabla^2 K(t) &= \frac{\rho'_{\sigma}(|t|)}{|t|} Id + \left[ \rho''_{\sigma}(|t|) - \frac{\rho'_{\sigma}(|t|)}{|t|} \right] \frac{tt^T}{|t|^2} \\ &= \frac{\rho'_{\sigma}(|t|)}{|t|} \left( Id - \frac{tt^T}{|t|^2} \right) + \rho''_{\sigma}(|t|) \frac{tt^T}{|t|^2}, \\ \partial_i \nabla^2 K(t) &= \left[ \rho''_{\sigma}(|t|) - \frac{\rho'_{\sigma}(|t|)}{|t|} \right] \frac{t_i}{|t|^2} Id + \left[ \rho'''_{\sigma}(|t|) \frac{t_i}{|t|} - \frac{t_i}{|t|^2} \left( \rho''_{\sigma}(|t|) - \frac{\rho'_{\sigma}(|t|)}{|t|} \right) \right] \frac{tt^T}{|t|^2} \\ & \quad + \left[ \rho''_{\sigma}(|t|) - \frac{\rho'_{\sigma}(|t|)}{|t|} \right] \frac{1}{|t|^2} \left( t e_i^T + e_i t^T - 2t_i \frac{tt^T}{|t|^2} \right) \\ &= \left[ \rho''_{\sigma}(|t|) - \frac{\rho'_{\sigma}(|t|)}{|t|} \right] \frac{1}{|t|^2} \left( t e_i^T + e_i t^T + t_i \left( Id - 3 \frac{tt^T}{|t|^2} \right) \right) + \rho'''_{\sigma}(|t|) \frac{t_i}{|t|} \frac{tt^T}{|t|^2}. \end{aligned}$$

Considering the definition of (21), we also obtain:

$$\|K^{(00)}(x, x')\| = |K(t)| = |\rho(\sigma^{-1}|t|)|, \qquad (23)$$

$$\|K^{(10)}(x, x')\|_x = \frac{\sigma}{\sqrt{-\rho''(0)}} \|\nabla K(t)\|_2 = \frac{1}{\sqrt{-\rho''(0)}} |\rho'(\sigma^{-1}|t|)|, \qquad (24)$$

$$\|K^{(11)}(x, x')\|_{x, x'} = \frac{\sigma^2}{(-\rho''(0))} \|\nabla^2 K(t)\|_2 \le \frac{1}{(-\rho''(0))} \left( \left| \frac{\rho'(\sigma^{-1}|t|)}{\sigma^{-1}|t|} \right| + \left| \rho''(\sigma^{-1}|t|) - \frac{\rho'(\sigma^{-1}|t|)}{\sigma^{-1}|t|} \right| \right), \qquad (25)$$

$$\|K^{(20)}(x, x')\|_x = \|K^{(11)}(x, x')\|_{x, x'}. \qquad (26)$$

## A.3 Uniform bounds

**Lemma A.3** For $(i, j) \in \{0, 1, 2\}$, we have $B_{ij} = \mathcal{O}(1)$.

*Proof.*
• $B_{00}$: considering (23), we have that $B_{00} \le \|\rho\|_{\infty}$, which is finite by Lemma A.2.

• $B_{10}$: considering (24), we have $B_{10} \le \sqrt{d+2} \|\rho'\|_{\infty} \le \sqrt{5} \|\rho'\|_{\infty}$. Using (22), we find

$$ 
\begin{aligned}
\rho'(s) &= 2^{d/2}\Gamma(d/2+1)(-(d/2)J_{d/2}(s)/s^{d/2+1} + J'_{d/2}(s)/s^{d/2}) \\
&= (2/s)^{d/2}\Gamma(d/2+1)(-(d/2)J_{d/2}(s)/s + (d/2)J_{d/2}(s)/s - J_{d/2+1}(s)) \\
&= (2/s)^{d/2}\Gamma(d/2+1)(-J_{d/2+1}(s)),
\end{aligned} 
\quad (27) $$

and we can use Lemma A.2 again to get $B_{10} = \mathcal{O}(1)$.

• $B_{11} = B_{20}$: considering (25) and (26), we have

$$ B_{11} = B_{20} \le (d+2) \left( \|\rho''\|_{\infty} + 2 \sup_{s>0} \left| \frac{\rho'(s)}{s} \right| \right). $$

Using (27) and (22) we get

$$ 
\begin{aligned}
\rho''(s) &= 2^{d/2}\Gamma(d/2+1)((d/2)J_{d/2+1}(s)/s^{d/2+1} - J'_{d/2+1}(s)/s^{d/2}) \\
&= (2/s)^{d/2}\Gamma(d/2+1)((d/2)J_{d/2+1}(s)/s - (d/2+1)J_{d/2+1}(s)/s + J_{d/2+2}(s)) \\
&= (2/s)^{d/2}\Gamma(d/2+1)(J_{d/2+2}(s) - J_{d/2+1}(s)/s)
\end{aligned} 
\quad (28) $$

and another application of Lemma A.2 yields that $B_{11} = B_{20} = \mathcal{O}(1)$.

• $B_{21}$: for any $x, x' \in \mathcal{X}$, defining $t = x - x'$, for every $q \in \mathbb{C}^d$, we have

$$ 
\begin{aligned}
[q]K^{(12)}(x, x') &= \sum_{i=1}^{d} q_i (\partial_i \nabla^2 K(t)) \\
&= \left[ \rho''_{\sigma}(|t|) - \frac{\rho'_{\sigma}(|t|)}{|t|} \right] \frac{1}{|t|^2} \left( tq^T + qt^T + \langle q, t \rangle \left( \mathrm{Id} - 3 \frac{tt^T}{|t|^2} \right) \right) + \rho'''_{\sigma}(|t|) \frac{\langle q, t \rangle}{|t|} \frac{tt^T}{|t|^2}.
\end{aligned} 
$$

Hence

$$ 
\begin{aligned}
\|K^{(12)}(x, x')\|_{x, x'} &= \sup_{|q|_x \le 1, \ |v_i|_x \le 1} \left| v_1^T \left( \sum_{i=1}^d q_i (\partial_i \nabla^2 K(t)) \right) v_2 \right| \\
&= \left( \frac{\sigma}{\sqrt{-\rho''(0)}} \right)^3 \sup_{|q| \le 1, \ |v_i| \le 1} \left| v_1^T \left( \sum_{i=1}^d q_i (\partial_i \nabla^2 K(t)) \right) v_2 \right| \\
&\le \left( \frac{\sigma}{\sqrt{-\rho''(0)}} \right)^3 \left( 6 \left\| \frac{1}{|\cdot|} \left( \rho''_{\sigma}(|\cdot|) - \frac{\rho'_{\sigma}(|\cdot|)}{|\cdot|} \right) \right\|_{\infty} + \|\rho'''_{\sigma}\|_{\infty} \right).
\end{aligned} 
$$

Now, we notice that

$$ 
\begin{aligned}
\left| \frac{1}{|t|} \left( \rho''_{\sigma}(|t|) - \frac{\rho'_{\sigma}(|t|)}{|t|} \right) \right| &= \sigma^{-3} \left| \frac{1}{\sigma^{-1}|t|} \left( \rho''(\sigma^{-1}|t|) - \frac{\rho'(\sigma^{-1}|t|)}{\sigma^{-1}|t|} \right) \right| \\
&\le \sigma^{-3} \sup_{s>0} \left| \frac{1}{s} \left( \rho''(s) - \frac{\rho'(s)}{s} \right) \right|
\end{aligned} 
$$

and $\|\rho'''_{\sigma}\|_{\infty} = \sigma^{-3}\|\rho'''_{\sigma}\|_{\infty}$, which yields

$$ \|K^{(12)}(x, x')\|_{x, x'} \le (-\rho''(0))^{-3/2} \left( 6 \sup_{s>0} \left| \frac{1}{s} \left( \rho''(s) - \frac{\rho'(s)}{s} \right) \right| + \|\rho'''_{\sigma}\|_{\infty} \right). $$
By (27) and (28) we get

$$ \rho''(s) - \rho'(s)/s = (2/s)^{d/2} \Gamma(d/2+1) J_{d/2+2}(s). \quad (29) $$

Moreover, using (22) and (28), we obtain

$$ \begin{aligned} \rho'''(s) &= (2/s)^{d/2} \Gamma(d/2+1) \left( -\frac{d}{2} \frac{J_{d/2+2}(s)}{s} + \frac{d}{2} \frac{J_{d/2+1}(s)}{s^2} + J'_{d/2+2}(s) - \frac{J'_{d/2+1}(s)}{s} + \frac{J_{d/2+1}(s)}{s^2} \right) \\ &= (2/s)^{d/2} \Gamma(d/2+1) \left( -\frac{d}{2} \frac{J_{d/2+2}(s)}{s} + \left(\frac{d}{2}+1\right) \frac{J_{d/2+1}(s)}{s^2} + \frac{d/2+2}{s} J_{d/2+2}(s) - J_{d/2+3}(s) \right) \\ &\quad - \left(\frac{d}{2}+1\right) \frac{J_{d/2+1}(s)}{s^2} + \frac{J_{d/2+2}(s)}{s} \\ &= (2/s)^{d/2} \Gamma(d/2+1) \left[ -J_{d/2+3}(s) + (3/s) J_{d/2+2}(s) \right]. \end{aligned} \quad (30) $$

Applying Lemma A.2 yields $B_{21} = \mathcal{O}(1)$.

• $B_{22}$: first, we have:

$$ \begin{aligned} \|K^{(22)}(x, x')\|_{x, x'} &= \sup_{|q_k|_x \le 1, \ |v_l|_x \le 1} \left| v_1^T \left( \sum_{i,j=1}^d q_{1,i} q_{2,j} \partial_{ij} \nabla^2 K(t) \right) v_2 \right| \\ &= \left( \frac{\sigma}{\sqrt{-\rho''(0)}} \right)^4 \sup_{|q_k|_2 \le 1} \left\| \sum_{i,j=1}^d q_{1,i} q_{2,j} \partial_{ij} \nabla^2 K(t) \right\|_2. \end{aligned} $$

Then, we recall that $\partial_i \nabla^2 K(t) = a_i(t)b_i(t) + c_i(t)$ with

$$ \begin{aligned} a_i(t) &\stackrel{\text{def.}}{=} \frac{1}{|t|^2} \left[ \rho''(|t|) - \frac{\rho'_{\sigma}(|t|)}{|t|} \right], \\ b_i(t) &\stackrel{\text{def.}}{=} \left( t e_i^T + e_i t^T + t_i \left( \mathrm{Id} - 3 \frac{tt^T}{|t|^2} \right) \right), \\ c_i(t) &\stackrel{\text{def.}}{=} \rho'''(|t|) \frac{t_i}{|t|} \frac{tt^T}{|t|^2}. \end{aligned} $$
Denoting by $E_{ij} \in \mathbb{R}^{d \times d}$ the matrix with entry 1 at row $i$ and column $j$ and zero otherwise, we have

$$
\begin{align*}
\partial_j a_i(t) &= \frac{1}{|t|^2} \left[ \rho_{\sigma}'''(|t|) \frac{t_j}{|t|} - \frac{t_j}{|t|^2} \left( \rho_{\sigma}''(|t|) - \frac{\rho_{\sigma}'(|t|)}{|t|} \right) \right] - 2 \frac{t_j}{|t|^4} \left[ \rho_{\sigma}''(|t|) - \frac{\rho_{\sigma}'(|t|)}{|t|} \right] \\
&= \frac{t_j}{|t|^3} \left[ \rho_{\sigma}'''(|t|) - \frac{1}{|t|} \left( \rho_{\sigma}''(|t|) - \frac{\rho_{\sigma}'(|t|)}{|t|} \right) - \frac{2}{|t|} \left( \rho_{\sigma}''(|t|) - \frac{\rho_{\sigma}'(|t|)}{|t|} \right) \right] \\
&= \frac{t_j}{|t|^3} \left[ \rho_{\sigma}'''(|t|) - \frac{3}{|t|} \left( \rho_{\sigma}''(|t|) - \frac{\rho_{\sigma}'(|t|)}{|t|} \right) \right],
\end{align*}
$$

$$
\partial_j b_i(t) = E_{ji} + E_{ij} + \delta_{ij} \left( \mathrm{Id} - 3 \frac{tt^T}{|t|^2} \right) - 3 \frac{t_i}{|t|^2} \left( t e_j^T + e_j t^T - 2 t_j \frac{tt^T}{|t|^2} \right),
$$

$$
\partial_j c_i(t) = \rho_{\sigma}^{(4)}(|t|) \frac{t_i t_j}{|t|^2} \frac{tt^T}{|t|^2} \\
+ \rho_{\sigma}^{\prime\prime\prime}(|t|) \left[ \left( \delta_{ij} - \frac{t_i t_j}{|t|^2} \right) \frac{1}{|t|} \frac{tt^T}{|t|^2} + \frac{t_i}{|t|} \frac{1}{|t|^2} \left( t e_j^T + e_j t^T + t_j \left( \mathrm{Id} - 3 \frac{tt^T}{|t|^2} \right) \right) \right],
$$

and hence

$$
\begin{align*}
\sum_{i,j=1}^{d} q_{1,i}q_{2,j}(\partial_j a_i(t))b_i(t) &= \frac{\langle q_2, t \rangle}{|t|^3} \left[ \rho_{\sigma}'''(|t|) - \frac{3}{|t|} \left( \rho_{\sigma}''(|t|) - \frac{\rho_{\sigma}'(|t|)}{|t|} \right) \right] \left( t q_1^T + q_1 t^T + \langle q_1, t \rangle \left( \mathrm{Id} - 3 \frac{tt^T}{|t|^2} \right) \right), \\
\sum_{i,j=1}^{d} q_{1,i}q_{2,j}a_i(t)(\partial_j b_i(t)) &= a_i(t) \left( q_1 q_2^T + q_2 q_1^T + \langle q_1, q_2 \rangle \left( \mathrm{Id} - 3 \frac{tt^T}{|t|^2} \right) - 3 \frac{\langle q_1, t \rangle}{|t|^2} \left( t q_2^T + q_2 t^T - 2 \langle q_2, t \rangle \frac{tt^T}{|t|^2} \right) \right), \\
\sum_{i,j=1}^{d} q_{1,i}q_{2,j}\partial_j c_i(t) &= \rho_{\sigma}^{(4)}(|t|) \frac{\langle q_1, t \rangle \langle q_2, t \rangle}{|t|^2} \frac{tt^T}{|t|^2} \\
&\quad + \rho_{\sigma}^{\prime\prime\prime}(|t|) \left[ \left( \langle q_1, q_2 \rangle - \frac{\langle q_1, t \rangle \langle q_2, t \rangle}{|t|^2} \right) \frac{1}{|t|} \frac{tt^T}{|t|^2} + \frac{\langle q_1, t \rangle}{|t|} \frac{1}{|t|^2} \left( t q_2^T + q_2 t^T + \langle q_2, t \rangle \left( \mathrm{Id} - 3 \frac{tt^T}{|t|^2} \right) \right) \right].
\end{align*}
$$

This ultimately yields

$$
\begin{align*}
\sup_{|q_k|_2 \le 1} \left\| \sum_{i,j=1}^d q_{1,i} q_{2,j} (\partial_j a_i(t)) b_i(t) \right\|_2 &\le 6 \left\| \frac{1}{|\cdot|} \left( \rho_{\sigma}'''(|\cdot|) - \frac{3}{|\cdot|} \left( \rho_{\sigma}''(|\cdot|) - \frac{\rho_{\sigma}'(|\cdot|)}{|\cdot|} \right) \right) \right\|_{\infty} \\
&\le 6 \sigma^{-4} \left( \sup_{s>0} \left| \frac{\rho'''(s)}{s} \right| + 3 \sup_{s>0} \left| \frac{1}{s^2} \left( \rho''(s) - \frac{\rho'(s)}{s} \right) \right| \right), \\
\sup_{|q_k|_2 \le 1} \left\| \sum_{i,j=1}^d q_{1,i} q_{2,j} a_i(t) \partial_j b_i(t) \right\|_2 &\le 18 \left\| \frac{1}{|\cdot|} \left( \rho_{\sigma}''(|\cdot|) - \frac{\rho_{\sigma}'(|\cdot|)}{|\cdot|} \right) \right\|_{\infty} \\
&= 18 \sigma^{-4} \sup_{s>0} \left| \frac{1}{s^2} \left( \rho''(s) - \frac{\rho'(s)}{s} \right) \right|, \\
\sup_{|q_k|_2 \le 1} \left\| \sum_{i,j=1}^d q_{1,i} q_{2,j} \partial_j c_i(t) \right\|_2 &\le \|\rho_{\sigma}^{(4)}\|_{\infty} + 8 \left\| \frac{\rho_{\sigma}'''(|\cdot|)}{|\cdot|} \right\|_{\infty} \\
&= \sigma^{-4} \left( \|\rho^{(4)}\|_{\infty} + \sup_{s>0} \left| \frac{\rho'''(s)}{s} \right| \right),
\end{align*}
$$

and we finally obtain

$$
B_{22} \le (-\rho''(0))^{-2} \left( \|\rho^{(4)}\|_{\infty} + 7 \sup_{s>0} \left| \frac{\rho'''(s)}{s} \right| + 36 \sup_{s>0} \left| \frac{1}{s^2} \left( \rho''(s) - \frac{\rho'(s)}{s} \right) \right| \right).
$$

In view of (29) and (30), another application of Lemma A.2 therefore yields $B_{22} = \mathcal{O}(1)$.

<figure><img src="image_7.png" /></figure>
## A.4 Near region bound

**Lemma A.4** For all $\mathfrak{D}_g(x, x') \le r_{\text{near}} \stackrel{\text{def.}}{=} 1/\sqrt{5}$ and all $v \in \mathbb{C}^d$ we have $-K^{(02)}(x, x') [v, v] \ge \bar{\epsilon}_2 |v|_x^2$ where $\bar{\epsilon}_2 \stackrel{\text{def.}}{=} 0.6 < r_{\text{near}}^{-2}$.

*Proof.* Let $r > 0$ and $x, x'$ be such that $\mathfrak{D}_g(x, x') \le r$. For any $v \in \mathbb{C}^d$, defining $t = x - x'$, we have

$$ 
\begin{aligned}
-K^{(02)}(x, x') [v, v] &= -v^T \nabla^2 K(t)v \\
&= -\frac{\rho'(\sigma^{-1}|t|)}{\sigma^{-1}|t|} \sigma^{-2}|v|^2 - \left[ \sigma^{-2} \rho''(\sigma^{-1}|t|) - \sigma^{-1} \frac{\rho'(\sigma^{-1}|t|)}{|t|} \right] \frac{(v \cdot t)^2}{|t|^2} \\
&\geq \left( -\frac{\rho'(\sigma^{-1}|t|)}{\sigma^{-1}|t|} - \left| \rho''(\sigma^{-1}|t|) - \frac{\rho'(\sigma^{-1}|t|)}{\sigma^{-1}|t|} \right| \right) \sigma^{-2}|v|^2 \\
&\geq \bar{\epsilon}_2(r) |v|_x^2,
\end{aligned}
$$

where

$$ 
\begin{aligned}
\bar{\epsilon}_2(r) &\stackrel{\text{def.}}{=} \frac{1}{(-\rho''(0))} \inf \left\{ -\frac{\rho'(s)}{s} - \left| \rho''(s) - \frac{\rho'(s)}{s} \right|, s\sqrt{-\rho''(0)} \in (-r, r) \right\} \\
&= \frac{1}{(-\rho''(0))} \inf \left\{ -\rho''(0) + \rho''(0) - \frac{\rho'(s)}{s} - \left| \rho''(s) - \frac{\rho'(s)}{s} \right|, s\sqrt{-\rho''(0)} \in (-r, r) \right\}.
\end{aligned}
$$

Proceeding as in the proof of Lemma A.2, we obtain:

$$ \left| \frac{\rho'(s)}{s} - \rho''(0) \right| \le (-\rho''(0)) \frac{2(e^{s^2/4} - 1)}{d + 4} \quad \text{and} \quad \left| \rho''(s) - \frac{\rho'(s)}{s} \right| \le (-\rho''(0)) \frac{s^2 e^{s^2/4}}{d + 4}. $$

If $s \le 1$ then

$$ \frac{2(e^{s^2/4} - 1)}{d + 4} \le \frac{e^{1/4} - 1}{3} \le 0.1 \quad \text{and} \quad \frac{s^2 e^{s^2/4}}{d + 4} \le \frac{e^{1/4}}{6} \le 0.3, $$

which shows $\bar{\epsilon}_2(r) \ge 0.6$ as soon as $r \le 1/\sqrt{5} \le 1/\sqrt{d+2}$. $\square$

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

## A.5 Far region bound

**Lemma A.5** For all $\mathfrak{D}_g(x, x') \ge r_{\text{near}}$ we have $|K(x, x')| \le 1 - \bar{\epsilon}_0$ where $\bar{\epsilon}_0 \stackrel{\text{def.}}{=} 0.07$.

*Proof.* If $\mathfrak{D}_g(x, x') = \sigma^{-1}|x - x'|/\sqrt{d+2} \ge r_{\text{near}} = 1/\sqrt{5}$ then $\sigma^{-1}|x - x'| \ge r_{\text{near}}\sqrt{d+2} \ge 2/\sqrt{5}$. We hence need to prove that $\rho(s) \le 1 - \bar{\epsilon}_0$ for every $s \ge 2/\sqrt{5}$. Using Lemma A.2, we obtain, for every $s > 0$:

$$ \rho(s) \le 0.8 \frac{2^{d/2} \Gamma(d/2 + 1)}{s^{1/3+d/2}} \le 0.8 \frac{2^{3/2} \Gamma(3/2 + 1)}{s^{4/3}}. $$

The right hand side is smaller than 0.9 on $(5/2, +\infty)$, but is not always smaller than 1 on $[2/\sqrt{5}, 5/2]$. On $[2/\sqrt{5}, 5/2]$, we use the other bound provided by Lemma A.2, namely:

$$ \rho(s) \le b(s) \stackrel{\text{def.}}{=} 1 - \frac{1}{d/2 + 1} \frac{s^2}{4} + \frac{1}{(d/2 + 1)(d/2 + 2)} \left( e^{s^2/4} - 1 - \frac{s^2}{4} \right). $$

We notice that $b$ is decreasing on $[2/\sqrt{5}, 2\sqrt{\log(3+d/2)}]$ and increasing on $[2\sqrt{\log(3+d/2)}, 5/2]$. Moreover, we have $b(2/\sqrt{5}) \le 0.93$ and $b(5/2) \le 0.63$. Hence, we finally obtain $\rho(s) \le 0.93$ for every $s \ge 2/\sqrt{5}$. $\square$

We now bound the kernel width $W(h, s)$ of $K$.

**Lemma A.6** If $h = \mathcal{O}(1)$ then $W(h, s) = \mathcal{O}(s^{2/(d+1)})$.
Proof. From the computations conducted in the proof of Lemma A.3 and the fact $J_{d/2}(t) = O(t^{-1/2})$ as $t \to +\infty$, we obtain that, for every $i, j \in \{0, 1, 2\}$ with $i + j \le 3$,

$$ \|K^{(ij)}(x, x')\|_{x,x'} \le \frac{1}{(\sigma^{-1}t)^{(d+1)/2}} $$

as soon as $\sigma^{-1}\|x - x'\| \gtrsim t$. As a result, provided $\Delta \gtrsim s^{2/(d+1)}$, for every $(x_k)_{k=1}^s \in S_\Delta$ we have

$$ \sum_{k=2}^{s} \| K^{(ij)}(x_1, x_k) \|_{x_1, x_k} \lesssim \sum_{k=2}^{s} \frac{1}{s} = O(1). $$

□

## A.6 Gradient bounds

For $j \in \{0, 1, 2\}$, we have $L_j(\omega) = \|D_j[\varphi_\omega](x)\|_x = (\sigma\|\omega\|_2/\sqrt{-\rho''(0)})^j$. Moreover, using the notations of [72, Section 5.1], we have

$$ \begin{aligned} \|D_2[\varphi_\omega](x) - D_2[\varphi_\omega](x') [\tau_{x \to x'}, \tau_{x \to x'}]\|_x & = \frac{\sigma^2}{(-\rho''(0))} \|\omega\|_2^2 |e^{-i\langle\omega, x\rangle} - e^{-i\langle\omega, x'\rangle}| \\ & \le \frac{\sigma^2}{(-\rho''(0))} \|\omega\|_2^2 |\langle\omega, x - x'\rangle| \\ & \le \frac{\sigma^2}{(-\rho''(0))} \|\omega\|_2^3 |x - x'| \\ & = (\sigma\|\omega\|_2/\sqrt{-\rho''(0)})^3 \mathfrak{D}_g(x, x'). \end{aligned} $$

This yields $L_3(\omega) \le (\sigma\|\omega\|_2/\sqrt{-\rho''(0)})^3$. Since $\omega \sim U(B(0, \sigma^{-1}))$, we have $\sigma\|\omega\|_2 \le 1$, and hence

$$ \|L_j\|_{\infty} \le (-\rho''(0))^{-j/2} $$

for every $j \in \{0, ..., 3\}$. We can therefore set $\bar{L}_j = (-\rho''(0))^{-j/2} = O(1)$ (see point (i) in the remark below Assumption 2 in [72, Section 5.2]).