Valentin Debarnot

Pierre Weiss

Blind inverse problems with isolated spikes

Assume that an unknown integral operator living in some known subspace is observed indirectly, by evaluating its action on a discrete measure containing a few isolated Dirac masses at an unknown location. Is this information enough to recover the impulse response location and the operator with a sub-pixel accuracy? We study this question and bring to light key geometrical quantities for exact and stable recovery. We also propose an in depth study of the presence of additive white Gaussian noise. We illustrate the well-foundedness of this theory on the challenging optical imaging problem of blind deconvolution and blind deblurring with unstationary operators.

To motivate this paper, let us start with a concrete problem in imaging. In Figure 1a, we simulated an image of fluorescent proteins observed with an optical microscope. Assume that an algorithm is able to recover the proteins locations at a sub-pixel accuracy from this image. By taking thousands of such images and stacking the protein locations, it is possible to break the diffraction limit and to construct an image with a resolution of the order of a nanometer. This principle was awarded the 2014 Nobel prize in chemistry [START_REF] Betzig | Imaging intracellular fluorescent proteins at nanometer resolution[END_REF][START_REF] William | Optical detection and spectroscopy of single molecules in a solid[END_REF].

From a mathematical viewpoint, this problem can be modelled as follows. Let M(R D ) denote the set of Radon measures, i.e. the dual of the set C 0 0 (R D ) of continuous functions vanishing at infinity. Let μ = N n=1 wn δ xn ∈ M(R D ) denote a Radon measure that encodes the protein locations ( xn ) and their intensities ( wn ). Assume that this measure is observed indirectly through a linear regularizing operator Ā : M(R D ) → C 0 0 (R D ):

y m = Āμ(z m ) + b m , (1)
where y ∈ R M is the observed data, Z = (z 1 , z 2 , . . . , z M ) denotes a set of sampling locations in R D and b = (b 1 , . . . , b M ) ∈ R M is some additive noise. The operator Ā is typically a convolution operator with some impulse response (e.g. an Airy pattern in Fourier optics) or a more complicated space-varying operator such as the one in Figure 1c. Numerous approaches have been developped over the years to recover the positions ( xn ) from the measurements y. We refer the interested reader to the summaries of the super-resolution challenges [START_REF] Sage | Quantitative evaluation of software packages for single-molecule localization microscopy[END_REF][START_REF] Sage | Super-resolution fight club: assessment of 2d and 3d single-molecule localization microscopy software[END_REF] for more insight on the possible approaches. The main hurdles to solve this problem are the following:

a) The number of measurements M can be large, making it essential to design computationally efficient methods.

b) The weights ( wn ) are usually unknown.

c) It is important to work off-the-grid to avoid biases in the location estimation.

d) The proteins can sometimes be aggregated in clusters, resulting in a difficult disentanglement of their individual locations.

e) Most importantly for this paper: the operator Ā is often only partially known, making it crucial to estimate both the positions and weights ( wn , xn ), but also the operator Ā itself.

The main objective of this work is to design certified methods, which are able to cope with the above difficulties. All the results will be stated under the following two assumptions.

Assumption 1.1 (The operators' structure). We assume that the operator Ā lives in a known finite dimensional subspace A of linear operators from M(R D ) to C 0 0 (R D ).

Throughout the paper, we let I ∈ N denote the subspace dimension and set A = span (A 1 , . . . , A I ). Any A ∈ A can therefore be parameterized by a vector γ = (γ i ) ∈ R I , with

Aµ = A(γ)µ def. = I i=1 γ i A i µ for any µ ∈ M(R D ).
(

) 2
The next assumption describes the sampling model considered in this work.

Assumption 1.2 (The observation model). Let (ν m ) 1≤m≤M in M(R D ) denote a collection of M linear forms on C 0 0 (R D ). Let μ ∈ M(R D ) denote a signal to reconstruct. We assume that we acquire the measurement vector y = (y 1 , . . . , y M ) with coordinates given by y m def.

= ν m , A(γ)μ + b m .

(

) 3
The observation model 1.2 allows to describe nearly any sampling device. For instance, the traditional pointwise sampling would consist in choosing ν m = δ zm , where (z m ) 1≤m≤M is a set of sampling locations. Fourier sampling can be modelled by setting ν m (ω) = exp(-i ω, z m ).

Finally, throughout the paper, we will assume that the measure μ is discrete.

with some sources isolated from the others. The meaning of isolated will be made precise in Theorem 3.4.

Our contribution

Our main contribution in this work is to propose a simple estimation method that strongly relies on Assumptions 1.1, 1.2 and 1.3 a) or b). Despite these restrictions, the proposed framework offers significant advantages:

• We can work with near arbitrary subspaces of operators A, beyond convolutions.

• We work under a general linear sampling model with arbitrary linear forms (ν m ).

• The proposed theory doesn't require a grid.

• It is rather simple and leads to recovery conditions that can be checked in advance (for some of them) or a posteriori (for some others).

• We provide strong stability results for the recovery of individual Dirac masses positions from simple correlation algorithms. Our results hold for noise levels of the order of the norm of the measurements, see Theorem 3.3. In comparison, some recent stability theories when sensing multiple source points (e.g. [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF]Theorem 2]) only provide asymptotic results when the noise level vanishes.

• We also provide explicit separation bounds, allowing the identification of isolated source in Theorem 3.4.

• We further refine the stability bounds, if b = (b m ) m ∈ R M is the realization of white Gaussian noise. This requires analyzing the suprema of continuous Gaussian processes and second order chaos. The resulting theory explains why correlation algorithms can perform extremely well even under large noise levels, see Theorem 3.5.

• We show how these results allow to stably estimate an operator Ā ∈ A in Theorems 3.7 and 3.8.

• Importantly, the proposed algorithms are simple to implement and efficient in practice.

The main limitations lie in Assumptions 1.1 and 1.3 above. The prior knowledge of a subspace of operators A is rather standard in the literature. It is actually the basis of bilinear inverse problems and of a popular trick called lifting. This assumption is quite realistic for the field of optics. When dealing with convolution operators, it is possible to use a principal component analysis to design a low dimensional orthogonal basis allowing to approximate efficiently any sufficiently smooth family of impulse responses. For more general spacevarying operators, we recently proposed an efficient calibration technique in [START_REF] Debarnot | A scalable estimator of sets of integral operators[END_REF][START_REF] Debarnot | Learning low-dimensional models of microscopes[END_REF], leading to low-dimensional subspaces of unstationnary operators.

Assumption 1.3 is -in a sense -more delicate to respect in practice. It is sometimes possible to sense some sources individually, one after the other using calibration techniques. This is the case for instance in PALM microscopy (photo-activable localization microscopy). This would lead to the model given in Assumption 1.3 a). Another favorable situation is depicted in Fig. 1, corresponding to model in Assumption 1.3 b). There, the impulse responses are not compactly supported but they possess a rather fast polynomial decay. Their intensity is soon dominated by the noise. We will see in Theorem 3.4 that their cross-interactions can be neglected under reasonable conditions. In both cases a) and b) though, we neglected important practical effects arising in applications such as quantization, auto-fluorescence background,... Though these application specific issues might be hard to handle, the proposed study still provides useful insights and guidelines.

Finally, let us stress out that the separated source points are used only to provide an estimate  of the operator Ā. Once it is estimated, we can use existing inverse problem solvers to recover any measure µ ∈ M(R D ) from its measurements Āµ. In the particular case of discrete measures, the independence or separation condition in Assumption 1.3 is not needed anymore, therefore taking care of problem d) above. This allows to address problems such as super-resolution or deblurring problems with the proposed formalism.

Related works

When the operator Ā is known, recovering μ is already a challenging problem, since the inverse problem is ill-posed and infinite dimensional. Specifying prior assumptions on the signal μ to certify its approximate recovery is essential [START_REF] Scherzer | Variational methods in imaging[END_REF]. A few mathematical breakthroughs were achieved in the recent past.

Off-the-grid total variation minimization with a known operator In [START_REF] De | Exact reconstruction using beurling minimal extrapolation[END_REF][START_REF] Emmanuel | Towards a mathematical theory of super-resolution[END_REF][START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF], the authors proposed to recover the individual source points by solving a generalization of the basis pursuit to an infinite dimensional setting. They showed that the recovery is stable given that the spikes are sufficiently separated. In [START_REF] Denoyelle | Support recovery for sparse super-resolution of positive measures[END_REF], the authors showed that the separation is not needed, provided that the weights ( wn ) are positive. However, the stability to noise deteriorates really fast when the distance between two source points gets lower than the Rayleigh criterion. From a numerical perspective, the solution of this problem can be found rather efficiently using techniques of semi-infinite programming [START_REF] Bredies | Inverse problems in spaces of measures[END_REF][START_REF] Denoyelle | The sliding frank-wolfe algorithm and its application to super-resolution microscopy[END_REF][START_REF] Flinth | On the linear convergence rates of exchange and continuous methods for total variation minimization[END_REF]. This type of approach is currently amongst the best competitors when a high density of proteins is used [START_REF] Sage | Super-resolution fight club: assessment of 2d and 3d single-molecule localization microscopy software[END_REF].

Gridded lifting for an unknown operator Assume that the operator Ā is unknown but lives in a known finite dimensional subspace A. Also assume that the positions ( xn ) are known, but that the weights ( wn ) are unknown. Under these hypotheses, an elegant solution to recover Ā and μ was proposed by Ahmed et al in [START_REF] Ahmed | Blind deconvolution using convex programming[END_REF] based on a trick called lifting. This approach allows to tackle bilinear problems of the form inf A∈A,u∈U

1 2 Au -y 2 2 , (6)
where U is a finite dimensional subspace of signals, by transforming the bilinear problem into a linear one restricted to rank-1 matrices. This nonconvex constraint can then be relaxed to a convex one by using the nuclear norm. This approach can be guaranteed to stably estimate ( wn ) and Ā under rather stringent assumptions. The assumptions were relaxed in a series of works [START_REF] Ahmed | Blind deconvolution using convex programming[END_REF][START_REF] Ling | Self-calibration and biconvex compressive sensing[END_REF][START_REF] Chi | Guaranteed blind sparse spikes deconvolution via lifting and convex optimization[END_REF][START_REF] Jung | Blind demixing and deconvolution at near-optimal rate[END_REF][START_REF] Ahmed | Leveraging diversity and sparsity in blind deconvolution[END_REF]. One important achievement was to allow to handle sparsity constraints over a fixed grid instead of subspace constraints. This is particularly relevant for the considered setting.

Off-the-grid lifting In [START_REF] Chi | Guaranteed blind sparse spikes deconvolution via lifting and convex optimization[END_REF], Y. Chi showed that the lifting trick could also be used when the unknown positions ( xn ) live off-the-grid, D = 1 and the operators in A are convolution operators. The approach was then extended to the 2 dimensional setting for convolution operators in [START_REF] Suliman | Blind two-dimensional super-resolution and its performance guarantee[END_REF]. In [START_REF] Chen | Vectorized hankel lift: A convex approach for blind super-resolution of point sources[END_REF], an alternative formulation was proposed based on the Hankel lifting for convolution operators in 1D. This approach is elegant but is currently restricted to convolution operators, while it is important in many applications to consider space variant systems. In addition, we will see that a convex relaxation may not be the most efficient approach from a practical viewpoint in the numerical experiments.

Multiple observations

The model 1.3 a), involves multiple observations when N > 1, for which specific theories have been developed in the discrete setting. In particular, a lot of attention was drawn to the case of multichannel blind deconvolution [START_REF] Amari | Multichannel blind deconvolution and equalization using the natural gradient[END_REF][START_REF] Sroubek | Multichannel blind iterative image restoration[END_REF], where multiple signals are convolved with a single filter. This leads to observations y n of the form y n = e u n , where e is an unknown filter and (u n ) are unknown inputs. Significant theoretical and numerical progress on this issue was made recently in the case where the inputs u n are sparse [START_REF] Wang | Blind deconvolution from multiple sparse inputs[END_REF][START_REF] Ahmed | Leveraging diversity and sparsity in blind deconvolution[END_REF][START_REF] Li | Multichannel sparse blind deconvolution on the sphere[END_REF][START_REF] Shi | Manifold gradient descent solves multi-channel sparse blind deconvolution provably and efficiently[END_REF]. The proposed setting is simpler since we assume that each input signal is the measurement of a single source. The added value of our result lies in the extension to operators beyond convolution models, off-the-grid estimation, and stronger stability results.

Preliminaries

All the proofs in this paper are postponed to the appendix.

Notation

Throughout the paper, we'll use the following conventions. A bar on top of a symbol indicates that it corresponds to the ground-truth (e.g. μ, wn , xn ). A hat indicates that it is an estimate (e.g. μ, ŵn , xn ). The symbols D, I, J, K, M, N denote cardinalities in N while d, i, j, k, m, n denote the associated indices. Bold fonts are used to design vectors and matrices and calligraphic fonts are used to design sets. For a vector x ∈ R D , x d denotes the d-th coefficient of x. The parentheses are used to evaluate functions. For instance, given f : R D → R, f (x) is the value of f at x. For µ ∈ M(R D ) and u ∈ C 0 0 (R D ), we let µ, u ∈ R denote the value of the linear form µ on u. We also let µ u denote the convolution product between µ and u defined for all

x ∈ R D by (µ u)(x) = µ, u(x -•) .
In all the paper, the notation •, • also refers to the usual scalar product on the vector space R N , where N ∈ N. For u ∈ R N , u 2 denotes the 2 -norm of u defined by u We let •, • L 2 (R D ) denote the usual scalar product of L 2 (R D ). For a compact and symmetric set Ω ⊂ R D , we let PW(Ω) denote the Paley-Wiener set of bandlimited functions on Ω, i.e. the set of functions in L 2 (R D ) that have a Fourier transform that vanishes outside Ω.

Further assumptions

Under Assumptions 1.1 and 1.2, the impulse response of an operator A(γ) at a location z ∈ R D is given by A(γ)δ z = I i=1 γ i A i δ z . This motivates introducing the matrix-valued function E : R D → R M ×I defined by

(E(z)) m,i def. = ν m , A i δ z . (7)
It will play an essential role in our analysis. Some of our results will depend on two additional hypotheses.

Assumption 2.1 (Identifiability of the operator from a single source). For all x ∈ D, the mapping E(x) : R I → R M is injective: we have

σ -Id E * (x)E(x) σ + Id (8) with 0 < σ -≤ σ + < +∞.
In what follows, we let κ def.

= σ+ σ-.

This assumption will be useful to guarantee that an operator can be stably estimated from a single source, once the location of a Dirac mass is known. It is not needed anymore if multiple sources are observed. Throughout the paper, we let R(x)

def.

= Ran(E(x)) [START_REF] De | Exact reconstruction using beurling minimal extrapolation[END_REF] denote the subspace of possible measurements for an impulse response located at x ∈ D and Π R(x) denote the orthogonal projection onto the range R(x).

Another important technical assumption is the following.

Assumption 2.2 (Identifiability of the Dirac masses location). The mapping E satisfies the following inequality for any pair x,

x ∈ D Π R(x) Π R( x) 2→2 ≤ 1 -φ( x -x 2 ) ( 10
)
for some non-decreasing function φ : R + → [0, 1] with φ(0) = 0 and φ(t) > 0 for t > 0.

This assumption allows to guarantee the stable recovery of the Dirac masses locations. This can be understood informally as follows. Take two locations x = x in D. Then, the two ranges R(x) and R( x) do not contain two identical elements. Hence, the knowledge of a measurement of the form E * ( x)γ should be enough to perfectly recover x. An additional implicit assumption is that the range R(x) cannot be the singleton {0} for any x. Indeed, taking x = x implies that Π R(x) 2→2 = 1.

Some intuition on Assumptions 2.1 and 2.2

Before stating our main results, we provide some intuition on the meaning of Assumption 2.1 and Assumption 2.2 and illustrate them through two examples.

An injectivity condition

Proposition 2.1. Under Assumptions 2.1 and 2.2, the mapping (x, γ) → E(x)γ is injective on (D × R I \{0}).

The injectivity of the mapping is a necessary condition to guarantee the identifiability of a position and an operator from a single measurement. For instance, it implies that -for any x -the subspace span(A i δ x , 1 ≤ i ≤ I) does not contain two elements that are shifted versions of each other. This hypothesis is essential to discard the standard ambiguity in blind deconvolution related to the fact that the signal and the convolution kernel can be shifted in opposite directions and still yield the same measurement vector, see e.g. [START_REF] Li | Identifiability and stability in blind deconvolution under minimal assumptions[END_REF].

A correlation condition Assumption 2.2 allows to control the correlation between measurements of an impulse response at x with an operator A(γ) and another at x with an operator A(γ). Indeed, we obtain using Cauchy-Schwarz inequality:

E(x)γ, E( x)γ = Π R(x) E(x)γ, Π R( x) E( x)γ = Π R( x) Π R(x) E(x)γ, E( x)γ ≤ Π R( x) Π R(x) 2→2 E(x)γ 2 E( x)γ 2 ≤ [1 -φ( x -x 2 )] E(x)γ 2 E( x)γ 2
A geometric condition The quantity Π R(x) Π R( x) 2→2 is related to the principal angle between the subspaces R(x) and R( x). To realize this, let us recall that the principal angle ∠ (U, V) between two subspaces U and V of a Hilbert space with norm • is defined by

cos (∠ (U, V)) = max u∈U ,v∈V u =0,v =0 u, v u • v . ( 11
)
We have:

Π R(x) Π R( x) 2→2 = sup u,v∈R M u 2= v 2=1 Π R(x) Π R( x) u, v = sup u∈R(x),v∈R(x) u 2= v 2=1 u, v = cos (∠ (R(x), R( x))) .

The case of convolution operators

In this paragraph, we aim at providing some insights on Assumptions 2.1 and 2.2 for the particular case of convolution operators. We work under the following assumption. Assumption 2.3. We assume that we are given an orthogonal1 family (e i ) 1≤i≤I of functions in PW(Ω) that vanish at infinity. The operators A i are convolutions with the filters e i , i.e. A i µ = e i µ for µ ∈ M(R D ).

The linear forms ν m describing the sampling device correspond to a Shannon sampler, i.e. ν m = δ zm , where the positions z m correspond to a Cartesian grid with a grid-size smaller than 2π diam(Ω) .

Under Assumption 2.3, any (u, v) ∈ PW(Ω) 2 satisfy

u, v L 2 (R D ) ∝ m∈N u(z m )v(z m ), (12)
which is a variant of the Shannon-Nyquist theorem, see e.g. [START_REF] Mallat | A wavelet tour of signal processing[END_REF]Thm 3.5].

Proposition 2.2 (Operator identifiability for convolution operators). Under Assumption (2.3), we have

E * (z)E(z) = Id, hence Assumption 2.1 is satisfied with σ -= σ + = 1.
Proposition 2.3 (Location identifiability for convolution operators). Let C : R D → R I×I denote the following cross-correlation matrix-valued function:

[C(x -x )] i,i def. = e i (• -x), e i (• -x ) L 2 (R D ) ( 13
)
Under Assumption (2.3), we have

Π R(x) Π R(x ) 2→2 = C(x -x ) 2→2 . (14)
Proposition 2.3 shows that the condition [START_REF] Emmanuel | Towards a mathematical theory of super-resolution[END_REF] characterizes the speed of decay of a cross-correlation matrix. For instance, consider the simplest case I = 1, corresponding to a convolution with a known filter e 1 . Then [START_REF] Emmanuel | Towards a mathematical theory of super-resolution[END_REF] simply measures how fast the auto-correlation function of e 1 decays away from 0. Intuitively, this information is central to derive stability results for algorithms that estimate the Dirac locations by finding correlation maxima. This statement is made precise in Theorem 3.3.

The case of product-convolution operators

To encode space varying operators, we now turn to product-convolution expansions [START_REF] Escande | Approximation of integral operators using product-convolution expansions[END_REF]. This decomposition allows to represent compactly most linear integral operators arising in applications. For the sake of the current paper, we work under the simplifying assumptions below. Assumption 2.4 (Product-convolution expansion). We assume that we are given an orthogonal family (e j ) 1≤j≤J of band-limited functions in PW(Ω), and another orthogonal family

(f k ) 1≤k≤K of functions in L 2 (R D ) ∩ C 0 0 (R D ). The family of observation operators A is a subspace of product-convolution expansions from M(R D ) to C 0 0 (R D
) defined as follows. For any A ∈ A, there exists a vector γ = (γ j,k ) ∈ R J×K such that for any µ ∈ M(R D ):

Aµ = A(γ)µ def. = J j=1 K k=1 γ j,k e j (f k µ). (15)
Similarly to Assumption 2.3, we assume that a Shannon sampler is used.

Letting i = (j, k), this implies that (E(z)) i,m = f k (z m )e j (z m -z).
Let us mention that the blurring operators appearing in optics can be represented very efficiently using this structure [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF]. In addition, we recently showed how a subspace of product-convolution operators A could be constructed in practice in optical imaging [START_REF] Bigot | Estimation of linear operators from scattered impulse responses[END_REF][START_REF] Debarnot | A scalable estimator of sets of integral operators[END_REF][START_REF] Debarnot | Learning low-dimensional models of microscopes[END_REF].

Proposition 2.4. Under Assumptions 1.1, 1.2 and 2.4, we have

Π R(x) Π R(x ) 2→2 = C(x -x ) 2→2 . (16)
However, for K ≥ 2, Assumption 2.1 is not valid: the mapping E(z) is not injective for any z.

As a consequence of this proposition, we will see that the identification of a product-convolution operator with K ≥ 2 is possible only under the condition N ≥ K, i.e. by observing multiple impulse responses.

Main results

Throughout this section, we will work under Assumptions 1.1, 1.2 and 1.3 a). For all 1 ≤ n ≤ N , the observation y n of a single Dirac mass at xn can be written as

y n = E( xn ) ᾱn + b n , (17)
where xn ∈ R D is an unknown location and ᾱn = wn γ ∈ R I is a vector co-linear to the unknown operator parameterization γ. The overall objective of this paper is to construct an estimate X = ( x1 , . . . , xN ) of X = ( x1 , . . . , xN ) and γ of γ and to certify their proximity, despite the perturbation term b n ∈ R M . The case of multiple sources in Model 1.3 b) will be treated as a particular case of [START_REF] Chi | Guaranteed blind sparse spikes deconvolution via lifting and convex optimization[END_REF], where b n coincides with the measurements of the N -1 sources different from n.

A simple two-step recovery algorithm

We aim at studying a simple recovery procedure.

Step 1: recovering the locations positions First we propose to estimate the positions xn by finding a global minimizer xn of the following problem

inf x∈D,α∈R I 1 2 E(x)α -y n 2 2 . (P x n )
Step 2: recovering the operator Second, depending on whether the weights w are known or unknown, we propose two different recovery strategies.

Case 1: known weights If w is known, we consider the following quadratic problem:

inf γ∈R I 1 2 N n=1 wn E( xn )γ -y n 2 2 , (P γ )
i.e. to find the operator's parameterization that provides the best fit to the observations y n , assuming that the Dirac mass locations are at xn .

Case 2: unknown weights If the true weights w are unknown, let

J(w, γ, X) def. = 1 2 N n=1 w n E( xn )γ -y n 2 2 .
We propose to solve the bilinear problem below inf

γ∈R I w∈R N J(w, γ, X). (P γ,w )
Additional regularization terms In the proposed formulation (P x n ), two implicit regularization terms are used: i) we look for a single Dirac mass located in D and ii) the operator lives in a known subspace. It is possible to add regularization terms to further constrain and stabilize the problem.

If the dimension I of the subspace of operators is large, stability issues should arise both for the estimation of the positions and of the operator. Indeed, multiple couples (x, γ) of positions and operators could lead to similar measurements. A possible solution to leverage this problem is to add a weighted 2 -regularization on γ of the form

1 2 I i=1 θ i γ 2 i
, where θ i are weights adapted to the problem at hand. Most of the theory developed in this paper could be extended to this setting as well. The main difference is that this regularization term introduces a bias in the operator estimate. This is why we prefer studying the unconstrained version above.

Second, notice that the problem (P γ,w ) suffers from the usual scaling ambiguity in bilinear inverse problem: if ( ŵ, γ) is a solution, then so is t ŵ, γ t for any t = 0. Without any further normalization, the weights w and parameterization γ can only be estimated up to a multiplicative factor. To avoid this problem a common solution is to constrain normalize the weight, e.g., with an affine constraint of the form γ, 1 = 1.

Estimating the locations

In this section, we establish a few results regarding the location estimate obtained by solving problem (P x n ). We first characterize and prove the existence of minimizers. We then study stability estimates, for generic noise, for white Gaussian noise and for perturbations with additional sources. These results constitute the most technical part of the paper.

A simple characterization of the locations

Let us start with an elementary observation. Proposition 3.1 (A simple characterization ). Assume that the problem (P x n ) admits a solution xn . Then

xn ∈ argmax x∈D 1 2 Π R(x) y n 2 2 . (18)
Notice that the existence of a minimizer is not automatic since the minimization domain D might be unbounded and the function x

→ 1 2 Π R(x) y n 2 2
might be discontinuous. In particular, singularities may appear if the dimension of the range R(x) varies over D. A sufficient condition for existence is given below.

Proposition 3.2 (A sufficient condition for existence ). Assume that

• The domain D is compact

• The mapping E is continuous over D and E(x) is of rank I for all x ∈ D.

Then Problem (P x n ) admits at least one solution xn .

Stability of the location estimates

In this section, we study the location estimation stability under various types of perturbations b n .

Generic perturbations

We start without assuming any specific struture to the noise term b n .

xn -xn 2 ≤ φ -1 + 2θ 2 + 4θ , (19)
where φ -1

+ (t) = inf {s s.t. φ(s) ≥ t} is the quantile function of φ (in particular φ -1 + = φ -1 if φ is bijective).
The above result is in sharp contrast with the asymptotic results available for multiple source points [START_REF] Emmanuel | Towards a mathematical theory of super-resolution[END_REF][START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF]. Indeed, we see that the noise level can be of the order of the signal's level and the bound [START_REF] Ahmed | Leveraging diversity and sparsity in blind deconvolution[END_REF] still yields a useful stability estimate. When looking at the proof, it actually gets clear that it cannot be improved for an arbitrary noise term b n , apart from the constant • Assumptions 2.2 is satisfied with

φ(t) = (t/a) b 1 + (t/a) b for some a > 0 and b > 0. (21)
• We have y n 2 ≤ c y 1 2 for all n ≥ 2 and some constant c > 0.

• The first source location x1 is isolated from the others:

x1 -xn 2 ≥ δ for n ≥ 2, δ ≥ 2a(τ (N -1)c) 1/b and τ ≥ 5.
Then any global maximizer of

H def. = 1 2 Π R(x) y 2 2 in the ball B δ/2 = {x, x - x1 2 ≤ δ/2} lies within the ball B r = {x, x -x1 2 ≤ r} with r ≤ 18.4 • a • 2 1/b τ (22)
The above theorem clarifies what is meant by "isolated". It depends on the speed of decay of the projection matrices, captured by the exponent b, a typical width a of the impulse responses and a relative amplitude of the Dirac mass to recover c. This result states that a source will be resolved with an accuracy a/τ if the distance of the surrounding molecules is higher than a(τ N c) 1/b . Notice that the function H may contain multiple local maxima. However, under a separation condition, the global maximizer within a sufficient large ball correspond to the actual location of the sources. Remark 3.1. Let us formulate a few additional remarks.

• The second condition is quite mild. It is for instance satisfied under Assumption 2.1 and if w1 ≥ c wn for some constant c > 0 and all n ≥ 2.

• Asymptotically in N , the factor 18.4 in (22) can be replaced by 1. The following bound emerges from the proof:

r ≤ a • 2 1/b τ • (τ 2 + τ )(N -1)c + 2τ 1 + (τ 2 -4τ -2)(N -1)c + (2τ -4)(N -1)c
.

• The behavior in N is tight in general, but quite pessimistic in average. Indeed, an adversarial situation corresponds to all the sources being located at the same place. In practice if we assume that the sources are well spread over the domain, better bounds can be obtained.

Additive white Gaussian noise Theorem 3.3 is tight for an arbitrary (possibly adversarial) noise term b n . In the case of white Gaussian noise though, the bound is pessimistic and can be improved significantly. The following theorem clarifies this aspect.

Theorem 3. [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF]. In what follows, c 1 , c 2 , c 3 , c 4 denote absolute constants (i.e., not depending on any parameter of the problem). Assume that:

• The noise is white and Gaussian b n ∼ N (0, σ2 Id).

• The minimization domain is the unit ball D = {x, x 2 ≤ 1}.

• The mapping x → Π R(x) is L-Lipschitz continuous with L ≥ c 1 : Π R(x) -Π R(x ) 2→2 ≤ L x -x 2 for x, x ∈ D. (23)
• The following inequality holds for x, x ∈ D:

Π R(x) Π R(x ) 2→2 ≤ 1 (β x -x 2 ) α with β > 0 and α > 1/2. ( 24
)
For any ρ > 0, set

= φ -1 c 2 • σ y 0,n 2 L β α α+1 • √ DI + ρ + σ 2 y 0,n 2 2 D log(L) + DI log(L) + ρ . (25)
Then

P ( xn -xn 2 ≥ ) ≤ exp (-c 3 ρ) + exp -c 4 ρ 2 . ( 26
)
Remark 3.2. Let us formulate a few remarks.

1. We did not keep track of the constants in the proof, but they are moderate, like

√ π
3. The condition L ≥ c 1 is not necessary. In fact, stronger statements can be obtained for L ≤ c 1 since we reach a different regime: for small L the projection matrices Π R(x) are similar on the unit disk and the noise b n generates less oscillations of the cost function. However, this condition seems less significant from a practical viewpoint and we decided to restrict the theorem for readability.

Application to a simple example

To illustrate Theorem 3.5, we provide a minimalist application below. We work in the 1D setting, i. To obtain a super-resolution effect, we want the precision to be smaller than Shannon's rate, i.e., a. This can be obtained by setting σ 1 √ τ a log(a) , ensuring that the term in the square root above is smaller than a constant. This result can be analyzed as follows:

• The scaling of σ as 1 √ a was to be expected, since this quantity corresponds to the amplitude of the signal.

• The multiplicative factor log(a) requires a bit more attention. By decreasing a, more sampling points are available on the interval D. Hence, the probability that the noise term b n correlates with the convolution kernel e gets higher. However, this probability only increases slowly with the number of sampling points. For instance, it is well known that the supremum of a random Gaussian vector in R M with mean 0 and covariance Id is of the order of log(M ). The multiplicative term log(a) reflects this phenomenon.

• The noise level can increase proportionally the square root of the oversampling factor τ -1 . This result illustrates how the oversampling factor allows to increase the localization accuracy for a fixed noise level, or on the contrary increase the noise level for a fixed localization accuracy.

• Under a white Gaussian noise assumption, Theorem 3.5 is significantly more powerful than Theorem 3.3. For this example, the2 -norm of the noise is not even bounded since there is an infinite number of samples. Therefore Theorem (3.3) cannot be applied. If we measured the 2 -norm of the noise on the samples in D only, it would scale in average as σ 2 a . The condition b n 2 2 y 0,n 2 2 would therefore translate to σ 1. In comparison, Theorem (3.5) allows for a scaling as σ 1 a 1/2 log(a) , which goes to infinity as a goes to 0!

Estimating the operator with known weights

In this section, we study the problem of estimating the operator parameterization γ under the assumption that the weights are wn are known. Once the positions xn have been estimated for every observation y n , the vector γ can be estimated by solving P γ . This is also equivalent to the following linear system:

N n=1 w2 n E * ( xn )E( xn ) γ = N n=1 wn E * ( xn )y n . (28)
Assumption 2.1 is sufficient to ensure that γ is unique with N = 1 observation. It is possible to guarantee the closeness between γ and γ under an additional Lipschitz regularity assumption on E. We start working with N = 1 observed impulse response. Theorem 3.7 (Stability of the operator estimate with a single observation ). Assume that N = 1 and that E is

√ σ + L E -Lipschitz continuous 2 : E(x) -E(x ) 2→2 ≤ √ σ + L E x -x 2 for all (x, x ) ∈ R D × R D . ( 29
)
Then, under Assumption 2.1, we have

γ -γ 2 γ 2 ≤ κ 3/2 b 1 2 y 0,1 2 + 2 ( x) (30)
with

2 ( x) = cκ 5/2 L E x -x 2 1 + b 1 2 y 0,1 2 + O x -x 2 2
for some absolute constant c.

Together with Theorem 3.3, this last result ensures that γ → γ when the noise level b 1 2 vanishes. This means that we can stably recover an operator when observing a single impulse response.

Unfortunately, Assumption 2.1 is not always met in practical situations of interest as outlined in Proposition 2.4. In that case, observing multiple impulse responses N > 1 can still make a stable estimation possible. In what follows, we let

Y 0 def. =    y 0,1 . . . y 0,N    and B def. =    b 1 . . . b N   
denote the stacked noiseless measurements and noise vectors. =

N n=1 w2 n E * (x n )E(x n ).
Let σ-= w2 -σ-and σ+ = w2 + σ+ and assume that σ-Id C( X) σ+ Id.

Similarly to Theorem 3.7, assume that E is σ+ L E -Lipschitz continuous and let κ = σ+ σ-. Then we have

γ -γ 2 γ 2 ≤ κ3/2 B 2 Y 0 2 + 2 ( X) (32)
with

2 ( X) ≤ cκ 5/2 L E X -X 2 1 + Y 0 2 + B 2 Y 0 + B 2 + O X -X 2 2 ,
for some absolute constant c.

Assumption ( 31) is a geometrical condition intertwining the locations of the Dirac masses and the observation mapping E. It can be hard to verify in advance. However it only requires computing the I × I matrix C( X), which can be achieved once X has been evaluated. The stable estimation of X on its side only depends on Assumption 2.2, which can be verified in advance and can be satisfied independently of Assumption 2.1. Hence, Theorem 3.8 actually yields a constructive result to guarantee the stable recovery of an operator with the following approach:

• If Assumption 2.2 is satisfied and the noise level is low, estimate X.

• Evaluate the condition number κ of C( X).

• If κ is sufficiently low, γ provides a good estimate of γ.

The case of unknown weights

Minimizing (P γ,w ) with respect to (w, γ) for a fixed X = X is a bilinear inverse problem. It received a considerable attention lately, with numerous progress both on the necessary and sufficient conditions to guarantee the recovery [START_REF] Ahmed | Blind deconvolution using convex programming[END_REF][START_REF] Jung | Blind demixing and deconvolution at near-optimal rate[END_REF][START_REF] Li | Identifiability and stability in blind deconvolution under minimal assumptions[END_REF][START_REF] Ahmed | Leveraging diversity and sparsity in blind deconvolution[END_REF][START_REF] Krahmer | On the convex geometry of blind deconvolution and matrix completion[END_REF], on the optimal stability to noise [START_REF] Chen | Convex and nonconvex optimization are both minimax-optimal for noisy blind deconvolution[END_REF], and on the numerical aspects through convex lifting [START_REF] Beinert | Tensor-free proximal methods for lifted bilinear/quadratic inverse problems with applications to phase retrieval[END_REF] or local optimization [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF][START_REF] Cambareri | Through the haze: a non-convex approach to blind gain calibration for linear random sensing models. Information and Inference[END_REF][START_REF] Traonmilin | The basins of attraction of the global minimizers of the non-convex sparse spike estimation problem[END_REF][START_REF] Zhu | Global optimality in low-rank matrix optimization[END_REF][START_REF] Li | Rapid, robust, and reliable blind deconvolution via nonconvex optimization[END_REF]. Our intention here is not to produce new results, but to discuss how the existing results apply to the current context. As an executive summary of what follows: the current theoretical results are still insufficient to guarantee a stable recovery in general, but we will see that some optimization methods perform well experimentally.

A low-dimensional bilinear problem

Let În

E( xn ) = Ûn V * n , (33)
where Ûn ∈ R M × În and Vn contain orthogonal columns and Û * n Ûn = Id. Hence, letting c n def.

= Û *

n y n , we obtain argmin

w∈R N ,γ∈R I J(w, γ, X) = argmin w∈R N ,γ∈R I 1 2 N n=1 Ûn V * n w n γ -y n 2 2 = argmin w∈R N ,γ∈R I 1 2 N n=1 V * n w n γ -c n 2 2 .
Letting B : R N × R I → R Î denote the following bilinear mapping:

B(w, γ) def. =    V * 1 w 1 γ . . . V * N w N γ    and c def. =    c 1 . . . c N    , (34)
we can rewrite J more compactly as J(w, γ, X) = 1 2 B(w, γ) -c 2 2 , and hence:

argmin w∈R N ,γ∈R I J(w, γ, X) = argmin w∈R N ,γ∈R I 1 2 B(w, γ) -c 2 2 . ( 35
)
Notice that the dimension M , which might be huge in applications, completely disappeared from this formulation.

A review of existing conditions for stable recovery

Recovering w and γ is possible only up to a multiplicative constant since J(tw, γ/t, X) = J(w, γ, X) for all t = 0. Now, consider the noiseless setting B = 0 and assume that the locations are perfectly recovered: X = X. In that situation, a necessary condition to recover ( w, γ) modulo the above scaling ambiguity is that there exists a unique pair (w, γ) with w 2 = 1 such that B(w, γ) = c. From our current understanding, deriving geometrical conditions to ensure this local injectivity condition still deserves some attention.

In [START_REF] Kech | Optimal injectivity conditions for bilinear inverse problems with applications to identifiability of deconvolution problems[END_REF][START_REF] Li | Identifiability and stability in blind deconvolution under minimal assumptions[END_REF], the authors study a more stringent global injectivity condition of the form

∀c ∈ R Î , ∃ a unique (w, γ) with w 2 = 1 s.t. B(w, γ) = c. ( 36
)
Their main result states that a necessary condition for B to be globally injective is that

Î ≥ 2(N + I) -4, (37)
which provides a rule on how to choose the number of measurements N . In addition, they prove that almost every bilinear mapping B with respect to the Lebesgue measure is globally injective provided that the inequality (37) holds. This is a beautiful contribution. Unfortunately, it suffers from two limitations in the current setting: First, the operator B that appears in our formulation possesses a peculiar structure which may fall in a set of 0 measure. Second, the result does not certify that a low complexity algorithm can actually recover the factors.

As for stability to noise, nearly all the existing results rely on some kind of randomness in the design of the bilinear mapping B. They do not apply to the current context where everything is deterministic. Overall, the results of the algorithms described below are therefore empirical.

Optimization methods

In this section, we review 3 algorithms to solve [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF].

Optimization of the factors Solving (35) can be achieved using local optimization over each factor w and γ [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF][START_REF] Zhu | Global optimality in low-rank matrix optimization[END_REF][START_REF] Li | Rapid, robust, and reliable blind deconvolution via nonconvex optimization[END_REF]. A simple approach consists in using an alternate minimization between the factors as outlined in Algorithm 1. Notice that every step of the algorithm can be performed efficiently since Algorithm 1 Alternating minimization Require: Initial guess:

w 1 ∈ R N . Require: Iteration number K. for all k = 1 → K -1 do γ k+1 = argmin γ∈R I 1 2 B(w k , γ) -c 2 2 .
w k+1 = argmin

w∈R N 1 2 B(w, γ k+1 ) -c 2 2 . end for return (w K , γ K ).
the dimensions of the problem are significantly reduced. This approach can be certified to recover a stable estimate ( ŵ, γ) of ( w, γ) provided that a clever initialization is used [START_REF] Zhu | Global optimality in low-rank matrix optimization[END_REF][START_REF] Li | Rapid, robust, and reliable blind deconvolution via nonconvex optimization[END_REF]. Sufficient recovery guarantees are for instance provided when the bilinear mapping B is chosen at random. This method also allows to easily incorporate constraints (e.g. nonnegativity) in the factors, which can sometimes allow a significantly improved reconstruction. In all our numerical experiments, we will use the spectral initialization from [START_REF] Li | Rapid, robust, and reliable blind deconvolution via nonconvex optimization[END_REF] as a starting point.

Optimization over rank-1 matrices

The bilinear mapping B(w, γ) can be rewritten as a linear mapping Λ on the rank-1 outer product T = wγ T : B(w, γ) = Λ(T ). Hence, we have:

inf w∈R N ,γ∈R I J(w, γ, X) = inf T ∈R N ×I ,rank(T )=1 1 2 Λ(T ) -c 2 2 . ( 38
)
The interest of the right-hand side in equation ( 38) compared to the left-hand side is that the scaling ambiguity is discarded. Letting T denote the set of rank-1 matrices, this alternative formulation can be solved using a projected gradient descent described in Algorithm 2. To the best of our knowledge, this Algorithm 2 Projected gradient descent Require: Initial guess: T ∈ R N ×I . Require:

Iteration number K. Compute τ = 1 Λ 2 2→2 using a power iteration. for all k = 1 → K -1 do T k+1 = Π T T k -τ Λ * ( Λ(T k ) -c) . end for Decompose T K = w K γ * K . return (w K , γ K ).
algorithm has not been analyzed so far. During the review process, we found a paper [START_REF] Eisenmann | Riemannian thresholding methods for row-sparse and low-rank matrix recovery[END_REF] describing a similar type of idea. Again, we will use the spectral initialization from [START_REF] Li | Rapid, robust, and reliable blind deconvolution via nonconvex optimization[END_REF] as a starting guess for this algorithm in the numerical experiments.

Convex relaxation using the nuclear norm Finally, a popular method [START_REF] Ahmed | Blind deconvolution using convex programming[END_REF][START_REF] Ahmed | Leveraging diversity and sparsity in blind deconvolution[END_REF][START_REF] Beinert | Tensor-free proximal methods for lifted bilinear/quadratic inverse problems with applications to phase retrieval[END_REF][START_REF] Chi | Guaranteed blind sparse spikes deconvolution via lifting and convex optimization[END_REF] is a convex relaxation using the nuclear norm. The usual convex relaxation of the nonconvex problem [START_REF] Traonmilin | The basins of attraction of the global minimizers of the non-convex sparse spike estimation problem[END_REF] is the following: inf

T ∈R N ×I , Λ(T )=c T * or inf T ∈R N ×I 1 2 Λ(T ) -c 2 2 + λ T * , (39)
where λ > 0 is a regularization parameter and • * is the nuclear norm, i.e. the sum of the singular values of T . This convex function over the space of matrices is well known to promote low-rank solutions since the extreme points of the associated unit ball are the rank-1 matrices [START_REF] Boyer | On representer theorems and convex regularization[END_REF]. The stable recovery of the tensor w γT has been established under rather stringent conditions based on random subspace assumptions [START_REF] Ahmed | Blind deconvolution using convex programming[END_REF][START_REF] Ahmed | Leveraging diversity and sparsity in blind deconvolution[END_REF]. Experimentally, the method seems to provide satisfactory results under much weaker conditions. From a numerical perspective, Problem (39) can be solved using a diversity of proximal algorithms, such as an accelerated proximal gradient descent or a Douglas-Rachford algorithm [START_REF] Patrick | Proximal splitting methods in signal processing[END_REF]. We do not further detail these algorithms, which are well documented in the literature.

Numerical experiments

The aim of this section is to illustrate the proposed theory using simple 1D examples and to explain the setting of the 2D experiment in Figure 1.

Convolution operators with known weights

We start with an illustration of Theorem 3.3 using convolution operators only. We focus on the case of pointwise sampling on [0, 1], by setting ν m = δ zm , with z m = m/M for m ∈ {1, . . . , M }. Notice that this case also covers the case of product-convolution operators since the ranges R(x) of convolution and product-convolution operators are identical.

The families of operators

We consider three families of convolution operators A 1 , A 2 and A 3 differing by the choice of the convolution filters.

Family A 1 is defined through a set of convolution operators A i with Gaussian filters (e i ) defined by:

e i (x) = exp(-x 2 /(2σ 2 i )) with σ i = 0.01 • i -1 I -1 + 0.03 • 1 - i -1 I -1 .
Using this family in a blind deconvolution problem allows to identify the variance of a Gaussian convolution filter. Gaussian convolution filters are amongst the most popular simplified point spread function models in microscopy.

Family A 2 is also defined using Gaussian convolution filters, but the standard deviation ranges in [0.03, 0.09] instead of [0.01, 0.03].

Family A 3 is defined with less regular convolution filters. Let ψ(x) = (1 -|x -1|) + denote the hat function.

e i = ψ(x • σ i ) where σ i = 0.02 • i -1 I -1 + 0.2 • 1 - i -1 I -1 .
In all settings we set I = 3. The filters corresponding to each family are displayed in Figure 2. We then orthogonalize the filters using a singular value decomposition on a very fine grid. This leads to a new family of orthogonal filters (e ⊥ i ) which will be used in all experiments to satisfy Assumption 2.3.

The inverse functions φ -1

As stated in Theorem 3.3, the critical element to guarantee a stable recovery of the locations xm is the function φ and its inverse, which characterizes the angle between the ranges R(x) and R(x ). To evaluate this function, we first sample the function Π R(0) Π R(k∆x) 2→2 on a fine grid. We store the result in the vector φ 0 (k)

def.

= 1 -Π R(0) Π R(k∆x) 2→2 with a sampling step ∆x. This function is not This problem is convex and can be solved using the CVX library [START_REF] Grant | Cvx: Matlab software for disciplined convex programming[END_REF] for instance. We use the solution φ of this problem in place of φ in Assumption 2.2. The inverse filters are displayed in Figure 3. The stability to noise is dependent on the speed of ascent of φ -1 + . As can be seen by comparing the two Gaussian families, the smallest the filter, the slower the ascent. Hence, very localized impulse responses should be easier to detect with a good accuracy than larger ones. Also notice that the regularity of the convolution kernels seem to have little importance since the inverses φ -1

1,+ and φ -1 3,+ behave roughly similarly in terms of speed of ascent.

Stability of the locations

Here, we study the robustness of the estimation to noise. To this end, we compute the empirical average of the error E(|x -x|) for various noise levels and realizations. The expectation is estimated by averaging 100 noise realizations. We fix γ once for all. We use white Gaussian noise, i.e. b n ∼ N (0, σ 2 Id), with σ = θ y 0 2 / √ M and θ ∈ [0, 2]. Figure 4 shows the resuting signals with M = 100 for the noise levels θ ∈ {0, 1, 2} and each family. Notice that θ = 1 corresponds to an expected norm of noise equal to the signal's norm. Hence, we consider rather extreme noise levels. We will see that the localization accuracy is surprisingly good in spite of this challenging setting.

The empirical estimate of E(|x-x|) is displayed with respect to the noise level σ in Figure 5. The error is displayed as a proportion of a pixel. For instance, an error of 0.1 means that the localization was accurate at a tenth of a pixel. Hence, we can expect a super-resolution effect for a precision below 0.5. This accuracy is obtained for all families up to the noise level θ = 1, corresponding to a measurement dominated by noise.

To end this experiment, we evaluate γ for all experiments and display the relative error γ-γ 2 γ 2 for all families of operators. Letting h = I i=1 γi e ⊥ i and ĥ = I i=1 γi e ⊥ i denote the true convolution filter and the estimated one, notice that we have

h-ĥ L 2 (R D ) h L 2 (R D ) = γ-γ 2
γ 2 , since the family (e ⊥ i ) is orthogonal. In Figure 6, we see that the reconstruction errors for any family of convolution filters behave really similarly. In particular, the errors using the family A 1 and the family A 2 are nearly identical. This might come as a surprise since the localization errors were significantly higher for the family A 2 , which is a scaled version of A 1 . This fact can be explained by the fact that the Lipschitz constant L E in Theorem 3.7 is inversely proportional to the scaling of the Gaussian, which compensates for the localization errors.

Product-convolution operators and unknown weights

The objective of this section is to compare the different algorithms described in Section 3.4 for the specific case of 1D product-convolution operators described in Assumption 2.4. In this experiment, we work on a grid and set X = X since the objective is not to assess the localization performance, but rather the ability to solve a bilinear inverse problem. We use the pointwise sampling model ν m = δ zm with z m = 10 • m/M and M = 1000. This corresponds to a uniform sampling of the interval [0, 10]. For the filters (e j ), we use the family of Gaussian convolution kernels A 1 with J = 3. We set the vectors f k as smooth random Gaussian processes by convolving a random vector with distribution N (0, Id) with a Gaussian filter of large variance. Once the family of operators is defined through the pairs of families (e j ) and (f k ), we can sample operators at random in this family by setting γ ∼ N (0, Id). In Figure 7, we visualize a set of operators indirectly by applying them to a Dirac comb with 4 spikes. As can be seen, the operators are space-varying. Their impulse responses belong to span(e j , j ∈ {1, . . . , J}).

To assess the performance of the different algorithms, we evaluate the percentage of perfect recovery results with various values of K and N . We run the algorithms 100 times with random locations for the N spikes xn , with random weights wn and with a random family (f k ). The proportions of perfectly recovered operators are displayed in Figure 8. For the considered families, the nuclear norm relaxation performs very poorly, suggesting that the relaxation approaches advocated both for discrete and gridless problems might not be the best competitors. In comparison, the alternate minimization (Algorithm 1) with the spectral initialization from [START_REF] Li | Rapid, robust, and reliable blind deconvolution via nonconvex optimization[END_REF] and the seemingly novel projected gradient descent (Algorithm 2) perform satisfactorily for a good range of values of K and N . Between those two, the projected gradient descent seems to provide better results for a wider range of parameters. A theoretical analysis of this idea might be worth an exploration. Unfortunately, no algorithm is able to succeed systematically. This might be related to the fact that the random positions (x n ) are badly located.

In this setting, the condition for global injectivity (37) reads:

N ≥ 2JK -4 J -2 ,
i.e. N ≥ 6K -4 for J = 3. We see the shortcomings of this rule in this experiment: perfect recovery does not always occur when this condition is satisfied, because the condition does not certify the success of an algorithm. And the algorithms manage to recover some operators when this condition is not met. However, it is clear that a necessary condition for identifiability is N ≥ K, since otherwise, even the problem with known weights cannot be identified.

A 2D experiment

To end this paper, we briefly describe how the experiment from Figure 1 was carried out. We generated a family of product-convolution operators with astigmatic impulse responses as follows. We set the family (e j ) as anisotropic Gaussian vectors with J = 8. We also set the family (f k ) as the monomials 1, x and y, resulting in K = 3 basis elements to describe the space variations. We launched the maximum correlation algorithm to locate the beads positions in Figure 1b. The average localization error is 0.015 pixel, despite a significant amount of additive Gaussian noise. We then discarded by hand the locations that were too close from each other (red stars). Notice that this part can be easily automatized by thresholding the minimal distance between adjacent locations. We kept the other locations (blue stars) to estimate the operator, resulting in a total of N = 27 impulse responses with a slightly inaccurate localization. We then used this information to recover the operator. Here we assumed that the weight ( wn ) were known and all equal to 1. The relative error between the operator estimated and the true one is 0.006. Here, we measured the distance between operators with the Hilbert-Schmidt norm. The whole process takes less than a second on a usual personal computer with Matlab.

Proof of Proposition 2.3

Proof. By Assumption 2.3, we have Π R(x) = E(x)E * (x). Then, by definition, we have

Π R(x) Π R(x ) 2→2 = E(x)E * (x)E(x )E * (x ) 2→2 = sup u∈R M , u 2=1 E(x)E * (x)E(x )E * (x )u 2 = sup v∈R I , v 2 =1 E(x)E * (x)E(x )v 2 = E(x)E * (x)E(x ) 2→2 = E * (x)E(x ) 2→2 = C(x -x ) 2→2 , using the fact that E * (x)E(x) = Id.
Proof of Proposition 2.4

Proof. The first part of the proof is trivial: the range R(x) of E(x) is unchanged. As for the second part, it suffices to realize that E contains columns which are colinear.

Proof of Proposition 3.1

Proof. We have

inf x∈D,α∈R I 1 2 E(x)α -y n 2 2 = inf x∈D inf α∈R I 1 2 E(x)α -y n 2 2 .
The minimum norm solution of the inner problem is given by α(x) = E + (x)y n , where E + (x) is the Moore-Penrose pseudo-inverse of E(x). Injecting this solution leads to the problem

inf x∈D 1 2 E(x)E + (x)y n -y n 2 2 = inf x∈D 1 2 Π R(x) y n -y n 2 2 = inf x∈D 1 2 y n 2 2 + 1 2 Π R(x) y n 2 2 -Π R(x) y n , y n = inf x∈D 1 2 y n 2 2 - 1 2 Π R(x) y n 2 2 .
Neglecting the constant term 1 2 y n 2 2 and changing the sign of the function yields the result.

Proof of Proposition 3.2

Proof. It is well known (see e.g., [START_REF] Stewart | On the continuity of the generalized inverse[END_REF]) that under the stated conditions, the pseudo-inverse E(x) + is a continuous mapping of x. This also implies the continuity of the projection mapping Π R(x) = E(x)E + (x) and of the whole function x → Π R(x) y n 2 2 . A continuous function over a compact domain admits at least one global minimizer.

Proof of Theorem 3.3

We start with a basic lemma. = f + and assume that the following sets are non-empty

X = argmax x∈S g(x) and X = argmax x∈S f (x).
Further assume that L ∞ (S) ≤ η for some η > 0 and that there exists an increasing function ϕ : R → R + such that

f (x) ≤ f ( x) -ϕ( x -x 2 ), ∀x ∈ S, x ∈ X . (40)
Then X = { x} is a singleton and any x ∈ X satisfies xx 2 ≤ ϕ -1 + (2η). Proof. By inequality [START_REF] Li | Rapid, robust, and reliable blind deconvolution via nonconvex optimization[END_REF] and strict monotonicity of φ, we have f ( x) > f (x) for all x = x. Hence x is the unique maximizer of f . We have

g( x) ≥ g( x) = f ( x) -( x) ≥ f ( x) -η. (41)
In addition, for any x ∈ S, we have

g(x) = f (x) + (x) ≤ f (x) + η ≤ f ( x) -ϕ( x -x 2 ) + η.
For any x ∈ S such that xx 2 > ϕ -1 + (2η), we have g(x) < f ( x) -η, and thus g(x) < g( x), which implies that x = x. The contraposition is that any x ∈ X satisfies xx 2 ≤ ϕ -1 + (2η). We are ready to prove Theorem 3.3.

Proof. i) Let y 0,n = wn E( xn )γ = y n -b n denote the noiseless measurements and let F 0,n (x, γ) = 1 2 wn E(x)γy 0,n 2 2 . We have F 0,n ( xn , γ) = 0. Let γ 0 (x) ∈ argmin γ∈R I F 0,n (x, γ) denote any minimizer (for instance the one given by the pseudo-inverse). Now, define G 0,n (x) def.

= F 0,n (x, γ 0 (x)). We have

G 0,n (x) = 1 2 y 0,n 2 2 -Π R(x) y 0,n 2 2 . By Assumption 2.2, R(x)∩R( x) = {0} for x = x. Hence, Π R(x) y 0,n 2 2 < y 0,n 2
2 for x = x and xn is the unique minimizer of G 0,n . Therefore, the function

H 0,n (x) def. = 1 2 Π R(x) y 0,n 2 2 = 1 2 Π R(x) y 0,n , y 0,n ,
also admits a unique maximizer in x. Overall, we see that under Assumption 2.2, G 0,n admits a unique minimizer equal to x. Under the additional Assumption 2.1, F 0,n admits a unique solution ( x, γ). ii) Now, let

F n (x, γ) = 1 2 wn E(x)γ -y 0,n 2
2 , γ(x) denote any minimizer of

F n w.r.t. γ, G n (x) = F n (x, γ(x)) and H n (x) def. = 1 2 Π R(x)
y n , y n . Let x denote any maximizer of H n and assume for now that we manage to obtain a bound of the form |H n (x) -H 0,n (x)| ≤ η for some η ≥ 0 and all x ∈ D. We have

H 0,n (x) = 1 2 Π R(x) y 0,n , y 0,n = 1 2 Π R(x) y 0,n , Π R( x) y 0,n = 1 2 Π R( x) Π R(x) y 0,n , y 0,n ≤ 1 2 Π R( x) Π R(x) 2→2 y 0,n 2 2 ≤ 1 2 (1 -φ( x -x 2 )) y 0,n 2 2 = H 0,n ( x) - 1 2 φ( x -x 2 ) y 0,n 2 2
by Assumption 2.2. Hence, we can use Lemma 5.1 with S = D, f (x) = H 0,n (x), g(x) = H n (x) and ϕ(r) = 1 2 φ(r) y 0,n 2 2 . This allows us to conclude that

xx 2 ≤ φ -1

+ 4η y 0,n 2 2 . ( 42
)
iii) The last remaining point is to control H 0,n -H n L ∞ (D) . We have

H n (x) = 1 2 Π R(x) (y 0,n + b n ) , y 0,n + b n = H 0,n (x) + Π R(x) y 0,n , b n + 1 2 Π R(x) (b n ) 2 2 .
Hence, using Cauchy-Schwarz inequality, we obtain for all x

|H n (x) -H 0,n (x)| ≤ Π R(x) y 0,n 2 Π R(x) b n 2 + 1 2 Π R(x) (b n ) 2 2 .
Using the facts that Π R(x) b n 2 ≤ b n 2 and that b n 2 ≤ θ y 0,n 2 , we obtain

H n -H 0,n L ∞ (D) ≤ y 0,n 2 2 θ + 1 2 θ 2 .
For the inequality [START_REF] Eisenmann | Riemannian thresholding methods for row-sparse and low-rank matrix recovery[END_REF] to make sense, we need 4(θ+

1 2 θ 2 ) ≤ 1, which is equivalent to θ < √ 6
2 -1 and Theorem 3.3 is proven.

We have

H(x) = H 1 (x) + Π R(x) y 1 , Π R(x) N n=2 y n 1 (x) + 1 2 Π R(x) N n=2 y n 2 2 2(x) , (44)
where

H 1 (x) = 1 2 Π R(x)
| 1 (x)| = Π R(x) y 1 , Π R(x) N n=2 y n ≤ N n=2 Π R(x) Π R( x1) y 1 , Π R(x) Π R( xn y n ≤ N n=2 (1 -φ( x -x1 2 ))(1 -φ( x -xn 2 )) y 1 2 y n 2 ≤ (N -1)c(1 -φ(δ)) y 1 2 2 ≤ (N -1)c(1 -φ(δ/2)) y 1 2 2 .
Similarly, using Jensen inequality, we obtain

| 2 (x)| = 1 2 Π R(x) N n=2 y n 2 2 ≤ N -1 2 N n=2 (1 -φ(δ/2)) 2 y n 2 2 ≤ (N -1) 2 c 2 2 (1 -φ(δ/2)) 2 y 1 2 2 . Hence L ∞ (B δ/2 ) ≤ (N -1) 2 c 2 2 (1 -φ(δ/2)) 2 y 1 2 2 + (N -1)c(1 -φ(δ/2)) y 1 2 2 . Set Z = (N -1)c. Lemma 5.1 states that any global maximizer x of H over B δ/2 satisfies x -x1 2 ≤ φ -1 4 Z 2 2 (1 -φ(δ/2)) 2 + Z(1 -φ(δ/2)) = a • 2 1/b • Z • (2 + Z + 2(δ/2a) b ) (1 + (δ/2a) b ) 2 -2Z 2 -4Z(1 + (δ/2a) b ) 1/b
, where the last line was obtained using a symbolic calculus software. Letting τ ≥ 5, the denominator above is positive. Setting δ = 2a(τ Z) 1/b , the above expression becomes

r ≥ a • 2 1/b τ • (τ 2 + τ )Z + 2τ 1 + (τ 2 -4τ -2)Z + (2τ -4)Z .
For τ ≥ 5, the above expression is bounded above by 18.4 • a•2 1/b τ , obtained for τ = 5 and Z → ∞.

Proof of Theorem 3.5

This theorem is the most technical and requires a few intermediate results. We start with an important observation. Proposition 5.2. Assume that b ∼ N (0, σ 2 Id) is white Gaussian noise of variance σ 2 . Define the following two random processes

∆ 1 (x) def. = Π R(x) y 0 , Π R(x) b and ∆ 2 (x) def. = 1 2 Π R(x) (b) 2 2 . ( 45
)
Then under Assumption 2.2

xx 2 ≤ φ -1

+ 2Ampl(∆ 1 + ∆ 2 ) y 0 2 2 , (46)
where Ampl(f )

def. = sup x∈D f (x) -inf x∈D f (x).
Proof. This proposition derives from point iii) in the previous proof. The inequality ( 42) can be improved as

x -x 2 ≤ φ -1 + 4 inf c∈R H -H 0 -c ∞ y 0 2 2 , (47)
since a constant term does not affect the location of a minimizer. We have

∆(x) def. = H(x) -H 0 (x) = Π R(x) y 0 , Π R(x) b + 1 2 Π R(x) (b) 2 2 .
It is therefore natural to set c = 1 2 (sup x∈R D ∆(x) -inf x∈R D ∆(x)), to minimize in the infinite norm in Problem [START_REF] Boucheron | Concentration inequalities: A nonasymptotic theory of independence[END_REF]. This proposition reveals that the critical quantity to control, to evaluate the localization error is the amplitude of the random process ∆ 1 + ∆ 2 . Obtaining tight analytical bounds for this is a difficult problem in general. Hopefully the following proposition shows that it can be evaluated efficiently using numerical procedures.

Proposition 5.3. We have Ampl(∆ 1 + ∆ 2 ) ≤ Ampl(∆ 1 ) + Ampl(∆ 2 ). In addition, the random variable Z 1 = Ampl(∆ 1 ) is sub-Gaussian:

P(|Z 1 -Z1 | ≥ t) ≤ 2 exp -t 2 /(8σ 2 y 0 2 2 )
and the random variable Z 2 = Ampl(∆ 2 ) is sub-exponential:

P(|Z 2 -Z2 | ≥ t) ≤ 2 exp -Ct/σ 2 ,
where Z1 and Z2 are the expectations of Z 1 and Z 2 and C is a universal constant.

Proof. First notice that

Ampl(∆ 1 + ∆ 2 ) = sup x∈D ∆ 1 (x) + ∆ 2 (x) -inf x∈D ∆ 1 (x) + ∆ 2 (x) ≤ sup x∈D ∆ 1 (x) + sup x∈D ∆ 2 (x) -inf x∈D ∆ 1 (x) + inf x∈D ∆ 2 (x) = Ampl(∆ 1 ) + Ampl(∆ 2 ).
We will treat the two random processes ∆ 1 and ∆ 2 separately. i) Consider the function f 1 (b) def.

= sup x∈D Π R(x) y 0 , b and define the random variable

V + 1 = f 1 (b) with mean V + 1 . Similarly, define V - 1 = inf x∈D Π R(x) y 0 , b . In addition, notice that Z 1 = Ampl(∆ 1 ) ≤ V + 1 -V - 1 .
We first show that f 1 is Lipschitz continuous. We have

f 1 (b + ) = sup x∈D Π R(x) y 0 , b + ≤ sup x∈D Π R(x) y 0 , b + y 0 2 2 f 1 (b + ) = sup x∈D Π R(x) y 0 , b + ≥ sup x∈D Π R(x) y 0 , b -y 0 2 2 30 Hence, |f 1 (b) -f 1 (b + )| ≤ y 0 2
2 and f 1 is y 0 2 -Lipschitz continuous. Using a Gaussian logarithmic Sobolev inequality [START_REF] Boucheron | Concentration inequalities: A nonasymptotic theory of independence[END_REF]Thm 5.6], we obtain that V + 1 is sub-Gaussian with

P |V + 1 -V + 1 | ≥ t ≤ 2 exp(-t 2 /(2σ 2 y 0 2 )).
The same result holds for V - 1 . Finally, the sum of two dependent sub-Gaussian variables with parameters σ 1 and σ 2 is sub-Gaussian with a sub-Gaussian parameter smaller than σ 1 + σ 2 , so that

P |Z 1 -Z1 | ≥ t ≤ 2 exp(-t 2 /(8σ 2 y 0 2 )).
ii) Now, define the random variable

Y + 2 def. = sup x∈D Π R(x) b 2 and V + 2 def. = 1 2 (Y + 2 ) 2 . The function b → sup x∈D Π R(x) b 2 is 1-Lipschitz continuous.
Hence using a Gaussian logarithmic Sobolev inequality again, we obtain that Y + 2 is sub-Gaussian with

P |Y + 2 -Ȳ + 2 | ≥ t ≤ 2 exp - t 2 2σ 2 .
Using [48, Lemma 2.7.6], we conclude that V + 2 is sub-exponential and satisfies

P |V + 2 -V + 2 | ≥ t ≤ 2 exp(-Ct/σ 2 ),
for an absolute constant C. We can make a similar proof for the random variable Y - 2 def.

= inf x∈D Π R(x) b 2 and conclude as in the previous proof.

Notice that the averages Z1 and Z2 have no reason to be defined in general since the amplitude could be infinite. We will see later that under suitable assumptions, they are well defined. Proposition (5.3) has two consequences. First, we see that the deviation around the mean of Z 1 scales as σ y 0 2 and the deviation around the mean of Z 2 scales as σ 2 . Second, Hoeffding [48, Thm 2.6.1] and Bernstein [48, 2.8.1] inequalities imply that computing an empirical average of Z 1 and Z 2 will converge rapidly to the true means Z1 and Z2 . Hence, it is possible to obtain a precise numerical estimate using an empirical average and we know that the probability that the variables deviate from the means by more than σ y 0 2 + σ 2 is negligible.

Unfortunately, the averages Z1 and Z2 are difficult to compute in general. Hence, the above proposition can only be used to estimate average deviations with a computer. We will now turn to bound Z1 and Z2 under additional regularity assumptions.

∈ R D Z1 ≤ c • L β α α+1 • σ √ DI y 0 2 , ( 48
)
where c is an absolute constant.

Proof. Here, we wish to control the supremum Z1 of the centered Gaussian process ∆ =

E ((∆ 1 (x) -∆ 1 (x )) 2 ). (49)
Let B(c, δ, dist) = {x ∈ R D , dist(c, x) ≤ δ} denote a ball of radius δ centered at c with respect to dist. Let S ⊆ D denote a set and define the covering number Θ(δ, S) as the minimum number of δ-balls needed to cover S. We recall the following simple results: for S ⊆ S , we have Θ(δ, S) ≤ Θ(δ, S ). In addition, if S can be partitioned as S = S 1 S 2 , we have Θ(δ, S) ≤ Θ(δ, S 1 ) + Θ(δ, S 2 ).

Dudley's theorem reads as follows

E sup x∈S ∆ 1 (x) ≤ 24 η/4 0 log(Θ(u, S)) du (50)
with

η def. = inf x∈S sup x ∈S dist(x, x ). ( 51
)
Bounding the pseudo-metric Using the fact that y 0 ∈ R( x) and that

y 0 = Π R( x) y 0 : dist(x, x ) 2 = E (∆ 1 (x) -∆ 1 (x )) 2 = E ( Π R(x) y 0 , Π R(x) b -Π R(x ) y 0 , Π R(x ) b ) 2 = E ( Π R(x) Π R( x) y 0 , b -Π R(x ) Π R( x) y 0 , b ) 2 = E Π R(x) -Π R(x ) Π R( x) y 0 , b 2 = E Π R(x) -Π R(x ) Π R( x) y 0 , Π R( x) b 2 ≤ Π R(x) -Π R(x ) Π R( x) 2 2→2 y 0 2 2 • E Π R( x) b 2 2 ≤ σ 2 I y 0 2 2 Π R(x) -Π R(x ) Π R( x) 2
2→2 .

In addition

Π R(x) -Π R(x ) Π R( x) 2→2 ≤ min 1, Π R(x) Π R( x) 2→2 + Π R(x ) Π R( x) 2→2 ≤ min 1, 1 β α x -x α 2 + 1 β α x -x α 2 Letting ζ def. = σ √ I y 0 2 and ψ(r) = 1 β α r α , we get dist(x, x ) ≤ ζ • min (1, ψ( x -x 2 ) + ψ( x -x 2 )) . (52)
In addition, the Lipschitz continuity of Π R(x) also implies that dist

(x, x ) 2 = E (∆ 1 (x) -∆ 1 (x )) 2 = E Π R(x) -Π R(x ) y 0 , Π R( x) b 2 (53) ≤ E Π R(x) -Π R(x ) 2 2→2 y 0 2 2 Π R( x) b 2 2
(54)

≤ σ 2 I y 0 2 2 L 2 x -x 2 2 . (55)
In what follows, we let dist(x, x )

def. = ζ • min (1, L x -x 2 , ψ( x -x 2 ) + ψ( x -x 2 )) (56)
and Θ denote the corresponding covering number. The inequality dist(x, x ) ≤ dist(x, x ) implies that Θ(δ, S) ≤ Θ(δ, S) for all S and δ.

Bounding the integral We have dist(x, x ) ≤ ζ for all x, x . Therefore:

E sup x∈D ∆ 1 (x) ≤ E sup x∈R D ∆ 1 (x) ≤ 24 ζ/4 0 log( Θ(u, R D )) du. (57)
By working over R D rather than D, we can assume that x = 0 without loss of generality. We have

ψ( x -x 2 ) + ψ( x -x 2 ) ≤ 2 β α min( x -x 2 , x -x 2 ) α .
Hence, if both x 2 and x 2 are larger than 2ζ

β α u 1/α , dist(x, x ) ≤ u. This implies that Θ u, x ∈ R D , x 2 ≥ 2ζ β α u 1/α = 1.
It remains to control the covering number of the central Euclidean ball of radius

u, x ∈ R D , x 2 ≤ 2ζ β α u 1/α ≤ 3 • 2ζ β α u 1 α • ζL u D . Letting a = 3 • 2 β α 1 α • L D and b = D α+1 α , we obtain: Θ u, R D ≤ a • u ζ -b + 1. (58)
The additional 1 in the expression above requires some attention since it cannot be easily integrated in the logarithm. To discard it, we can integrate the "metric entropy" log( Θ(u, R D )) only up to the smallest u, say u s such that Θ(u, R D ) ≥ 2. Up to this value, a • (u/ζ) -b is necessarily larger than 1, so that u s ≤ ζa 1/b and we obtain

E sup x∈R D ∆ 1 (x) ≤ 24 us 0 log(a • u -b + 1) du ≤ 24 ζa 1/b 0 log(2a • (u/ζ) -b ) du = 24ζ • √ b • (2a) 1/b 2 -1/b 0 log(v -1 ) dv ≤ 24ζ • √ b • (2a) 1/b 1 0 log(v -1 ) dv = 24ζ • √ b • (2a) 1/b • √ π 2 .
Letting c denote an absolute constant, we finally obtain the bound

E sup x∈R D ∆ 1 (x) ≤ c • D α + 1 α • β -α/(α+1) • L α α+1 • σ √ I y 0 2 .
The term α+1 α can be absorbed in the constant since α ≥ 1/2. To conclude, we use the fact that Z1 ≤ 2E (sup x∈R D ∆ 1 (x)). Proposition 5.5 (Control of Z2 ). Under the assumptions of Theorem 3.5, we have, for all x ∈ D

P Z 2 ≥ c 1 σ 2 (D log(L) + D log(L)I) + t ≤ 2 exp -c 2 min t 2 σ 4 (D log(L) + I) , t σ 2
for some absolute constants c, c 1 , c 2 . This implies

Z2 ≤ c • σ 2 D log(L) + D log(L)I . (59)
Proof. Controlling Z2 amounts to studying the supremum of a so-called Gaussian chaos of order 2. This problem arises in different fields and has been addressed using tools of generic chaining in [START_REF] La | Estimates of moments and tails of gaussian chaoses[END_REF][START_REF] Talagrand | The generic chaining: upper and lower bounds of stochastic processes[END_REF][START_REF] Krahmer | Suprema of chaos processes and the restricted isometry property[END_REF][START_REF] Boucheron | Concentration inequalities: A nonasymptotic theory of independence[END_REF]. We will make use of the following theorem, which has been rewritten using our notation. where the covering number is evaluated w.r.t. the spectral distance • 2→2 . Then where we used the inequality erfc(z) < exp(-z 2 )

P sup x∈D |∆ 2 (x) -E(∆ 2 (x))| ≥ σ 2 (c 2 E + t) ≤ 2 exp -c 3 min t 2 σ 4 V 2 ,
√ πz to obtain the last result. If L is larger than exp(1)/L, the term Recall that for a nonnegative random variable X, E(X) = ∞ 0 P(X > t) dt. Hence, the concentration inequality P Ampl(∆ 2 ) ≥ σ 2 (c 2 E + t) ≤ 2 exp -c 3 min We now have all the ingredients to prove the main Theorem 3.5. We have for any real t 1 , t 2

t 2 σ 4 V 2 ,
P Z 1 + Z 2 -Z1 -Z2 ≥ t 1 + t 2 ≤ P ([Z 1 -Z1 ≤ t 1 ] ∩ [Z 2 -Z2 ≤ t 2 ]) c ≤ P [Z 1 -Z1 ≥ t 1 ] ∪ [Z 2 -Z2 ≥ t 2 ] ≤ P Z 1 -Z1 ≥ t 1 + P Z 2 -Z2 ≥ t 2 ≤ 2 exp - t 2 1 8σ 2 y 0 2 2 + 2 exp - ct 2 σ 2 .
where we used Proposition 5.3 to obtain the last inequality. Let ρ > 0 denote an arbitrary number and set t 1 = cρσ y 0 2 , t 2 = cσ 2 ρ. Injecting these values in the above inequality, we obtain

P Z 1 + Z 2 -Z1 -Z2 ≥ t 1 + t 2 ≤ exp -cρ 2 + exp (-cρ) (60)
for some universal constant c. In words, the probability that Z 1 + Z 2 deviates from Z1 + Z2 by more than ρ(σ y 0 2 + σ 2 ) decays exponentially fast in ρ. Now let ∈ R + and let t = φ( ) ∈ [0, 1]. Using the fact that both φ and φ -1 are increasing, we obtain

P ( x -x 2 ≥ ) = P x -x 2 ≥ φ -1 (t) (61)

Proof of Proposition 3.6

Proof. We take the convention of [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. Define sinc(x)

def.

= sin(πx)/(πx) and let F denote the Fourier transform on L 2 (R) defined by Notice that the normalization by 1/ √ a ensures that e L 2 (R) = 1 for all a > 0. The Fourier transform of e is given by: Using the Taylor expansion of the sinc, we obtain for all x ∈ R 1 -sinc(x) ≤ x 2 π 2 6 .

This allows to take φ(r) = r 2 π 2 3a 2 in Assumption 2.2, φ -1 (r) = a √ 3r π and

Π R(x) -Π R(x ) 2 2→2 ≤ 2(x -x ) 2 π 2 6a 2 implying L ≤ π √ 3a .
The fact that we observe a single source located at x1 ∈ D with weight w1 = 1 implies that y 0,1 = E(x 1 ) and that y 0

Now, to correctly sample the signal, we can set b = τ a for some τ ≤ 1. We obtain using -log(t) -log(t) for small t:

We have E( x) = E( x) + ∆ with ∆ 2→2 ≤ √ σ + L E xx 2 . Hence P ( x) = P ( x) + ∆ with

∆ 2→2 ≤ 2L E x -x 2 √ σ + E( x) 2→2 + σ + L 2 E x -x 2 2 .
(66)

Figure 1 :

 1 Figure 1: A sketch of the contribution: (a) a noisy image of the action of an unknown operator on a few Dirac masses, (b) detection of isolated spikes (in green) and of clusters (in red), (c) an operator estimate using only the green impulse responses, applied to a Dirac comb. The results in (b) and (c) were obtained using the algorithms proposed in this paper, see Section 4.3 for the technical details.

Assumption 1 . 3 (

 13 The signal structure). Let D ⊆ R D denote a domain, ( xn ) 1≤n≤N denote a collection of N points in D and ( wn ) 1≤n≤N denote N nonzero weights. We will work using one of the two models below. a) Independent sources We observe N independent sources (μ n ) of the form μn = wn δ xn . (4) b) Multiple sources We observe a single measure μ of the form μ = N n=1 wn δ xn ,

2 2 = 2 F

 22 u, u . For two matrices M 1 , M 2 in R M ×N , the notation M * 1 stands for the transpose (or trans-conjugate) of M 1 , M 1 , M 2 F = Tr(M * 1 M 2 ) denotes the Frobenius scalar product and M 1 = M 1 , M 1 F the Frobenius norm. The notation σ min (M 1 ) stands for smallest non-zero singular value of M 1 .

Theorem 3 . 3 (

 33 Stability of the Dirac locations). Let y 0,n = E( xn )γ denote the noiseless measurements and assume that b n 2 ≤ θ y 0,n 2 with θ < . Then, under Assumption 2.2, the following inequality holds:

√ 6 2 - 1 .

 61 Additional sources In this paragraph, we study what happens for a superposition of sources, i.e., if we replace the model 1.3 a) by 1.3 b). This leads to the following measurement vector y m def. = N n=1 wn ν m , A(γ)δ xn (20) with xn ∈ D for all n. Can we still recover a location -say x1 -if it is sufficiently distant from the others? In what follows, we let y = (y 1 , . . . , y m ) denote the complete measurement vector and y n = wn E( xn )α denote the measurement associated to the n-th source. Theorem 3.4 (Stability with spurious sources ). Let y = (y 1 , . . . , y m ) be generated according to the model (20) with N ≥ 2. Under the following hypotheses:

  e., set D = 1 and D = [-1, 1]. Let us define the normalized sinc function as sinc(x) def. = sin(πx)/(πx). Consider the following convolution kernel e(x) def. = 1 √ a sinc(x/a), for some scale parameter a > 0. (27) We assume that the family A from Assumption 1.1 is of dimension I = 1. The operator A 1 is a convolution with the filter e. The sampling model from Assumption 1.2 corresponds the usual point-wise sampling on a grid satisfying the Nyquist rate, i.e., ν m = δ zm with z m = bm and b ≤ a. Letting τ = b/a, the oversampling factor can therefore be defined as τ -1 . Assume that we observe a single source located at x1 ∈ D with weight w1 = 1. Proposition 3.6 (Stability for a single band-limited convolution kernel ). Under the above assumptions, set ρ > 1 and let = ca σ √ τ a(1 + ρ) + τ aσ 2 (-log(a) + ρ) . Then P(|x -x| ≥ ) ≤ exp (-c 1 ρ) .

Theorem 3 . 8 (

 38 Stability of the operator estimate with multiple observations ). Given X = (x 1 , . . . , x n ), let w-= min n | wn |, w+ = max n | wn | and C(X) def.

  def. = dim(R( xn )) and Î = N n=1 În . Using a singular value decomposition, we can decompose E( xn ) as

Figure 2 :

 2 Figure 2: The different families of convolution filters used in Section 4.1.

Figure 3 : 1 2 φ -φ 0 2 2

 3122 Figure 3: The corresponding inverse functions φ -1 i,+ (in red) for i ∈ {1, 2, 3}, for the different convolution systems.

Figure 4 :

 4 Figure 4: Examples of measurements vectors y for different noise levels and each family.

Figure 5 :Figure 6 :

 56 Figure 5: Average localization error E(|x-x|) as a fraction of a pixel for different noise levels θ ∈ [0, 1] and the three families A 1 , A 2 and A 3 .

Figure 7 :

 7 Figure 7: Examples of random product-convolution operators for a fixed family (e j ) and (f k ) and 3 random realizations of γ.

Figure 8 :

 8 Figure 8: Percentages of perfect recovery results for the different algorithms.

Lemma 5 . 1 .

 51 Let D ∈ N, S ⊆ R D , f : S → R and : S → R denote two mappings. Define g def.

5. 5 4 Proof.

 54 Proof of Theorem 3.Assume that μ = N n=1 wn δ xn . Then, we get

y 1 2 2 . 2 .

 122 Under Assumption 2.2, H 1 possesses a unique global maximizer at x1 . Let δ > 0 denote some radius and B δ/2 = {x ∈ R D , x -x1 2 ≤ δ/2}. Similarly to the previous proof, we now apply Lemma 5.1 with S = B δ/2 , f = H 1 , = 1 + 2 and ϕ(t) = 1 2 φ(t) y 1 2 Consider a point x ∈ B(δ/2). Assumption 2.2 and the hypothesis x1 -xn 2 ≥ δ for all n ≥ 2 allows to derive the following inequalities

Proposition 5 . 4 (

 54 Control of Z1 ). Under the assumptions of Theorem 3.5, we have, for all x

Theorem 5 . 6 (= c 1 d2→2

 561 Theorem 3.1 in [51]). Let S = {Π R(x) , x ∈ D} denote the family of projections. Let d 2→2 (S) = sup A∈S A 2→2 and d F (S) = sup A∈S A F . Let us define the following Dudley integral γ u, S) du,

t σ 2 = c 1 d2→2

 21 U , with E = γ 2 +γd F (S)+d F (S)d 2→2 (S), U = d 2 2→2 (S), V = d 2→2 (S) (γ + d F (S)) .and where c 1 , c 2 , c 3 are absolute constants. In what follows, the constants c, c 1 , c 2 , c 3 may change from line to line, but are always absolute. Notice that Ampl(∆ 2 ) ≤ 2 sup x∈D |∆ 2 (x) -E(∆ 2 (x))|. We have d F (S) = √ I since dim(R(x)) = I and d 2→2 (S) = 1 since the matrices Π R(x) are projections. The Lipschitz continuity assumption (23) also yields Θ(u, S) ≤ c L u D .To evaluate Dudley's integral, we can separate two cases. Let us first assume that cL > 1. In that case, we get γ

1 2 √ 1 √= c 1 d2→2V ≤ c 1 •

 12111 log(cL)is smaller than log(cL) and can be discarded, up to changing the multiplicative constant. Therefore, we getγ ≤ c D log(cL).Now, let us propose another bound covering the case cL ≤ 1 to justify the remark following the main theorem. Similarly to the proof of Proposition 5.4, let u s denote the smallest u > 0 such that Θ(u, S) ≥ 2. The previous inequality yields c L us D ≥ 2 and hence u s ≤ cL/2 1/D ≤ cL. Therefore, we get that γ L is larger than a universal constant c, we obtain: E ≤ c 1 • D log(L) + D log(L)I + √ I . D log(L) + √ I and U = 1.

=

  E (Ampl(∆ 2 )) ≤ c • (σ 2 E + σ 2 V + σU ) ≤ cσ 2 (D log(L) + D log(L)I + √ I).

y 0 2 2= c • σ y 0 2 + ρ + σ 2 D

 222 It remains to use Propositions 5.4 and 5.5 to obtainZ1 + Z2 + cρ(σ y 0 2 + σ 2 ) ≤ c σ 2 (D log(L) + DI log(L)) + L β α α+1• σ √ DI y 0 2 + ρ(σ y 0 2 + σ 2 ) log(L) + DI log(L) + ρ .

F

  exp(-ixω) dx.We recall thatF -1 (g)(x) = 1 2π R g(ω) exp(ixω) dω,and the Parseval identity readsu, v L 2 (R) = 1 2π F(u), F(v) L 2 (R) .

F

  (e)(ω) = R e(x) exp(-ixω) dx = √ a • 1 Ωa (ω),withΩ a = [-π/a, π/a].We recall the Shannon-Whittaker interpolation formula. For a given b ≤ a, let ψ m (t)

2 = 2 L 2 | 2 -

 2222 m ) ∈Z is an orthogonal basis of PW(Ω a ) for the usual scalar product on L 2 (R). In addition, for u ∈ PW(Ω a ),u, ψ m L 2 (R) = √ b • u(bm).Hence for any u ∈ PW(Ω a ), we get the Shannon interpolation formula:u(x) = m∈Z u, ψ m L 2 (R) ψ m (x) = √ b m∈Z u(mb)ψ m (x).Therefore, for any u, v ∈ PW(Ω a ), we obtainu, v L 2 (R) = b m∈Z u(mb)v(mb).The assumptions imply that E(x) is the infinite vector in 2 defined by E(x) = (e(x -z m )) m∈Z . Using Parseval identity, we get for all x, x ∈ R thatΠ R(x) Π R(x ) 2→2 ≤ | E(x), E(x ) 2 | E(x) 2 E(x ) 2 = | R e(t -x)e(t -x ) dt| e(• -x) L 2 (R) e(• -x ) L 2 (R) = 1 a R sinc((t -x)/a)sinc((t -x )/a) dt = 1 a • 2π R F(sinc((• -x)/a))(ω)F(sinc((• -x )/a))(ω) dω x |/a) .Hence, the Assumption of Theorem 3.5 holds with α = 1 and β = 1/a. In addition, we haveΠ R(x) -Π R(x ) e(• -x) -e(• -x ) -ix•) -exp(-ix •)]F(e) exp(-ixω) -exp(-ix ω)| 2 dω 2 cos((x -x)ω) dω = a π ω -sin((x -x)ω) x -x π/a -π/a = 2 [1 -sinc((x -x)/a)]

2 2 = 1 b

 21 . Plugging these expressions in and ignoring the constants, we get that = φ -1 c σ √ b(1 + ρ) + bσ 2 (-log(a) + -log(a) + ρ) .

= φ - 1

 1 c σ √ τ a(1 + ρ) + τ aσ 2 (-log(a) + ρ) = ca σ √ τ a(1 + ρ) + τ aσ 2 (-log(a) + ρ) .

5. 8

 8 Proof of Theorem 3.7 Proof. Let P (x) def. = E * (x)E(x). By definition, we have γ = P ( x) -1 E * ( x)(y 0,1 + b 1 )

  = P (φ( xx 2 ) ≥ t)As previously, set t 1 = cρσ y 0 2 , t 2 = cσ 2 ρ and

		P rop.5.2 ≤ P	2Ampl(∆ 1 + ∆ 2 ) 2 y 0 2	≥ t
		≤ P Z 1 + Z 2 ≥	t y 0 2	2 2	.	(62)
	t y 0 2	2 2	= Z1 + Z2 + t

1 + t 2 .

Combining (60) and (62), we obtain

P ( xx 2 ≥ ) ≤ exp -cρ 2 + exp (-cρ) for = φ -1 2 Z1 + Z2 + t 1 + t 2

The orthogonality is not a strong assumption, since any family can be orthogonalized.

, 2. The highest multiplicative constant is 24, which appears in Dudley's inequality.2. In this statement, we work over the unit ball. However, any compact domain would yield a similar result, up to multiplicative factors of the diameter. The important thing is to work over a compact domain to restrict the family of possible ranges {R(x), x ∈ D}.

The scaling in √ σ + is natural considering Assumption 2.1.

Acknowledgments

The authors wish to thank Kristian Bredies for his careful analysis of the paper, which allowed to correct some flaws in the concentration results. They thank Emmanuel Soubies and Axel Flinth for their comments and advices on the manuscript. They also wish to thank Landry Duguet for his initial exploration of the problem during a one month internship. P. Weiss thanks T. Rezk for fruitful discussions.

This work was supported by Fondation pour la Recherche Médicale (FRM grant number ECO20170637521 to V. Debarnot), by the ANR project Optimization on Measures Spaces ANR-17-CE23-0013-01, from ANR-3IA Artificial and Natural Intelligence Toulouse Institute and from the ANR Micro-blind.

Proof. We show that for all (x, γ), (x , γ ) ∈ R D × R I \{0} E(x)γ = E(x )γ ⇒ γ = γ .

If x = x : by Assumption 2.1 the mapping γ → E(x)γ is injective, meaning that for all γ, γ ∈ R I \{0} we have

If x = x : by Assumption 2.2 we have R(x) ∩ R(x ) = {0}, which implies that for any γ, γ ∈ R I \{0}, we have E(x)γ = E(x )γ . We prove it by contradiction. Assume that R(x) ∩ R(x ) = {0}. Then, there ex-

where we used [START_REF] Bredies | Inverse problems in spaces of measures[END_REF] to obtain the last equality.

The linear system to recover γ can be rewritten as

with δ = ∆(y 0,1 +b 1 )+E * ( x)b 1 . Under Assumption 2.1, the unique solution of P ( x)γ = E * ( x)y 0,1 is γ. If x is sufficiently close to x, we have ∆ 2→2 < σ - and P ( x) -1 ∆ 2→2 < 1. We can now use standard results of linear algebra, see e.g. [52, 3.6], to obtain that

.

Letting x → x, we see that the noise term ccontains an error term 1 defined by

which does not vanish as x → x and an additional term 2 ( x) which does vanish.

Using the fact that

2 ), we obtain

for some absolute constant c.

Proof of Theorem 3.8

The proof is nearly identical to the previous one and we omit it for brevity.