Maxime Ferreira Da Costa

Yuejie Chi
email: yuejiechi@cmu.edu

Local Geometry of Nonconvex Spike Deconvolution from Low-Pass Measurements

Keywords: nonconvex spike deconvolution, preconditioned gradient descent, local geometry Contents

Spike deconvolution is the problem of recovering the point sources from their convolution with a known point spread function, which plays a fundamental role in many sensing and imaging applications. In this paper, we investigate the local geometry of recovering the parameters of point sources-including both amplitudes and locations-by minimizing a natural nonconvex least-squares loss function measuring the observation residuals. We propose preconditioned variants of gradient descent (GD), where the search direction is scaled via some carefully designed preconditioning matrices. We begin with a simple fixed preconditioner design, which adjusts the learning rates of the locations at a different scale from those of the amplitudes, and show it achieves a linear rate of convergence-in terms of entrywise errors-when initialized close to the ground truth, as long as the separation between the true spikes is sufficiently large. However, the convergence rate slows down significantly when the dynamic range of the source amplitudes is large. To bridge this issue, we introduce an adaptive preconditioner design, which compensates for the learning rates of different sources in an iteration-varying manner based on the current estimate. The adaptive design provably leads to an accelerated convergence rate that is independent of the dynamic range, highlighting the benefit of adaptive preconditioning in nonconvex spike deconvolution. Numerical experiments are provided to corroborate the theoretical findings.

Introduction

Spike deconvolution, also known as super-resolution [START_REF] Donoho | Superresolution via sparsity constraints[END_REF], is the task of recovering a stream of point sources from their convolution with a point spread function (PSF). This classical problem is at the core of many sensing and imaging modalities, including but not limited to radar, sonar, optical imaging, neuroimaging, and communication systems [START_REF] Potter | Sparsity and compressed sensing in radar imaging[END_REF][START_REF] Zhu | Super-resolution in SAR imaging: Analysis with the atomic norm[END_REF][START_REF] Zhu | Faster STORM using compressed sensing[END_REF][START_REF] Berger | Sparse channel estimation for multicarrier underwater acoustic communication: From subspace methods to compressed sensing[END_REF]. The PSF, which models the physical limitations of the imaging device involved in the experimental process, is commonly assumed to act as a band-limited and shift-invariant low-pass filter on the point sources [START_REF] Lindberg | Mathematical concepts of optical superresolution[END_REF][START_REF] Born | Principles of optics: electromagnetic theory of propagation, interference and diffraction of light[END_REF]. The sharpness of the original sources, modeled by Dirac impulses, is degraded through the convolution process, introducing undesirable ambiguity on the complex amplitudes and locations of the sources. The spike deconvolution task amounts to inverting the low-passing effects of the PSF, and to recovering the original sources as precisely as possible.

There is a rich literature on algorithmic investigations of the spike deconvolution problem, ranging from classical root-finding methods such as Prony's method, subspace methods such as MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF][START_REF] Liao | MUSIC for single-snapshot spectral estimation: Stability and super-resolution[END_REF], ESPRIT [START_REF] Roy | ESPRIT-estimation of signal parameters via rotational invariance techniques[END_REF] and matrix pencil [START_REF] Moitra | Super-resolution, extremal functions and the condition number of Vandermonde matrices[END_REF], to more recent optimization methods such as atomic norm minimization (a.k.a. total variation minimization) [START_REF] Castro | Exact reconstruction using Beurling minimal extrapolation[END_REF][START_REF] Tang | Compressed sensing off the grid[END_REF][START_REF] Candès | Towards a mathematical theory of super-resolution[END_REF][START_REF] Chi | Harnessing sparsity over the continuum: Atomic norm minimization for superresolution[END_REF][START_REF] Heckel | Super-resolution radar[END_REF] and basis pursuit [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF][START_REF] Herman | High-resolution radar via compressed sensing[END_REF]. While classical approaches harness the algebraic properties of complex exponentials by mapping the observations onto a low-dimensional linear subspace to recover the parameters of interest, optimization methods, on the other hand, attempt to recover the parameters via minimizing some carefully-designed loss function. As such, optimization methods tend to be more versatile in adapting to different imaging modalities, as well as amenable to modern advances in large-scale optimization. Inspired by the development of compressive sensing [START_REF] Candès | Near-Optimal Signal Recovery From Random Projections: Universal Encoding Strategies?[END_REF][START_REF] Donoho | Compressed sensing[END_REF], initial approaches for spike deconvolution relies on a discretization of the spike locations, and then attempts to recover a sparse solution using sparsity-promoting convex relaxations such as the LASSO [START_REF] Zhu | Faster storm using compressed sensing[END_REF][START_REF] Tang | Sparse recovery over continuous dictionaries-just discretize[END_REF]. However, the fundamental issue of basis mismatch [START_REF] Chi | Sensitivity to basis mismatch of compressed sensing for spectrum analysis and beamforming[END_REF] inherent to the discretization process may significantly hinder the localization performance, and increasing the grid size to reach finer precision levels leads to higher computational cost. Therefore, there has been a surge of interest in developing provably correct convex programs-such as the atomic norm minimization framework mentioned earlier-for spike deconvolution over the continuum in recent years, with strong performance guarantees developed under sufficient separations between the point sources [START_REF] Fernandez-Granda | Super-resolution of point sources via convex programming[END_REF][START_REF] Da | A tight converse to the spectral resolution limit via convex programming[END_REF][START_REF] Li | Approximate support recovery of atomic line spectral estimation: A tale of resolution and precision[END_REF][START_REF] Da | On the stable resolution limit of total variation regularization for spike deconvolution[END_REF]. Nonetheless, the atomic norm framework requires solving a semidefinite program whose complexity scales at least cubically with respect to the signal length; and therefore is computationally expensive and memory inefficient. In addition, although it is in principle possible to examine the so-called dual polynomial to localize the sources [START_REF] Candès | Towards a mathematical theory of super-resolution[END_REF], it often boils down to an additional post-processing step on the output of the convex program to recover the source parameters, which may hamper the guarantee of the overall procedure.

Motivated by the recent success of nonconvex methods, especially simple first-order methods, in various signal estimation and machine learning tasks [START_REF] Chi | Nonconvex optimization meets low-rank matrix factorization: An overview[END_REF][START_REF] Chizat | On the global convergence of gradient descent for over-parameterized models using optimal transport[END_REF], we are interested in understanding the efficacy of first-order methods in nonconvex spike deconvolution. In fact, first-order methods have already been popular empirically for spike deconvolution, but little is known about their theoretical underpinnings [START_REF] Huang | Super-resolution image reconstruction for high-density threedimensional single-molecule microscopy[END_REF]. As a first step, this paper focuses on the local geometry and performance guarantees of recovering the parameters of point sources-including both amplitudes and locations-by minimizing a natural nonconvex least-squares loss function measuring the observation residuals.

Observation model

Formally, we formulate the spike deconvolution problem as follows. Consider a vector of 2r parameters θ ⋆ = [a ⋆ 1 , . . . , a ⋆ r , τ ⋆ 1 , . . . , τ ⋆ r ] T , where a ⋆ ℓ ∈ C and τ ⋆ ℓ ∈ R correspond to the complex amplitude and location of the ℓ-th spike, respectively, ℓ = 1, . . . , r. Denoting by M the set of Radon measures over the reals, we assume that the point source signal µ ⋆ ∈ M to resolve is of the form

µ ⋆ = µ(θ ⋆ ) = r ℓ=1 a ⋆ ℓ δ τ ⋆ ℓ , (1)
where δ τ stands for the Dirac function located at τ ∈ R. Let us further denote the largest and the smallest amplitude of the spikes as

a ⋆ max = max 1≤ℓ≤r |a ⋆ ℓ | = ∥a ⋆ ∥ ∞ , a ⋆ min = min 1≤ℓ≤r |a ⋆ ℓ | .
The dynamic range of the measure µ ⋆ , an important quantity that will be used repetitively later, is thus defined as a ⋆ max /a ⋆ min . Denoting by g ∈ L 1 (R) the PSF, the temporal signal x ∈ L 1 (R) resulting from the convolution of the point source signal µ ⋆ and the PSF g reads

x(τ ) = (g * µ ⋆ ) (τ ) = r ℓ=1 a ⋆ ℓ g(τ -τ ⋆ ℓ ), (2)
where * denotes the convolution product.

A versatile observation model commonly encountered in practice considers the measurements to be taken from a uniform sampling of the Fourier transform of the temporal signal x. Denote by F(•) the Fourier transform of a measure µ lying in M, given by

F(µ)(f ) = R e -i2πf τ dµ(τ ), ∀µ ∈ M, ∀f ∈ R. (3)
Denote by G = F(g) and X = F(x) the Fourier transform of g and x, respectively. We assume that the PSF is band-limited within the bandwidth B = 1 so that it constrains no frequency greater than 1/2, i.e. G(f ) = 0 for |f | >1 2 . 1 For convenience, we assume that an odd number N = 2n + 1 of measurements are taken in the Fourier domain, uniformly spaced in the bandwidth

[-1 2 , 1 2 ]. The sampled signal x ∈ C N writes x = Φ(µ ⋆ ) = diag(g) F(µ) - n N , . . . , F(µ) n N ⊤ , (4)
where Φ : M → C N represents the observation operator and g ∈ C N is a vector with generic term

g k = G k N for k = -n, . . . , n.
Up to a scaling, we assume g to have unit Euclidean norm, i.e. ∥g∥ 2 = 1. The goal of spike deconvolution is thus to recover the measure µ ⋆ , or equivalently, the parameter θ ⋆ , from x.

Our contributions

In the rest of this paper, we assume that the model order r is known, and consider a natural nonconvex loss function, which aims to minimize the quadratic loss of the parameters θ = [a 1 , . . . , a r , τ 1 , . . . , τ r ] T of the Radon measures, given by

min θ L(θ) = 1 2 ∥Φ(µ(θ)) -x∥ 2 2 . (5)
Due to the nonlinear form of the parameters, the loss function L(θ) is clearly nonconvex. As a first step towards nonconvex spike deconvolution, we are interested in understanding the local geometry of the loss function [START_REF] Berger | Sparse channel estimation for multicarrier underwater acoustic communication: From subspace methods to compressed sensing[END_REF] and its implications on the computation efficacy of first-order methods. Without loss of generality, we assume the autocorrelation function h of the PSF g to be a triangular low pass function, i.e.

h(τ ) = ∞ -∞ f (u + τ )f (u)du = C sin 2 πτ 2 πτ 2 2 , (6)
where C > 0 is a constant, and that their respective Fourier transform H(f ) and G(f ) are linked through the relation

G(f ) = H(f ) for all f ∈ R. It comes that g k = G k N = H k N = 1 n + 1 1 - |k| n + 1 , (7)
for all k = -n, . . . n after rescaling with the constraint ∥g∥ 2 = 1. Note that the triangular low-pass function and the Fejér kernel-its discrete counterpart-play an important role in the deconvolution literature and have been extensively proposed as a convolution kernel to evaluate the norm and distance between Radon measures [START_REF] Candès | Towards a mathematical theory of super-resolution[END_REF]. Additionally, our main results can be re-derived with any other bandlimited PSF g by following analogous reasoning, as long as it autocorrelation h is an absolutely integrable function.

Concretely, we propose and analyze preconditioned variants of gradient descent (GD), where the search direction is scaled via some carefully designed preconditioning matrices. Our contributions are summarized as follows.

• We begin with a simple fixed preconditioner design, which adjusts the learning rates of the locations at a different scale from those of the amplitudes, and show it achieves a linear rate of convergence-in terms of entrywise errors-when initialized close to the ground truth, as long as the separation between the true spikes is sufficiently large. However, the convergence rate slows down significantly when the dynamic range of the source amplitudes is large.

• To bridge this issue, we introduce an adaptive preconditioner design, which compensates the learning rates of different sources in an iteration-varying manner based on the current estimate. The adaptive design provably leads to an accelerated convergence rate that is independent of the dynamic range, highlighting the benefit of adaptive preconditioning in nonconvex spike deconvolution.

Our result is based on understanding the geometric properties of scaled Hessian matrices via a set of novel summation bounds on the absolute sums of sampled Fejér kernels and higher-order derivatives, which might be of independent interest in other contexts.

Related work

The closest work to ours on recovering spike signals from low-pass observations using nonconvex optimization is [START_REF] Traonmilin | The basins of attraction of the global minimizers of the non-convex sparse spike estimation problem[END_REF]. The radius of the basin of attraction for gradient descent is characterized whenever the observation operator satisfies the restricted isometry property over the set of well-separated sparse measures. Although the problem setup is versatile, specializing this result in our context of low-pass measurements yields a convergence region whose size scales inversely with the number of sources, which is pessimistic when the number of sources is large. Moreover, the analysis in [START_REF] Traonmilin | The basins of attraction of the global minimizers of the non-convex sparse spike estimation problem[END_REF] focuses on the Euclidean error of the parameters, while we focus on the entrywise error, which is more meaningful for gauging the recovery quality of the point sources. Projected gradient methods, which merge pairs of colliding spikes at each iteration, have been proposed in [START_REF] Traonmilin | Projected gradient descent for non-convex sparse spike estimation[END_REF][START_REF] Bénard | Fast off-the-grid sparse recovery with over-parametrized projected gradient descent[END_REF] but without theoretical convergence guarantees.

Our work can be viewed as falling into a growing line of research on developing provably efficient nonconvex methods-especially first-order methods-for high-dimensional signal estimation, examples including phase retrieval [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF][START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF], low-rank matrix estimation [START_REF] Sun | Guaranteed matrix completion via non-convex factorization[END_REF][START_REF] Chi | Nonconvex optimization meets low-rank matrix factorization: An overview[END_REF], blind (sparse) deconvolution [START_REF] Li | Rapid, robust, and reliable blind deconvolution via nonconvex optimization[END_REF][START_REF] Chen | Convex and nonconvex optimization are both minimax-optimal for noisy blind deconvolution under random designs[END_REF][START_REF] Lee | Blind recovery of sparse signals from subsampled convolution[END_REF], dictionary learning [START_REF] Sun | Complete dictionary recovery over the sphere I: Overview and the geometric picture[END_REF][START_REF] Gilboa | Efficient dictionary learning with gradient descent[END_REF], multi-channel sparse deconvolution [START_REF] Shi | Manifold gradient descent solves multi-channel sparse blind deconvolution provably and efficiently[END_REF][START_REF] Qu | Exact recovery of multichannel sparse blind deconvolution via gradient descent[END_REF], and so on. In particular, the preconditioned gradient methods considered in this paper are motivated by [START_REF] Tong | Accelerating ill-conditioned low-rank matrix estimation via scaled gradient descent[END_REF][START_REF] Tong | Scaling and scalability: Provable nonconvex low-rank tensor estimation from incomplete measurements[END_REF][START_REF] Tong | Low-rank matrix recovery with scaled subgradient methods: Fast and robust convergence without the condition number[END_REF], which demonstrated that preconditioning can efficiently accelerate the convergence of gradient descent in ill-conditioned low-rank estimation.

Notation and paper organization

Vectors and matrices are denoted by boldface and capital boldface letters, respectively. Vectors x ∈ C N with odd dimension N = 2n + 1 are indexed between -n and n, so that x = [x -n , . . . , x n ] ⊤ for convenience. Transpose and Hermitian transpose of a vector or a matrix A are denoted by A ⊤ and A H , respectively. Furthermore, the adjoint of the operator Φ is denoted by Φ * . We write 1 d and 0 d the all-one and null vector (or matrix) in dimension d, respectively. With a slight abuse of notation, we denote by |a|, |a| 2 , a -1 the vector with entries equal to the modulus, the squared modulus, and the inverse of the entries of a, respectively. The element-wise product between two vectors a and a ′ is written as a ⊙ a ′ . We denote by D ℓ the space of ℓ-times differentiable functions of the real variable. For any function h ∈ D ℓ , we write its ℓth derivative h (ℓ) . We denote by ⟨• , •⟩ and ⟨• , •⟩ R = ℜ (⟨• , •⟩) the usual inner product and real inner product between Radon measure, respectively. Additionally, we let δ (ℓ) ∈ M be the functional which satisfies

g (ℓ) (τ ) = δ (ℓ) τ , g , ∀g ∈ D ℓ , ∀τ ∈ R. ( 8
)
Fejér kernel We denote by F N (•) the normalized Fejér kernel of order N = 2n + 1 defined by

F N (t) = 1 n + 1 n k=-n 1 - |k| n + 1 e i2πkt = sin 2 (π(n+1)t) (n+1) 2 sin 2 (πt) if t / ∈ Z 1 otherwise. ( 9
)
The Fejér kernel is a trigonometric polynomial, hence it is infinitely differentiable. We point out that the second derivative of F N (•) at the origin satisfies

F ′′ N (0) = 1 n + 1 n k=-n -4π 2 k 2 1 - |k| n + 1 = - 2 3 π 2 n(n + 2) < 0. (10)
Some of its properties, key to this paper, are derived and discussed in Appendix A. The Fejér kernel plays an important role in the sequel as the Gramian Φ * Φ : M → M of the observation operator Φ is a convolution product with F N , i.e.

Φ * Φ(µ) = F N * µ, ∀µ ∈ M. (11)
Wrap-around distance For any set of r points τ = {τ 1 , . . . , τ r } ⊂ T, we denote by ∆(τ ) is minimal wrap-around distance, defined by

∆(τ ) ≜ min ℓ̸ =ℓ ′ inf p∈Z |τ ℓ -τ ℓ ′ + p| . (12)
The rest of this paper is organized as follows. Section 2 starts by defining the preconditioned gradient methods and two of its designs with provable local convergence guarantees, using a fixed preconditioner and an adaptive preconditioner in Section 2.2 and Section 2.3, respectively. Section 3 provides the analysis of the main theorems by controlling the conditioning of the scaled Hessian matrix of the loss function in a neighborhood of the ground truth. Numerical experiments are provided in Section 4 to corroborate our findings. Finally, a brief conclusion is drawn in Section 5.

2 How does preconditioning help local convergence?

Preconditioned gradient descent

Recognizing that the parameters corresponding to the amplitudes and locations may require different treatments, we consider iterates of preconditioned gradient descent (GD) to recover the ground truth parameters, where the preconditioner can possibly be iteration-varying. Given an initialization point θ 0 ∈ C 2r , the update sequence of preconditioned GD is obtained by successively moving oppositely along the direction of a linear transform of the gradient. More specifically, the update rule reads

θ k+1 = θ k -P k ∇L(θ k ), (13)
where

θ k = [a ⊤ k , τ ⊤ k ]
⊤ is the k-th iterate, P k ∈ C 2r×2r is a preconditioning matrix (also called preconditioner)

that can vary at each iteration; the choice of P k will be detailed momentarily. Here, it is worth noticing that there are no additional learning rates in [START_REF] Tang | Compressed sensing off the grid[END_REF], which can be thought of as already absorbed and set within the preconditioner P k . By analogy with the celebrated Newton-Raphson method, which selects

P k = ∇ 2 L(θ k ) -1
(which might however be computationally expensive), the role of the preconditioning matrix P k is to balance the local optimization landscape towards a quadratic function to improve the convergence rate towards a local minimum over the vanilla gradient method. By basic calculation, the gradient ∇L(θ) at point θ = [a 1 , . . . , a r , τ 1 , . . . , τ r ] T is given by

dL(θ) da j = Φ(δ τj ), Φ(µ(θ) -µ(θ ⋆ )) = δ τj , Φ * Φ(µ(θ) -µ(θ ⋆ )) = δ τj , F N * (µ(θ) -µ(θ ⋆ )) = r ℓ=1 a ℓ F N (τ j -τ ℓ ) - r ℓ=1 a ⋆ ℓ F N (τ j -τ ⋆ ℓ ) (14a)
for j = 1, . . . , r, and similarly,

dL(θ) dτ j = Φ(a j δ ′ τj ), Φ(µ(θ) -µ(θ ⋆ )) R = a j δ ′ τj , Φ * Φ(µ(θ) -µ(θ ⋆ )) R = a j δ ′ τj , F N * (µ(θ) -µ(θ ⋆ )) R = ℜ a j r ℓ=1 a ℓ F ′ N (τ j -τ ℓ ) - r ℓ=1 a ⋆ ℓ F ′ N (τ j -τ ⋆ ℓ ) (14b)
for j = 1, . . . , r.

In this paper, we are particularly interested in preconditioning matrices P k that are diagonally structured, so they do not add computation overhead compared with vanilla gradient methods. In the sequel, we study the basin of attraction and the convergence rate of preconditioned GD for two different preconditioning strategies. The first consists of selecting a time-invariant, diagonal preconditioning matrix P = P k whose role is to judiciously renormalize the learning rates between the amplitudes a ∈ C r and the locations τ ∈ R r which are of different units. The second strategy seeks to dynamically update the preconditioning matrix P k based on the current iterate to better approximate the inverse of the Hessian matrix around the point θ k and accelerate convergence.

Invariant preconditioning

In this section, we seek to recover the ground truth parameter θ ⋆ from an instance of the preconditioned GD algorithm [START_REF] Tang | Compressed sensing off the grid[END_REF] where the sequence of preconditioning matrices is constant, i.e. P = P k for all k ∈ N. We fix

P = diag 1 r -F ′′ N (0) -1 A -2 1 r , (15)
where A > 0 is an input parameter that controls the ratio between the learning rate applied to the amplitudes a k of the sources and that applied to the locations τ k of the sources throughout the iterative process.

Performance metric To gauge the performance, we define by S ∈ C 2r×2r the weighting matrix

S = diag a ⋆-1 -F ′′ N (0)1 r , , (16)
where F ′′ N (0) is given in [START_REF] Roy | ESPRIT-estimation of signal parameters via rotational invariance techniques[END_REF], and study the convergence properties of preconditioned GD in terms of the infinity norm weighted by the matrix S, i.e.

∥S (θ

k -θ ⋆ )∥ ∞ = max j a k,j -a ⋆ j a ⋆ j , -F ′′ N (0) τ k,j -τ ⋆ j . (17)
Intuitively, the role of this scaling is to analyze a unitless metric that decorrelates the error with the dynamic range of the sources and with the problem dimension, as we have -F ′′ N (0) = O(n) and that the error on the source locations is expected to be inversely proportional to the number of observation:

∥τ k -τ ⋆ ∥ ∞ = O n -1 .
The following theorem establishes the linear convergence of preconditioned GD with a fixed preconditioning matrix P whenever the input parameter A is properly set, and the initial point θ 0 is close enough to the ground truth θ ⋆ , as long as the true spikes are sufficiently separated.

Theorem 1 (Linear convergence with invariant preconditioner). Suppose that n ≥ 2 and that the input parameter A satisfies ∥a ⋆ ∥ ∞ ≤ 3 2 A. Moreover, assume that

η := 276.21 A 2 ∥a ⋆ ∥ ∞ (a ⋆ min ) 3 ((n + 1) ∆(τ ⋆ )) -2 < 1, (18)
then if the initial point θ 0 = [a ⊤ 0 , τ ⊤ 0 ] ⊤ satisfies ∥S (θ 0 -θ ⋆ )∥ ∞ ≤ 1 2 , (19)
the iterates {θ k } of preconditioned GD (13) with a fixed preconditioner (15) converge towards θ ⋆ according to

∥S (θ k -θ ⋆ )∥ ∞ ≤ 1 - 1 4 (a ⋆ min ) 2 A 2 (1 -η) k ∥S (θ 0 -θ ⋆ )∥ ∞ ( 20
)
for all k ∈ N.

Theorem 1 indicates that preconditioned GD admits a linear rate of convergence as long as the separation condition ∆(τ ⋆ ) is sufficiently large with respect to the dynamic range, i.e.

(n + 1) ∆(τ ⋆ ) ≳ ∥a ⋆ ∥ ∞ a ⋆ min 3/2 . ( 21
)
Additionally, our finding is independent of the number of sources r of the input measure, both in terms of the size of the basin of attraction (cf. ( 19)) and the convergence rate. Faster convergence rate are achieved for smaller values of the parameter 0 < η < 1, when the separation ∆(τ ⋆ ) of the true spikes is larger or the dynamic range

∥a ⋆ ∥ ∞ a ⋆
min of the amplitudes is smaller. However, even for small values of η, the convergence rate ρ predicted by Theorem 1 is lower bounded by

ρ ≥ 1 - 1 4 a ⋆ min A 2 ≥ 1 - 9 16 a ⋆ min ∥a ⋆ ∥ ∞ 2 . ( 22
)
This suggests that a high dynamic range will lead to a slow convergence rate, independently of the separation ∆(τ ⋆ ). Additionally, the convergence guarantees established in Theorem 1 demand to adjust the input parameter A as a function of ∥a ⋆ ∥ ∞ , which can be impractical in scenarios with no postulate on the norm of the source amplitudes.

Adaptive preconditioning

In order to mitigate the limitations of the fixed preconditioning strategy presented in Section 2.2, we propose to study an instance of preconditioned GD where the preconditioner P k varies at each iteration and is selected as a function of the current iterate θ k . In particular, we fix

P k = diag 1 r -F ′′ N (0) -1 |a k | -2 . ( 23
)
Similar to Theorem 1, the next theorem guarantees a linear convergence rate of the iterates towards the ground truth θ ⋆ , provided a good enough initialization point θ 0 , as long as the true spikes are sufficiently separated.

Theorem 2 (Linear convergence with adaptive preconditioner). Suppose that n ≥ 2, and assume that

γ := 11.60 ∥a∥ ∞ a ⋆ min ((n + 1)∆(τ ⋆ )) -2 < 1 2 , ( 24
)
then if the initial point θ 0 = [a ⊤ 0 , τ ⊤ 0 ] ⊤ satisfies ∥S (θ 0 -θ ⋆ )∥ ∞ ≤ 1 - 2 3 , (25)
the iterates {θ k } of preconditioned GD (13) with an adaptive preconditioner (23) converge towards θ ⋆ according to

∥S (θ k -θ ⋆ )∥ ∞ ≤ 1 2 + γ k ∥S (θ 0 -θ ⋆ )∥ ∞ (26)
for all k ∈ N.

Theorem 2 guarantees that preconditioned GD with an adaptive preconditioner achieves a constant linear rate of convergence in a similar basin of attraction, provided that the separation condition ∆(τ ⋆ ) is sufficiently large with respect to the dynamic range, i.e.

(n + 1) ∆(τ ⋆ ) ≳ ∥a ⋆ ∥ ∞ a ⋆ min 1/2 , ( 27
)
which is much weaker than the requirement for the case using a fixed preconditioner, as indicated in Theorem 1. Consequently, this highlights the benefit of adaptive preconditioning in accelerating the convergence in the presence of high dynamic ranges for nonconvex spike deconvolution.

Remark. We have not attempted to fully optimize the constants in the above theorems. Therefore, their values are set in a quite pessimistic fashion; see Section 4 for numerical experiments.

Analysis

This section is devoted to proving the two main results of this paper comprised in Theorem 1 and Theorem 2. Before entering the core of the proofs, we first provide some warm-up analysis that will be required in the latter proofs.

Preliminaries

Contraction of entrywise errors

The convergence analysis of the preconditioned gradient method presented in Theorem 1 and Theorem 2 calls for understanding of the contraction properties of the sequence {∥S (θ k -θ ⋆ )∥ ∞ }. Starting from the update rule [START_REF] Tang | Compressed sensing off the grid[END_REF], leveraging ∇L(θ ⋆ ) = 0, and applying the fundamental theorem of calculus, we have

S(θ k+1 -θ ⋆ ) = S (θ k -P k ∇L(θ k ) -θ ⋆ ) = S(θ k -θ ⋆ ) -SP k (∇L(θ k ) -∇L(θ ⋆ )) = S(θ k -θ ⋆ ) -SP k 1 0 ∇ 2 L (θ ⋆ + u (θ k -θ ⋆ )) du (θ k -θ ⋆ ) = I -SP k 1 0 ∇ 2 L (θ ⋆ + u (θ k -θ ⋆ )) du S -1 S (θ k -θ ⋆ ) . (28)
Let

S k = {θ ⋆ + u (θ k -θ ⋆ ) |u ∈ [0, 1
]} be the line segment that connects θ ⋆ and θ k in C 2r , and denote by ρ k the quantity

ρ k ≜ max θ∈S k I -SP k H(θ)S -1 ∞ , (29)
where

H(θ) = ∇ 2 L(θ)
is the Hessian of the loss function L at point θ. Continuing to bound [START_REF] Chi | Nonconvex optimization meets low-rank matrix factorization: An overview[END_REF] yields

∥S (θ k+1 -θ ⋆ )∥ ∞ ≤ I -SP k 1 0 H (θ ⋆ + u (θ k -θ ⋆ )) du S -1 ∞ ∥S (θ k -θ ⋆ )∥ ∞ = 1 0 I -SP k H (θ ⋆ + u (θ k -θ ⋆ )) S -1 du ∞ ∥S (θ k -θ ⋆ )∥ ∞ ≤ 1 0 I -SP k H (θ ⋆ + u (θ k -θ ⋆ )) S -1 ∞ du ∥S (θ k -θ ⋆ )∥ ∞ ≤ max θ∈S k I -SP k H(θ)S -1 ∞ ∥S (θ k -θ ⋆ )∥ ∞ = ρ k ∥S (θ k -θ ⋆ )∥ ∞ . (30)
Hence, the crux of the convergence analysis is to show that ρ k < 1 (and control the size of ρ k ) uniformly over the segment S k whenever the point θ k lies in an appropriate region centered around the ground truth θ ⋆ . Further analysis towards that goal requires an explicit derivation of the Hessian matrix H(θ), which is done next.

Hessian decomposition

Recall that H(θ) = ∇ 2 L(θ) ∈ C 2r×2r is the Hessian matrix of the loss function L in (5) at the point θ = [a ⊤ , τ ⊤ ] ⊤ . We decompose H(θ) as

H(θ) = H a,a (θ) H a,τ (θ) H H a,τ (θ) H τ,τ (θ) , (31)
where each block is of size r × r, with generic terms

[H a,a (θ)] (i,j) = d 2 f (θ) da i da j , [H a,τ (θ)] (i,j) = d 2 f (θ) da i dτ j , [H τ,τ (θ)] (i,j) = d 2 f (θ) dτ i dτ j
for i, j = 1 . . . , r. A direct calculation of the Hessian matrix H(θ) (see, e.g. [START_REF] Traonmilin | The basins of attraction of the global minimizers of the non-convex sparse spike estimation problem[END_REF]) yields a decomposition of the form

H(θ) = G(θ) + E(θ), (32)
where the terms G(θ) ∈ C 2r×2r and E(θ) ∈ C 2r×2r are described in the sequel.

Structure of G(θ)

The matrix G(θ) can be written as

G(θ) = diag 1 r -F ′′ N (0)a H D(τ ) diag 1 r -F ′′ N (0)a , (33)
with D(τ ) ∈ C 2r×2r given with a block structure

D(τ ) = D 0 (τ ) D 1 (τ ) D 1 (τ ) H D 2 (τ ) . ( 34
)
The entries of the blocks D 0 (τ ), D 1 (τ ), D 2 (τ ) ∈ C r×r are composed of

[D 0 (τ )] (i,j) = F N (τ i -τ j ), (35a)
[D 1 (τ )] (i,j) = -F ′ N (τ i -τ j )/ -F ′′ N (0), (35b)
[D 2 (τ )] (i,j) = F ′′ N (τ i -τ j )/F ′′ N (0) (35c)
for all i, j = 1, . . . , r. As shall be seen, the matrix G(θ) is a relatively well-conditioned matrix whose spectrum can be controlled as a function of the separation parameter between the spikes (n + 1)∆(τ ⋆ ), the dynamic range of the amplitudes

∥a∥ ∞ a ⋆ min
, and the distance of θ to the ground truth parameter θ ⋆ .

Structure of E(θ)

The matrix E(θ) is given by the block structure decomposition

E(θ) = 0 r×r E 1 (θ) E 1 (θ) H E 2 (θ) , (36)
where the entries of E 1 (θ), E 2 (θ) ∈ C r×r are given as

[E 1 (θ)] (i,j) = δ ′ τj , F N * (µ(θ) -µ(θ ⋆ )) i = j; 0 i ̸ = j; (37a) [E 2 (θ)] (i,j) = a j δ ′′ τj , F N * (µ(θ) -µ(θ ⋆ )) i = j; 0 i ̸ = j. (37b)
It is worth noting that E(θ ⋆ ) = 0, hence E can be interpreted as a perturbation term that grows as θ deviates from θ ⋆ . With the above preliminaries, we are now ready to prove the main results of this paper. The following two sections present our convergence analyses for preconditioned GD using a fixed preconditioner and an adaptive preconditioner, respectively.

Proof of Theorem 1

We recall that

P k = P = diag 1 r -F ′′ N (0) -1 A -2 1 r
for all k ∈ N when the preconditioner is fixed. The contraction analysis [START_REF] Huang | Super-resolution image reconstruction for high-density threedimensional single-molecule microscopy[END_REF] presented in Section 3.1.1 suggests that the contraction rate in Theorem 1 is controlled by the quantity SP H(θ)S -1 -I ∞ in neighborhood around the ground truth θ ⋆ . The next theorem, whose proof is deferred to Appendix B, provides a uniform bound on this quantity on a neighborhood of the ground truth.

Theorem 3 (Uniform bound of the Hessian). Suppose that n ≥ 2, (n + 1)∆(τ ⋆ ) ≥ 16.5 and ∥τ -

τ ⋆ ∥ ∞ ≤ 1 4 ∆(τ ⋆ ). Let θ k = [a ⊤ k , τ ⊤ k ] ⊤ and assume that A ≥ max 3 2 ∥a ⋆ ∥ ∞ , ∥a k ∥ ∞ then the exist two positive constants K ∆ = 2.13, (38a)
K θ = 44.42, (38b)
such that for all θ = [a ⊤ , τ ⊤ ] ⊤ ∈ S k with ∥S (θ k -θ ⋆ )∥ ∞ < 1, we have that SP H(θ)S -1 -I ∞ ≤ 1 - a ⋆ min A 2 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 + (4K ∆ + K θ ∥S (θ k -θ ⋆ )∥ ∞ ) ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ ⋆ )) -2 (1 + ∥S (θ k -θ ⋆ )∥ ∞ ) 2 . ( 39
)
Theorem 3 provides a bound on SP H(θ)S -1 -I ∞ depending on the quantity ∥a ⋆ ∥∞ a ⋆ min

((n + 1)∆(τ ⋆ )) -2 ,
which can be made small enough under the hypothesis of Theorem 1. We proceed with the rest of the proof by induction. For the base case, it is trivial that (20) holds for k = 0. We start the induction by assuming that

∥S (θ k -θ ⋆ )∥ ∞ ≤ 1 - 1 4 a ⋆ min A 2 (1 -η) k ∥S (θ 0 -θ ⋆ )∥ ∞ (40)
holds for some k ∈ N. We begin by verifying the assumptions of Theorem 3.

• First, as the dynamic range

∥a ⋆ ∥ ∞ a ⋆ min ≥ 1
, it is easy to see that the hypothesis [START_REF] Herman | High-resolution radar via compressed sensing[END_REF] immediately implies that (n + 1)∆(τ ⋆ ) ≥ 16.6.

• By the definition [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF] and the induction hypothesis [START_REF] Lee | Blind recovery of sparse signals from subsampled convolution[END_REF], it follows that

∥τ k -τ ⋆ ∥ ∞ ≤ 1 -F ′′ N (0) ∥S (θ k -θ ⋆ )∥ ∞ ≤ 1 -F ′′ N (0) ∥S (θ 0 -θ ⋆ )∥ ∞ ≤ √ 3 2π 1 n + 1 ≤ 1 n + 1 ≤ ∆(τ ⋆ ) 4 . (41)
• Furthermore, we have that

∥a k ∥ ∞ ≤ ∥a ⋆ ∥ ∞ + ∥a k -a ⋆ ∥ ∞ ≤ ∥a ⋆ ∥ ∞ (1 + ∥S (θ k -θ ⋆ )∥ ∞ ) < ∥a ⋆ ∥ ∞ (1 + ∥S (θ 0 -θ ⋆ )∥ ∞ ) ≤ 3 2 ∥a ⋆ ∥ ∞ ≤ A. (42)
Hence, the assumptions of Theorem 3 hold, which yields

ρ k = max θ∈S k SP H(θ)S -1 -I ∞ ≤ 1 - (a ⋆ min ) 2 4A 2 + 9 4 4K ∆ + 1 2 K θ ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ ⋆ )) -2 ≤ 1 - (a ⋆ min ) 2 4A 2 1 -9 4K ∆ + 1 2 K θ A 2 ∥a ⋆ ∥ ∞ a ⋆ min 3 ((n + 1)∆(τ ⋆ )) -2 ≤ 1 - (a ⋆ min ) 2 4A 2 1 -276.21 A 2 ∥a ⋆ ∥ ∞ a ⋆ min 3 ((n + 1)∆(τ ⋆ )) -2 ≤ 1 - (a ⋆ min ) 2 4A 2 (1 -η) , (43)
where we substituted the definition (18) of η in the last line. It results from the iterative analysis (30) that the next update θ k+1 obeys

∥S (θ k+1 -θ ⋆ )∥ ∞ ≤ ρ k ∥S (θ k -θ ⋆ )∥ ∞ ≤ 1 - (a ⋆ min ) 2 4A 2 (1 -η) ∥S (θ k -θ ⋆ )∥ ∞ ≤ 1 - (a ⋆ min ) 2 4A 2 (1 -η) k+1 ∥S (θ 0 -θ ⋆ )∥ ∞ , (44)
which concludes the proof of Theorem 1.

Proof of Theorem 2

We proceed with the proof of Theorem 2 analogously to the proof of Theorem 1 presented in Section 3.2. First, we establish the following intermediate theorem that controls the conditioning of the scaled Hessian matrix SP k (θ)HS -1 uniformly over the segment S k as a function of the weighted infinity-norm distance ∥S (θ k -θ ⋆ )∥ ∞ . The proof of Theorem 4 is deferred to Appendix B.3.

Theorem 4 (Uniform bound of the Hessian). Suppose that n ≥ 2, (n + 1)∆(τ ⋆ ) ≥ 4.7 and ∥τ -τ ⋆ ∥ ∞ ≤ 1 4 ∆(τ ⋆ ) then the exists two positive constants

K ∆ ≤ 2.32, (45a)
K θ ≤ 75.80, (45b)
such that for all θ = [a ⊤ , τ ⊤ ] ⊤ ∈ S k satisfying ∥S (θ k -θ ⋆ )∥ ∞ < 1 we have that SP k H(θ)S -1 -I ∞ ≤ 1 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 -1 + (4K ∆ + K θ ∥S (θ k -θ ⋆ )∥ ∞ ) ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ ⋆ )) -2 (1 -∥S (θ k -θ ⋆ )∥ ∞ ) 2 . ( 46
)
The rest of the proof follows similarly by induction. For the base case, it is trivial that the initial point θ 0 verifies [START_REF] Li | Approximate support recovery of atomic line spectral estimation: A tale of resolution and precision[END_REF]. We now assume that θ k satisfies

∥S (θ k -θ ⋆ )∥ ∞ ≤ 1 2 + γ k ∥S (θ 0 -θ ⋆ )∥ ∞ (47)
for some k ∈ N. Let's verify the assumptions of Theorem 4.

• First, as the dynamic range

∥a ⋆ ∥ ∞ a ⋆ min ≥ 1
, the assumption [START_REF] Fernandez-Granda | Super-resolution of point sources via convex programming[END_REF] easily implies that (n + 1)∆(τ ⋆ ) ≥ 4.7.

• By the definition [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF] and the induction hypothesis [START_REF] Tong | Low-rank matrix recovery with scaled subgradient methods: Fast and robust convergence without the condition number[END_REF], it follows that

∥τ k -τ ⋆ ∥ ∞ ≤ 1 -F ′′ N (0) ∥S (θ k -θ ⋆ )∥ ∞ ≤ 1 -F ′′ N (0) ∥S (θ 0 -θ ⋆ )∥ ∞ ≤ √ 3 π 1 - 2 3 1 n + 1 ≤ 1 n + 1 ≤ ∆(τ ⋆ ) 4 . (48)
Hence, the assumptions of Theorem 4 hold. Noticing that the function f (u) := 1 (1-u) 2 is increasing over [0, 1), we have that

1 (1-∥S(θk-θ ⋆ ) 2 ∥ ∞ ) ≤ 1 2 .
Together with Theorem 4, this yields the bound

ρ k = max θ∈S k SP k H(θ)S -1 -I ∞ ≤ 1 2 + 4K ∆ + K θ 1 - 2 3 ∥a ⋆ ∥ ∞ 2a ⋆ min ((n + 1)∆(τ )) -2 ≤ 1 2 + 11.60 ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ ⋆ )) -2 ≤ 1 2 + γ < 1, (49)
where we substituted the definition [START_REF] Fernandez-Granda | Super-resolution of point sources via convex programming[END_REF] of γ in the third inequality. It results from the iterative analysis (30) that the next update θ k+1 satisfies

∥S (θ k+1 -θ ⋆ )∥ ∞ ≤ρ k ∥S (θ k -θ ⋆ )∥ ∞ ≤ 1 2 + γ ∥S (θ k -θ ⋆ )∥ ∞ ≤ 1 2 + γ k+1 ∥S (θ 0 -θ ⋆ )∥ ∞ , (50)
which concludes the proof of Theorem 2.

Numerical experiments

This section provides a numerical validation of Theorem 1 and Theorem 2. In the following experiments, the signal length is set to N = 65 (i.e. n = 32). The ground truth signal is composed of r = 6 sources placed in the interval [-1 2 , 1 2 ) while ensuring that (n + 1)∆(τ ⋆ ) ≥ 2, which is a more optimistic separation condition than what the theorems' statements suggest. Additionally, the dynamic range is denoted by κ = ∥a ⋆ ∥ ∞ a ⋆ min . The complex amplitudes a ⋆ ∈ C r are selected independently and uniformly at random in a complex annulus with bounds 1 ≤ |a ⋆ ℓ | ≤ κ. The input parameter of the invariant preconditioning scheme is set at

A = 3 2 ∥a ⋆ ∥ ∞ .
Size of the basin of attraction Of critical importance in the analysis of Theorem 1 and Theorem 2 is the distance between the initial parameter θ 0 and the ground truth θ ⋆ . We start by comparing the success rates of both preconditioning schemes on reconstructing the ground truth as a function of the initialization distance ∥S (θ 0 -θ ⋆ )∥ ∞ . In each experiment, the starting point θ 0 is drawn uniformly over the set of points equidistant to θ ⋆ . An experiment is labeled as a success if ∥S (θ 200 -θ ⋆ )∥ ∞ ≤ 10 -2 after 200 iterations. Figure 1 suggests that, for both schemes, the size of the basin of attraction is independent of the dynamic range κ, and is around the order of magnitude ∥S (θ 0 -θ ⋆ )∥ ∞ ≃ 1. This suggests that the numerical constants (1/2 and 1 -2/3 ∼ 0.184, respectively) set forth in Theorem 1 and Theorem 2, respectively, are pessimistic and nonconvex spike deconvolution performs in a much more benign manner than predicted by our theory, indicating room for further refinements. Linear convergence using a spectral initialization In practice, several ad hoc initialization methods could be envisaged to produce an initial point θ 0 that falls in the basin of attraction of the preconditioned gradient descent methods. Herein, we proceed by uniformly discretizing the spectral domain over N elements. Given the knowledge of the ground truth model order r, the initial locations τ 0 are selected as the r elements of the discrete grid whose weighted Fourier transform best describes the observation x. Mathematically, consider the following optimization problem

u 0 = arg min u 1 2 ∥x -diag(g)F N u∥ 2 2 s.t. ∥u∥ 0 ≤ r, (51)
where F N ∈ C N ×N is a discrete Fourier transform matrix, and the ℓ 0 -norm denotes the cardinality of the support. Writing I 0 = supp (u 0 ) ⊂ [-n, . . . , n] the support of the solution u 0 of (51), the parameter θ 0 is constructed in a second stage by selecting τ 0 = [ k1 N , . . . , kr N ] ⊤ where k ℓ ∈ I 0 , ℓ = 1, . . . , r and a 0 = u |I0 as the restriction of u to the elements in I 0 . The program (51) is itself a non-convex sparse reconstruction problem, which we approximate the solution using the orthogonal matching pursuit algorithm [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF]. The proposed initialization procedure offers several benefits over more classical methods: It is highly scalable, robust to high dynamic range, and does not involve any polynomial root finding subroutine.

Figure 2 pictures the convergence rate of preconditioned GD under the invariant and adaptive preconditioning schemes, respectively. For both schemes, θ 0 is selected according to the previously described initialization procedure. It can be seen that, although both preconditioning schemes ensure a linear converge of the iterate sequence, the convergence rate with a fixed preconditioner degrades as the dynamic range of the sources increases. In contrast, the one with an adaptive preconditioner remains unchanged. Additionally, the adaptive preconditioning scheme benefits from faster convergence rates for a given dynamic range. These experimental results corroborate the theoretical findings presented in Section 2. Noisy recovery We next examine the performance of preconditioned GD in the presence of noise. We assume observations of the form x = Φ(µ ⋆ ) + w, where w is white Gaussian noise, and estimate µ ⋆ by minimizing (5) starting from an initial point θ 0 obtained by the spectral initialization procedure described above. Figure 3 draws the statistical error ∥S (θ 200 -θ ⋆ )∥ ∞ of both preconditioning schemes after 200 iterations -when convergence is reached -as a function of the signal-to-noise ratio (SNR), defined as

SNR = ∥Φ(µ ⋆ )∥ 2 2 / ∥w∥ 2 2 .
The results are benchmarked against the Cramér-Rao bound (CRB) [START_REF] Scharf | Geometry of the Cramér-Rao bound[END_REF]. Both statistical errors remain close to the CRB under a sufficiently large SNR, providing an empirical validation of the robustness of the proposed algorithms.

Conclusion

This work proposed and analyzed preconditioned gradient methods for nonconvex spike deconvolution using both fixed and adaptive preconditioners, and demonstrated that for ground truth with sufficiently separated spikes, the proposed methods achieve a linear rate of convergence that is independent of the number of spikes, as long as a close enough initialization is provided near the ground truth. In particular, by designing the preconditioner to compensate adaptively for the amplitude profile of the spikes, it is possible to accelerate the convergence rate to be dimension-free and independent of the dynamic range, while the convergence using a fixed preconditioner slows down when the dynamic range is large. Our work thus highlights the importance of preconditioning in accelerating convergence in nonconvex spike deconvolution.

As a first step towards understanding the efficacy of first-order methods for spike deconvolution, this works opens up several interesting directions for further investigation.

• Initialization schemes. One immediate direction is to analyze initialization schemes that produce initial estimates that fall into the basin of attraction, which we suspect the procedure described in Section 4 is a good candidate.

• Model order. For simplicity, it is assumed that the model order r is known perfectly, which might not hold in practice. It is of great interest to develop modified algorithms when the model order is overspecified, which has recently been examined comprehensively in [START_REF] Xu | The power of preconditioning in overparameterized low-rank matrix sensing[END_REF] for low-rank estimation from small random initializations.

• General observations. Another direction is to extend the analysis to more general observation operators, possibly including random sampling, missing data, as well as corruptions. This may necessarily require a reformulation of the loss function, such as a nonsmooth and nonconvex formulation using the least absolute deviation [START_REF] Tong | Low-rank matrix recovery with scaled subgradient methods: Fast and robust convergence without the condition number[END_REF] to improve robustness.

• Separation condition. Last but not least, it is of great importance to study to what extent it is possible to relax the success condition in terms of the separation condition, possibly with additional positive constraints of the source amplitudes.

A Summation bounds of the Fejér kernel

The purpose of this section is to present Lemma 5, which delivers fundamental bounds on the absolute sum of the Fejér kernel and its derivatives at sampled points of interest. Although specific to the Fejér kernel, Lemma 5 could be adapted to any other absolutely integrable point spread function without a significant change in the proof structure.

Lemma 5 (Uniform bounds on the Fejér kernel). Suppose that n ≥ 2. Let τ = {τ 1 , . . . , τ r } ⊂ T, and let α > 0 be such that (n + 1)∆(τ ) ≥ α. Let u = {u i,j } i̸ =j ⊂ R be a set of r(r-1) 2 real numbers that are absolutely bounded by β such that (n + 1) max

i̸ =j {|u i,j |} ≜ β < α 2 . ( 52
)
Then the inequalities

max i j̸ =i F (ℓ) N (τ j -τ i + u i,j ) ≤ C ℓ (n + 1) ℓ ((n + 1)∆(τ )) -2 (53)
hold for ℓ = 0, 1, 2, 3, where the constants C ℓ only depend on α and β and are given by

C 0 = 4 π 2 (α -2β) -1 α, (54a)
C 1 = 4 π (α -2β) -1 + 8 π 2 (α -2β) -2 α, (54b)
C 2 = 80 9 (α -2β) -1 + 16 π (α -2β) -2 + 64 3π 2 (α -2β) -3 α, (54c)
C 3 = 16 3 (α -2β) -1 π + 1488 27 (α -2β) -2 + 192 π (α -2β) -3 + 192 π 2 (α -2β) -4 α. ( 54d
)
Proof. First, the function F N is a trigonometric polynomial, hence infinitely differentiable. Let us begin by examining the general expression of F N (t) and its derivatives up to the third order. Assuming t / ∈ Z, by basic calculation, we have that

F N (t) = sin 2 (π(n + 1)t) (n + 1) 2 sin 2 (πt) , (55a)
F ′ N (t) = 2π n + 1 cos((n + 1)πt) sin((n + 1)πt) csc 2 (πt) - 2π (n + 1) 2 sin((n + 1)πt) cot(πt) csc 2 (πt), (55b)
F ′′ N (t) = 2π 2 cos 2 ((n + 1)πt) -sin 2 ((n + 1)πt) csc 2 (πt) - 8π 2 n + 1 cos((n + 1)πt) sin((n + 1)πt) cot(πt) csc 2 (πt) + 2π 2 (n + 1) 2 sin 2 (π(n + 1)t) 2 cot 2 (πt) + 1 csc 2 (πt) (55c)
F ′′′ N (t) = -(n + 1)8π 3 cos ((n + 1)πt)sin ((n + 1)πt) csc 2 (πt) -12π 3 cos 2 ((n + 1)πt) -sin 2 ((n + 1)πt) cot(πt) csc 2 (πt) + 12π 3 n + 1 cos ((n + 1)πt)sin ((n + 1)πt) 3 cot 2 (πt) + 1 csc 2 (πt) - 8π 3 (n + 1) 2 sin 2 ((n + 1)πt) 3 cot 3 (πt) + 2 cot(πt) csc 2 (πt). (55d)
|F N (t)| ≤ 1 (n + 1) 2 csc 2 (πt), ( 56a
)
|F ′ N (t)| ≤ π n + 1 csc 2 (πt) + 2π (n + 1) 2 |cot(πt)| csc 2 (πt), ( 56b
)
|F ′′ N (t)| ≤ 2π 2 csc 2 (πt) + 4π 2 n + 1 |cot(πt)| csc 2 (πt) + 2π 2 (n + 1) 2 2 cot 2 (πt) + 1 csc 2 (πt) (56c)
|F ′′′ N (t)| ≤ (n + 1)4π 3 csc 2 (πt) + 12π 3 |cot(t)| csc 2 (t) + 12π 3 n + 1 3 cot 2 (πt) + 1 csc 2 (πt) + 8π 3 (n + 1) 2 3 cot 3 (πt) + 2 |cot(πt)| csc 2 (πt). (56d)
To continue, observe that a basic building block in the bound (56) takes the following function form, which is denoted by

h ℓ (t) = cot ℓ (πt) csc 2 (πt), ℓ = 0, 1, 2. (57)
The functions h ℓ , ℓ = 0, 1, 2 are even, 1-periodic, and continuous, non-negative, decreasing and convex over the interval (0, 1 2 ]. Define the three sums S ℓ (n, τ , u, i), ℓ = 0, 1, 2 by

S ℓ (n, τ , u, i) = j̸ =i h ℓ (τ j -τ i + u i,j ) . (58)
It is straightforward to see that the terms of interest can be bounded in terms of S ℓ (n, τ , u, i) in view of (56). We shall claim the following bound of S ℓ (n, τ , u, i) holds, which will be proven at the end of proof:

S ℓ (n, τ , u, i) ≤ 2 π∆(τ ) ℓ+2 1 (ℓ + 1) (1 -2α -1 β) ℓ+1 , ℓ = 0, 1, 2. (59)
Interestingly, the above bound on S ℓ (n, τ , u, i) does not depend on n or i. We are now ready to establish the bound (53) of interest. First, we have

max i j̸ =i |F N (τ j -τ i + u i,j )| ≤ j̸ =i 1 (n + 1) 2 h 0 (τ j -τ i + u i,j ) = 1 (n + 1) 2 S 0 (n, τ , u, i) ≤ 4 π 2 ∆(τ ) -2 (n + 1) 2 (1 -2α -1 β) = 4 π 2 (α -2β) -1 α ((n + 1)∆(τ )) -2 . ( 60
)
Similarly, for the case of the first derivative, we have

max i j̸ =i |F ′ N (τ j -τ i + u i,j )| ≤ j̸ =i π n + 1 h 0 (τ j -τ i + u i,j ) + 2π (n + 1) 2 h 1 (τ j -τ i + u i,j ) = π n + 1 S 0 (n, τ , u, i) + 2π (n + 1) 2 S 1 (n, τ , u, i) ≤ 4∆(τ ) -2 π(n + 1) (1 -2α -1 β) + 8∆(τ ) -3 π 2 (n + 1) 2 (1 -2α -1 β) 2 ≤ 4 π (α -2β) -1 + 8 π 2 (α -2β) -2 α(n + 1) ((n + 1)∆(τ )) -2 , (61)
where the last line uses (n + 1)∆(τ ) ≥ α. Moving onto the second derivative, it follows

max i j̸ =i |F ′′ N (τ j -τ i + u i,j )| ≤ j̸ =i 2π 2 h 0 (τ j -τ i + u i,j ) + 4π 2 n + 1 h 1 (τ j -τ i + u i,j ) + 2π 2 (n + 1) 2 (2h 2 (τ j -τ i + u i,j ) + h 0 (τ j -τ i + u i,j )) = 2π 2 S 0 (n, τ , u, i) + 4π 2 n + 1 S 1 (n, τ , u, i) + 2π 2 (n + 1) 2 (2S 2 (n, τ , u, i) + S 0 (n, τ , u, i)) ≤ 2π 2 1 + 1 (n + 1) 2 S 0 (n, τ , u, i) + 4π 2 n + 1 S 1 (n, τ , u, i) + 4π 2 (n + 1) 2 S 2 (n, τ , u, i) ≤ 20π 2 9 S 0 (n, τ , u, i) + 4π 2 n + 1 S 1 (n, τ , u, i) + 4π 2 (n + 1) 2 S 2 (n, τ , u, i) ≤ 80∆(τ ) -2 9 (1 -2α -1 β) + 16∆(τ ) -3 π(n + 1) (1 -2α -1 β) 2 + 64∆(τ ) -4 3π 2 (n + 1) 2 (1 -2α -1 β) 3 ≤ 80 9 (α -2β) -1 + 16 π (α -2β) -2 + 64 3π 2 (α -2β) -3 α(n + 1) 2 ((n + 1)∆(τ )) -2 . ( 62
)
Finally, for the third derivative, it holds

max i j̸ =i |F ′′′ N (τ j -τ i + u i,j )| ≤ j̸ =i (n + 1)4π 3 h 0 (τ j -τ i + u i,j ) + 12π 3 h 1 (τ j -τ i + u i,j ) + 12π 3 n + 1 (3h 2 (τ j -τ i + u i,j ) + h 0 (τ j -τ i + u i,j )) + 8π 3 (n + 1) 2 (3h 3 (τ j -τ i + u i,j ) + 2h 1 (τ j -τ i + u i,j )) = (n + 1)4π 3 S 0 (n, τ , u, i) + 12π 3 S 1 (n, τ , u, i) + 12π 3 n + 1 (3S 2 (n, τ , u, i) + S 0 (n, τ , u, i)) + 8π 3 (n + 1) 2 (3S 3 (n, τ , u, i) + 2S 1 (n, τ , u, i)) ≤ (n + 1)4π 3 1 + 3 (n + 1) 2 S 0 (n, τ , u, i) + 12π 3 1 + 4 3 (n + 1) 2 S 1 (n, τ , u, i) + 36π 3 n + 1 S 2 (n, τ , u, i) + 24π 3 (n + 1) 2 S 3 (n, τ , u, i) ≤ (n + 1) 64 3 π ∆(τ ) -2 1 -2α -1 β + 2488 27 ∆(τ ) -3 (1 -2α -1 β) 2 + 192 π(n + 1) ∆(τ ) -4 (1 -2α -1 β) 3 + 192 π 2 (n + 1) 2 ∆(τ ) -5 (1 -2α -1 β) 4 ≤ 16 3 (α -2β) -1 π + 1488 27 (α -2β) -2 + 192 π (α -2β) -3 + 192 π 2 (α -2β) -4 α(n + 1) 3 ((n + 1)∆(τ )) -2 . ( 63
)
The proof is thus completed if we can prove (59), which is the focus of the rest of the proof.

Proof of (59) We fix the index i and take the convention u i,i = 0. As the functions {h ℓ } are 1-periodic, the quantity S ℓ (n, τ , β, i) is invariant by integer translations of {τ j }'s. Therefore, one can make the assumption, up to a modulo considerations and a reordering of the indices that the sequence {τ j -τ i + u i,j } j is within the range [-1 2 , 1 2 ) and in an ascending order, so that

- 1 2 ≤ τ 1 -τ i + u i,1 < τ 2 -τ i + u i,2 < • • • < τ r -τ i + u i,r < 1 2 .
Denote by r + and r -the number of positive and negative elements in the set {τ j -τ i + u i,j } j , respectively. As τ j -τ i + u i,j = 0 if and only if i = j, we have that r + + r -= r -1. Using the separation condition, and as h ℓ is decreasing over (0, 1 2 ] and even, we have that

0 ≤ h ℓ (τ i+j -τ i + u i,j ) ≤ h ℓ j∆(τ ) - β n + 1 , j = 1, . . . , r + ; 0 ≤ h ℓ (τ i-j -τ i + u i,j ) ≤ h ℓ -j∆(τ ) + β n + 1 = h ℓ j∆(τ ) - β n + 1 , j = 1, . . . , r -.
We can subsequently bound the sum (58) as

S ℓ (n, τ , u, i) = r+ j=1 h ℓ τ i+j -τ i + u i,(i+j) + r- j=1 h ℓ τ i-j -τ i + u i,(i-j) ≤ r+ j=1 h ℓ j∆(τ ) - β n + 1 + r- j=1 h ℓ j∆(τ ) - β n + 1 ≤ 2 ⌈ r-1 2 ⌉ j=1 h ℓ j∆(τ ) - β n + 1 . ( 64
)
Identifying and associating the right-hand side of (64) to a Riemann sum with a mid-point rule and recalling that h ℓ is decreasing and convex over (0, 1 2 ] leads to the majorant

S ℓ (n, τ , u, i) ≤ 2 (⌈ r-1 2 ⌉+ 1 2 )∆(τ)-β n+1 ∆(τ ) 2 -β n+1 h ℓ (u)du ≤ 2 1 2 ∆(τ ) 2 -β n+1 h ℓ (u)du = 2 H ℓ 1 2 -H ℓ ∆(τ ) 2 - β n + 1 , (65)
where we used the non-negativity of h ℓ and the inequality r-1 2 + 1 2 ∆(τ ) -β n+1 ≤ 1 2 in the second inequality. Here, {H ℓ }'s are primitives of the functions {h ℓ } over the interval (0, 1 2 ]. Moreover, we have

H ℓ (t) = - 1 π (ℓ + 1) cot ℓ+1 (πt) + c (66)
for ℓ = 0, 1, 2, where c is an arbitrary constant. This yields, with the inequality 0 ≤ cot(πu) ≤ 1 πu for all 0 < u < 1 2 , a further simplification of (65):

∆(τ )S ℓ (n, τ , u, i) ≤ 2 π(ℓ + 1) cot ℓ+1 π ∆(τ ) 2 - β n + 1 ≤ 2 π ℓ+2 ∆(τ ) -ℓ-1 1 1 -2β∆(τ ) -1 n+1 ℓ+1 ≤ 2 π ℓ+2 ∆(τ ) -ℓ-1 1 (1 -2α -1 β) ℓ+1 , (67)
which immediately leads to the claimed bound (59).

B Proof of the uniform Hessian bounds

This section is dedicated to establish Theorem 3 and Theorem 4.

B.1 Technical lemmas

Lemma 6 and Lemma 7 provide bounds on core quantities that are involved in the decomposition of the Hessian matrix H(θ) characterized in Section 3.1.2. Their proofs are presented in Appendix B.4 and Appendix B.5, respectively.

Lemma 6. Suppose that n ≥ 2 and let τ ⊂ T be such that (n + 1)∆(τ ) ≥ α for some α > 0, then there exists a constant K ∆ with

K ∆ := max C 0 + 3 √ 3 4π C 1 , 3 √ 3 4π C 1 + 27 16π 2 C 2 (68)
where the constants C 0 , C 1 , C 2 > 0 are defined in in (54) with parameters α and β = 0 such that

∥D(τ ) -I∥ ∞ ≤ K ∆ ((n + 1)∆(τ )) -2 . ( 69
)
Lemma 7. Suppose that n ≥ 2 and let τ ⋆ = [τ ⋆ 1 , . . . , τ ⋆ r ] ⊤ and τ = [τ 1 , . . . , τ r ] ⊤ be two vectors of points around the torus. Assume that (n + 1)∆(τ ⋆ ) ≥ α. As long as

(n + 1) ∥τ -τ ⋆ ∥ ∞ ≤ β < α 2 , we have Φ(δ ′ τj ), Φ (µ(θ)) -µ(θ ⋆ )) ≤ (C 1 ∥a -a ⋆ ∥ ∞ + C 2 ∥a ⋆ ∥ ∞ (n + 1) ∥τ -τ ⋆ ∥ ∞ ) (n + 1) ((n + 1)∆(τ )) -2 , (70a)
Φ(δ ′′ τj ), Φ (µ(θ)) -µ(θ ⋆ )) ≤ (C 2 ∥a -a ⋆ ∥ ∞ + C 3 ∥a ⋆ ∥ ∞ (n + 1) ∥τ -τ ⋆ ∥ ∞ ) (n + 1) 2 ((n + 1)∆(τ )) -2 (70b)
for all j = 1, . . . , r, where the constants C 1 , C 2 , C 3 > 0 are defined in (54) with parameters (α, β).

B.2 Proof of Theorem 3

Recalling the expression of the Hessian in [START_REF] Traonmilin | The basins of attraction of the global minimizers of the non-convex sparse spike estimation problem[END_REF], it follows that

SP H(θ)S -1 -I = SP G(θ)S -1 -I + SP E(θ)S -1 . (71)
We proceed to bound SP G(θ)S -1 -I ∞ and SP E(θ)S -1 ∞ separately, and then combine them via the triangle inequality.

Step 1: bound SP G(θ)S -1 -I ∞ From ( 16) and [START_REF] Traonmilin | Projected gradient descent for non-convex sparse spike estimation[END_REF], we have that

SP G(θ)S -1 -I = SP diag 1 r -F ′′ N (0)a H D(τ ) diag 1 r -F ′′ N (0)a S -1 -I = diag a ⋆-1 A -2 a H D(τ ) diag a ⋆ a -I = diag a ⋆-1 A -2 a H (D(τ ) -I) diag a ⋆ a + diag 1 r A -2 |a| 2 -I. ( 72
)
This immediately yields from the triangle inequality that

SP G(θ)S -1 -I ∞ ≤ diag a ⋆-1 A -2 ⊙ a H (D(τ ) -I) diag a ⋆ a ∞ + diag 1 r A -2 ⊙ |a| 2 -I ∞ ≤ diag a ⋆-1 A -2 ⊙ a H ∞ ∥D(τ ) -I∥ ∞ diag a ⋆ a ∞ + diag 1 r A -2 ⊙ |a| 2 -I ∞ ≤ max j 1 a ⋆ j , |a j | A 2 max j a ⋆ j , |a j | ∥D(τ ) -I∥ ∞ + max j 1 - |a j | 2 A 2 . ( 73
)
The three maxima in (73) can be controlled using the basic relation

1 -∥S(θ k -θ ⋆ )∥ ∞ ≤ a ⋆ j -a k,j -a ⋆ j a ⋆ j ≤ |a j | a ⋆ j ≤ a ⋆ j + a k,j -a ⋆ j a ⋆ j ≤ 1 + ∥S(θ k -θ ⋆ )∥ ∞ , (74)
for j = 1, . . . , r and by exploiting the assumption ∥a ⋆ ∥ ∞ ≤ A as follows

max j 1 a ⋆ j , |a j | A 2 ≤ max j 1 a ⋆ j , a ⋆ j A 2 (1 + ∥S(θ k -θ ⋆ )∥ ∞ ) ≤ 1 a ⋆ min (1 + ∥S(θ k -θ ⋆ )∥ ∞ ) , (75a)
max j a ⋆ j , |a j | ≤ max j a ⋆ j , a ⋆ j (1 + ∥S(θ k -θ ⋆ )∥ ∞ ) ≤ ∥a ⋆ ∥ ∞ (1 + ∥S(θ k -θ ⋆ )∥ ∞ ) , (75b)
max j 1 - |a j | 2 A 2 ≤ 1 -min j |a j | 2 A 2 ≤ 1 - (a ⋆ min ) 2 A 2 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) . (75c)
The bounds (75) and Lemma 6 imply with (73) that

SP G(θ)S -1 -I ∞ ≤ ∥a ⋆ ∥ ∞ a ⋆ min (1 + ∥S (θ k -θ ⋆ )∥ ∞ ) 2 ∥D(τ ) -I∥ ∞ + 1 - (a ⋆ min ) 2 (1 -∥S (θ k -θ ⋆ )∥ ∞ ) 2 A 2 ≤ ∥a ⋆ ∥ ∞ a ⋆ min (1 + ∥S (θ k -θ ⋆ )∥ ∞ ) 2 K ∆ ((n + 1)∆(τ )) -2 + 1 - (a ⋆ min ) 2 (1 -∥S (θ k -θ ⋆ )∥ ∞ ) 2 A 2 , ( 76
)
where the constant K ∆ defined in Lemma 6 is numerically evaluated to K ∆ ≤ 2.32 by setting the parameter α = 16.5.

Step 2: bound SP E(θ)S -1 ∞ Using the block diagonal structure of the matrix E(θ) defined in [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF], we have that

SP E(θ)S -1 =   0 1 √ -F ′′ N (0) diag a ⋆-1 E 1 (θ) 1 √ -F ′′ N (0) diag A -2 a ⋆ H E 1 (θ) -1 F ′′ N (0) A -2 E 2 (θ)   . (77)
Therefore, it follows

SP E(θ)S -1 ∞ ≤ 1 -F ′′ N (0) max j max 1 a ⋆ j , a ⋆ j A 2 δ ′ τj , F N * (µ(θ) -µ(θ ⋆ )) - 1 F ′′ N (0) max j |a j | A 2 δ ′′ τj , F N * (µ(θ) -µ(θ ⋆ )) ≤ (a ⋆ min ) -1 (1 + ∥S(θ k -θ ⋆ )∥ ∞ ) 1 -F ′′ N (0) max j δ ′ τj , F N * (µ(θ) -µ (θ ⋆ )) - 1 F ′′ N (0) max j δ ′′ τj , F N * (µ(θ) -µ(θ ⋆ )) , (78)
where we used the bounds (75) in the second line. From (78), it can be seen that bounding the quantity of interest amounts to controlling δ ′ τj , F N * (µ(θ) -µ(θ ⋆ )) and δ ′′ τj , F N * (µ(θ) -µ(θ ⋆ )) , which can be achieved by applying Lemma 7. Substituting the expression provided by Lemma 7 into (78) leads to

SP E(θ)S -1 ∞ ≤ C 1 (n + 1) -F ′′ N (0) + C 2 (n + 1) 2 -F ′′ N (0) ∥a -a ⋆ ∥ ∞ ∥a ⋆ ∥ ∞ + C 2 (n + 1) 2 -F ′′ N (0) + C 3 (n + 1) 3 (-F ′′ N (0)) 3/2 -F ′′ N (0) ∥τ k -τ ⋆ ∥ ∞ • ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ )) -2 (1 + ∥S(θ k -θ ⋆ )∥ ∞ ) . ( 79
E(θ)S -1 ∞ ≤ K a ∥a -a ⋆ ∥ ∞ ∥a ⋆ ∥ ∞ + K τ -F ′′ N (0) ∥τ k -τ ⋆ ∥ ∞ • ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ )) -2 (1 + ∥S(θ k -θ ⋆ )∥ ∞ ) , (80)
where the constants K a and K τ are given by

K a = C 1 27 16π 2 + C 2 27 16π 2 ≤ 4.40, (81a)
K τ = C 2 27 16π 2 + C 3 27 16π 2 3/2 ≤ 6.71. ( 81b
)
Step 3: combine the bounds The scaled Hessian matrix SP H(θ)S -1 can be controlled over S k by the triangle inequality as follows

SP H(θ)S -1 -I ∞ ≤ SP (G(θ) + E(θ)) S -1 -I ∞ ≤ SP G(θ)S -1 -I ∞ + SP E(θ)S -1 ∞ . (82)
Furthermore, we note that the inequality |∆(τ ) -∆(τ ⋆ )| ≤ 2 ∥τ -τ ⋆ ∥ ∞ holds for every τ , τ ⋆ ⊂ R. This implies, with the assumption ∥τ -τ ⋆ ∥ ∞ ≤ 1 4 ∆(τ ⋆ ), that ∆(τ ) ≥ 1 2 ∆(τ ⋆ ). Substituting the bounds given in (76) and (80) yields

SP H(θ)S -1 -I ∞ ≤ 1 - a ⋆ min 2 A 2 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 + K ∆ + K a ∥a -a ⋆ ∥ ∞ ∥a ⋆ ∥ ∞ + K τ -F ′′ N (0) ∥τ k -τ ⋆ ∥ ∞ ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ )) -2 (1 + ∥S (θ k -θ ⋆ )∥ ∞ ) 2 ≤ 1 - a ⋆ min 2 A 2 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 + (K ∆ + (K a + K τ ) ∥S (θ k -θ ⋆ )∥ ∞ ) ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ )) -2 (1 + ∥S (θ k -θ ⋆ )∥ ∞ ) 2 ≤ 1 - a ⋆ min 2 A 2 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 + (K ∆ + (K a + K τ ) ∥S (θ k -θ ⋆ )∥ ∞ ) ∥a ⋆ ∥ ∞ a ⋆ min 4 ((n + 1)∆(τ ⋆ )) -2 (1 + ∥S (θ k -θ ⋆ )∥ ∞ ) 2 ≤ 1 - a ⋆ min 2 A 2 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 + (4K ∆ + K θ ∥S (θ k -θ ⋆ )∥ ∞ ) ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ ⋆ )) -2 (1 + ∥S (θ k -θ ⋆ )∥ ∞ ) 2 , (83)
where we defined in the last line the constant K θ as

K θ = 4 (K a + K τ ) ≤ 44.42. (84)
This concludes the proof of the theorem.

B.3 Proof of Theorem 4

We proceed analogously to the proof of Theorem 3 presented in Appendix B.2. We start from the expansion [START_REF] Traonmilin | The basins of attraction of the global minimizers of the non-convex sparse spike estimation problem[END_REF] of the Hessian matrix to get

SP k H(θ)S -1 -I = SP k G(θ)S -1 -I + SP k E(θ)S -1 , (85)
and proceed to bound SP k G(θ)S -1 -I ∞ and SP k E(θ)S -1 ∞ individually before recombining them via the triangle inequality.

Step 1: bound SP k G(θ)S -1 -I ∞ From ( 16) and ( 33), we have that

SP k G(θ)S -1 -I = SP k diag 1 r -F ′′ N (0)a H D(τ ) diag 1 r -F ′′ N (0)a H S -1 -I = diag a ⋆-1 |a k | -2 ⊙ a H D(τ ) diag a ⋆ a -I = diag a ⋆-1 |a k | -2 ⊙ a H (D(τ ) -I) diag a ⋆ a + diag 1 r |a k | -2 ⊙ |a| 2 -I. ( 86
)
This immediately yields from the triangle inequality

SP k G(θ)S -1 -I ∞ ≤ diag a ⋆-1 |a k | -2 ⊙ a H (D(τ ) -I) diag a ⋆ a ∞ + diag 1 r |a k | -2 ⊙ |a| 2 -I ∞ ≤ diag a ⋆-1 |a k | -2 ⊙ a H ∞ ∥D(τ ) -I∥ ∞ diag a ⋆ a ∞ + diag 1 r |a k | -2 ⊙ |a| 2 -I ∞ ≤ max j 1 a ⋆ j , a ⋆ j |a k,j | 2 max j a ⋆ j , |a j | ∥D(τ ) -I∥ ∞ + max j a ⋆ j 2 |a k,j | 2 -1 (87)
The three maxima in (87) can be controlled using the basic relation

a ⋆ j |a k,j | ≤ a ⋆ j a ⋆ j -a k,j -a ⋆ j ≤ 1 1 -∥S(θ k -θ ⋆ )∥ ∞ , (88)
as follows

max j 1 a ⋆ j , a ⋆ j |a k,j | 2 ≤ max j 1 a ⋆ j , 1 a ⋆ j 1 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 ≤ 1 a ⋆ min 1 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 (89a) max j a ⋆ j , |a j | ≤ max j a ⋆ j , a ⋆ j (1 + ∥S(θ k -θ ⋆ )∥ ∞ ) ≤ ∥a ⋆ ∥ ∞ (1 + ∥S(θ k -θ ⋆ )∥ ∞ ) (89b) max j a ⋆ j 2 |a k,j | 2 -1 ≤ max j 1 - a ⋆ j 2 |a k,j | 2 , a ⋆ j 2 |a k,j | 2 -1 ≤ max j 1, a ⋆ j 2 |a k,j | 2 -1 ≤ 1 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 -1. (89c)
The bounds (89) and Lemma 6 conclude with (87) on

SP k G(θ)S -1 -I ∞ ≤ ∥a ⋆ ∥ ∞ a ⋆ min 1 + ∥S (θ k -θ ⋆ )∥ ∞ (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 ∥D(τ ) -I∥ ∞ + 1 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 -1 ≤ ∥a ⋆ ∥ ∞ a ⋆ min 1 + ∥S (θ k -θ ⋆ )∥ ∞ (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 K ∆ ((n + 1)∆(τ )) -2 + 1 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 -1, (90)
where the constant K ∆ defined in Lemma 6 is numerically evaluated to K ∆ ≤ 2.32 by setting the parameter α = 4.7.

Step 2: bound SP k E(θ)S -1 ∞ Again, using the block diagonal structure of the matrix E(θ) defined in [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF], we similarly have that

SP k E(θ)S -1 =   0 1 √ -F ′′ N (0) diag a ⋆-1 E 1 (θ) 1 √ -F ′′ N (0) diag |a k | -2 ⊙ a ⋆ H E 1 (θ) -1 F ′′ N (0) |a k | -2 E 2 (θ)   . (91)
Therefore, the quantity of interest can be bounded as follows

SP k E(θ)S -1 ∞ ≤ 1 -F ′′ N (0) max j max 1 a ⋆ j , a ⋆ j |a j,k | 2 δ ′ τj , F N * (µ(θ) -µ(θ ⋆ )) - 1 F ′′ N (0) max j |a j | |a j,k | 2 δ ′′ τj , F N * (µ(θ) -µ(θ ⋆ )) ≤ (a ⋆ min ) -1 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 1 -F ′′ N (0) max j δ ′ τj , F N * (µ(θ) -µ (θ ⋆ )) - 1 F ′′ N (0) max j δ ′′ τj , F N * (µ(θ) -µ(θ ⋆ )) , (92)
where we used the inequalities (89) on the second line. Substituting the expression provided by Lemma 7 into (92) leads to

SP k E(θ)S -1 ∞ ≤ C 1 n + 1 -F ′′ N (0) + C 2 (n + 1) 2 -F ′′ N (0) ∥a -a ⋆ ∥ ∞ ∥a ⋆ ∥ ∞ + C 2 (n + 1)
SP k E(θ)S -1 ∞ ≤ K a ∥a -a ⋆ ∥ ∞ ∥a ⋆ ∥ ∞ + K τ -F ′′ N (0) ∥τ k -τ ⋆ ∥ ∞ ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ )) -2 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 , ( 94
)
where the constants K a and K τ are given by

K a = C 1 27 16π 2 + C 2 27 16π 2 ≤ 5.26, (95a)
K τ = C 2 27 16π 2 + C 3
Finally, the inequality |∆(τ ) -∆(τ ⋆ )| ≤ 2 ∥τ -τ ⋆ ∥ ∞ holds for every τ , τ ⋆ ⊂ R. This implies, with the assumption ∥τ -τ ⋆ ∥ ∞ ≤ 1 4 ∆(τ ⋆ ), that ∆(τ ) ≥ 1 2 ∆(τ ⋆ ). Substituting the bounds given in (90) and (94) yields

SP k H(θ)S -1 -I ∞ ≤ 1 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 -1 + K ∆ (1 + ∥S (θ k -θ ⋆ )∥ ∞ ) + K a ∥a -a ⋆ ∥ ∞ ∥a ⋆ ∥ ∞ + K τ -F ′′ N (0) ∥τ k -τ ⋆ ∥ ∞ ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ )) -2 (1 -∥S (θ k -θ ⋆ )∥ ∞ ) 2 ≤ 1 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 -1 + (K ∆ (1 + ∥S (θ k -θ ⋆ )∥ ∞ ) + (K a + K τ ) ∥S (θ k -θ ⋆ )∥ ∞ ) ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ )) -2 (1 -∥S (θ k -θ ⋆ )∥ ∞ ) 2 ≤ 1 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 -1 + (K ∆ + (K ∆ + K a + K τ ) ∥S (θ k -θ ⋆ )∥ ∞ ) ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ )) -2 (1 -∥S (θ k -θ ⋆ )∥ ∞ ) 2 ≤ 1 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 -1 + (K ∆ + (K ∆ + K a + K τ ) ∥S (θ k -θ ⋆ )∥ ∞ ) ∥a ⋆ ∥ ∞ a ⋆ min 4 ((n + 1)∆(τ ⋆ )) -2 (1 -∥S (θ k -θ ⋆ )∥ ∞ ) 2 ≤ 1 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 -1 + (4K ∆ + K θ ∥S (θ k -θ ⋆ )∥ ∞ ) ∥a ⋆ ∥ ∞ a ⋆ min ((n + 1)∆(τ ⋆ )) -2 (1 -∥S(θ k -θ ⋆ )∥ ∞ ) 2 , (97)
where defined the constant K θ as

K θ = 4 (K ∆ + K a + K τ ) ≤ 75.80. (98)
This concludes the proof of the theorem.

B.4 Proof of Lemma 6

Leveraging the block structure (34) of the matrix D(τ ), it boils down to controlling

∥D(τ ) -I∥ ∞ ≤ max {∥D 0 (τ ) -I∥ ∞ + ∥D 1 (τ )∥ ∞ , ∥D 2 (τ ) -I∥ ∞ + ∥D 1 (τ )∥ ∞ } , (99)
which can be accomplished with the aid of Lemma 5. Specifically, recalling the expressions in [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF], applying Lemma 5 with parameters α and β = 0, and noticing that the inequality (n+1) 2 -F ′′ N (0) ≤ 27 16π 2 holds whenever n ≥ 2, the quantities of interest in (99) can be controlled as follows

∥D 0 (τ ) -I∥ ∞ = max i j̸ =i
|F N (τ j -τ i )| ≤ C 0 ((n + 1)∆(τ )) -2 , (100a)

∥D 1 (τ )∥ ∞ = 1 -F ′′ N (0) max i j |F ′ N (τ j -τ i )| ≤ C 1 n + 1 -F ′′ N (0) ((n + 1)∆(τ )) -2 ≤ 3 √ 3 4π C 1 ((n + 1)∆(τ )) -2 , (100b)

B.5 Proof of Lemma 7

For the first inequality, following the definition, we have that

Φ(δ ′ τj ), Φ (µ(θ) -µ(θ ⋆ )) = δ ′ τj , Φ * (Φ (µ(θ) -µ(θ ⋆ ))) = r ℓ=1 a ℓ F ′ N (τ j -τ ℓ ) - r ⋆ ℓ=1 a ⋆ ℓ F ′ N (τ j -τ ⋆ ℓ ) = r ℓ=1 (a ℓ -a ⋆ ℓ )F ′ N (τ j -τ ℓ ) + r ℓ=1 a ⋆ ℓ (F ′ N (τ j -τ ℓ ) -F ′ N (τ j -τ ⋆ ℓ )) ≤ r ℓ=1 (a ℓ -a ⋆ ℓ )F ′ N (τ j -τ ℓ ) + r ℓ=1 a ⋆ ℓ (F ′ N (τ j -τ ℓ ) -F ′ N (τ j -τ ⋆ ℓ )) , (102)
where the last line used the triangle inequality. To proceed, we control the two terms separately. Using Hölder's inequality and Lemma 5, we obtain Next, the second term can be bounded by applying Hölder's inequality, the mean-value theorem and Lemma 5 with parameter β ≥ (n + 1) ∥τ -τ ⋆ ∥ ∞ as follows

r ℓ=1 a ⋆ ℓ (F ′ N (τ j -τ ℓ ) -F ′ N (τ j -τ ⋆ ℓ )) ≤ ∥a ⋆ ∥ ∞ r ℓ=1 |F ′ N (τ j -τ ℓ ) -F ′ N (τ j -τ ⋆ ℓ )| ≤ ∥a ⋆ ∥ ∞ ∥τ -τ ⋆ ∥ ∞ r ℓ=1 sup |u ℓ |≤∥τ -τ ⋆ ∥ ∞ |F ′′ N (τ j -τ ℓ + u ℓ )| ≤ ∥a ⋆ ∥ ∞ ∥τ -τ ⋆ ∥ ∞ C 2 (n + 1) 2 ((n + 1)∆(τ )) -2 . ( 104
)
Plugging the previous two bounds into the inequality (102) reduces to Φ(δ ′ τj ), Φ (µ(θ) -µ(θ ⋆ ))

≤ (C 1 ∥a -a ⋆ ∥ ∞ + C 2 ∥a ⋆ ∥ ∞ (n + 1) ∥τ -τ ⋆ ∥ ∞ ) (n + 1) ((n + 1)∆(τ )) -2 . ( 105
)
Moreover, we may show that through analogous reasoning that Φ(δ ′′ τj ), Φ (µ(θ) -µ(θ ⋆ ))

≤ (C 2 ∥a -a ⋆ ∥ ∞ + C 3 ∥a ⋆ ∥ ∞ (n + 1) ∥τ -τ ⋆ ∥ ∞ ) (n + 1) 2 ((n + 1)∆(τ )) -2 , (106) which concludes the proof.

1 1 . 2 1 . 3 1 . 4 5 2 2 . 1 2 . 2 2 . 3

 11213145212223 Observation model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Our contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Notation and paper organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . How does preconditioning help local convergence? 5 Preconditioned gradient descent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Invariant preconditioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Adaptive preconditioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

= 6 Figure 1 :

 61 Figure 1: Success rate of the invariant and adaptive preconditioning schemes on reconstructing the ground truth θ ⋆ as a function of the initialization distance ∥S (θ 0 -θ ⋆ )∥ ∞ for different dynamic ranges κ. The results are averaged over 1000 randomized trials.

Figure 2 :

 2 Figure 2: Convergence rates of the iterate sequence of preconditioned GD towards the ground truth as the dynamic range κ varies for: (a) the invariant preconditioning scheme; (b) the adaptive preconditioning scheme.

Figure 3 :

 3 Figure 3: Statistical error ∥S (θ 200 -θ ⋆ )∥ ∞ of both preconditioned GD schemes as a function of the SNR. The dynamic range is set to κ = 3, and the results are averaged over 1000 randomized trials.

  Using the four trigonometric bounds |sin(a)| ≤ 1, |cos(a)| ≤ 1, |cos(a) sin(a)| ≤ 1 2 , and cos 2 (a) -sin 2 (a) ≤ 1 for a ∈ R, and the triangle inequality, (55) can be further absolutely bounded as

2 ( 1 -

 21 ∥S(θ k -θ ⋆ )∥ ∞ ) 2 . (93) Evaluating the constants C 1 ≤ 3.24, C 2 ≤ 22.90 and C 3 ≤ 105.55 defined in (54) with parameters α = 4.7 and β = α 4 = 1.2, altogether with the inequality (n+1) 2 -F ′′ N (0) ≤ 27 16π 2 under the assumption n ≥ 2 yields

Step 3 :

 3 combine the bounds The scaled Hessian matrix SP k H(θ)S -1 can be controlled over S k by the triangle inequality as followsSP k H(θ)S -1 -I ∞ ≤ SP k (G(θ) + E(θ)) S -1 -I ∞ ≤ SP k G(θ)S -1 -I ∞ + SP k E(θ)S -1∞ .

∥D 2 ( 2 (

 22 j -τ i )| ≤ C substituting (100) into (99) leads to ∥D(τ ) -I∥ ∞ ≤ max C 0 the desired statement.

  -a ⋆ ℓ )F ′ N (τ j -τ ℓ ) ≤ ∥a -a ⋆ ∥ ∞ r ℓ=1 |F ′ N (τ j -τ ℓ )| ≤ C 1 (n + 1) ∥a -a ⋆ ∥ ∞ ((n + 1)∆(τ ))

)

  Evaluating the constants C 1 ≤ 2.75, C 2 ≤ 19.08 and C 3 ≤ 48.74 defined in (54) with parameters α = 16.5 and β = α 4 = 4.125, altogether with the inequality (n+1) 2

	-F ′′ N (0) ≤ 27 16π 2 under the assumption n ≥ 2 yields
	SP

The normalization B = 1 is made without loss of generality, up to a rescaling of the source locations τ ⋆ ℓ . Assuming that τ ⋆ ℓ ∈ [-T

,T 2 ], where T is the length of the observation window, the deconvolution problem only depends on the time-bandwidth product T • B[START_REF] Chi | Guaranteed Blind Sparse Spikes Deconvolution via Lifting and Convex Optimization[END_REF].

† YThe work of Y. Chi is supported in part by Office of Naval Research under N00014-19-1-2404, and by National Science Foundation under CAREER ECCS-1818571 and ECCS-2126634.