Quentin Denoyelle
email: denoyelle@ceremade.dauphine.fr

Vincent Duval
email: vincent.duval@inria.fr

Gabriel Peyré
email: gabriel.peyre@ens.fr

Emmanuel Soubies
email: emmanuel.soubies@epfl.ch

The Sliding Frank-Wolfe Algorithm and its Application to Super-Resolution Microscopy

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

1.1. Super-Resolution using the BLASSO Super-resolution consists in retrieving the fine scale details of a possibly noisy signal from coarse scale information. The importance of recovering the high frequencies of a signal comes from the fact that there is often a physical blur in the acquisition process, such as diffraction in optical systems, wave reflection in seismic imaging or spikes recording from neuronal activity.

In resolution theory [START_REF] Aj | Resolution: a survey[END_REF], the two-point resolution criterion defines the ability of a system to resolve two points of equal intensities. It is defined as a distance, namely the Rayleigh criterion, which only depends on the system. In the case of the ideal low-pass filter (i.e. convolution with the Dirichlet kernel) with cutoff frequency f c , the Rayleigh criterion is 1/f c . Then, super-resolution in signal processing consists in developing techniques which enable to retrieve information below the Rayleigh criterion.

Let us introduce in a more formal way the problem which will be the framework of this article. Let X be a connected subset of R d with non-empty interior or the d-dimensional torus T d (d ∈ N * ) and M(X) the Banach space of bounded Radon measures on X. The latter can be seen as the topological dual of C 0 (X, R), the space of continuous functions on X that vanish at infinity. We consider a given integral operator Φ : M(X) → H, where H is a separable Hilbert space, whose kernel ϕ is supposed to be a smooth function (see Definition 4 for the technical assumptions made on ϕ), i.e. ∀m ∈ M(X), Φm def. = X ϕ(x)dm(x).

(

) 1
The operator Φ models the acquisition process. It includes translation-invariant operators such as convolutions (i.e. ϕ(x) = φ(•-x)) as well as non-translation invariant operators such as the Laplace transform (X = R * + and ϕ(x) = (t → e -tx ) ∈ L 2 (R + )) considered in the present paper.

The sparse spikes super-resolution problem aims at recovering an approximation of an unknown input discrete measure m a0,x0 def.

= N i=1 a 0,i δ x0,i from noisy measurements y def.

= y 0 + w where y 0 def.

= Φm a0,x0 and w ∈ H models the acquisition noise. Here a 0,i ∈ R are the amplitudes of the Dirac masses at positions x 0,i ∈ X. This is an ill-posed inverse problem and the BLASSO is a way to solve it in a stable way by introducing a sparsity-enforcing convex regularization.

1.1.1. From the LASSO to the BLASSO The common practice in sparse spike recovery relies on 1 regularization which is known as LASSO in statistic [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] or basis pursuit in the signal processing community [START_REF] Scott Shaobing Chen | Atomic decomposition by basis pursuit[END_REF]. Given a grid of possible positions, the reconstruction problem is addressed as the minimization of a quadratic error subject to an 1 penalization. The 1 prior provides solutions with few nonzero coefficients and can be computed efficiently with convex optimization methods. Moreover, recovery guarantees have been proved under certain assumptions [START_REF] David | Super-resolution via sparsity constraints[END_REF].

Following recent works (see for instance [START_REF] Badri Narayan Bhaskar | Atomic norm denoising with applications to line spectral estimation[END_REF][START_REF] Bredies | Inverse problems in spaces of measures[END_REF][START_REF] Candès | Towards a mathematical theory of superresolution[END_REF][START_REF] De | Exact reconstruction using Beurling minimal extrapolation[END_REF][START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF]), we consider instead sparse spike estimation methods which operate over a continuous domain, i.e. without resorting to some sort of discretization on a grid. The inverse problem is solved over the space of Radon measures which is a non-reflexive Banach space. This continuous "grid-free" setting makes the mathematical analysis easier and allows us to make precise statement about the location of the recovered spikes.

The technique that we are considering in this paper consists in solving a convex optimization problem that uses the total variation norm, which is the counterpart of the 1 norm for measures. It favors the emergence of spikes in the solution and is defined by

∀m ∈ M(X), |m|(X) def. = sup ψ∈C0 X ψdm ; ψ ∞,X 1 . (2)
In particular, for m a0,x0

def.

= N i=1 a 0,i δ x0,i ,

|m a0,x0 |(X) = a 0 1 ,
which shows in a way that the total variation norm generalizes the 1 norm to the continuous setting of measures (i.e. no discretization grid is required). When no noise is contaminating the data, one considers the classical basis pursuit, defined originally in [START_REF] Scott Shaobing Chen | Atomic decomposition by basis pursuit[END_REF] in a finite dimensional setting, written here over the space of Radon measures min m∈M(X) |m|(X) s.t. Φm = y 0 .

(P 0 (y 0 ))

This problem is studied in [START_REF] Candès | Towards a mathematical theory of superresolution[END_REF], in the case where Φ is an ideal low-pass filter on the torus X = T.

When the signal is noisy, i.e. when one observes y = y 0 + w, with w ∈ H, we may rather consider the problem min m∈M(X)

1 2
Φm -y 2 H + λ|m|(X).

(P λ (y))

Here λ > 0 is a parameter that should be adapted to the noise level w H . This problem is coined as BLASSO [START_REF] De | Exact reconstruction using Beurling minimal extrapolation[END_REF].

1.1.2. BLASSO performance analysis In order to quantify the recovery performance of the methods P 0 (y 0 ) and P λ (y), the following two questions arise:

(i) Does the solutions of P 0 (y 0 ) recover the input measure m a0,x0 ?

(ii) How close is the solution of P λ (y) to the solution of P 0 (y 0 ) ?

When the amplitudes of the spikes are arbitrary complex numbers, the answers to the above questions require a large enough minimum separation distance ∆(m a0,x0 ) between the spikes, where ∆(m a0,x0 ) def.

= min

i =j d X (x 0,i , x 0,j ).

(

When X = T, d X is the geodesic distance on the circle ∀x, y ∈ R, d X (x + Z, y + Z) = min k∈Z |x -y + k|.

In [START_REF] Candès | Towards a mathematical theory of superresolution[END_REF], the authors shows that for the ideal low-pass filter, m a0,x0 is the unique solution of P 0 (y 0 ) provided that ∆(m a0,x0 )

C fc where C > 0 is a universal constant and f c the cutoff frequency of the ideal low-pass filter. In the same paper, it is shown that C 2 when a 0 ∈ C N and C 1.87 when a 0 ∈ R N . In [START_REF] Fernandez-Granda | Super-resolution of point sources via convex programming[END_REF], the constant C is further refined to C

1.26 when a 0 ∈ R N . Suboptimal lower bounds on C were given in [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF][START_REF] Tang | Resolution limits for atomic decompositions via markov-bernstein type inequalities[END_REF]. Moreover, it was recently shown in [START_REF] Ferreira | A tight converse to the spectral resolution limit via convex programming[END_REF] that necessarily C 1 in the sense that for all ε > 0, and for f c large enough, there exist measures with ∆(m a0,x0 ) (1 -ε)/f c which are not identifiable using (P 0 (y 0 )).

The second question receives partial answers in [START_REF] Azaïs | Spike detection from inaccurate samplings[END_REF][START_REF] Bredies | Inverse problems in spaces of measures[END_REF][START_REF] Candès | Super-resolution from noisy data[END_REF][START_REF] Fernandez-Granda | Support detection in super-resolution[END_REF]. In [START_REF] Bredies | Inverse problems in spaces of measures[END_REF], it is shown that if the solution of P 0 (y 0 ) is unique then the measures recovered by P λ (y) converge in the weak-* sense to the solution of P 0 (y 0 ) when λ → 0 and w H /λ → 0. In [START_REF] Candès | Super-resolution from noisy data[END_REF], the authors measure the reconstruction error using the L 2 norm of an ideal low-pass filtered version of the recovered measures. In [START_REF] Azaïs | Spike detection from inaccurate samplings[END_REF][START_REF] Fernandez-Granda | Support detection in super-resolution[END_REF], error bounds are given on the locations of the recovered spikes with respect to those of the input measure m a0,x0 . However, those works provide little information about the geometrical structure of the measures recovered by P λ (y). That point is addressed in [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF] where the authors show that under the Non Degenerate Source Condition, there exists a unique solution to P λ (y) with the exact same number of spikes as the original measure provided that λ and w H /λ are small enough. Moreover in that regime, this solution converges to the original measure when the noise drops to zero.

BLASSO for positive spikes. For positive spikes (i.e. a 0,i > 0), the picture is radically different. Exact recovery of m a0,x0 without noise (i.e. (w, λ) = (0, 0)) holds whatever the distance between the spikes [START_REF] De | Exact reconstruction using Beurling minimal extrapolation[END_REF], but stability constants explode as ∆(m a0,x0 ) → 0. However, the authors in [START_REF] Veniamin | Super-resolution of positive sources: the discrete setup[END_REF] show that stable recovery is obtained if the signal-to-noise ratio grows faster than O(1/∆ 2N ). This closely matches the optimal lower bounds of O(1/∆ 2N -1 ) obtained by combinatorial methods [START_REF] Demanet | The recoverability limit for superresolution via sparsity[END_REF].

Finally, provided a certain nondegeneracy condition, it was recently shown in [START_REF] Denoyelle | Support recovery for sparse superresolution of positive measures[END_REF] that support recovery is guaranteed in the presence of noise if the signal-to-noise ratio grows faster than O(1/∆ 2N -1 ).

Solving the BLASSO

As the BLASSO is an optimization problem over the infinite dimensional space of Radon measures M(X), its resolution is challenging. We review in this section the existing approaches to tackle this problem. They can be roughly divided into three main families although there exists a flurry of generalizations and extensions that must be considered separately.

Fixed spatial discretization.

A common approach consists in constraining the measure to be supported on a grid. This leads to a finite dimensional convex optimization problem-known as LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] or basis pursuit [START_REF] Scott Shaobing Chen | Atomic decomposition by basis pursuit[END_REF]-for which there exist numerous solvers. These include the block-coordinate descent (BCD) algorithm [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF][START_REF] Tong | Coordinate descent algorithms for lasso penalized regression[END_REF], the homotopy/LARS algorithm [START_REF] Efron | Least angle regression[END_REF][START_REF] Soussen | Homotopy based algorithms for 0 -regularized least-squares[END_REF], or proximal forward-backward splitting algorithms [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] such as the Iterative Soft Thresholding (IST) [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]. Although simples to implement, the latters are in general slow to converge (the error in the objective function is typically of the order of O(1/k), where k is the number of iterations) [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF][START_REF] Mário | An EM algorithm for wavelet-based image restoration[END_REF]. However, there exist accelerated versions such as FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], which benefit from a better non-asymptotic rate of convergence (O(1/k 2 )). Finally, it is noteworthy that these proximal methods enjoy a linear asymptotic rate (see for instance [START_REF] Liang | Activity identification and local linear convergence of forward-backward-type methods[END_REF]), but this regime can be slow to reach.

The main limitation of these grid-based methods is that, in order to go below the Rayleigh limit and perform super-resolution, the grid must be thin enough. This leads to theoretical and practical issues. Indeed, refining the grid not only increases the computational cost of each iteration, but it also deteriorates the conditioning of the linear operator to invert. Hence, in practice, these methods provide solutions which are composed of small clusters of non-zero coefficents around each "true" spike. A way to mitigate this issue is to perform a post processing by replacing each cluster of spikes by its center of mass, as proposed in [START_REF] Tang | Sparse recovery over continuous dictionaries-just discretize[END_REF][START_REF] Flinth | Exact solutions of infinite dimensional total-variation regularized problems[END_REF]. This drastically reduces the number of false positive spikes although it is hard to analyze theoretically and can be unstable. Instead, one can also consider methods based on safe rules [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF] which perform a progressive pruning of the grid and keep only active sets of weights [START_REF] Massias | From safe screening rules to working sets for faster lasso-type solvers[END_REF]. Finally, it has been shown in [START_REF] Duval | Sparse spikes super-resolution on thin grids I: the Lasso[END_REF][START_REF] Duval | Sparse spikes super-resolution on thin grids II: the continuous basis pursuit[END_REF] that the solution of the LASSO, in a small noise regime and when the step size tends to zero, contains pairs of spikes around the true ones.

Fixed spectral discretization and semidefinite programming (SDP) formulation. In [START_REF] Candès | Towards a mathematical theory of superresolution[END_REF], the authors propose a reformulation of the Basis Pursuit for measures into an equivalent finite dimensional SDP for which solvers exist. Similarly, one can get an SDP formulation of the BLASSO. However, these equivalences are only true in a 1-dimensional setting. In higher dimensions (d 2), one needs to use the so-called Lasserre's hierarchy [START_REF] Bernard | Moments, positive polynomials and their applications[END_REF][START_REF] Bernard | Global optimization with polynomials and the problem of moments[END_REF]. This principle has been used for the super-resolution problem in [START_REF] De Castro | Exact solutions to super resolution on semi-algebraic domains in higher dimensions[END_REF].

The resolution of SDPs can be tackled through proximal splitting methods [START_REF] Toh | An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems[END_REF] as well as interior point methods [START_REF] Boyd | Convex optimization[END_REF]. However, the overall complexity of the latter is polynomial in O(f 2d c ), where d is the dimension of the domain X, which restricts its application to small dimensional problems. This limitation has led to recent developments [START_REF] Catala | A low-rank approach to off-the-grid sparse deconvolution[END_REF] where the authors proposed a relaxed low rank SDP formulation of the BLASSO in order to use a Frank-Wolfe-type method (see below). The resulting method enjoys the better overall complexity of O(f d c log(f c )) per iteration. Finally, note that these SDP-based approaches are restricted to certain type of forward operators (typically Fourier measurements). In contrast, grid-based proximal methods as well as Frank-Wolfe (directly on the BLASSO, see below) can be used for a larger class of operators Φ.

Optimization over the space of measures. In order to directly solve the BLASSO, one needs to design algorithms that do not use any Hilbertian structure and can instead deal with measures. The benefit is the fact that one can exploit advantageously the continuous setting of the problem (typically moving continuously spikes over the domain). In contrast to fixed spatial or spectral discretization methods, these algorithms proceed by iteratively adding new spikes, i.e. Dirac masses, to the recovered measure.

The Frank-Wolfe (FW) algorithm [START_REF] Frank | An algorithm for quadratic programming[END_REF] (see Section 4), also called the Conditional Gradient Method (CGM) [START_REF] Evgenii | Constrained minimization methods[END_REF], solves optimization problems of the form min m∈C f (m), where C is a weakly compact convex set of a topological vector space and f is a differentiable convex function (in the case of the BLASSO, m is a Radon measure). It proceeds by iteratively minimizing a linearized version of f . No Hilbertian structure is used which makes it well suited to work on the space of Radon measures. It has been proven under a curvature condition on f (which holds on a Banach space for smooth functions having a Lipschitz gradient) that the rate of convergence of this algorithm in the objective function is O(1/k) (see for instance [START_REF] Fedorovich | Approximate methods in optimization problems[END_REF]). However, it is possible to improve the convergence speed of FW by replacing the current iterate by any "better" candidate m ∈ C that further decreases the objective function f . This simple idea has led to several successful variations of the standard FW algorithm. For instance, the authors of [START_REF] Bredies | Inverse problems in spaces of measures[END_REF] proposed a modified Frank-Wolfe algorithm for the BLASSO where the final step updates the amplitudes and positions of spikes by a gradient descent on a non-convex optimization problem. Moving the spikes positions takes advantage of the continuous framework of the problem (the domain X is not discretized) which is the main ingredient that leads to a typical N -step convergence observed empirically. Finally, this approach has later been used in [START_REF] Boyd | The alternating descent conditional gradient method for sparse inverse problems[END_REF] and provides state of the art results in many sparse inverse problems such as matrix completion or Single Molecule Localization Microscopy (SMLM) [START_REF] Holden | Super-resolution fight club[END_REF][START_REF] Sage | Quantitative evaluation of software packages for single-molecule localization microscopy[END_REF].

Other methods for super-resolution

The Prony method [START_REF] De Prony | Essai expérimental et analytique sur les lois de la dilatabilité des fluides élastiques et sur celles de la force expansive de la vapeur de l'eau et de la vapeur de l'alcool à différentes températures[END_REF] and its successors such as MUSIC (MUltiple SIgnal Classification) [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], ESPRIT (Estimation of Signal Parameters by Rotational Invariance Techniques) [START_REF] Kailath | ESPRIT-estimation of signal parameters via rotational invariance techniques[END_REF], or Matrix Pencil [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF], are spectral methods which perform spikes localization from low frequency measurements. They do not need any discretization and enable to recover exactly the initial signal in the noiseless case as long as there are enough observations compared to the number of distinct frequencies [START_REF] Liao | MUSIC for single-snapshot spectral estimation: stability and super-resolution[END_REF]. Extensions to deal with noise have been developped in [START_REF] Cadzow | Signal enhancement-a composite property mapping algorithm[END_REF][START_REF] Condat | Cadzow denoising upgraded: a new projection method for the recovery of Dirac pulses from noisy linear measurements[END_REF] and stability is known under a minimum separation distance [START_REF] Liao | MUSIC for single-snapshot spectral estimation: stability and super-resolution[END_REF]. Greedy algorithms constitute another class of popular methods for sparse super-resolution. The Matching Pursuit (MP) [START_REF] Mallat | Matching pursuit with time-frequency dictionaries[END_REF] adds new spikes by finding the ones that best correlate with the residual. The Orthogonal Matching Pursuit (OMP) [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF][START_REF] Soussen | Joint k-step analysis of orthogonal matching pursuit and orthogonal least squares[END_REF][START_REF] Herzet | Relaxed recovery conditions for OMP/OLS by exploiting both coherence and decay[END_REF] is similar to MP but imposes that the current estimate of the observations, i.e. Φ( k i=1 a i δ xi ), is always orthogonal to the residual. Hence, the amplitudes of the Dirac masses are updated by an orthogonal projection after every support update (i.e. addition of a new spike). It is noteworthy that there exist many generalizations/variants of OMP. For instance, the results of OMP can be improved with a backtracking step at each iteration, allowing to remove non reliable spikes from the support of the reconstructed measure [START_REF] Huang | Backtracking-based matching pursuit method for sparse signal reconstruction[END_REF].

These greedy pursuit algorithms can be applied without grid discretization [START_REF] Jacques | A geometrical study of matching pursuit parametrization[END_REF] which enables the use of local optimizations over the spikes' positions [START_REF] Eftekhari | Greed is super: A fast algorithm for super-resolution[END_REF]. Finally, let us mention the class of nonconvex optimization methods which include the well known Iterative Hard Thresholding (IHT) [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF][START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF]

Contributions

Our first set of contributions, detailed in Section 3, studies the BLASSO performance in the special case of several types of Laplace transforms. This theoretical study is motivated by the use of these Laplace transform for certain types of fluorescence microscopy imaging devices. Our main finding is that for positive spikes, these operators can be stably inverted without minimum separation distance. This study makes use of the theoretical tools developed in our previous work [START_REF] Denoyelle | Support recovery for sparse superresolution of positive measures[END_REF].

Our algorithmic contributions are detailed in Section 4, where we introduce the Sliding Frank-Wolfe, which is an extension of the initial FW solver proposed in [START_REF] Bredies | Inverse problems in spaces of measures[END_REF]. Proposition 5 shows that this algorithm, used to solve the BLASSO, enjoys the same convergence property as the classical Frank-Wolfe algorithm (weak-* convergence with a rate in the objective function of O(1/k)). Our main theoretical contribution is Theorem 3 which proves that our algorithm converges towards the unique solution of the BLASSO in a finite number of iterations.

Section 5 makes the connection between these two sets of contributions, by showcasing the SFW algorithm for 3-D PALM/STORM super-resolution fluorescence microscopy. We study its performance for several imaging operators, among which some relies on the inversion of a Laplace transform along the depth axis.

The code to reproduce the numerical illustrations of this article can be found online at https://github.com/qdenoyelle.

Notations and Definitions

This section gathers some useful notations and definitions.

Ground space and measures. We frame our theoretical and numerical analysis of the BLASSO on the space of Radon measure over a set X.

Definition 1 (Set X of positions of spikes). The set of positions of spikes, denoted X, is supposed to be a subset of R d with non-empty interior Definition 2 (Continuous functions on X). Let (Y, • Y ) be a normed space. We denote by C c (X, Y ) the space of Y -valued continuous functions with compact support, by C 0 (X, Y ) the set of continuous functions that vanish at infinity i.e.

∀ε > 0, ∃K ⊂ X compact, sup x∈X\K ϕ(x) Y ε,
and by C k (X, Y ) the set of k-times differentiable functions on X. Note that when X is compact, C c (X, Y ) and C 0 (X, Y ) are simply the set C (X, Y ) of continuous functions on X.

Now we can define rigorously the space of real bounded Radon measures on X.

Definition 3 (Set M(X) of Radon measures). We denote by M(X) the set of real bounded Radon measures on X which is the topological dual of C 0 (X, R) endowed with • ∞,X (the supremum norm for functions defined on X).

By the Riesz representation theorem, M(X) is also the set of regular real Borel measures with finite total mass on X. See [START_REF] Rudin | Real and Complex Analysis, 3rd Ed[END_REF] for more details on Radon measures. Kernels. This paragraph details the assumptions that we use in the following on the kernel ϕ. We recall that the operator Φ : M(X) → H, which models the acquisition process of the source signal, has the form:

∀m ∈ M(X), Φm def. = X ϕ(x)dm(x). (5)
The above quantity is well-defined (as a Bochner integral) as soon as ϕ is continuous and bounded. In order to apply some results of [START_REF] Denoyelle | Support recovery for sparse superresolution of positive measures[END_REF], we add the hypotheses that are summarized below.

Definition 4 (Admissible kernels ϕ). We denote by KER (k) , the set of admissible kernels of order k. A function ϕ : X → H belongs to KER (k) if:

• ϕ ∈ C k (X, H),
• For all p ∈ H, x ∈ X → ϕ(x), p H vanishes at infinity,

• for all 0 i k, sup

x∈X D i ϕ(x) H < +∞.
where D i ϕ is the i-th differential of ϕ.

Operators. Given x = (x 1 , . . . , x N ) ∈

• X N , we denote by Φ x : R N → H the linear operator such that:

∀a ∈ R N , Φ x (a) def. = N i=1 a i ϕ(x i ), (6)
and by Γ

x : (R N × R N × • • • × R N d
) → H the linear operator defined by:

∀(a, b 1 , . . . , b d ) ∈ R N × (R N ) d , Γ x       a b 1 . . . b d       def. = N i=1   a i ϕ(x i ) + d j=1 b j,i ∂ j ϕ(x i )   . ( 7
)
We may also write Γ x = Φ x (Φ x ) (1) , where (Φ x ) (1) (sometimes denoted by Φ x ) stacks all the first order derivatives of ϕ for the different positions x i . Similarly we define (Φ x ) (k) for k 1 by stacking all the derivatives of order k. Finally, Γ + x refers to the pseudo-inverse of Γ x .

When d = 1, given x c ∈ • X, we denote by ϕ k ∈ H the k th derivative of ϕ at x c , i.e. ϕ k def. = ϕ (k) (x c ). (8)
In particular, ϕ 0 = ϕ(x c ). Given k ∈ N, we then define:

Ψ k def. = ϕ 0 ϕ 1 . . . ϕ k . (9)
Injectivity Assumption. In order to avoid degeneracy issues we sometimes assume the following injectivity assumption of the operator when restricted to discrete spikes.

Definition 5. Let ϕ : X → H. For all k ∈ N, we say that the hypothesis I k holds at = Φm a0,x0 + w are the noisy observations of a measure composed of a sum of Dirac masses. The optimality of a measure m λ for P λ (y) is characterized by the fact that the function

x c ∈ • X if
η λ def. = Φ * p λ where p λ def. = 1 λ (y -m λ ) (10)
satisfies η λ ∞,X 1. The function η λ is then called a dual certificate. When one is interested in the recovery of the support, i.e. finding a solution m a,x of P λ (y) composed of exactly the same number of Dirac masses as the initial measure m a0,x0 , in a small noise regime, an important object is the so-called vanishing derivatives precertificate introduced in [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF].

Definition 6 (Vanishing Derivatives Precertificate, [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF]). If Γ x0 has full column rank, there is a unique solution to the problem

inf { p H ; ∀i = 1, . . . , N, (Φ * p)(x 0,i ) = sign(a 0,i ), (Φ * p) (x 0,i ) = 0 R d } .
Its solution p V is given by

p V = (Γ + x0 ) * sign(a 0 ) 0 (R d ) N , (11)
and we define the vanishing derivatives precertificate as η V def.

= Φ * p V (Γ x0 is defined in Equation 7).

One can show that is if η V ∞,X
1 then η V is a a so-called valid certificate, which assures that m a0,x0 is a solution to the constrained problem (corresponding to setting w = 0 and λ → 0 in P λ (y)) min Φm=y0 |m|(X) where y 0 def.

= Φm a0,x0 .

(P 0 (y 0 ))

More importantly, if it satisfies a stronger nondegeneracy condition detailed in Definition 7 below, then η V also ensures the stable recovery of the support in a small noise regime when solving the BLASSO. This result proved in [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF] is stated in Theorem 1.

Definition 7 (Nondegeneracy of η V , [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF]). We say that η V is nondegenerate if

       ∀ x ∈ X \ N i=1 {x 0,i }, |η V (x)| < 1, ∀ i ∈ {1, . . . , N }, det(D 2 η V (x 0,i )) = 0. ( 12
)
Theorem 1 (Exact Support Recovery, [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF]). Assume that ϕ ∈ KER (2) , Γ x0 has full column rank and η V is nondegenerate. Then there exists

C > 0 such that if (λ, w) ∈ R * + × H satisfies: max (λ, w H /λ) C,
then there is a unique solution m a,x to P λ (y) composed of N Dirac masses such that (a, x) = g(λ, w) where g is C 1 . In particular, by taking the regularization parameter λ = w H /C proportional to the noise level, one obtains

|(a, x) -(a 0 , x 0 )| ∞ = O( w H ),
where |•| ∞ is the ∞ norm for vectors. Figure 9 displays some example of η V associated to several Φ operators for 3-D super-resolution fluorescence microscopy. This shows that for these inverse problems, the BLASSO stably recovers the support of the input measure if the noise level is not too high.

The Super-Resolution Problem

In this section, X is considered to be 1-dimensional and we now tackle the superresolution problem in presence of noise using the BLASSO. In this setting, we assume that the Dirac masses of the initial measure have positive amplitudes and cluster at some point x c ∈ • X. We parametrize this cluster as

m a0,tz0 def. = N i=1 a 0,i δ xc+tz0,i where a 0,i > 0, z 0,i ∈ R,
and where the parameter t > 0 controls the separation distance between the spikes of the input measure.

In [START_REF] Denoyelle | Support recovery for sparse superresolution of positive measures[END_REF], the authors proved that the recovery of the support in presence of noise in the limit t → 0 is controlled by the 2N -1 vanishing derivatives precertificate.

Proposition 1 (2N -1 Vanishing Derivatives Precertificate, [START_REF] Denoyelle | Support recovery for sparse superresolution of positive measures[END_REF]). If I 2N -1 holds at x c (see Definiton 5), there is a unique solution to the problem

inf p H ; (Φ * p)(x c ) = 1, (Φ * p) (x c ) = 0, . . . , (Φ * p) (2N -1) (x c ) = 0 .
We denote by p W its solution, given by

p W = (Ψ + 2N -1 ) * δ 2N where δ 2N def. = (1, 0, . . . , 0) T ∈ R 2N , (13)
and we define the 2N -1 vanishing derivatives precertificate as η W def.

= Φ * p W (see Equation 9for the definition of Ψ 2N -1 ).

Figure 1 shows η W in the case of a Gaussian convolution kernel. Remark 1. From Proposition 1, one can easily see that η W can equivalently be written as

∀x ∈ X, η W (x) = 2N -1 k=0 α k ∂ (k) 2 c ϕ (x, x c ), ( 14
)
where c ϕ is the correlation kernel associated to the correlation operator Φ * Φ, namely c ϕ (x, x ) = ϕ(x), ϕ(x ) H , and the coefficients α k are defined by the equations

∀k ∈ {0, . . . , 2N -1}, η (k) W (x c ) = δ k 0 . (15)
If η W satisfies some nondegeneracy property (see Definition 8) then one can prove that the recovery of the support in a small noise regime when t → 0 is possible. Theorem 2 (see [START_REF] Denoyelle | Support recovery for sparse superresolution of positive measures[END_REF]) makes this statement precise by quantifying the scaling between the noise level and the separation t to ensure the recovery. Definition 8. Assume that I 2N -1 holds at x c and ϕ ∈ KER (2N ) . We say that η W is

(2N -1)-nondegenerate if η (2N ) W (x c ) = 0 and for all x ∈ X \ {x c }, |η W (x)| < 1. Theorem 2. Suppose that ϕ ∈ KER (2N +1) and that η W is (2N -1)-nondegenerate.
Then there exist positive constants t 0 , C, M (depending only on ϕ, a 0 and z 0 ) such that for all 0 < t < t 0 , for all (λ, w) ∈ B 0, Ct 2N -1 with w H /λ C,

• the BLASSO has a unique solution,

• that solution has exactly N spikes, and it is of the form m a,xc+tz , with (a, z) = g(λ, w) (where g is a C 2N function), • the following inequality holds

|(a, z) -(a 0 , z 0 )| ∞ M |λ| t 2N -1 + w H t 2N -1 .
In the next section, we prove that the main assumption of Theorem 2 (the nondegeneracy of η W ) is satisfied for some operators Φ associated to Laplace measurements.

BLASSO for Laplace Inversion

Most existing theoretical studies of super-resolution are focussed on translationinvariant operator Φ (convolution or Fourier measurements), see Section 1.1. In contrast, this section presents new results for one of the most fundamental nontranslation invariant operator: the Laplace transform (and variants).

The behavior of the Laplace transform is radically different from the one of the Fourier transform, and understanding the impact of the lack of translation invariance on super-resolution is relevant for many applications in imaging, including those considered in Section 5. A first argument in favor of the BLASSO for the Laplace transform is the study provided in [START_REF] Duval | A characterization of the Non-Degenerate Source Condition in Super-Resolution[END_REF]. It essentially shows that the recovery of N positive spikes with stability of the support is possible using at least 2N measurements, regardless of the spacing of the spikes (and the spacings of the samples). The stability is asserted by showing that η V,t and η W are nondegenerate, using abstract T-systems arguments.

Our strategy here is different, as we provide closed form expressions for η W for these operators in order to show its nondegeneracy. The results presented here are thus complementary to those of [START_REF] Duval | A characterization of the Non-Degenerate Source Condition in Super-Resolution[END_REF], providing additional theoretical guarantees which backup our numerical observations. The main differences are

• we provide closed form expressions to η W ,

• some of the impulse responses we consider are L 2 -normalized, a case which is not covered by the theory of [START_REF] Duval | A characterization of the Non-Degenerate Source Condition in Super-Resolution[END_REF],

• we cannot deal with arbitrary samplings µ, contrary to [START_REF] Duval | A characterization of the Non-Degenerate Source Condition in Super-Resolution[END_REF].

In this section, we suppose that N spikes are clustered at the position x c ∈ • X (which appears in the following results because of the non translation invariance of the kernel).

In the next section, we first detail the different continuous operators considered. Then, Section 3.3 gives explicit formulas for η W in two different setups and shows that η W is (2N -1)-nondegenerate. Finally, Section 3.4 provides some numerical material concerning η W when the continuous kernels are approximated by a sampling.

Laplace Operators

We suppose in this section that X = [x min , x max ] ⊂ R * + is a compact interval, and that H = L 2 (R + , µ) for some Radon measure µ on R + . A generic Laplace measurement kernel is defined as

∀x ∈ X, ϕ(x) def. = s → ξ(x)e -sx ∈ H. ( 16
)
This choice ensures that ϕ defines a valid operator Φ for all the the Laplace-like transform models presented below, provided e -xmins dµ(s) has sufficiently many finite moments (in the following we require finite moments of order 4N -1). The kernel is parametrized by a positive Radon measure µ on X (which models the sampling pattern) and a non-negative weighting function ξ ∈ C (X, R) (which takes care of the normalization of the measurement). The adjoint operator is thus defined as

(Φ * p)(x) = ξ(x) R+ e -sx p(s)dµ(s).
The choice of µ is let to the experimentalist and corresponds to the way samples are chosen. A discrete measure µ = K k=1 µ k δ s k corresponds to using a finite set of samples values s k . In this case, one can equivalently consider finite-dimensional observations H = R K and define ϕ(x) def.

= (ξ(x)µ k e -s k x ) K k=1 ∈ H. A continuous measure dµ(s) = h µ (s)ds is a mathematical idealization, where a high value of h µ (s) indicates that a high number of measurements have been taken for the index s (or equivalently that there is less noise for this measurement). On contrast, a value µ(s) = 0 indicates that this measurement is not available.

In contrast, ξ can be freely chosen but strongly impacts the BLASSO problem by weighting the contribution of each position. The design of such a spatially-varying weighting is crucial (and non trivial) here because the operator Φ is not translationinvariant. The most frequent normalization for LASSO-type problems is

ξ(x) 2 = 1 R+ e -2sx dµ(s) , (17)
which guarantees that ϕ(x) H = 1 for all x ∈ X. See Section 3.3.2 for more details for this normalization. Note that both µ and ξ can be independently chosen, since they operate separately on the input and output variables x and s.

Correlation kernel. The properties of the BLASSO problem (and also the implementation of BLASSO solvers) only depend on the correlation operator Φ * Φ (rather than on the operator Φ itself).

This operator reads (Φ * Φm)(x) = X c ϕ (x, x )dm(x ) where c ϕ is a symmetric positive kernel. For Laplace-type operators, it reads

∀x, x ∈ X, c ϕ (x, x ) = ξ(x)ξ(x ) R+ e -(x+x )s dµ(s).
The choice of normalization [START_REF] Scott Shaobing Chen | Atomic decomposition by basis pursuit[END_REF] ensures that c ϕ (x, x) = 1.

We now detail in the following sections several particular cases covered by Equation 16and study the associated η W .

Preliminaries Results

This section gathers preliminary results useful for the computation of η W . One begins with two elementary lemmas. Their proofs are left to the reader. The first one is a simple consequence of the Faa di Bruno formula.

(x c ) = 1, η (x c ) = 0, . . . , η (2N -1) (x c ) = 0, (18)
if and only if ν

def. = η • h satisfies ν(t c ) = 1, ν (t c ) = 0, . . . , ν (2N -1) (t c ) = 0. (19)
Moreover, in that case,

ν (2N ) (t c ) = η (2N ) (x c )(h (t c )) 2N .
The next one follows from the general Leibniz rule.

η(x c ) = 1, η (x c ) = 0, . . . , η (2N -1) (x c ) = 0, (20) then P def.
= η × g satisfies:

P (x c ) = g(x c ), P (x c ) = g (x c ), . . . , P (2N -1) (x c ) = g (2N -1) (x c ). (21)
In particular, if

P ∈ R 2N -1 [T ]
, then P is the Taylor expansion of g at x c of order 2N -1, and η

(2N ) W (x c ) = -g (2N ) (x c )/g(x c ) provided that g(x c ) = 0.
c ϕ (x, x ) = 1 x + x . ( 22
)
The following Proposition provides a formula for η W in this unnormalized continuous setting and proves that it is nondegenerate. In Figure 2, one sees that when the position x c where the spikes cluster increases, the curvature of η W at x c decreases. This means that it is harder in this situation to perform the recovery. It reflects the exponential decay of the kernel ϕ.

Proposition 2. η W is (2N -1)-nondegenerate. More precisely, we have ∀x ∈ X, η W (x) = 1 - x -x c x + x c 2N . (23
Proof of Proposition 2. From Equations ( 14) and [START_REF] De | Exact reconstruction using Beurling minimal extrapolation[END_REF], one sees that η W has the form

η W (x) = 2N k=1 β k (x + x c ) k , where β k ∈ R. We set h : t → (1/t -x c ), ν def. = η • h so that ν(t) = 2N k=1 β k t k , is a polynomial with degree at most 2N with ν(0) = 0. By Lemma 1, ν satisfies (19) at t c def. = 1 2xc . As a result, ν(t) = 1 + β 2N (t -t c ) 2N . The constant β 2N is fixed by the condition ν(0) = 0, so that ν(t) = 1 -t-tc tc 2N
, and η W is given by [START_REF] De Castro | Exact solutions to super resolution on semi-algebraic domains in higher dimensions[END_REF].

The 2N derivative is ν (2N ) (t c ) = -(2N )! (tc) 2N , so that η W (x c ) = -(2N )! (2xc) 2N < 0.

L 2 -Normalized Laplace Operator

We choose µ = L, where L is the Lebesgue measure on R + , and

∀x ∈ X, ξ(x) = 1 R+ e -2sx ds = √ 2x, so that for all x ∈ X, ϕ(x) : s → √ 2xe -sx and ϕ(x) H = 1. One gets ∀x, x ∈ X, c ϕ (x, x ) def. = ϕ(x), ϕ(x ) H = 2 √ xx x + x . (24)
The following Proposition provides a formula for η W in this normalized setting and proves that it is nondegenerate. Proposition 3. η W is (2N -1)-nondegenerate. More precisely, we have the following formula: In Figure 3, one sees that when the position x c where the spikes cluster increases then the curvature of η W at x c decreases. The interpretation is the same as in the previous paragraph.

∀x ∈ X, η W (x) = 2 √ xx c x + x c N -1 k=0 (2k)! 2 2k (k!) 2 x -x c x + x c 2k . ( 25
Proof of Proposition 3. From the general Leibniz rule, we have for all n ∈ {0, . . . , 2N -1} and for all x, x ∈ X:

d n dx n (c ϕ (x, x )) = 2 √ x n k=0 n k d n-k dx n-k √ x d k dx k 1 x + x
Evaluating this expression at x = x c , one gets that:

∂ n 2 c ϕ (x, x c ) = √ x n k=0 α k (x + x c ) k+1 ,
for some coefficients α k ∈ R. As a result, η W is the unique function the form

η W (x) = √ x 2N -1 k=0 β k (x + x c ) k+1
for some coefficients β k ∈ R, which satisfies [START_REF] Candès | Towards a mathematical theory of superresolution[END_REF]. As before, we set t

= 1 x+xc , that is x = h(t) def. = 1 t -x c
, and h is a diffeomorphism of (0, 1/x c ) onto (0, +∞). Then:

η W • h(t) = 1 t -x c tP (t) = t -t 2 x c P (t),
where

P (T ) = 2N -1 k=0 β k T k ∈ R 2N -1 [T ]
. By Lemma 1 and Lemma 2, P is the Taylor expansion of order 2N -1 of

g : t → 1 √ t-t 2 xc at t c = h -1 (x c ) = 1 2xc . Setting t = u + 1 2xc
, we note that:

1 √ t -t 2 x c = 2 √ x c 1 -(2ux c ) 2 and 1 √ 1 -z 2 = N -1 k=0 (2k)! 2 2k (k!) 2 z 2k + o(z 2N -1 ). One deduces that 1 √ t -t 2 x c = 2 √ x c N -1 k=0 (2k)! 2 2k (k!) 2 [2x c (t -t c )] 2k + o((t -t c ) 2N -1 ).
As a result, P is given by

P (t) = 2 √ x c N -1 k=0 (2k)! 2 2k (k!) 2 [2x c (t -t c )] 2k and η W • h(t) = t -t 2 x c P (t) (26) = 1 - +∞ k=M (2k)! 2 2k (k!) 2 [2x c (t -t c )] 2k +∞ k=0 (2k)! 2 2k (k!) 2 [2x c (t -t c )] 2k . (27)
One sees that |η W • h(t)| < 1 for all t ∈ (0, 1 xc ) \ { 1 2xc }, and by Lemma 2,

(η W • h) (2N ) (t c ) = -g (2N ) (t c )/g(t c ) = - ((2N )!) 2 (N !) 2 x 2N c < 0 (28) so that η W • h (hence η W ) is (2N -1)-nondegenerate. One recovers η W by composing with h -1 , noting that 2x c (t -t c ) = xc-x x+xc .

Sampled Approximations

The previous two cases (normalized and unnormalized versions of the Laplace transform) correspond to mathematical idealizations. In practice, one needs to restrict the sampling patterns by limiting their ranges and considering discrete samples. The following two setups are involved in the application of Section 5.

Discretized Unnormalized Laplace.

We assume that µ = K-1 k=0 δ s k and ξ = 1. Then ϕ(x) = (e -s k x ) K-1 k=0 ∈ R K and:

c ϕ (x, x ) = K-1 k=0 e -s k (x+x ) . Discretized L 2 -normalized Laplace. We let µ = K-1 k=0 δ s k and ξ(x) = K-1 k=0 e -2s k x -1/2
. Then ϕ(x) = ξ(x)(e -s k x ) K-1 k=0 ∈ R K , ϕ(x) H = 1 and:

c ϕ (x, x ) = ξ(x)ξ(x ) K-1 k=0 e -s k (x+x ) .
In contrast to the continuous setups of Section 3.3, we do not have closed-form expressions for η W . However, if a sequence of measures, e.g. µ n = Kn-1 k=0 µ n,k δ s n,k converges in a suitable sense towards the Lebesgue measure µ = L, the following proposition shows that the corresponding η W must be nondegenerate for n large enough. We consider both the unnormalized and L 2 -normalized setups, corresponding respectively to

c ϕn (x, x ) = R+ e -(x+x
) 29
Then, both in the unnormalized and the L 2 -normalized case, for n large enough, the 2N -1 vanishing derivatives precertificate η W,n is (2N -1)-nondegenerate.

Proof. Let us denote by Ψ

[n]

2N -1 = (ϕ 0 , . . . , ϕ 2N -1 ) (resp. Ψ 2N -1 ) the impulse response derivatives corresponding to µ n (resp. µ = L), and by η W the 2N -1 vanishing derivatives precertificate for µ = L. First, in view of Sections 3.3.1 and 3.3.2, we observe that the result follows immediately if we prove that

lim n→+∞ Ψ [n] * 2N -1 Ψ [n] 2N -1 = Ψ * 2N -1 Ψ 2N -1 , (30)
(as it implies the linear independence of (ϕ 0 , . . . , ϕ 2N -1 ) for n large enough), and that

∀i ∈ {0, 1, . . . , 2N }, lim n→+∞ η (i) W,n -η (i) W ∞,X = 0, (31)
(as it implies |η W,n (x)| < 1 for x = x c and η

(2N )
W,n (x c ) < 0 for n large enough). We recall from ( 14) that η W,n is given by η

W,n (x) = 2N -1 i=0 α [n] i ∂ (i) 2 c ϕn (x, x c ) where α [n] = (Ψ [n] * 2N -1 Ψ [n]
2N -1 ) -1 δ 2N (provided the matrix is invertible), and the (i, j)entry of (Ψ

[n] * 2N -1 Ψ [n] 2N -1 ) is ∂ (i) 1 ∂ (j)
2 c ϕn (x c , x c ). As a consequence, both [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF] and [START_REF] Cardoso | Topography of cells revealed by variable-angle total internal reflection fluorescence microscopy[END_REF] are established if we can prove that

lim n→+∞ sup x,x ∈[xmin,xmax] ∂ (i) 1 ∂ (j) 2 c ϕn (x, x ) -∂ (i) 1 ∂ (j) 2 c ϕ (x, x ) = 0, (32)
for all i ∈ {0, . . . , 2N }, j ∈ {0, . . . , 2N -1}. First, we prove (32) in the unnormalized case, i.e. c ϕn (x, x ) = R+ e -(x+x )s dµ n (s). The dominated convergence theorem ensures that

∂ (i) 1 ∂ (j) 2 c ϕn (x, x ) =
R+ s i+j e -(x+x )s dµ n (s) (and similarly for c ϕ and µ).

Let (x, x ) ∈ [x min , x max ] 2 and let ψ ∈ C c (R + ) such that ψ(s) = 1 for s ∈ [0, 1], ψ (
s) = 0 for s 2, and 0 ψ 1 on R + . We denote by C the supremum in [START_REF] David | Super-resolution via sparsity constraints[END_REF].

Let ε > 0 and A > 0. Then,

R+ s i+j e -(x+x )s dµ n (s) - R+ s i+j e -(x+x )s ds R+ s i+j e -(x+x )s ψ s A dµ n (s) - R+ s i+j e -(x+x )s ψ s A ds def. = a + R+ s i+j e -(x+x )s (1 -ψ s A )dµ n (s) =b + R+ s i+j e -(x+x )s (1 -ψ s A )ds =c .
We have We choose A > 0 sufficiently large so that +∞ A

c +∞ A (1 + s 4N -1 )e -2xmins
(1 + s 4N -1 )e -2xmins ds ε and e -xminA C ε, hence max(b, c) ε. Now, to prove that a is uniformly small for (x, x ) ∈ [x min , x max ] 2 as n → +∞, we apply Lemma 3 to ((x, x ), s) → s i+j e -(x+x )s ψ s A defined on [x min , x max ] 2 × [0, 2A]. This yields the desired result.

The proof for the normalized case readily follows from the uniform convergence of the unnormalized case and the fact that the normalization factors

ξ n (x) = R+ e -2sx dµ n (s) -1/2 R+ e -2sxmax dµ n (s) -1/2
are upper bounded by some positive constant independent of n. Proof. We note that the mapping (η, ν) → S ηdν is continuous on C (S) × M(S). Since x → ψ(x, •) is continuous from X to C (S), the mapping

F : (x, ν) -→ S ψ(x, s)dν(s) (33)
is continuous on X × M(S). Now, since S is compact, M(S) is the dual of the Banach space C (S), and the Banach-Steinhaus theorem implies that there exists R > 0 such that sup n |µ n |(S) R (and |µ|(S) R).

The subspace B R def.

= {ν ∈ M(S) ; |ν|(S) R} is metrizable for the weak-* topology and compact. As a result, the mapping F is uniformly continuous on the compact X × B R . In particular, as

µ n → µ in B R , sup x∈X S ψ(x, s)dµ n (s) - S ψ(x, s)dµ(s) → 0.
Figure 4 illustrates this convergence between the precertificates in the unnormalized case.

The SFW Algorithm

In this section, we present the Sliding Frank-Wolfe (see Algorithm 2), a new version of the modified Frank-Wolfe algorithm introduced in [START_REF] Bredies | Inverse problems in spaces of measures[END_REF]. Moreover, we prove in Theorem 3 that it converges in a finite number of steps under mild assumptions. The code can be found in https://github.com/qdenoyelle.

We suppose in this section that X ⊂ R d is compact, or X = T d with d ∈ N * and ϕ ∈ KER (2) (see Definition 4).

The Algorithm

Frank-Wolfe Algorithm. The Frank-Wolfe (FW) algorithm [START_REF] Frank | An algorithm for quadratic programming[END_REF], also called the Conditional Gradient Method (CGM) [START_REF] Evgenii | Constrained minimization methods[END_REF] solves the following optimization problem

min m∈C f (m), ( 34
)
where C is a weakly compact convex set of a Banach space, and f is a differentiable convex function. For instance, in the case of sparse recovery problems, m is a measure and C is a subset of M(X). A chief advantage of FW with respect to most first order optimization scheme (such as gradient descent or proximal splitting method) is that it does not rely on any underlying Hilbertian structure and only makes use of directional derivatives. It is thus particularly well adapted to optimize over the space of Radon measures. The algorithm is detailed in Algorithm 1.

Algorithm 1 Frank-Wolfe Algorithm 1: for k = 0, . . . , n do 2:

Minimize:

s [k] argmin s∈C f (m [k] ) + df (m [k] )[s -m [k] ]. 3: if df (m [k] )[s [k] -m [k] ] = 0 then 4:
m [k] solution of [START_REF] Duval | Sparse spikes super-resolution on thin grids I: the Lasso[END_REF]. Stop.

5:

else 6:

Step research:

γ [k] ← 2 k+2 or γ [k] argmin γ∈[0,1] f (m [k] + γ(s [k] -m [k]
)).

7:

Update:

m [k+1] ← m [k] + γ [k] (s [k] -m [k]
).

8:

end if 9: end for Let us note that the FW algorithm is naturally endowed with a stopping criterion in Step 3 (see for instance [26, Ch. 3, Sec.1.2]) which is equivalent to the standard optimality condition for constrained convex problems

∀s ∈ C, df (m [k] )[s -m [k] ] 0. ( 35
)
Frank-Wolfe for the BLASSO. The FW algorithm cannot be applied directly to the BLASSO because it is an optimization problem over M(X) which is not bounded and the objective function

∀m ∈ M(X), T λ (m) def. = 1 2 Φm -y 2 H + λ|m|(X), (36)
is not differentiable. Instead, we propose to consider an equivalent problem to the BLASSO, using an epigraphical lift (following an idea of [START_REF] Harchaoui | Conditional gradient algorithms for norm-regularized smooth convex optimization[END_REF]), which is presented in Lemma 4.

Lemma 4. The BLASSO

min m∈M(X) T λ (m) def. = 1 2 Φm -y 2 H + λ|m|(X). (P λ (y)) is equivalent to min (t,m)∈C Tλ (m, t) def. = 1 2 Φm -y 2 H + λt, ( Pλ (y))
where we defined

C def. = {(t, m) ∈ R + × M(X) ; |m|(X) t M } and M def. = y 2

H

2λ . The equivalence stated in Lemma 4 is to be understood in the following sense: m is a solution to (P λ (y)) if and only if (t, m) is a solution to ( Pλ (y)) for some t 0. Moreover, in that case t = |m|(X) and Tλ (m, t) = T λ (m). As a result, one can directly translate the FW algorithm (see Algorithm 1) to P λ (y).

Proof. Let m * be a minimizer of T λ on M(X), then we have

T λ (m * ) T λ (0) = λM. (37)
Hence, one can restrict the BLASSO to the set of measures m ∈ M(X) such that |m|(X) M and Pλ (y) is obtained using an epigraphical representation.

The next two remarks discuss the applicability of standard results on FW to the BLASSO. Remark 2 (Well-posedness). The FW algorithm is well defined for Pλ (y). Indeed, Tλ is a differentiable functional on the Banach space R × M(X), with differential

d Tλ (t, m) : (t , m ) -→ X Φ * (Φm -y)dm + λt . ( 38
)
Although C is not weakly compact (otherwise, by the Eberlein-Shmulyan theorem, M(X) would be reflexive), it is compact for the weak-* topology: as d Tλ (t, m) is represented by (λ, Φ * (Φm -y)) ∈ R × C 0 (X), it does reach its minimum on C. Remark 3 (Rate of convergence). Let us note that d Tλ is Lipschitz continuous (because ϕ ∈ KER (2) ), hence by classical results for the study of the convergence of the FW algorithm, one obtains the O(1/k) rate of convergence in the objective function for any minimizing sequence for the BLASSO.

Lemma 5 ([26, Th. 3.1.7]). Let (t k , m [k] ) k∈N be a sequence generated by Algorithm 1 applied to Pλ (y). Then, there exists C 1 > 0 such that for any m * solution of P λ (y)

we have ∀k ∈ N * , T λ (m [k] ) -T λ (m * ) C 1 k . (39)
Next, we discuss how the minimization step yields a greedy approach and a natural stopping criterion. The following two remarks also crucially relate the algorithm to the dual certificate of (P λ (y)). Remark 4 (Greedy approach). Obviously, the FW algorithm is only interesting if, in step 2 of Algorithm 1, one is able to minimize the linear form s → d Tλ (t [k] , m [k] )[s] on C. That linear form reaches its minimum at least at one extreme point of C, i.e. s = (0, 0) or points of the form s = (M, ±M δ x ) for x ∈ X. Finding a minimizer among those points amounts to finding a point x in

argmin x∈X ± 1 λ Φ * (y -Φm [k] ) (x) + 1 λM, or equivalently in argmax x∈X η [k] (x) -1 where η [k] def. = 1 λ Φ * (y -Φm [k] )
(note the similarity of η [k] with the dual certificate defined in [START_REF] Boyd | The alternating descent conditional gradient method for sparse inverse problems[END_REF]). As a consequence, at each Step 7 of Algorithm 1, a new spike is created at some point in argmax X η [k] (unless s = (0, 0) is optimal, which means that η [k] ∞,X 1). This spike creation step is at the core of the algorithms in [START_REF] Bredies | Inverse problems in spaces of measures[END_REF] and [START_REF] Boyd | The alternating descent conditional gradient method for sparse inverse problems[END_REF]. Remark 5 (Stopping criterion). It is interesting to relate the stopping criterion

(t [k] , m [k] ) ∈ argmin s∈C d Tλ (t [k] , m [k] )[s],
with the dual certificates for (P λ (y)). As noted above (see Equation 35), the stopping criterion is equivalent to (t [k] , m [k] ) being a solution, hence t

[k] = |m [k] |(X). If m [k] = 0, without loss of generality we write m [k] = N [k] i=1 a [k] i δ x [k] i where the x [k] i 's are distinct, so that t [k] = |m [k] |(X) = i a [k] i . We also set ε [k] i def.
= sign(a = d Tλ (t [k] , m [k] ).

Assume first that |m [k] |(X) < M , so that the smallest face of C which contains (t [k] , m [k] ) is

F def. = conv (0, 0), (M, M ε [k] 1 δ x [k] 1 ), . . . , (M, M ε [k] N [k] δ x [k] N [k]
) .

Since argmin s∈C L is a face of C containing (t [k] , m [k] ) (see [START_REF] Tyrell | Convex analysis[END_REF]Sec. 18]), it must contain F . Hence

L(0, 0) = L(M, M ε [k] 1 δ x [k] 1 ) = • • • = L(M, M ε [k] N [k] δ x [k] N [k] ) = min C L. (40)
Now, if |m [k] |(X) = M , it means that Tλ (t [k] , m [k] ) = Tλ (0, 0), so that by convexity of Tλ and optimality of (t [k] , m [k] ) one has

L(t [k] , m [k] ) = d Tλ (t [k] , m [k] )[t [k]
, m [k] ] = 0 = L(0, 0). As the smallest face which contains (t [k] , m [k] ) is

F def. = conv (M, M ε [k] 1 δ x [k] 1 ), . . . , (M, M ε [k] N [k] δ x [k] N [k]
) , we deduce as above that (40) holds.

In particular, L(0, 0) inf x∈X L(M, ±M δ x ) yields

0 inf x∈X -η [k] (x) + 1 , (41)
that is η [k] ∞,X 1. Moreover L(t [k] , m [k] ) = N [k] j=1 a [k] j L M, M ε [k] j δ x [k] j N j=1 a [k] j L(M, ±M δ x [k] j ) , yields - N [k] j=1 a [k] j η [k] (x [k] j ) - N [k] j=1 a [k] j η [k] (x [k] j ) , from which we deduce η [k] (x [k] j ) = sign(a [k] j
). As a result, when the FW algorithm stops (if it does), we observe that the quantity η [k] it has constructed is the dual certificate for (P λ (y)). If m [k] = 0, the argument is similar (as (41) must hold).

The Sliding Frank-Wolfe algorithm. Applying directly Algorithm 1 yields a sequence of measures (m [k] ) k∈N which weakly-* converges towards some solution m * in a greedy way. But the generated measures m [k] are not very sparse compared to m * , each Dirac mass of m * being approximated by a multitude of Dirac masses of m [k] with inexact positions. It is therefore suggested in [START_REF] Bredies | Inverse problems in spaces of measures[END_REF], and strongly advocated in [START_REF] Boyd | The alternating descent conditional gradient method for sparse inverse problems[END_REF], to modify the Frank-Wolfe algorithm for the resolution of the BLASSO and to let the Dirac positions move.

m [k] = N [k] i=1 a [k] i δ x [k] i , a [k] i ∈ R, x [k] i pairwise distincts, find x [k] * ∈ X s.t.: x [k] * ∈ arg max x∈X |η [k] (x)| where η [k] def. = 1 λ Φ * (y -Φm [k] ), 4: if |η [k] (x [k] * )| 1 then 5:
m [k] is a solution of P λ (y). Stop.

6:

else 7:
Obtain

m [k+1/2] = N [k]
i=1 a

[k+1/2] i δ x [k] i + a [k+1/2] N [k] +1 δ x [k] *
, s.t.:

a [k+1/2] ∈ arg min a∈R N [k] +1 1 2 Φ x [k+1/2] a -y 2 H + λ a 1 where x [k+1/2] = (x [k] 1 , . . . , x [k] N [k] , x [k] * ) 8:
Obtain

m [k+1] = N [k] +1 i=1 a [k+1] i δ x [k+1] i , s.t.: (a [k+1] , x [k+1] ) ∈ arg min (a,x)∈R N [k] +1 ×X N [k] +1 1 2 Φ x a -y 2 H + λ a 1 ,
using a non-convex solver initialized with (a [k+1/2] , x [k+1/2] ).

9:

Eventually remove zero amplitudes Dirac masses from m [k+1] .

10:

end if 11: end for One important feature of the FW algorithm, as noted in [START_REF] Jaggi | Revisiting frank-wolfe: Projection-free sparse convex optimization[END_REF][START_REF] Boyd | The alternating descent conditional gradient method for sparse inverse problems[END_REF], is that in the update step 7, the point m [k+1] may be replaced with any point m ∈ C which has lower energy, without breaking the convergence property and the convergence rate. The Frank-Wolfe algorithm with our modified update step is described in Algorithm 2, we call it the Sliding Frank-Wolfe (SFW) algorithm. Since the t variable is only auxiliary in ( Pλ (y)), we omit it and we formulate directly Algorithm 2 in terms of m only.

As we detail below, the algorithm slightly (but crucially) differs from the one in [START_REF] Boyd | The alternating descent conditional gradient method for sparse inverse problems[END_REF]. The main ingredient is to replace the final update with the minimization of a non-convex minimization problem updating both the positions and the amplitudes of the spikes (whereas [START_REF] Boyd | The alternating descent conditional gradient method for sparse inverse problems[END_REF] update successively the amplitudes and the positions). Remark 6 (Links between FW applied to Pλ (y) and the SFW). Algorithm 2 is a valid variant of FW, as the update step decreases more the energy than the standard convex combination using γ [k] . Indeed,

T λ (m [k+1] ) T λ (m [k+1/2] ) T λ (m [k] + γ [k] (sign(η [k] (x [k] * ))M δ x [k] * -m [k] )).
It is noteworthy that other forms were previously used in [START_REF] Bredies | Inverse problems in spaces of measures[END_REF][START_REF] Boyd | The alternating descent conditional gradient method for sparse inverse problems[END_REF], but, to our knowledge, the update procedure (Steps 7 and 8) described in the present paper is new. As we show in Theorem 3, optimizing over both the amplitudes and the positions is essential to prove the convergence of the algorithm in a finite number of iterations.

Remark 7 (Stopping criterion of the SFW). One may observe that the condition

df (m [k] )[s [k] -m [k] ] = 0 of Algorithm 1 (or equivalently m [k] ∈ argmin s∈C df (m [k] )[s]) has been replaced with |η [k] (x [k] * )| 1.
In fact the optimality conditions for the non-convex local descent (Step 8) at iteration k -1 imply

∀i ∈ {1, . . . , N [k] }, η [k] (x [k] i ) = sign(a [k] i ), whereas |η [k] (x [k] * )| 1 implies η [k] ∞,X 1, hence η [k] is a valid dual certificate.
With the words of Remark 5, Step 8 implies that L M, M ε

[k] j δ x [k] j = 0 for 1 j N [k] , whereas the condition |η [k] (x [k] * )| 1 means 0 = L(0, 0) = min C L.
As m is a convex combination of those points, we deduce that (|m [k] |(X), m [k] ) ∈ argmin C L, that is the optimality condition [START_REF] Duval | Sparse spikes super-resolution on thin grids II: the continuous basis pursuit[END_REF]. Remark 8 (Adaptation for the positive BLASSO). In many applications, one is often interested in recovering positive spikes (see for example in Section 5). As a result, in these cases it is better to add a positivity constraint m 0 to the BLASSO. This leads to several changes in Algorithm 2

• the stopping condition |η [k] (x [k] * )| 1 becomes η [k] (x [k] * ) 1, • the LASSO is solved on R N [k] +1 + , • the optimization problem of Step 8 is solved on R N [k] +1 + × X N [k] +1 .

Implementation details.

• A Newton method, initialized by a grid search, is used to to find the maximum of |η [k] | over the compact domain X in step 3. The size of the grid depends on the operator Φ. For example, when Φ is the convolution by the Dirichlet kernel with cutoff frequency f c , we choose a number of points proportional to f c .

• The LASSO problem at step 7 is solved using the fast iterative shrinkage thresholding algorithm (FISTA) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF].

• To solve the non-convex optimization problem at step 8, we deploy a bounded BFGS. It allows to enforce the positions x i to be in the compact domain X and to preserve the sign of the amplitudes a i . These constraints ensure the differentiability of the objective function which is required by BFGS.

Study of the Convergence of the SFW Algorithm

We now study the convergence properties of the Sliding Frank-Wolfe algorithm presented last section (see Algorithm 2). Our main result is Theorem 3 where one shows that if η λ = 1 λ Φ * (y -Φm a,x ), where m a,x = N i=1 a i δ xi , is the unique solution of P λ (y) and is nondegenerate (see Equation ( 42)), then Algorithm 2 recovers m a,x in a finite number of iterations. But, first, one shows that our algorithm produces a sequence of measures (m [k] ) k∈N that converges towards m * (if m * ∈ M(X) is the unique solution of the BLASSO) for the weak-* topology on M(X).

Proposition 5. Let (m [k] ) k∈N be the sequence obtained from the Sliding Frank-Wolfe algorithm. Then it has an accumulation point for the weak-* topology on M(X), and that point is a solution to (P λ (y)).

Proof. By Remark 6, we know that (m [k] ) k∈N is a sequence obtained by applying Algorithm 1 to Pλ (y) where the final update is step 7 and 8 of the SFW. As a result, using Lemma 5, one gets that for any m * solution of P λ (y),

∀k ∈ N, T λ (m [k] ) -T λ (m * ) C 1 k .
Hence (m [k] ) is a bounded minimizing sequence. One can extract from it a subsequence that converges towards some m ∈ M(X) (with |m|(X) M ) for the weak-* topology. Since T λ is convex and l.s.c., it is also weak-* l.s.c. so that one obtains:

T λ (m) = T λ (m * ).
Hence m is a solution of P λ (y).

From this Proposition, one easily deduces the following Corollary.

Corollary 1. If m * ∈ M(X)
is the unique solution of P λ (y) then (m [k] ) k∈N weak-* converges towards m * .

In fact, under mild assumptions, our algorithm even converges towards the solution of the BLASSO in a finite number of iterations, thanks to the displacement of the spikes over the continuous domain X. For the sake of clarity, we state and prove this Theorem in the case of d = 1 but the changes for d ∈ N * can be easily done.

Theorem 3. Suppose that ϕ ∈ KER (2) , that m a,x = N i=1 a i δ xi is the unique solution of P λ (y), and that

η λ = 1 λ Φ * (y -Φm a,x ) is nondegenerate, i.e. ∀x ∈ X \ N i=1 {x i }, |η λ (x)| < 1 and ∀i ∈ {1, . . . , N }, η λ (x i ) = 0. ( 42
)
Then Algorithm 2 recovers m a,x after a finite number of steps (i.e. there exists k ∈ N such that m [k] = m a,x ).

Proof. Since m a,x is the unique solution of P λ (y), one knows by Corollary 1 that the sequence (m [k] ) k∈N produced by Algorithm 2 converges for the weak-* topology towards m a,x . As Φ is weak-* to weak continuous and by defining p [k] def. = 1 λ (y -Φm [k] ), one gets that (p [k] ) k∈N converges towards p λ in the weak topology of H and that η [k] def.

= Φ * p [k] converges pointwise towards η λ . Then one can show that Φ * is a compact operator. Indeed, for any bounded subset A ∈ H, one can check easily that Φ * A is equicontinuous and pointwise relatively compact so that by Ascoli theorem Φ * A is relatively compact for the strong topology of C 0 (X, R). As a result one can extract a subsequence of (η [k] ) k∈N that converges towards η λ in uniform norm. η λ is then the unique accumulation point in the uniform norm of the bounded sequence (η [k] ) k∈N hence its convergence towards η λ in uniform norm. One can repeat this argument for (η [k] ) k∈N and (η [k] ) k∈N (since ϕ ∈ KER (2) ), obtaining for all j ∈ {0, 1, 2}

(η [k] ) (j) • ∞,X -→ k→+∞ η (j) λ . (43)
Because η λ is nondegenerate, there exists a small neighborhood around each x i on which η λ = 0. Hence, we deduce from Equation ( 43) that there exist ε > 0 and k 1 ∈ N such that:

∀k k 1 , ∀i ∈ {1, . . . , N }, ∀x ∈]x i -ε, x i + ε[, η [k] (x) = 0.
We denote in the following

I xi,ε def. = ]x i -ε, x i + ε[, ∀i ∈ {1, . . . , N }.
Since m [k] converges towards m a,x in the weak-* topology and |m a,x | does not charge the boundary of

I xi,ε , we have ∀i ∈ {1, . . . , N }, m [k] (I xi,ε ) → m a,x (I xi,ε ) = a i = 0,
so that there exists k 2 ∈ N such that for all k k 2 , m [k] has at least one spike in each I xi,ε . In particular m [k] has at least N spikes.

Again, from Equation ( 43), since (η [k] ) k∈N converges uniformly towards η λ , one deduces that there exists k 3 ∈ N such that for all k k 3 :

Sat ± (η [k] ) ⊂ Sat ± (η λ ) ⊕ (] -ε, ε[×{0}) ,
where the set of saturation points of a given η ∈ C 0 (X, R) is defined as:

Sat ± (η) def. = {(x, v) ∈ X × {-1, 1}; η(x) = v} . Moreover, ∀x ∈ X \ N i=1 I xi,ε , |η [k] (x)| < 1.
In particular for k k 3 , m [k] has no spikes in X \ N i=1 I xi,ε because it would contradict the optimality conditions of Step 8 of Algorithm 2: for all i ∈ {1, . . . , N

[k] }, η [k] (x [k] i ) = sign(a [k] i ). Suppose now that k max(k 1 , k 2 , k 3 ).
Then m [k] has at least one spike in each neighborhood of x i and no spikes outside. Moreover |η [k] | < 1 outside the neighborhoods and η [k] = 0 inside. Let i ∈ {1, . . . , N } and denote x

[k] j ∈ I xi,ε a position of a spike of m [k] . From the optimality conditions of Step 8, one has also that η

[k] (x [k] j ) = 0. This combined with η [k] = 0 in I xi,ε implies that |η [k] | < 1, except at x [k]
j . Hence, m [k] has exactly one spike in this neighborhood. As a consequence, we proved that m [k] has exactly N spikes (one inside each neighborhood) and:

∀x ∈ X \ N i=1 {x [k] i }, |η [k] (x)| < 1.
Hence m [k] , composed of N spikes, is a solution of P λ (y). Since m a,x is supposed to be the unique solution of P λ (y), one concludes that:

m [k] = m a,x ,
i.e. the algorithm recovers m a,x in a finite number of iterations.

Note that one proved the convergence in a finite number of iterations but not exactly N iterations if m a,x is composed of N spikes. However in practice this is exactly what we observe.

Illustration of the N -Steps Convergence of the SFW

We now illustrate how the algorithm works and we shows that it converges in exactly N iterations in practice (when the noise level and the regularization parameter are appropriate, i.e. max(λ, w H /λ) is low enough).

We consider X = [0, 1] and a convolution operator with a sampled Gaussian kernel for Φ

Φ : m ∈ M(X) → [0,1] ϕdm ∈ R K where ϕ(x) = 1 √ 2πσ 2 e - ( i-1 K-1 -x) 2 2σ 2 1 i K .
We set σ = 0.05 and K = 100. The initial measure used is m a0,x0 = 1.3δ 0.3 +0.8δ 0.37 + 1.4δ 0.7 and the noise is small (w = 10 -4 w 0 where w 0 = randn(K)). Figure 5 shows η V for this configuration. One can see that it is nondegenerate. Hence, in a small noise now regime, with the appropriate choice of λ, there is a unique measure solution of P + λ (y) which is composed of the same number of spikes as m a0,x0 . Moreover, by Theorem 3, the SFW algorithm recovers it in a finite number of iterations. The decrease of the objective function throughout the algorithm iterations (cumulative iterations of BFGS) is presented in Figure 6. As indicated by the two vertical black lines, which show the intermediate iterations, the algorithm converges in exactly 3 iterations. One can observe an important decrease of the objective function each time a spike is added. Also, it is noteworthy that BFGS converges with very few iterations when k = 0 and k = 1 (first two spikes added) and that the main computational load for the non-convex step occurs for k = 2 (more iterations of BFGS). Figure 7 shows m [k] and η [k] at different times of the algorithm. More precisely, for k ∈ {0, 1, 2}, we display the initial measure m a0,x0 , the recovered measure, and the associated η. Moreover, we present them after the LASSO step (i.e. m [k+1/2] and η [k+1/2] ) as well as after the BFGS step (i.e. m [k+1] and η [k+1] ) .

One remarks, as expected, that for all i, η [k+1/2] (x i ) = 1, η [k+1] (x i ) = 1 and η [k+1] (x i ) = 0. In the first two main iterations, the spikes are almost not moved by the BFGS. However, at the last iteration, the displacement of the positions and amplitudes of the spikes is crucial to obtain η [k+1] ∈ ∂|m [k+1] |(X), and thus recover the solution of P + λ (y) in three steps.

Single Molecule Localization Microscopy

The field of fluorescent microscopy has experienced an important revolution during the past two decades with the emergence of super-resolution techniques. These modalities, such as structured illumination microscopy (SIM) [START_REF] Gustafsson | Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy[END_REF], stimulated emission depletion (STED) [START_REF] Hell | Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy[END_REF], or single molecule localization microscopy (SMLM)-which includes photoactivated localization microscopy (PALM) [START_REF] Betzig | Imaging intracellular fluorescent proteins at nanometer resolution[END_REF][START_REF] Hess | Ultra-high resolution imaging by fluorescence photoactivation localization microscopy[END_REF] and stochastic optical reconstruction microscopy (STORM) [START_REF] Michael | Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (storm)[END_REF]-bypass the diffraction limit so as to reach unprecedented nanoscale resolution. The main principle behind these methods relies on a combined use of optics and numerical processing, which is commonly called computational imaging. The resolution improvement is thus directly related to the performance of the reconstruction algorithms employed to process the acquired data.

SMLM techniques use photoactivables fluorescent probes to sequentially image a subset of activated molecules. Then, dedicated algorithms are deployed to precisely extract the position of these molecules. While the difficulty of the localization problem increases with the density of activated molecules per acquisitions, low density activations drastically reduce the temporal resolution of the system which makes the method limited for live imaging. Hence, current trends in SMLM concern the development of efficient algorithms dealing with high density data for which classical point-spread function (PSF) fitting or centroid localization methods [START_REF] Henriques | Quickpalm: 3d real-time photoactivation nanoscopy image processing in imagej[END_REF] fail. In particular, off-the-grid sparse regularized methods have shown their efficiency for high density settings [START_REF] Huang | Super-resolution image reconstruction for high-density three-dimensional single-molecule microscopy[END_REF][START_REF] Boyd | The alternating descent conditional gradient method for sparse inverse problems[END_REF]. For a complete review and comparisons of existing methods, we refer the reader to the two recent SMLM challenges [START_REF] Sage | Quantitative evaluation of software packages for single-molecule localization microscopy[END_REF][START_REF] Sage | Super-resolution fight club: A broad assessment of 2d & 3d single-molecule localization microscopy software[END_REF].

Initially introduced for two-dimensional imaging, SMLM has been extended to 3D thanks to PSF engineering. The principle relies on the design of PSFs which vary in the axial direction (i.e. z) in order to encode an information about the depth of molecules. Conventional PSF models include astigmatism [START_REF] Huang | Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy[END_REF] and double-helix [START_REF] Rama | Three-dimensional, single-molecule fluorescence imaging beyond the diffraction limit by using a double-helix point spread function[END_REF]. An alternative to PSF engineering is to record simultaneously multiple focal planes, as in the biplane modality [START_REF] Juette | Three-dimensional sub-100 nm resolution fluorescence microscopy of thick samples[END_REF]. It is noteworthy that these two approaches can also be combined as in [START_REF] Huang | 3d multifocus astigmatism and compressed sensing (3d macs) based superresolution reconstruction[END_REF] where the authors use both an astigmatism PSF and multi-focal acquisitions.

In this section, we study the performance of the SFW algorithm on both astigmatism and double-helix modalities with various number of focal planes (typically from 1 to 4). We emphasize that conventional astigmatism and double-helix SMLM devices-in particular commercial ones-use a single focal plane. As opposed to singlefocal acquisitions, multi-focal acquisitions require to mount and synchronize several cameras in parallel. To the best of our knowledge, such a setting has only been reported by Huang et al [START_REF] Huang | 3d multifocus astigmatism and compressed sensing (3d macs) based superresolution reconstruction[END_REF] for the astigmatism SMLM. Moreover, we propose to compare these two modalities to an alternative approach where depth information is extracted from multi-angle total internal reflection fluorescence (MA-TIRF) microscopy acquisitions. Such an approach has never been reported yet and we expect our numerical simulations to serve as a proof of concept for further developments. One of the main interest in combining SMLM with MA-TIRF is that classical PSFs, which are better localized laterally than astigmatism or double-helix, can be used. This would reduce the difficulty of lateral molecule localization for high density settings while recovering the depth through the MA-TIRF acquisitions.

Forward Operators

In this section, we define the forward operator Φ for the three modalities considered in this paper. The first two correspond to conventional three-dimensional SMLM with astigmatism or double-helix PSFs. The third one, on the contrary, uses a MA-TIRF excitation in order to get an information about the depth of molecules. The operator Φ : M(X) → R N1N2K maps the Radon measures m ∈ M(X) to the discrete noiseless measurements Φm

∈ R N1N2K , Φm = X ϕ(x)dm(x). (44)
It is fully characterized by the function ϕ : X → R N1N2K . Hence, for each modality, we only have to define ϕ. In the following, X def.

= [0, b 1 ] × [0, b 2 ] × [0, b 3
] is a subset of R 3 , and we write x = (x 1 , x 2 , x 3 ) ∈ X. Then, we consider a camera containing N 1 × N 2 pixels and we denote the center of the ith pixel by (c i,1 , c i,2 ). Finally, we provide expressions of ϕ which enclose the integration over camera pixels

Ω i def. = (c i,1 , c i,2 ) + - b 1 2N 1 , b 1 2N 1 × - b 2 2N 2 , b 2 2N 2 ⊂ Ω def. = [0, b 1 ] × [0, b 2 ].
Astigmatism model. This modality provides depth information using an astigmatism deformation of the PSF with respect to the axial direction z. It is customary to model the latter with a Gaussian function whose variances σ 1 and σ 2 vary with z according to [START_REF] Huang | Super-resolution image reconstruction for high-density three-dimensional single-molecule microscopy[END_REF][START_REF] Kirshner | Can localization microscopy benefit from approximation theory?[END_REF] σ 1 (z)

def. = σ 0 1 + αz -β d 2 and σ 2 (z) def. = σ 1 (-z). ( 45
)
The constants involved in [START_REF] Harchaoui | Conditional gradient algorithms for norm-regularized smooth convex optimization[END_REF] can be calibrated from real data [START_REF] Huang | Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy[END_REF][START_REF] Kirshner | Can localization microscopy benefit from approximation theory?[END_REF]. Then, integrating this Gaussian model over camera pixels, we have for all i ∈ {1, . . . , N 1 N 2 } and k ∈ {1, . . . , K}

[ϕ(x)] i,k def. = 1 2πσ 1 (x 3 -z k )σ 2 (x 3 -z k ) Ωi e - (x 1 -s 1 ) 2 2σ 2 1 (x 3 -z k ) + (x 2 -s 2 ) 2 2σ 2 2 (x 3 -z k ) ds 1 ds 2 ,
where (z k ) K k=1 are the positions of the considered focal planes.

Double-helix model. Here, depth information is obtained by using a PSF formed out of two lobes which coil around each other along z to form a double-helix shape. In this paper, we model these lobes by two Gaussian functions with fixed variances σ 1 = σ 2 , and with a center whose lateral position (r 1 , r 2 ) (respectively, (-r 1 , -r 2 )) varies with z according to r 1 (z)

def.

= ω 2 cos(θ(z)) and r 2 (z)

def.

= -ω 2 sin(θ(z)) where θ(z) = θ speed z.

Parameters ω > 0 and θ speed > 0 correspond to the distance between the two Gaussian and the rotation speed of the double-helix (rad/nm), respectively. Then, integrating this model over camera pixels, we have for all i ∈ {1, . . . , N 1 N 2 } and k ∈ {1, . . . , K}

[ϕ(x)] i,k def. = 1 2πσ 1 σ 2 u∈{-1,1} Ωi e - (x 1 +ur 1 (x 3 -z k )-s 1 ) 2 2σ 2 1 + (x 2 +ur 2 (x 3 -z k )-s 2 ) 2 2σ 2 2 ds 1 ds 2 ,
where (z k ) K k=1 are the positions of the considered focal planes.

MA-TIRF model. With this modality, each activated set of molecules is imaged using K ∈ N TIRF illuminations with incident angles (α k ) K k=1 . Let n i > 0 and n t > 0 be the refractive indices of the incident (i.e. glass coverslip) and the transmitted (i.e. sample) medium, respectively. A TIRF excitation is obtained when the incident angle α is greater than the critical angle α c = arcsin(n t /n i ) for which we have total internal reflection of the light within the incident medium. This phenomenon produces an evanescent wave which decays in the transmitted medium as exp(-sx 3 ), where s = (4πn i )/λ sin 2 (α) -sin 2 (α c ) is the penetration depth and λ is the wavelength of the incident laser beam [START_REF] Axelrod | Cell-substrate contacts illuminated by total internal reflection fluorescence[END_REF][START_REF] Axelrod | Total internal reflection fluorescence microscopy[END_REF]. Because the decay of this evanescent excitation vary with the incident angle, the depth of biological structures can be recovered with a nanometric precision from multi-angle acquisitions [START_REF] Boulanger | Fast high-resolution 3d total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF][START_REF] Cardoso | Topography of cells revealed by variable-angle total internal reflection fluorescence microscopy[END_REF][START_REF] Zheng | Three-dimensional super-resolved live cell imaging through polarized multi-angle tirf[END_REF]. Combining this principle with SMLM techniques lead to a forward model Φ defined, for all i ∈ {1, . . . , N 1 N 2 } and k ∈ {1, . . . , K}, by

[ϕ(x)] i,k def. = ξ(x 3 )e -s k x3 2πσ 1 σ 2 Ωi e - (x 1 -s 1 ) 2 2σ 2 1 + (x 2 -s 2 ) 2 2σ 2 2 ds 1 ds 2 , (47)
where ξ(z)

= K k=1 e -2s k z -1/2
. This model comes from the combination of a lateral convolution with the axial TIRF excitation. Here the PSF of the system is assumed to be a Gaussian with variances σ 1 = σ 2 , and to be constant along x 3 (because only a thin layer of few hundred nanometers is excited by the evanescent wave). The values (s k ) K k=1 correspond to the penetration depths associated to the incident angles (α k ) K k=1 . Remark 9. One particularity of the MA-TIRF modality is that the kernel ϕ in [START_REF] Henriques | Quickpalm: 3d real-time photoactivation nanoscopy image processing in imagej[END_REF] is separable. This can be exploited numerically to reduce the overall algorithm complexity.

Illustrations and numerical computation of η V . Examples of noiseless measurements y 0 = Φm a0,x0 with m a0,x0 = δ (1.5,2.5,0.1) + δ (1.5,3,0.5) + δ (2,5,0.7) + δ (4.5,3.5,0.4) + δ (5,1,0.2) .

are presented in Figure 8 for the three modalities. The parameters used for these simulations are provided in Table 1. One can observe the effect of the three modalities on molecules at different depths. For the astigmatism modality, the orientation along which the PSF is defocuced indicates the position of the molecule with respect to the focal plane (above/below). Moreover, the larger is this defocucing, the deeper is the molecule. In the case of the double-helix modality, we can clearly see the rotation of the PSF with depth. Finally, for the MA-TIRF modality, we can observe that the recorded intensities for deep molecules decrease, with the incident angle, faster than the intensity for molecules which are close to the glass coverslip (i.e. x 3 = 0). Although, for these three-dimensional models, an explicit expression of η V seems challenging to come by, the latter can be computed numerically for specific points x ∈ X. A representation of η V for the measure given in (48) at x 3 = 0.1 and x 3 = 0.5 is depicted in Figure 9. For the three modalities, we have that η V (1.5, 2.5, 0.1) = η V (1.5, 3, 0.5) = 1 and otherwise η V is smaller than 1. Hence, η V seems nondegenerate and a measure composed of the same number of Dirac masses as m a0,x0 can be recovered by the SFW algorithm.

Simulation setting

Imaged Structure. Simulations were performed using the microtubules-like structure depicted in Figure 10. It has been generated within the volume x 1 (µm)

X = [0, b 1 ] × [0, b 2 ] × [0, b 3 ] ⊂ R
Figure 8: Noiseless acquisitions y0 for the measure ma 0 ,x 0 given in ( 48) and K = 4. The parameters used for these simulations are given in Table 1. The color of the molecules represent their depths: 0 (red) -0.8µm (blue).

small steps). Then, in order to give a width to the filaments, each point x ∈ X randomly chosen on one of the curves is replaced by a point randomly chosen in a ball centered at x with radius 10 nm. Thus, simulated filaments have a diameter of 20 nm.

Simulation of noiseless acquisitions. The N tot ∈ N * molecules of the simulated structure are divided into n ∈ N * sparse set of N ∈ N * molecules using a random permutation (i.e. N tot = n × N ). This models the sequential stochastic activation of fluorophores used in SMLM. For each of the n subsets of molecules, we define a Radon measure composed of a sum of Dirac masses-located at the position of the molecules-with positive amplitudes

m a0,x0 = N i=1
a 0,i δ x0,i where a 0,i > 0 and x 0,i ∈ X.

The amplitudes are randomly generated within [1, 1.5]. An example of a set of activated molecules is shown in Figure 10 (black crosses). Now let (N 1 × N 2 ) be the size of the grid of pixels on the detector plane, and K be the number of focal planes (or the number of TIRF "angles", see Section 5.1) which are recorded. Then, the noiseless measurements y 0 for an activated measure m a0,x0 follow the model

y 0 = Φm a0,x0 , (50)
where Φ is defined in (44). Finally, it is noteworthy that in practice the number of activated molecules varies from one activation to another around an average value (which depends on the power of the excitation laser beam). However, fixing this number to N for each activated set of molecules allows us to better control the density of spikes in order to study the behaviour of the algorithm when the latter increases.

Noise Model. There are two predominant sources of noise in microscopy data.

• The shot noise which is inherent to the quantum nature of light (random emissions of photons). It is well modeled by a Poisson distribution whose intensity is the number of photon collected at each pixel. Given the noiseless acquisition y 0 , we normalize it such that max i∈{1,...,N1N2} K k=1

[y 0 ] i,k = n photon , (51)
where n photon > 0 denotes the maximal photon budget per pixel and controls the noise level. Then, each entry of y 0 is replaced by a realization of a Poisson distribution P with parameter [y 0 ] i,k . It is noteworthy from ( 51) that the level of noise not only increases as n photon decreases, but it also increases with K. • The readout noise w G of the camera. It is usually modeled by a Gaussian distribution with variance σ 2 . Finally the noisy data are given by y = P(y 0 ) + w G .

(52)

Results

For each of the three modalities presented in Section 5.1 (Double-Helix, Astigmatism, MA-TIRF), acquisitions where simulated using the optical parameters gathered in Table 1. These parameters have been tuned according to the experimental PSF used in the SMLM challenge [START_REF] Sage | Super-resolution fight club: A broad assessment of 2d & 3d single-molecule localization microscopy software[END_REF]. Finally, we generated different experiments by varying the density of molecules N ∈ {5, 10, 15} as well as the number of focal planes (or angles for the TIRF model) K ∈ {1, 2, 3, 4}.

Metrics for evaluation

In order to assess the quality of the reconstructed volumes, we consider standard metrics which reflect both the detection rate and the localization error [START_REF] Sage | Super-resolution fight club: A broad assessment of 2d & 3d single-molecule localization microscopy software[END_REF][START_REF] Sage | Quantitative evaluation of software packages for single-molecule localization microscopy[END_REF]. Given a recovered frame and a tolerance radius r > 0, we pair estimated molecules and ground truth (GT) molecules when the distance between them is lower than r. Paired estimated molecules are then referred as true positive (TP) while unpaired ones as false positive (FP). Finally, the unpaired GT molecules are identified as false negative (FN). These quantities being determined for each frame, we can compute the Jaccard index (Jac), the Recall (Rec) and the Precision (Pre) metrics,

Jac = #TP #TP + #FP + #FN Rec = #TP #TP + #FN Pre = #TP #TP + #FP . (53)
The Jaccard index measures the overall performance of detection by giving a measure of similarity between the two sets of points. The Recall and Precision metrics can then be used to measure the ability of an algorithm to minimize FN and FP detection, respectively. Finally, the TP molecules are used to compute the root mean squared error (RMSE) along each dimension

RMSE x1 = 1 #TP i∈TP ([x i ] 1 -[x 0,i ] 1 ) 2 , (54)
and similarly for RMSE x2 and RMSE x3 . Note that, by construction, the RMSE is bounded by the radius r. Hence, in the following, we use different values for r depending on the metric of interest. α c + αmax-αc K-1 (k -1) Incident angles α max sin -1 (NA/n i ) Maximal incident angle 5.3.2. Choice of the regularization parameter λ For each experiment (i.e. N ∈ {5, 10, 15} and K ∈ {1, 2, 3, 4}), we choose the value of the regularization parameter λ which maximizes the Jaccard index for a radius of r = 0.02 (i.e. 20 nm). This training step was performed over a small subset of initial measures m a0,x0 (i.e. frames). Then the recovery was done on the complete dataset using the optimal λ found.

Discussion

The evolution of Jaccard, Recall, and Precision metrics with respect to K are depicted in Figure 11. As expected, they all increase with K. However, although the improvement is significant from K = 1 to K = 2, higher values only provide marginal gains. This can be explained by the fact that the photon budget n photon is distributed over the K acquisitions (see equation ( 51)). Hence, the additional axial information brought by increasing the number of acquisitions per activation should be balanced by the higher noise corrupting the data. Another observation from these plots concerns the degradation of the performance as the density (i.e. the number of molecules N ) increases. These results also bring useful information in order to improve existing systems. Let us recall that current commercial systems includes Astigmatism and Double-Helix modalities with one focal plane (i.e. K = 1). Hence, it can be inferred from our simulations that recording an image at two focal planes for each activation of molecules would not only improve significantly the reconstruction quality but make the reconstructions more robust when the density of molecules increases. These observations corroborate the study in [START_REF] Huang | 3d multifocus astigmatism and compressed sensing (3d macs) based superresolution reconstruction[END_REF] where the authors use a multi-focus astigmatism system. However, to preserve a reasonable temporal resolution, multifocal acquisitions require to synchronize several cameras [START_REF] Huang | 3d multifocus astigmatism and compressed sensing (3d macs) based superresolution reconstruction[END_REF] which can be expensive and lead to delicate calibration procedures (e.g. alignment and PSF aberrations for each camera). In that respect, the proposed combination of SMLM with MA-TIRF offers an interesting alternative to improve existing systems. First, it has the potential to provide reconstructions whose quality compares favorably with the Double-Helix model while improving over the Astigmatism modality. Second, it only requires the use of galvanometric mirrors to control the incident angle [START_REF] Boulanger | Fast high-resolution 3d total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF]. It is noteworthy that commercial SMLM systems generally use a single TIRF illumination to limit the illumination depth. Finally, as for the multi-focus strategy, MA-TIRF requires some calibrations (e.g. incident angles) for which there exist dedicated procedures [START_REF] Boulanger | Fast high-resolution 3d total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF][START_REF] Soubies | A framework for multi-angle tirf microscope calibration[END_REF]. Remark 10. Although the PSFs used for these simulations have been adjusted using experimental PSFs, they remain idealistic. This is particularly the case for the Double-Helix which in practice deviates from two Gaussian lobes that coil around each other along z [START_REF] Sage | Super-resolution fight club: A broad assessment of 2d & 3d single-molecule localization microscopy software[END_REF]. In contrast, the Gaussian model yields a precise approximation of the MA-TIRF (i.e. widefield) PSF [START_REF] Zhang | Gaussian approximations of fluorescence microscope point-spread function models[END_REF]. The main simplification for the latter lies in the fact that each molecule is activated only during one set of multi-angle acquisitions. This would not be the case with a real implementation of the system and the model should be improved by considering the temporal aspect of the acquisition. However, the present study constitutes a first proof-of-concept and future developments will consider a more sophisticated model.

The results in terms of RMSE presented in Figure 12 lead to similar interpretations. First, the detection accuracy is increasing with K while decreasing with N . Second, we can observe that the differences between the Double-Helix and the MA-TIRF models mainly come from the precision in x 3 . Indeed, they both lead to the same lateral RMSE (around 5nm when N = 5 and 12nm at the highest density N = 15), but the Double-Helix enjoys a better axial RMSE. This reflects the challenging problem that constitues the inversion of the Laplace transform, which is related to the MA-TIRF model. Nevertheless, the SWF algorithm performs quite well at this task (see also Figures 13 and 14). Another observation concerns the fact that the Double-Helix can reach a better axial than lateral RMSE. This fact, which was also observed in the recent SMLM challenge [START_REF] Sage | Super-resolution fight club: A broad assessment of 2d & 3d single-molecule localization microscopy software[END_REF], can be explained by the large lateral support of the Double-Helix PSF as well as its good axial discrimination.

Finally, three-dimensional representations of the recovered structures are presented in Figures 13 and 14 for a fixed K = 4 and N = 10, respectively. These figures complete and illustrate the observations made with the computed metrics.

Conclusion

This paper demonstrated from both theoretical and practical perspectives the Sliding Frank-Wolfe Algorithm, in particular when facing challenging non-translation invariant operator such as the Laplace kernels. Such operators lead to difficulties in estimating the spikes positions which is efficiently addressed by non-convex update step of the grid location. The BLASSO method, coupled with this Sliding Frank-Wolfe solver, is well adapted to these non-convolutive operators because it does not rely on spectral (Fourier) methods and can be analyzed theoretically through the prism of convex duality and vanishing certificates.

MA-TIRF

Astigmatism Double-Helix

K = 1 K = 2 K = 3

•X,

  or T d with d ∈ N * . Definition 1 covers the particular case of X = R d , X = T d or any compact subset with non-empty interior of R d .

N = 7 Figure 1 :-(•-x) 2 2σ 2 )

 7122 Figure 1: ηW for a Gaussian convolution (x ∈ R, ϕ(x) = e -(•-x) 2 2σ 2 ) for several numbers of spikes and σ = 1.

Lemma 1 .

 1 Let I, I ⊂ R be open intervals, and h : I → I be a smooth diffeomorphism. Let x c ∈ I, t c := h -1 (x c ) ∈ I , and let η : I → R be a smooth function. Then η satisfies η

Lemma 2 .

 2 Let I be an open interval, t c ∈ I and let g : I → R, η : I → R be two smooth functions. If η satisfies:

3. 3 .

 3 Explicit Formulas for η W in Continuous Settings 3.3.1. Classical Laplace Operator We suppose that µ = L, where L is the Lebesgue measure on R + , and ξ = 1. Then one has

  xc = 1

Figure 2 :

 2 Figure 2: ηW for the unnormalized Laplace model for a varying xc with fixed N = 2 and a fixed xc = 1 with varying N ∈ {2, 4, 6}.

Figure 3 :

 3 Figure 3: ηW for the normalized Laplace model for a varying xc with fixed N = 2 and a fixed xc = 1 with varying N ∈ {2, 4, 6}.

  )s dµ n (s), andc ϕn (x, x ) = ξ n (x)ξ n (x ) R+ e -(x+x)s dµ n (s) where ξ n (x) = R+ e -2xs dµ n (s), and similarly for c ϕ and µ = L. Proposition 4. Let (µ n ) n∈N be a sequence of positive measures which converges towards the Lebesgue measure µ in the local weak-* topology, i.e. ∀ψ ∈ C c (R + ), lim n→+∞ R+ ψ(s)dµ n (s) = R+ ψ(s)ds, and such that sup n∈N R+ (1 + s 4N -1 )e -xmins dµ n (s) < +∞. (

  ds, and b e -xminA R+ (1 + s 4N -1 )e -xmins dµ n (s) e -xminA C.

Lemma 3 .

 3 Let X and S be two compact metric spaces, and ψ ∈ C (X × S). If {µ n } n∈N and µ are Radon measures such that µ n * µ in the weak-* convergence of M(S), then lim n→+∞ S ψ(x, s)dµ n (s) = S ψ(x, s)dµ(s), uniformly in x ∈ X.

Figure 4 :

 4 Figure 4: Approximation of ηW for the unnormalized continuous Laplace operator (see Proposition 2) by the ηW obtained for discretized unnormalized Laplace operators.

Algorithm 2 1 :

 21 Sliding Frank-Wolfe Algorithm Initialize with m [0] = 0 and n = 0. 2: for k = 0, . . . , n do 3:

Figure 5 :

 5 Figure 5: ηV for ma 0 ,x 0 = 1.3δ0.3 + 0.8δ0.37 + 1.4δ0.7.

3 TFigure 6 :

 36 Figure 6: Values of the objective function throughout the SFW algorithm (cumulative iterations of the BFGS). The vertical black lines separate the main outer iterations of the algorithm.

k

  = 0. Start of the loop.

1 k = 1 .

 11 End of the loop.

1 k = 2 .

 12 End of the loop.

Figure 7 :

 7 Figure 7: Main steps of the SFW algorithm.

3

 3 where b 1 = b 2 = 6.4 µm and b 3 = 0.8 µm. (49)The filaments were obtained by randomly sampling many points along four curves defined by polynomial equations. To ensure a uniform distribution of the points along the curves, we first parametrized each curve by a piecewise linear function (with very 0

Figure 9 : 1 Figure 10 :

 9110 Figure 9: Numerical computation of ηV at x3 = 0.1 (top) and x3 = 0.5 (bottom) for the three models and the measure ma 0 ,x 0 given in (48). The colormap ranges from 0 (blue) to 1

Figure 11 :

 11 Figure 11: Evolution of Jaccard, Recall and Precision metrics with respect to K, for a radius of detection r = 0.02 (20nm).

Figure 12 :

 12 Figure 12: Evolution of the RMSE (nm) with respect to K, for a radius of detection r = 0.1 (100nm).

Figure 14 :

 14 Figure 14: Recovered structures for N = 10.

  We use the ∞ norm, |•| ∞ , for vectors of R N or R 2N , whereas the notation • refers to an operator norm (on matrices, or bounded linear operators). • H is the norm on H associated to the inner product •, • H . • ∞,X denotes the L ∞ norm for functions defined on X.

	where y	def.
		and only if
		ϕ ∈ KER (k) and (ϕ 0 , . . . , ϕ k ) are linearly independent in H.	(I k )
	Norms. 2. Reminders on the BLASSO
	2.1. Recovery of the Support in Presence of Noise
	Let x 0 ∈	• X N , a 0 ∈ (R \ {0}) N and m a0,x0 =	N i=1 a 0,i δ x0,i . The BLASSO is the
	variational problem
			min m∈M(X)	1 2	Φm -y	2 H + λ|m|(X),	(P λ (y))

Table 1 :

 1 Parameters used for data simulation.

		Parameter	Value	Description
		b 1 = b 2	6.4µm	Region of interest
		b 3	0.8µm	Maximal depth of molecules
	All modalities	N 1 = N 2 NA n i n t λ	64 1.49 1.515 1.333 0.66µm	Detector grid size Objective numerical aperture Refractive index incident medium Refractive index transmitted medium Excitation wavelength
		n photon	1000	Photon budget
		σ	10 -4	Variance of Gaussian noise
	Astigmatism	σ 0 β d α (z k ) K k=1	0.42λ /NA 0.2µm λ n i /(2NA 2 ) -0.79 kb 3 /(K + 1)	PSF variance at focus Depth for which the variance is minimal Parameter related to the depth-of-field Scaling constant Focal planes
	Double-Helix	σ 1 = σ 2 ω θ speed (z k ) K k=1	0.42λ /NA 1µm 0.3846π rad/µm kb 3 /(K + 1)	PSF variance Distance between the two PSF lobes Rotation speed of the PSF Focal planes
	MA-TIRF	σ 1 = σ 2 (α k ) K k=1	0.42λ /NA	PSF variance

Acknowledgement

The authors would like to thank Laure Blanc-Féraud for initiating this collaboration and for stimulating discussions. The work of Gabriel Peyré has been supported by the European Research Council (ERC project NORIA). The work of

Emmanuel Soubies has been supported by the European Research Council (ERC project GlobalBioIm).