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

HAL
open science

# Blind inverse problems with isolated spikes

Valentin Debarnot, Pierre Weiss

► To cite this version:

Valentin Debarnot, Pierre Weiss. Blind inverse problems with isolated spikes. Information and Inference, 2022, 10.1093/imaiai/iaac015. hal-03272066v2

HAL Id: hal-03272066

https://inria.hal.science/hal-03272066v2

Submitted on 2 Nov 2021

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

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

Valentin Debarnot & Pierre Weiss

November 2, 2021

## Abstract

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.

## 1 Introduction

<figure><img src="image_2.png" /><figcaption>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.</figcaption></figure>

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 [1, 2].
From a mathematical viewpoint, this problem can be modelled as follows. Let $\mathcal{M}(\mathbb{R}^D)$ denote the set of Radon measures, i.e. the dual of the set $C_0^0(\mathbb{R}^D)$ of continuous functions vanishing at infinity. Let $\bar{\mu} = \sum_{n=1}^N \bar{w}_n \delta_{\bar{x}_n} \in \mathcal{M}(\mathbb{R}^D)$ denote a Radon measure that encodes the protein locations $(\bar{x}_n)$ and their intensities $(\bar{w}_n)$. Assume that this measure is observed indirectly through a linear regularizing operator $\bar{A} : \mathcal{M}(\mathbb{R}^D) \to C_0^0(\mathbb{R}^D)$:

$$y_m = \bar{A}\bar{\mu}(z_m) + b_m, \qquad (1)$$

where $y \in \mathbb{R}^M$ is the observed data, $Z = (z_1, z_2, \dots, z_M)$ denotes a set of sampling locations in $\mathbb{R}^D$ and $b = (b_1, \dots, b_M) \in \mathbb{R}^M$ is some additive noise. The operator $\bar{A}$ 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 developed over the years to recover the positions $(\bar{x}_n)$ from the measurements $y$. We refer the interested reader to the summaries of the super-resolution challenges [3, 4] 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 $(\bar{w}_n)$ 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 $\bar{A}$ is often only partially known, making it crucial to estimate both the positions and weights $(\bar{w}_n, \bar{x}_n)$, but also the operator $\bar{A}$ 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 $\bar{A}$ lives in a known finite dimensional subspace $\mathcal{A}$ of linear operators from $\mathcal{M}(\mathbb{R}^D)$ to $C_0^0(\mathbb{R}^D)$.

Throughout the paper, we let $I \in \mathbb{N}$ denote the subspace dimension and set $\mathcal{A} = \text{span}(A_1, \dots, A_I)$. Any $A \in \mathcal{A}$ can therefore be parameterized by a vector $\gamma = (\gamma_i) \in \mathbb{R}^I$, with

$$A\mu = A(\gamma)\mu \stackrel{\text{def.}}{=} \sum_{i=1}^{I} \gamma_i A_i \mu \quad \text{for any } \mu \in \mathcal{M}(\mathbb{R}^D). \qquad (2)$$

The next assumption describes the sampling model considered in this work.
**Assumption 1.2** (The observation model). Let $(\nu_m)_{1\le m\le M}$ in $\mathcal{M}(\mathbb{R}^D)$ denote a collection of $M$ linear forms on $C_0^0(\mathbb{R}^D)$. Let $\bar{\mu} \in \mathcal{M}(\mathbb{R}^D)$ denote a signal to re-construct. We assume that we acquire the measurement vector $\mathbf{y} = (y_1, \dots, y_M)$ with coordinates given by

$$y_m \stackrel{\text{def.}}{=} \langle \nu_m, A(\bar{\gamma})\bar{\mu} \rangle + b_m. \quad (3)$$

The observation model 1.2 allows to describe nearly any sampling device. For instance, the traditional pointwise sampling would consist in choosing $\nu_m = \delta_{z_m}$, where $(z_m)_{1\le m\le M}$ is a set of sampling locations. Fourier sampling can be modelled by setting $\nu_m(\omega) = \exp(-i\langle \omega, z_m \rangle)$.

Finally, throughout the paper, we will assume that the measure $\bar{\mu}$ is discrete.

**Assumption 1.3** (The signal structure). Let $\mathcal{D} \subseteq \mathbb{R}^D$ denote a domain, $(\bar{x}_n)_{1\le n\le N}$ denote a collection of $N$ points in $\mathcal{D}$ and $(\bar{w}_n)_{1\le n\le N}$ denote $N$ nonzero weights. We will work using one of the two models below.

a) **Independent sources** We observe $N$ independent sources $(\bar{\mu}_n)$ of the form

$$\bar{\mu}_n = \bar{w}_n \delta_{\bar{x}_n}. \quad (4)$$

b) **Multiple sources** We observe a single measure $\bar{\mu}$ of the form

$$\bar{\mu} = \sum_{n=1}^{N} \bar{w}_n \delta_{\bar{x}_n}, \quad (5)$$

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

## 1.1 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 $\mathcal{A}$, beyond convolutions.

* We work under a general linear sampling model with arbitrary linear forms $(\nu_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. [5, 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 \in \mathbb{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 $\tilde{A} \in \mathcal{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 $\mathcal{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 space-varying operators, we recently proposed an efficient calibration technique in [6, 7], leading to low-dimensional subspaces of unstationary 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 $\tilde{A}$ of the operator $\tilde{A}$. Once it is estimated, we can use existing inverse problem solvers to recover any measure $\mu \in \mathcal{M}(\mathbb{R}^D)$ from its measurements $\tilde{A}\mu$. 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.

## 1.2 Related works

When the operator $\tilde{A}$ is known, recovering $\tilde{\mu}$ is already a challenging problem, since the inverse problem is ill-posed and infinite dimensional. Specifying prior assumptions on the signal $\tilde{\mu}$ to certify its approximate recovery is essential [8]. A few mathematical breakthroughs were achieved in the recent past.
**Off-the-grid total variation minimization with a known operator** In [9, 10, 5], 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 [11], the authors showed that the separation is not needed, provided that the weights ($\bar{w}_n$) 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 [12, 13, 14]. This type of approach is currently amongst the best competitors when a high density of proteins is used [4].

**Gridded lifting for an unknown operator** Assume that the operator $\bar{A}$ is unknown but lives in a known finite dimensional subspace $\mathcal{A}$. Also assume that the positions ($\bar{x}_n$) are known, but that the weights ($\bar{w}_n$) are unknown. Under these hypotheses, an elegant solution to recover $\bar{A}$ and $\bar{\mu}$ was proposed by Ahmed et al in [15] based on a trick called lifting. This approach allows to tackle bilinear problems of the form

$$ \inf_{A \in \mathcal{A}, u \in \mathcal{U}} \frac{1}{2} \|Au - y\|_2^2, \qquad (6) $$

where $\mathcal{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 ($\bar{w}_n$) and $\bar{A}$ under rather stringent assumptions. The assumptions were relaxed in a series of works [15, 16, 17, 18, 19]. 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 [17], Y. Chi showed that the lifting trick could also be used when the unknown positions ($\bar{x}_n$) live off-the-grid, $D = 1$ and the operators in $\mathcal{A}$ are convolution operators. The approach was then extended to the 2 dimensional setting for convolution operators in [20]. In [21], 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 [22, 23], where multiple signals are convolved with a single filter. This leads to observations $y_n$ of the form $y_n = e \star 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 [24, 19, 25, 26]. 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.

# 2 Preliminaries

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

## 2.1 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. $\bar{\mu}$, $\bar{w}_n$, $\bar{x}_n$). A hat indicates that it is an estimate (e.g. $\hat{\mu}$, $\hat{w}_n$, $\hat{x}_n$). The symbols $D, I, J, K, M, N$ denote cardinalities in $\mathbb{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 $\mathbf{x} \in \mathbb{R}^D$, $x_d$ denotes the $d$-th coefficient of $\mathbf{x}$. The parentheses are used to evaluate functions. For instance, given $f: \mathbb{R}^D \to \mathbb{R}$, $f(\mathbf{x})$ is the value of $f$ at $\mathbf{x}$.

For $\mu \in \mathcal{M}(\mathbb{R}^D)$ and $u \in C_0^0(\mathbb{R}^D)$, we let $\langle \mu, u \rangle \in \mathbb{R}$ denote the value of the linear form $\mu$ on $u$. We also let $\mu \star u$ denote the convolution product between $\mu$ and $u$ defined for all $\mathbf{x} \in \mathbb{R}^D$ by $(\mu \star u)(\mathbf{x}) = \langle \mu, u(\mathbf{x} - \cdot) \rangle$.

In all the paper, the notation $\langle \cdot, \cdot \rangle$ also refers to the usual scalar product on the vector space $\mathbb{R}^N$, where $N \in \mathbb{N}$. For $\mathbf{u} \in \mathbb{R}^N$, $\|\mathbf{u}\|_2$ denotes the $l^2$-norm of $\mathbf{u}$ defined by $\|\mathbf{u}\|_2^2 = \langle \mathbf{u}, \mathbf{u} \rangle$. For two matrices $M_1, M_2$ in $\mathbb{R}^{M \times N}$, the notation $M_1^*$ stands for the transpose (or trans-conjugate) of $M_1$, $\langle M_1, M_2 \rangle_F = \text{Tr}(M_1^* M_2)$ denotes the Frobenius scalar product and $\|M_1\|_F^2 = \langle M_1, M_1 \rangle_F$ the Frobenius norm. The notation $\sigma_{\min}(M_1)$ stands for smallest non-zero singular value of $M_1$.

We let $\langle \cdot, \cdot \rangle_{L^2(\mathbb{R}^D)}$ denote the usual scalar product of $L^2(\mathbb{R}^D)$. For a compact and symmetric set $\Omega \subset \mathbb{R}^D$, we let $\mathcal{P}W(\Omega)$ denote the Paley-Wiener set of band-limited functions on $\Omega$, i.e. the set of functions in $L^2(\mathbb{R}^D)$ that have a Fourier transform that vanishes outside $\Omega$.

## 2.2 Further assumptions

Under Assumptions 1.1 and 1.2, the impulse response of an operator $A(\gamma)$ at a location $z \in \mathbb{R}^D$ is given by $A(\gamma)\delta_z = \sum_{i=1}^I \gamma_i A_i \delta_z$. This motivates introducing the matrix-valued function $\mathbf{E}: \mathbb{R}^D \to \mathbb{R}^{M \times I}$ defined by

$$ (\mathbf{E}(z))_{m,i} \stackrel{\text{def.}}{=} \langle \nu_m, A_i \delta_z \rangle. \qquad (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 $\mathbf{x} \in \mathcal{D}$, the mapping $\mathbf{E}(\mathbf{x}): \mathbb{R}^I \to \mathbb{R}^M$ is injective: we have

$$ \sigma_{-} \mathbf{Id} \preceq \mathbf{E}^*(\mathbf{x}) \mathbf{E}(\mathbf{x}) \preceq \sigma_{+} \mathbf{Id} \qquad (8) $$

with $0 < \sigma_- \le \sigma_+ < +\infty$. In what follows, we let $\kappa \stackrel{\text{def.}}{=} \frac{\sigma_+}{\sigma_-}$.
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

$$ \mathcal{R}(x) \stackrel{\text{def.}}{=} \text{Ran}(\mathcal{E}(x)) \qquad (9) $$

denote the subspace of possible measurements for an impulse response located at $x \in \mathcal{D}$ and $\Pi_{\mathcal{R}(x)}$ denote the orthogonal projection onto the range $\mathcal{R}(x)$. Another important technical assumption is the following.

**Assumption 2.2** (Identifiability of the Dirac masses location). The mapping $\mathcal{E}$ satisfies the following inequality for any pair $x, \bar{x} \in \mathcal{D}$

$$ \|\Pi_{\mathcal{R}(x)} \Pi_{\mathcal{R}(\bar{x})}\|_{2 \to 2} \le 1 - \phi(\|x - \bar{x}\|_2) \qquad (10) $$

for some non-decreasing function $\phi: \mathbb{R}_+ \to [0, 1]$ with $\phi(0) = 0$ and $\phi(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 \ne \bar{x}$ in $\mathcal{D}$. Then, the two ranges $\mathcal{R}(x)$ and $\mathcal{R}(\bar{x})$ do not contain two identical elements. Hence, the knowledge of a measurement of the form $\mathcal{E}^*(\bar{x})\bar{\gamma}$ should be enough to perfectly recover $\bar{x}$. An additional implicit assumption is that the range $\mathcal{R}(x)$ cannot be the singleton $\{0\}$ for any $x$. Indeed, taking $x = \bar{x}$ implies that $\|\Pi_{\mathcal{R}(x)}\|_{2 \to 2} = 1$.

## 2.3 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, \gamma) \mapsto \mathcal{E}(x)\gamma$ is injective on $(\mathcal{D} \times \mathbb{R}^I \setminus \{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\delta_x, 1 \le i \le 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. [27].

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

$$ \begin{aligned} \langle \mathcal{E}(x)\gamma, \mathcal{E}(\bar{x})\bar{\gamma} \rangle &= \langle \Pi_{\mathcal{R}(x)} \mathcal{E}(x)\gamma, \Pi_{\mathcal{R}(\bar{x})} \mathcal{E}(\bar{x})\bar{\gamma} \rangle \\ &= \langle \Pi_{\mathcal{R}(\bar{x})} \Pi_{\mathcal{R}(x)} \mathcal{E}(x)\gamma, \mathcal{E}(\bar{x})\bar{\gamma} \rangle \\ &\le \|\Pi_{\mathcal{R}(\bar{x})} \Pi_{\mathcal{R}(x)}\|_{2 \to 2} \|\mathcal{E}(x)\gamma\|_2 \|\mathcal{E}(\bar{x})\bar{\gamma}\|_2 \\ &\le [1 - \phi(\|x - \bar{x}\|_2)] \|\mathcal{E}(x)\gamma\|_2 \|\mathcal{E}(\bar{x})\bar{\gamma}\|_2 \end{aligned} $$
**A geometric condition** The quantity $||\Pi_{\mathcal{R}(x)}\Pi_{\mathcal{R}(\bar{x})}||_{2\to2}$ is related to the principal angle between the subspaces $\mathcal{R}(x)$ and $\mathcal{R}(\bar{x})$. To realize this, let us recall that the principal angle $\angle(\mathcal{U}, \mathcal{V})$ between two subspaces $\mathcal{U}$ and $\mathcal{V}$ of a Hilbert space with norm $||\cdot||$ is defined by

$$ \cos(\angle(\mathcal{U}, \mathcal{V})) = \max_{\substack{u \in \mathcal{U}, v \in \mathcal{V} \\ u \neq 0, v \neq 0}} \frac{\langle u, v \rangle}{||u|| \cdot ||v||}. \quad (11) $$

We have:

$$ \begin{aligned} ||\Pi_{\mathcal{R}(x)}\Pi_{\mathcal{R}(\bar{x})}||_{2\to2} &= \sup_{\substack{\mathbf{u}, \mathbf{v} \in \mathbb{R}^M \\ ||\mathbf{u}||_2 = ||\mathbf{v}||_2 = 1}} \langle \Pi_{\mathcal{R}(x)}\Pi_{\mathcal{R}(\bar{x})}\mathbf{u}, \mathbf{v} \rangle \\ &= \sup_{\substack{\mathbf{u} \in \mathcal{R}(x), \mathbf{v} \in \mathcal{R}(x) \\ ||\mathbf{u}||_2 = ||\mathbf{v}||_2 = 1}} \langle \mathbf{u}, \mathbf{v} \rangle = \cos(\angle(\mathcal{R}(x), \mathcal{R}(\bar{x}))). \end{aligned} $$

## 2.4 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 orthogonal <fnref n="1" /> family $(e_i)_{1\le i\le I}$ of functions in $\mathcal{PW}(\Omega)$ that vanish at infinity. The operators $A_i$ are convolutions with the filters $e_i$, i.e. $A_i\mu = e_i \star \mu$ for $\mu \in \mathcal{M}(\mathbb{R}^D)$.

The linear forms $\nu_m$ describing the sampling device correspond to a Shannon sampler, i.e. $\nu_m = \delta_{z_m}$, where the positions $z_m$ correspond to a Cartesian grid with a grid-size smaller than $\frac{2\pi}{\text{diam}(\Omega)}$.

Under Assumption 2.3, any $(u, v) \in \mathcal{PW}(\Omega)^2$ satisfy

$$ \langle u, v \rangle_{L^2(\mathbb{R}^D)} \propto \sum_{m \in \mathbb{N}} u(z_m)v(z_m), \quad (12) $$

which is a variant of the Shannon-Nyquist theorem, see e.g. [28, Thm 3.5].

**Proposition 2.2** (Operator identifiability for convolution operators). Under Assumption (2.3), we have $E^*(z)E(z) = \text{Id}$, hence Assumption 2.1 is satisfied with $\sigma_- = \sigma_+ = 1$.

**Proposition 2.3** (Location identifiability for convolution operators). Let $\mathcal{C}: \mathbb{R}^D \to \mathbb{R}^{I \times I}$ denote the following cross-correlation matrix-valued function:

$$ [\mathcal{C}(x - x')]_{i,i'} \stackrel{\text{def.}}{=} \langle e_i(\cdot - x), e_{i'}(\cdot - x') \rangle_{L^2(\mathbb{R}^D)} \quad (13) $$

Under Assumption (2.3), we have

$$ ||\Pi_{\mathcal{R}(x)}\Pi_{\mathcal{R}(\bar{x})}||_{2\to2} = ||\mathcal{C}(x - \bar{x})||_{2\to2}. \quad (14) $$

<footnote n="1">¹The orthogonality is not a strong assumption, since any family can be orthogonalized.</footnote>
Proposition 2.3 shows that the condition (10) 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 (10) 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.

## 2.5 The case of product-convolution operators

To encode space varying operators, we now turn to product-convolution expansions [29]. 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 \le j \le J}$ of band-limited functions in $PW(\Omega)$, and another orthogonal family $(f_k)_{1 \le k \le K}$ of functions in $L^2(\mathbb{R}^D) \cap C_0^0(\mathbb{R}^D)$.

The family of observation operators $\mathcal{A}$ is a subspace of product-convolution expansions from $\mathcal{M}(\mathbb{R}^D)$ to $C_0^0(\mathbb{R}^D)$ defined as follows. For any $A \in \mathcal{A}$, there exists a vector $\gamma = (\gamma_{j,k}) \in \mathbb{R}^{J \times K}$ such that for any $\mu \in \mathcal{M}(\mathbb{R}^D)$:

$$A\mu = A(\gamma)\mu \stackrel{\text{def.}}{=} \sum_{j=1}^{J} \sum_{k=1}^{K} \gamma_{j,k} e_j \star (f_k \odot \mu). \quad (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 [30]. In addition, we recently showed how a subspace of product-convolution operators $\mathcal{A}$ could be constructed in practice in optical imaging [31, 6, 7].

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

$$\| \Pi_{\mathcal{R}(x)} \Pi_{\mathcal{R}(x')} \|_{2 \to 2} = \| C(x - x') \|_{2 \to 2}. \quad (16)$$

However, for $K \ge 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 \ge 2$ is possible only under the condition $N \ge K$, i.e. by observing multiple impulse responses.

## 3 Main results

Throughout this section, we will work under Assumptions 1.1, 1.2 and 1.3 a). For all $1 \le n \le N$, the observation $y_n$ of a single Dirac mass at $\bar{x}_n$ can be written as

$$y_n = E(\bar{x}_n)\bar{\alpha}_n + b_n, \quad (17)$$
where $\bar{x}_n \in \mathbb{R}^D$ is an unknown location and $\bar{\alpha}_n = \bar{w}_n \bar{\gamma} \in \mathbb{R}^I$ is a vector co-linear to the unknown operator parameterization $\bar{\gamma}$. The overall objective of this paper is to construct an estimate $\hat{X} = (\hat{x}_1, \dots, \hat{x}_N)$ of $\hat{X} = (\bar{x}_1, \dots, \bar{x}_N)$ and $\hat{\gamma}$ of $\bar{\gamma}$ and to certify their proximity, despite the perturbation term $b_n \in \mathbb{R}^M$. The case of multiple sources in Model 1.3 b) will be treated as a particular case of (17), where $b_n$ coincides with the measurements of the $N-1$ sources different from $n$.

## 3.1 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 $\bar{x}_n$ by finding a global minimizer $\hat{x}_n$ of the following problem

$$ \inf_{x \in \mathcal{D}, \alpha \in \mathbb{R}^I} \frac{1}{2} \| \mathcal{E}(x)\alpha - y_n \|_2^2. \quad (\mathcal{P}_n^x) $$

**Step 2: recovering the operator** Second, depending on whether the weights $\bar{w}$ are known or unknown, we propose two different recovery strategies.

**Case 1: known weights** If $\bar{w}$ is known, we consider the following quadratic problem:

$$ \inf_{\gamma \in \mathbb{R}^I} \frac{1}{2} \sum_{n=1}^{N} \| \bar{w}_n \mathcal{E}(\hat{x}_n) \gamma - y_n \|_2^2, \quad (\mathcal{P}^\gamma) $$

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 $\hat{x}_n$.

**Case 2: unknown weights** If the true weights $\bar{w}$ are unknown, let

$$ J(w, \gamma, \hat{X}) \stackrel{\text{def.}}{=} \frac{1}{2} \sum_{n=1}^{N} \| w_n \mathcal{E}(\hat{x}_n) \gamma - y_n \|_2^2. $$

We propose to solve the bilinear problem below

$$ \inf_{\substack{\gamma \in \mathbb{R}^I \\ w \in \mathbb{R}^N}} J(w, \gamma, \hat{X}). \quad (\mathcal{P}^{\gamma, w}) $$

**Additional regularization terms** In the proposed formulation $(\mathcal{P}_n^x)$, two implicit regularization terms are used: i) we look for a single Dirac mass located in $\mathcal{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, \gamma)$ of positions and operators could lead to similar measurements. A possible solution to leverage this problem is to add a weighted $l^2$-regularization on $\gamma$ of the form $\frac{1}{2} \sum_{i=1}^I \theta_i \gamma_i^2$, where $\theta_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 $(\mathcal{P}^{\gamma,w})$ suffers from the usual scaling ambiguity in bilinear inverse problem: if $(\hat{w}, \hat{\gamma})$ is a solution, then so is $(t\hat{w}, \frac{\hat{\gamma}}{t})$ for any $t \neq 0$. Without any further normalization, the weights $\bar{w}$ and parameterization $\bar{\gamma}$ 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 $\langle \gamma, 1 \rangle = 1$.

## 3.2 Estimating the locations

In this section, we establish a few results regarding the location estimate obtained by solving problem $(\mathcal{P}_n^x)$. 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.

### 3.2.1 A simple characterization of the locations

Let us start with an elementary observation.

**Proposition 3.1 (A simple characterization).** Assume that the problem $(\mathcal{P}_n^x)$ admits a solution $\hat{x}_n$. Then

$$ \hat{x}_n \in \underset{x \in \mathcal{D}}{\operatorname{argmax}} \frac{1}{2} \| \Pi_{\mathcal{R}(x)} y_n \|_2^2. \quad (18) $$

Notice that the existence of a minimizer is not automatic since the minimization domain $\mathcal{D}$ might be unbounded and the function $x \mapsto \frac{1}{2}\|\Pi_{\mathcal{R}(x)}y_n\|_2^2$ might be discontinuous. In particular, singularities may appear if the dimension of the range $\mathcal{R}(x)$ varies over $\mathcal{D}$. A sufficient condition for existence is given below.

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

* The domain $\mathcal{D}$ is compact

* The mapping $\mathcal{E}$ is continuous over $\mathcal{D}$ and $\mathcal{E}(x)$ is of rank 1 for all $x \in \mathcal{D}$. Then Problem $(\mathcal{P}_n^x)$ admits at least one solution $\hat{x}_n$.

### 3.2.2 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 structure to the noise term $b_n$.

**Theorem 3.3 (Stability of the Dirac locations).** Let $y_{0,n} = \mathcal{E}(\bar{x}_n)\bar{\gamma}$ denote the noiseless measurements and assume that $\|b_n\|_2 \le \theta\|y_{0,n}\|_2$ with $\theta < \frac{\sqrt{6}}{2} - 1 \simeq 0.225$. Then, under Assumption 2.2, the following inequality holds:

$$ \|\dot{x}_n - \bar{x}_n\|_2 \le \phi_+^{-1}(2\theta^2 + 4\theta), \quad (19) $$

where $\phi_+^{-1}(t) = \inf\{s \text{ s.t. } \phi(s) \ge t\}$ is the quantile function of $\phi$ (in particular $\phi_+^{-1} = \phi^{-1}$ if $\phi$ is bijective).
The above result is in sharp contrast with the asymptotic results available for multiple source points [10, 5]. Indeed, we see that the noise level can be of the order of the signal's level and the bound (19) 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 $\frac{\sqrt{6}}{2} - 1$.

**Additional sources** In this paragraph, we study what happens for a super-position 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 \stackrel{\text{def.}}{=} \sum_{n=1}^{N} \bar{w}_n \langle \nu_m, A(\bar{\gamma}) \delta_{\bar{x}_n} \rangle \quad (20)$$

with $\bar{x}_n \in \mathcal{D}$ for all $n$. Can we still recover a location - say $\bar{x}_1$ - if it is sufficiently distant from the others? In what follows, we let $\mathbf{y} = (y_1, \dots, y_m)$ denote the complete measurement vector and $\mathbf{y}_n = \bar{w}_n E(\bar{x}_n) \alpha$ denote the measurement associated to the $n$-th source.

**Theorem 3.4** (Stability with spurious sources). Let $\mathbf{y} = (y_1, \dots, y_m)$ be generated according to the model (20) with $N \ge 2$. Under the following hypotheses:

* Assumptions 2.2 is satisfied with

$$\phi(t) = \frac{(t/a)^b}{1 + (t/a)^b} \quad \text{for some } a > 0 \text{ and } b > 0. \quad (21)$$

* We have $\|\mathbf{y}_n\|_2 \le c \|\mathbf{y}_1\|_2$ for all $n \ge 2$ and some constant $c > 0$.

* The first source location $\bar{x}_1$ is isolated from the others:

$$\|\bar{x}_1 - \bar{x}_n\|_2 \ge \delta \text{ for } n \ge 2, \delta \ge 2a(\tau(N-1)c)^{1/b} \text{ and } \tau \ge 5.$$

Then any global maximizer of $H \stackrel{\text{def.}}{=} \frac{1}{2} \|\Pi_{\mathcal{R}(x)} \mathbf{y}\|_2^2$ in the ball $\mathcal{B}_{\delta/2} = \{\mathbf{x}, \|\mathbf{x} - \bar{x}_1\|_2 \le \delta/2\}$ lies within the ball $\mathcal{B}_r = \{\mathbf{x}, \|\mathbf{x} - \bar{x}_1\|_2 \le r\}$ with

$$r \le 18.4 \cdot \frac{a \cdot 2^{1/b}}{\tau} \quad (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/\tau$ if the distance of the surrounding molecules is higher than $a(\tau Nc)^{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 $\bar{w}_1 \ge c' \bar{w}_n$ for some constant $c' > 0$ and all $n \ge 2$.
• Asymptotically in $N$, the factor 18.4 in (22) can be replaced by 1. The following bound emerges from the proof:

$$r \leq \frac{a \cdot 2^{1/b}}{\tau} \cdot \left[ \frac{(\tau^2 + \tau)(N-1)c + 2\tau}{1 + (\tau^2 - 4\tau - 2)(N-1)c + (2\tau - 4)(N-1)c} \right].$$

• 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.5.** 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 \sim N(0, \sigma^2 \mathbf{I}_d)$.

• The minimization domain is the unit ball $\mathcal{D} = \{\mathbf{x}, \|\mathbf{x}\|_2 \le 1\}$.

• The mapping $\mathbf{x} \mapsto \Pi_{\mathcal{R}(\mathbf{x})}$ is $L$-Lipschitz continuous with $L \ge c_1$:

$$\|\Pi_{\mathcal{R}(\mathbf{x})} - \Pi_{\mathcal{R}(\mathbf{x}')}\|_{2 \to 2} \le L \|\mathbf{x} - \mathbf{x}'\|_2 \text{ for } \mathbf{x}, \mathbf{x}' \in \mathcal{D}. \quad (23)$$

• The following inequality holds for $\mathbf{x}, \mathbf{x}' \in \mathcal{D}$:

$$\|\Pi_{\mathcal{R}(\mathbf{x})} \Pi_{\mathcal{R}(\mathbf{x}')}\|_{2 \to 2} \le \frac{1}{(\beta \|\mathbf{x} - \mathbf{x}'\|_2)^\alpha} \quad \text{with } \beta > 0 \text{ and } \alpha > 1/2. \quad (24)$$

For any $\rho > 0$, set

$$\epsilon = \phi^{-1} \left( c_2 \cdot \frac{\sigma}{\|\mathbf{y}_{0,n}\|_2} \left[ \left( \frac{L}{\beta} \right)^{\frac{\alpha}{\alpha+1}} \cdot \sqrt{DI} + \rho \right] + \frac{\sigma^2}{\|\mathbf{y}_{0,n}\|_2^2} \left[ D \log(L) + \sqrt{DI \log(L)} + \rho \right] \right). \quad (25)$$

Then

$$\mathbb{P}(\|\hat{\mathbf{x}}_n - \bar{\mathbf{x}}_n\|_2 \ge \epsilon) \le \exp(-c_3\rho) + \exp(-c_4\rho^2). \quad (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 $\frac{\sqrt{\pi}}{2}$, 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 $\{\mathcal{R}(\mathbf{x}), \mathbf{x} \in \mathcal{D}\}$.
3. The condition $L \ge c_1$ is not necessary. In fact, stronger statements can be obtained for $L \le c_1$ since we reach a different regime: for small $L$ the projection matrices $\Pi_{\mathbb{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.

### 3.2.3 Application to a simple example

To illustrate Theorem 3.5, we provide a minimalist application below. We work in the 1D setting, i.e., set $D = 1$ and $\mathcal{D} = [-1, 1]$. Let us define the normalized sinc function as $\text{sinc}(x) \stackrel{\text{def.}}{=} \sin(\pi x)/(\pi x)$. Consider the following convolution kernel

$$e(x) \stackrel{\text{def.}}{=} \frac{1}{\sqrt{a}} \text{sinc}(x/a), \text{ for some scale parameter } a > 0. \quad (27)$$

We assume that the family $\mathcal{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., $\nu_m = \delta_{z_m}$ with $z_m = bm$ and $b \le a$. Letting $\tau = b/a$, the oversampling factor can therefore be defined as $\tau^{-1}$. Assume that we observe a single source located at $\bar{x}_1 \in \mathcal{D}$ with weight $\bar{w}_1 = 1$.

**Proposition 3.6 (Stability for a single band-limited convolution kernel).** *Under the above assumptions, set $\rho > 1$ and let*

$$\epsilon = ca\sqrt{(\sigma\sqrt{\tau a}(1 + \rho) + \tau a\sigma^2(-\log(a) + \rho))}.$$

Then

$$\mathbb{P}(|\bar{x} - \bar{x}_1| \ge \epsilon) \le \exp(-c_1\rho).$$

To obtain a super-resolution effect, we want the precision $\epsilon$ to be smaller than Shannon's rate, i.e., $\epsilon \lesssim a$. This can be obtained by setting $\sigma \lesssim \frac{1}{\sqrt{\tau 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 $\sigma$ as $\frac{1}{\sqrt{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 $\mathcal{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 $\mathbb{R}^M$ with mean 0 and covariance $\mathbf{I}$d 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 $\tau^{-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, the $ℓ^2$-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 $\frac{\sigma^2}{a}$. The condition $\|b_n\|_2^2 \lesssim \|y_{0,n}\|_2^2$ would therefore translate to $\sigma \lesssim 1$. In comparison, Theorem (3.5) allows for a scaling as $\sigma \lesssim \frac{1}{a^{1/2}\log(a)}$, which goes to infinity as $a$ goes to 0!

### 3.3 Estimating the operator with known weights

In this section, we study the problem of estimating the operator parameterization $\tilde{\gamma}$ under the assumption that the weights are $\tilde{w}_n$ are known. Once the positions $\hat{x}_n$ have been estimated for every observation $y_n$, the vector $\hat{\gamma}$ can be estimated by solving $\mathcal{P}\gamma$. This is also equivalent to the following linear system:

$$ \left( \sum_{n=1}^{N} \tilde{w}_n^2 \mathbf{E}^*(\hat{x}_n) \mathbf{E}(\hat{x}_n) \right) \gamma = \sum_{n=1}^{N} \tilde{w}_n \mathbf{E}^*(\hat{x}_n) y_n. \quad (28) $$

Assumption 2.1 is sufficient to ensure that $\hat{\gamma}$ is unique with $N=1$ observation.

It is possible to guarantee the closeness between $\tilde{\gamma}$ and $\hat{\gamma}$ under an additional Lipschitz regularity assumption on $\mathbf{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 $\mathbf{E}$ is $\sqrt{\sigma_+}L_E$-Lipschitz continuous<fnref n="2" />:

$$ \|\mathbf{E}(\mathbf{x}) - \mathbf{E}(\mathbf{x}')\|_{2 \to 2} \le \sqrt{\sigma_+} L_E \|\mathbf{x} - \mathbf{x}'\|_2 \quad \text{for all } (\mathbf{x}, \mathbf{x}') \in \mathbb{R}^D \times \mathbb{R}^D. \quad (29) $$

Then, under Assumption 2.1, we have

$$ \frac{\|\hat{\gamma} - \tilde{\gamma}\|_2}{\|\tilde{\gamma}\|_2} \le \kappa^{3/2} \frac{\|b_1\|_2}{\|y_{0,1}\|_2} + \epsilon_2(\hat{x}) \quad (30) $$

with

$$ \epsilon_2(\hat{x}) = c\kappa^{5/2} L_E \|\hat{x} - \bar{x}\|_2 \left( 1 + \frac{\|b_1\|_2}{\|y_{0,1}\|_2} \right) + O \left( \|\hat{x} - \bar{x}\|_2^2 \right) $$

for some absolute constant $c$.

Together with Theorem 3.3, this last result ensures that $\hat{\gamma} \to \tilde{\gamma}$ 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 \stackrel{\text{def.}}{=} \begin{pmatrix} y_{0,1} \\ \vdots \\ y_{0,N} \end{pmatrix} \quad \text{and} \quad B \stackrel{\text{def.}}{=} \begin{pmatrix} b_1 \\ \vdots \\ b_N \end{pmatrix} $$

denote the stacked noiseless measurements and noise vectors.

<footnote n="2">²The scaling in $\sqrt{\sigma_+}$ is natural considering Assumption 2.1.</footnote>
**Theorem 3.8** (Stability of the operator estimate with multiple observations). Given $\mathbf{X} = (\mathbf{x}_1, \dots, \mathbf{x}_n)$, let $\bar{w}_- = \min_n |\bar{w}_n|$, $\bar{w}_+ = \max_n |\bar{w}_n|$ and

$$C(\mathbf{X}) \stackrel{\text{def.}}{=} \sum_{n=1}^{N} \bar{w}_n^2 E^*(\mathbf{x}_n) E(\mathbf{x}_n).$$

Let $\tilde{\sigma}_- = \bar{w}_-^2 \hat{\sigma}_-$ and $\tilde{\sigma}_+ = \bar{w}_+^2 \hat{\sigma}_+$ and assume that

$$\tilde{\sigma}_{-} \mathrm{Id} \preceq C(\hat{\mathbf{X}}) \preceq \tilde{\sigma}_{+} \mathrm{Id}. \qquad (31)$$

Similarly to Theorem 3.7, assume that $\mathbf{E}$ is $\sqrt{\hat{\sigma}_+} L_E$-Lipschitz continuous and let $\tilde{\kappa} = \frac{\hat{\sigma}_+}{\hat{\sigma}_-}$. Then we have

$$\frac{\|\hat{\gamma} - \tilde{\gamma}\|_2}{\|\tilde{\gamma}\|_2} \le \tilde{\kappa}^{3/2} \frac{\|\mathbf{B}\|_2}{\|\mathbf{Y}_0\|_2} + \epsilon_2(\hat{\mathbf{X}}) \qquad (32)$$

with

$$\epsilon_2(\hat{\mathbf{X}}) \le c \tilde{\kappa}^{5/2} L_E \|\bar{\mathbf{X}} - \hat{\mathbf{X}}\|_2 \left( 1 + \frac{\|\mathbf{Y}_0\|_2 + \|\mathbf{B}\|_2}{\|\mathbf{Y}_0 + \mathbf{B}\|_2} \right) + O \left( \|\bar{\mathbf{X}} - \hat{\mathbf{X}}\|_2^2 \right),$$

for some absolute constant $c$.

Assumption (31) is a geometrical condition intertwining the locations of the Dirac masses and the observation mapping $\mathbf{E}$. It can be hard to verify in advance. However it only requires computing the $I \times I$ matrix $C(\hat{\mathbf{X}})$, which can be achieved once $\hat{\mathbf{X}}$ has been evaluated. The stable estimation of $\hat{\mathbf{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 $\hat{\mathbf{X}}$.

* Evaluate the condition number $\tilde{\kappa}$ of $C(\hat{\mathbf{X}})$.

* If $\tilde{\kappa}$ is sufficiently low, $\hat{\gamma}$ provides a good estimate of $\tilde{\gamma}$.

## 3.4 The case of unknown weights

Minimizing $(\mathcal{P}^{\gamma,w})$ with respect to $(w, \gamma)$ for a fixed $\mathbf{X} = \hat{\mathbf{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 [15, 18, 27, 19, 32], on the optimal stability to noise [33], and on the numerical aspects through convex lifting [34] or local optimization [35, 36, 37, 38, 39, 40]. 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.
### 3.4.1 A low-dimensional bilinear problem

Let $\hat{I}_n \stackrel{\text{def.}}{=} \dim(\mathcal{R}(\hat{x}_n))$ and $\hat{I} = \sum_{n=1}^N \hat{I}_n$. Using a singular value decomposition, we can decompose $\mathcal{E}(\hat{x}_n)$ as

$$ \mathcal{E}(\hat{x}_n) = \hat{U}_n \hat{V}_n^*, \quad (33) $$

where $\hat{U}_n \in \mathbb{R}^{M \times \hat{I}_n}$ and $\hat{V}_n$ contain orthogonal columns and $\hat{U}_n^* \hat{U}_n = \text{Id}$. Hence, letting $c_n \stackrel{\text{def.}}{=} \hat{U}_n^* y_n$, we obtain

$$ \begin{aligned} \underset{\mathbf{w} \in \mathbb{R}^N, \gamma \in \mathbb{R}^I}{\operatorname{argmin}} \quad & J(\mathbf{w}, \gamma, \hat{\mathbf{X}}) = \underset{\mathbf{w} \in \mathbb{R}^N, \gamma \in \mathbb{R}^I}{\operatorname{argmin}} \frac{1}{2} \sum_{n=1}^N \| \hat{U}_n \hat{V}_n^* w_n \gamma - \mathbf{y}_n \|_2^2 \\ & = \underset{\mathbf{w} \in \mathbb{R}^N, \gamma \in \mathbb{R}^I}{\operatorname{argmin}} \frac{1}{2} \sum_{n=1}^N \| \hat{V}_n^* w_n \gamma - c_n \|_2^2. \end{aligned} $$

Letting $\hat{\mathcal{B}}: \mathbb{R}^N \times \mathbb{R}^I \to \mathbb{R}^I$ denote the following bilinear mapping:

$$ \hat{\mathcal{B}}(\mathbf{w}, \gamma) \stackrel{\text{def.}}{=} \begin{pmatrix} \hat{V}_1^* w_1 \gamma \\ \vdots \\ \hat{V}_N^* w_N \gamma \end{pmatrix} \quad \text{and} \quad \mathbf{c} \stackrel{\text{def.}}{=} \begin{pmatrix} c_1 \\ \vdots \\ c_N \end{pmatrix}, \quad (34) $$

we can rewrite $J$ more compactly as $J(\mathbf{w}, \gamma, \hat{\mathbf{X}}) = \frac{1}{2} \|\hat{\mathcal{B}}(\mathbf{w}, \gamma) - \mathbf{c}\|_2^2$, and hence:

$$ \underset{\mathbf{w} \in \mathbb{R}^N, \gamma \in \mathbb{R}^I}{\operatorname{argmin}} \quad J(\mathbf{w}, \gamma, \hat{\mathbf{X}}) = \underset{\mathbf{w} \in \mathbb{R}^N, \gamma \in \mathbb{R}^I}{\operatorname{argmin}} \frac{1}{2} \|\hat{\mathcal{B}}(\mathbf{w}, \gamma) - \mathbf{c}\|_2^2. \quad (35) $$

Notice that the dimension $M$, which might be huge in applications, completely disappeared from this formulation.

### 3.4.2 A review of existing conditions for stable recovery

Recovering $\mathbf{w}$ and $\gamma$ is possible only up to a multiplicative constant since

$$ J(t\mathbf{w}, \gamma/t, \hat{\mathbf{X}}) = J(\mathbf{w}, \gamma, \hat{\mathbf{X}}) \text{ for all } t \neq 0. $$

Now, consider the noiseless setting $\mathcal{B} = 0$ and assume that the locations are perfectly recovered: $\hat{\mathbf{X}} = \hat{\mathbf{X}}$. In that situation, a necessary condition to recover $(\bar{\mathbf{w}}, \bar{\gamma})$ modulo the above scaling ambiguity is that there exists a unique pair $(\mathbf{w}, \gamma)$ with $\|\mathbf{w}\|_2 = 1$ such that $\hat{\mathcal{B}}(\mathbf{w}, \gamma) = \mathbf{c}$. From our current understanding, deriving geometrical conditions to ensure this local injectivity condition still deserves some attention.

In [41, 27], the authors study a more stringent global injectivity condition of the form

$$ \forall \mathbf{c} \in \mathbb{R}^I, \exists \text{ a unique } (\mathbf{w}, \gamma) \text{ with } \|\mathbf{w}\|_2 = 1 \text{ s.t. } \hat{\mathcal{B}}(\mathbf{w}, \gamma) = \mathbf{c}. \quad (36) $$

Their main result states that a necessary condition for $\hat{\mathcal{B}}$ to be globally injective is that

$$ \hat{I} \ge 2(N+I) - 4, \quad (37) $$
which provides a rule on how to choose the number of measurements $N$. In addition, they prove that almost every bilinear mapping $\hat{\mathcal{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 $\hat{\mathcal{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 $\hat{\mathcal{B}}$. They do not apply to the current context where everything is deterministic. Overall, the results of the algorithms described below are therefore empirical.

### 3.4.3 Optimization methods

In this section, we review 3 algorithms to solve (35).

**Optimization of the factors** Solving (35) can be achieved using local optimization over each factor $w$ and $\gamma$ [36, 39, 40]. 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

<caption>**Algorithm 1 Alternating minimization**</caption>

**Require:** Initial guess: $w_1 \in \mathbb{R}^N$.
**Require:** Iteration number $K$.
**for all** $k = 1 \rightarrow K - 1$ **do**
  $\gamma_{k+1} = \underset{\gamma \in \mathbb{R}^I}{\operatorname{argmin}} \frac{1}{2} \|\hat{\mathcal{B}}(w_k, \gamma) - c\|_2^2$.
  $w_{k+1} = \underset{w \in \mathbb{R}^N}{\operatorname{argmin}} \frac{1}{2} \|\hat{\mathcal{B}}(w, \gamma_{k+1}) - c\|_2^2$.
**end for**
**return** $(w_K, \gamma_K)$.

the dimensions of the problem are significantly reduced. This approach can be certified to recover a stable estimate $(\tilde{w}, \tilde{\gamma})$ of $(\tilde{w}, \tilde{\gamma})$ provided that a clever initialization is used [39, 40]. Sufficient recovery guarantees are for instance provided when the bilinear mapping $\hat{\mathcal{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 [40] as a starting point.

**Optimization over rank-1 matrices** The bilinear mapping $\hat{\mathcal{B}}(w, \gamma)$ can be rewritten as a linear mapping $\hat{\Lambda}$ on the rank-1 outer product $T = w\gamma^T$: $\hat{\mathcal{B}}(w, \gamma) = \hat{\Lambda}(T)$. Hence, we have:

$$ \inf_{w \in \mathbb{R}^N, \gamma \in \mathbb{R}^I} J(w, \gamma, \hat{X}) = \inf_{T \in \mathbb{R}^{N \times I}, \mathrm{rank}(T)=1} \frac{1}{2} \|\hat{\Lambda}(T) - c\|_2^2. \quad (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 $\mathcal{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

<caption>**Algorithm 2 Projected gradient descent**</caption>

**Require:** Initial guess: $T \in \mathbb{R}^{N \times I}$.
**Require:** Iteration number $K$.

Compute $\tau = \frac{1}{\|\hat{\Lambda}\|_2^2}$ using a power iteration.

**for** all $k = 1 \rightarrow K-1$ **do**
    $T_{k+1} = \Pi_{\tau}(T_k - \tau\hat{\Lambda}^*(\hat{\Lambda}(T_k) - c))$.

**end for**

Decompose $T_K = w_K\gamma_K^*$.

**return** $(w_K, \gamma_K)$.

algorithm has not been analyzed so far. During the review process, we found a paper [42] describing a similar type of idea. Again, we will use the spectral initialization from [40] as a starting guess for this algorithm in the numerical experiments.

**Convex relaxation using the nuclear norm** Finally, a popular method [15, 19, 34, 17] is a convex relaxation using the nuclear norm. The usual convex relaxation of the nonconvex problem (38) is the following:

$$ \inf_{T \in \mathbb{R}^{N \times I}, \hat{\Lambda}(T)=c} \|T\|_* \quad \text{or} \quad \inf_{T \in \mathbb{R}^{N \times I}} \frac{1}{2} \|\hat{\Lambda}(T) - c\|_2^2 + \lambda \|T\|_*, \quad (39) $$

where $\lambda > 0$ is a regularization parameter and $\| \cdot \|_*$ 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 [43]. The stable recovery of the tensor $\bar{w}\bar{\gamma}^T$ has been established under rather stringent conditions based on random subspace assumptions [15, 19]. 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 [44]. We do not further detail these algorithms, which are well documented in the literature.

# 4 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.

## 4.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 $\nu_m = \delta_{z_m}$,
<figure><img src="image_3.png" /><figcaption>Figure 2: The different families of convolution filters used in Section 4.1.</figcaption></figure>

with $z_m = m/M$ for $m \in \{1, \dots, M\}$. Notice that this case also covers the case of product-convolution operators since the ranges $\mathcal{R}(x)$ of convolution and product-convolution operators are identical.

### 4.1.1 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\sigma_i^2)) \text{ with } \sigma_i = 0.01 \cdot \frac{i-1}{I-1} + 0.03 \cdot \left(1 - \frac{i-1}{I-1}\right).$$

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 $\psi(x) = (1 - |x - 1|)_+$ denote the hat function.

$$e_i = \psi(x \cdot \sigma_i) \text{ where } \sigma_i = 0.02 \cdot \frac{i-1}{I-1} + 0.2 \cdot \left(1 - \frac{i-1}{I-1}\right).$$

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^\perp$) which will be used in all experiments to satisfy Assumption 2.3.

### 4.1.2 The inverse functions $\phi^{-1}$

As stated in Theorem 3.3, the critical element to guarantee a stable recovery of the locations $\bar{x}_m$ is the function $\phi$ and its inverse, which characterizes the angle between the ranges $\mathcal{R}(x)$ and $\mathcal{R}(x')$. To evaluate this function, we first sample the function $\|\Pi_{\mathcal{R}(0)}\Pi_{\mathcal{R}(k\Delta x)}\|_{2 \to 2}$ on a fine grid. We store the result in the vector $\phi_0(k) \stackrel{\text{def.}}{=} 1 - \|\Pi_{\mathcal{R}(0)}\Pi_{\mathcal{R}(k\Delta x)}\|_{2 \to 2}$ with a sampling step $\Delta x$. This function is not
<figure><img src="image_4.png" /><figcaption>Figure 3: The corresponding inverse functions $\phi_{i,+}^{-1}$ (in red) for $i \in \{1, 2, 3\}$, for the different convolution systems.</figcaption></figure>

necessarily nondecreasing. Hence, we find the closest nondecreasing function by solving an isotonic regression problem of the form:

$$ \inf_{\phi} \frac{1}{2} \|\phi - \phi_0\|_2^2 \text{ with } \phi_{k+1} - \phi_k \ge 0 \text{ and } \phi \ge \phi_0. $$

This problem is convex and can be solved using the CVX library [45] for instance. We use the solution $\hat{\phi}$ of this problem in place of $\phi$ in Assumption 2.2. The inverse filters are displayed in Figure 3. The stability to noise is dependent on the speed of ascent of $\phi_+^{-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 $\phi_{1,+}^{-1}$ and $\phi_{3,+}^{-1}$ behave roughly similarly in terms of speed of ascent.

### 4.1.3 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(|\hat{x} - \bar{x}|)$ for various noise levels and realizations. The expectation is estimated by averaging 100 noise realizations. We fix $\tilde{\gamma}$ once for all. We use white Gaussian noise, i.e. $b_n \sim N(0, \sigma^2 \mathbf{Id})$, with $\sigma = \theta \|y_0\|_2 / \sqrt{M}$ and $\theta \in [0, 2]$. Figure 4 shows the resulting signals with $M = 100$ for the noise levels $\theta \in \{0, 1, 2\}$ and each family. Notice that $\theta = 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(|\hat{x} - \bar{x}|)$ is displayed with respect to the noise level $\sigma$ 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 $\theta = 1$, corresponding to a measurement dominated by noise.

To end this experiment, we evaluate $\hat{\gamma}$ for all experiments and display the relative error $\frac{\|\hat{\gamma} - \tilde{\gamma}\|_2}{\|\tilde{\gamma}\|_2}$ for all families of operators. Letting $\bar{h} = \sum_{i=1}^I \tilde{\gamma}_i e_i^\perp$ and $\hat{h} = \sum_{i=1}^I \hat{\gamma}_i e_i^\perp$ denote the true convolution filter and the estimated one, notice that we have $\frac{\|\bar{h} - \hat{h}\|_{L^2(\mathbb{R}^D)}}{\|\bar{h}\|_{L^2(\mathbb{R}^D)}} = \frac{\|\hat{\gamma} - \tilde{\gamma}\|_2}{\|\tilde{\gamma}\|_2}$, since the family $(e_i^\perp)$ is orthogonal. In
<figure><img src="image_5.png" /><figcaption>(a) $\theta = 0, A_1$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) $\theta = 0, A_2$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) $\theta = 0, A_3$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) $\theta = 1, A_1$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(e) $\theta = 1, A_2$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(f) $\theta = 1, A_3$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(g) $\theta = 2, A_1$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(h) $\theta = 2, A_2$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(i) $\theta = 2, A_3$</figcaption></figure>

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

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.

## 4.2 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 $\tilde{X} = \bar{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 $\nu_m = \delta_{z_m}$ with $z_m = 10 \cdot 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, \mathbf{I}_d)$ with a Gaussian filter of large variance.
<figure><img src="image_14.png" /><figcaption>Figure 5: Average localization error $E(|\hat{x} - \bar{x}|)$ as a fraction of a pixel for different noise levels $\theta \in [0, 1]$ and the three families $\mathcal{A}_1, \mathcal{A}_2$ and $\mathcal{A}_3$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 6: Boxplots of the relative errors $\frac{\|\hat{\gamma} - \bar{\gamma}\|_2}{\|\bar{\gamma}\|_2}$ using a single measurement for various noise levels.</figcaption></figure>

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 $\gamma \sim \mathcal{N}(0, \mathbf{I}_d)$. 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 $\text{span}(e_j, j \in \{1, \dots, 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 $\bar{x}_n$, with random weights $\bar{w}_n$ 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 [40] 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 $(\bar{x}_n)$ are badly located.

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

$$N \geq \frac{2JK - 4}{J - 2},$$

i.e. $N \geq 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
<figure><img src="image_16.png" /><figcaption>Figure 7: Examples of random product-convolution operators for a fixed family ($e_j$) and ($f_k$) and 3 random realizations of $\gamma$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 8: Percentages of perfect recovery results for the different algorithms.</figcaption></figure>

algorithms manage to recover some operators when this condition is not met. However, it is clear that a necessary condition for identifiability is $N \ge K$, since otherwise, even the problem with known weights cannot be identified.

## 4.3 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 astig- matic impulse responses as follows. We set the family ($e_j$) as anisotropic Gaus- sian 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 po- sitions 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 in- accurate localization. We then used this information to recover the operator. Here we assumed that the weight ($\tilde{w}_n$) 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.

# 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.

# 5 Proofs

## 5.1 Proofs of the propositions from Section 2.3

### Proof of Proposition 2.1

*Proof.* We show that for all $(x, \gamma), (x', \gamma') \in \mathbb{R}^D \times \mathbb{R}^I \setminus \{0\}$

$$E(x)\gamma = E(x')\gamma' \Rightarrow \gamma = \gamma'$$

If $x = x'$: by Assumption 2.1 the mapping $\gamma \mapsto E(x)\gamma$ is injective, meaning that for all $\gamma, \gamma' \in \mathbb{R}^I \setminus \{0\}$ we have

$$E(x)\gamma = E(x)\gamma' \Rightarrow (x, \gamma) = (x, \gamma').$$

If $x \neq x'$: by Assumption 2.2 we have $\mathcal{R}(x) \cap \mathcal{R}(x') = \{0\}$, which implies that for any $\gamma, \gamma' \in \mathbb{R}^I \setminus \{0\}$, we have $E(x)\gamma \neq E(x')\gamma'$. We prove it by contradiction. Assume that $\mathcal{R}(x) \cap \mathcal{R}(x') \neq \{0\}$. Then, there exists $\gamma \in \mathbb{R}^I \setminus \{0\}$ such that $E(x)\gamma \in \mathcal{R}(x')$. Let $v \stackrel{\text{def.}}{=} \frac{E(x)\gamma}{\|E(x)\gamma\|_2}$. Then $\|\Pi_{\mathcal{R}(x')}(\Pi_{\mathcal{R}(x)}\|_{2 \to 2} \ge \|\Pi_{\mathcal{R}(x')}(v)\|_2 = \|v\|_2 = 1$. Nevertheless, by Assumption 2.2 we have $\|\Pi_{\mathcal{R}(x')}(\Pi_{\mathcal{R}(x)}\|_{2 \to 2} < 1$ for $x \neq x'$. This concludes the proof. $\square$

### Proof of Proposition 2.2

*Proof.* We have $[E(z)]_{m,i} = \langle e_i(\cdot - z), \delta_{z_m} \rangle = e_i(z_m - z)$. Hence

$$[E^*(z)E(z)]_{i,i'} = \sum_{m \in \mathbb{N}} e_i(z_m - z)e_{i'}(z_m - z) = \begin{cases} 1 & \text{if } i = i' \\ 0 & \text{otherwise,} \end{cases}$$

where we used (12) to obtain the last equality. $\square$
### Proof of Proposition 2.3

*Proof.* By Assumption 2.3, we have $\Pi_{\mathcal{R}(x)} = E(x)E^*(x)$. Then, by definition, we have

$$ 
\begin{aligned}
\|\Pi_{\mathcal{R}(x)}\Pi_{\mathcal{R}(x')}'\|_{2\to2} &= \|E(x)E^*(x)E(x')E^*(x')\|_{2\to2} \\
&= \sup_{u\in\mathbb{R}^M, \|u\|_2=1} \|E(x)E^*(x)E(x')E^*(x')u\|_2 \\
&= \sup_{v\in\mathbb{R}^I, \|v\|_2=1} \|E(x)E^*(x)E(x')v\|_2 \\
&= \|E(x)E^*(x)E(x')\|_{2\to2} = \|E^*(x)E(x')\|_{2\to2} \\
&= \|C(x-x')\|_{2\to2},
\end{aligned}
$$

using the fact that $E^*(x)E(x) = \text{Id}$.

□

### Proof of Proposition 2.4

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

## 5.2 Proof of Proposition 3.1

*Proof.* We have

$$ \inf_{x \in \mathcal{D}, \alpha \in \mathbb{R}^I} \frac{1}{2} \|E(x)\alpha - y_n\|_2^2 = \inf_{x \in \mathcal{D}} \inf_{\alpha \in \mathbb{R}^I} \frac{1}{2} \|E(x)\alpha - y_n\|_2^2. $$

The minimum norm solution of the inner problem is given by $\alpha(x) = E^+(x)y_n$, where $E^+(x)$ is the Moore-Penrose pseudo-inverse of $E(x)$. Injecting this solution leads to the problem

$$ 
\begin{aligned}
\inf_{x \in \mathcal{D}} \frac{1}{2} \|E(x)E^+(x)y_n - y_n\|_2^2 &= \inf_{x \in \mathcal{D}} \frac{1}{2} \|\Pi_{\mathcal{R}(x)}y_n - y_n\|_2^2 \\
&= \inf_{x \in \mathcal{D}} \frac{1}{2} \|y_n\|_2^2 + \frac{1}{2} \|\Pi_{\mathcal{R}(x)}y_n\|_2^2 - \langle \Pi_{\mathcal{R}(x)}y_n, y_n \rangle = \inf_{x \in \mathcal{D}} \frac{1}{2} \|y_n\|_2^2 - \frac{1}{2} \|\Pi_{\mathcal{R}(x)}y_n\|_2^2.
\end{aligned}
$$

Neglecting the constant term $\frac{1}{2}\|y_n\|_2^2$ and changing the sign of the function yields the result. □

## 5.3 Proof of Proposition 3.2

*Proof.* It is well known (see e.g., [46]) 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 $\Pi_{\mathcal{R}(x)} = E(x)E^+(x)$ and of the whole function $x \mapsto \|\Pi_{\mathcal{R}(x)}y_n\|_2^2$. A continuous function over a compact domain admits at least one global minimizer. □
## 5.4 Proof of Theorem 3.3

We start with a basic lemma.

**Lemma 5.1.** Let $D \in \mathbb{N}$, $S \subseteq \mathbb{R}^D$, $f: S \to \mathbb{R}$ and $\epsilon: S \to \mathbb{R}$ denote two mappings. Define $g \stackrel{\text{def.}}{=} f + \epsilon$ and assume that the following sets are non-empty

$$ \hat{\mathcal{X}} = \underset{\mathbf{x} \in S}{\operatorname{argmax}} g(\mathbf{x}) \quad \text{and} \quad \bar{\mathcal{X}} = \underset{\mathbf{x} \in S}{\operatorname{argmax}} f(\mathbf{x}). $$

Further assume that $\|\epsilon\|_{L^\infty(S)} \le \eta$ for some $\eta > 0$ and that there exists an increasing function $\varphi: \mathbb{R} \to \mathbb{R}_+$ such that

$$ f(\mathbf{x}) \le f(\bar{\mathbf{x}}) - \varphi(\|\mathbf{x} - \bar{\mathbf{x}}\|_2), \quad \forall \mathbf{x} \in S, \bar{\mathbf{x}} \in \bar{\mathcal{X}}. \quad (40) $$

Then $\bar{\mathcal{X}} = \{\bar{\mathbf{x}}\}$ is a singleton and any $\hat{\mathbf{x}} \in \hat{\mathcal{X}}$ satisfies $\|\hat{\mathbf{x}} - \bar{\mathbf{x}}\|_2 \le \varphi_+^{-1}(2\eta)$.

*Proof.* By inequality (40) and strict monotonicity of $\phi$, we have $f(\bar{\mathbf{x}}) > f(\mathbf{x})$ for all $\mathbf{x} \ne \bar{\mathbf{x}}$. Hence $\bar{\mathbf{x}}$ is the unique maximizer of $f$. We have

$$ g(\hat{\mathbf{x}}) \ge g(\bar{\mathbf{x}}) = f(\bar{\mathbf{x}}) - \epsilon(\bar{\mathbf{x}}) \ge f(\bar{\mathbf{x}}) - \eta. \quad (41) $$

In addition, for any $\mathbf{x} \in S$, we have

$$ \begin{aligned} g(\mathbf{x}) = f(\mathbf{x}) + \epsilon(\mathbf{x}) &\le f(\mathbf{x}) + \eta \\ &\le f(\bar{\mathbf{x}}) - \varphi(\|\mathbf{x} - \bar{\mathbf{x}}\|_2) + \eta. \end{aligned} $$

For any $\mathbf{x} \in S$ such that $\|\mathbf{x} - \bar{\mathbf{x}}\|_2 > \varphi_+^{-1}(2\eta)$, we have $g(\mathbf{x}) < f(\bar{\mathbf{x}}) - \eta$, and thus $g(\mathbf{x}) < g(\bar{\mathbf{x}})$, which implies that $\mathbf{x} \ne \hat{\mathbf{x}}$. The contraposition is that any $\hat{\mathbf{x}} \in \hat{\mathcal{X}}$ satisfies $\|\hat{\mathbf{x}} - \bar{\mathbf{x}}\|_2 \le \varphi_+^{-1}(2\eta)$. $\square$

We are ready to prove Theorem 3.3.

*Proof.* i) Let $\mathbf{y}_{0,n} = \bar{w}_n \mathbf{E}(\bar{\mathbf{x}}_n)\bar{\gamma} = \mathbf{y}_n - \mathbf{b}_n$ denote the noiseless measurements and let $F_{0,n}(\mathbf{x}, \gamma) = \frac{1}{2}\|\bar{w}_n \mathbf{E}(\mathbf{x})\gamma - \mathbf{y}_{0,n}\|_2^2$. We have $F_{0,n}(\bar{\mathbf{x}}_n, \bar{\gamma}) = 0$. Let $\gamma_0(\mathbf{x}) \in \operatorname{argmin}_{\gamma \in \mathbb{R}^l} F_{0,n}(\mathbf{x}, \gamma)$ denote any minimizer (for instance the one given by the pseudo-inverse). Now, define $G_{0,n}(\mathbf{x}) \stackrel{\text{def.}}{=} F_{0,n}(\mathbf{x}, \gamma_0(\mathbf{x}))$. We have $G_{0,n}(\mathbf{x}) = \frac{1}{2}(\|\mathbf{y}_{0,n}\|_2^2 - \|\mathbf{\Pi}_{\mathcal{R}(\mathbf{x})}\mathbf{y}_{0,n}\|_2^2)$. By Assumption 2.2, $\mathcal{R}(\mathbf{x}) \cap \mathcal{R}(\bar{\mathbf{x}}) = \{0\}$ for $\mathbf{x} \ne \bar{\mathbf{x}}$. Hence, $\|\mathbf{\Pi}_{\mathcal{R}(\mathbf{x})}\mathbf{y}_{0,n}\|_2^2 < \|\mathbf{y}_{0,n}\|_2^2$ for $\mathbf{x} \ne \bar{\mathbf{x}}$ and $\bar{\mathbf{x}}_n$ is the unique minimizer of $G_{0,n}$. Therefore, the function

$$ H_{0,n}(\mathbf{x}) \stackrel{\text{def.}}{=} \frac{1}{2} \|\mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} \mathbf{y}_{0,n}\|_2^2 = \frac{1}{2} \langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} \mathbf{y}_{0,n}, \mathbf{y}_{0,n} \rangle, $$

also admits a unique maximizer in $\bar{\mathbf{x}}$. Overall, we see that under Assumption 2.2, $G_{0,n}$ admits a unique minimizer equal to $\bar{\mathbf{x}}$. Under the additional Assumption 2.1, $F_{0,n}$ admits a unique solution $(\bar{\mathbf{x}}, \bar{\gamma})$.

ii) Now, let $F_n(\mathbf{x}, \gamma) = \frac{1}{2}\|\bar{w}_n \mathbf{E}(\mathbf{x})\gamma - \mathbf{y}_{0,n}\|_2^2$, $\gamma(\mathbf{x})$ denote any minimizer of $F_n$ w.r.t. $\gamma$, $G_n(\mathbf{x}) = F_n(\mathbf{x}, \gamma(\mathbf{x}))$ and $H_n(\mathbf{x}) \stackrel{\text{def.}}{=} \frac{1}{2}\langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})}\mathbf{y}_n, \mathbf{y}_n \rangle$. Let $\hat{\mathbf{x}}$ denote any maximizer of $H_n$ and assume for now that we manage to obtain a bound of the form $|H_n(\mathbf{x}) - H_{0,n}(\mathbf{x})| \le \eta$ for some $\eta \ge 0$ and all $\mathbf{x} \in \mathcal{D}$. We have

$$ \begin{aligned} H_{0,n}(\mathbf{x}) &= \frac{1}{2} \langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} \mathbf{y}_{0,n}, \mathbf{y}_{0,n} \rangle = \frac{1}{2} \langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} \mathbf{y}_{0,n}, \mathbf{\Pi}_{\mathcal{R}(\bar{\mathbf{x}})} \mathbf{y}_{0,n} \rangle \\ &= \frac{1}{2} \langle \mathbf{\Pi}_{\mathcal{R}(\bar{\mathbf{x}})} \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} \mathbf{y}_{0,n}, \mathbf{y}_{0,n} \rangle \le \frac{1}{2} \| \mathbf{\Pi}_{\mathcal{R}(\bar{\mathbf{x}})} \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} \|_{2 \to 2} \| \mathbf{y}_{0,n} \|_2^2 \\ &\le \frac{1}{2} (1 - \phi(\|\mathbf{x} - \bar{\mathbf{x}}\|_2)) \| \mathbf{y}_{0,n} \|_2^2 = H_{0,n}(\bar{\mathbf{x}}) - \frac{1}{2} \phi(\|\mathbf{x} - \bar{\mathbf{x}}\|_2) \| \mathbf{y}_{0,n} \|_2^2 \end{aligned} $$
by Assumption 2.2. Hence, we can use Lemma 5.1 with $S = \mathcal{D}$, $f(x) = H_{0,n}(x)$, $g(x) = H_n(x)$ and $\varphi(r) = \frac{1}{2}\phi(r)\|y_{0,n}\|_2^2$. This allows us to conclude that

$$ \|\hat{x} - \bar{x}\|_2 \le \phi_+^{-1} \left( \frac{4\eta}{\|y_{0,n}\|_2^2} \right). \quad (42) $$

iii) The last remaining point is to control $\|H_{0,n} - H_n\|_{L^\infty(\mathcal{D})}$. We have

$$ \begin{aligned} H_n(x) &= \frac{1}{2} \langle \Pi_{\mathcal{R}(x)} (y_{0,n} + b_n), y_{0,n} + b_n \rangle \\ &= H_{0,n}(x) + \langle \Pi_{\mathcal{R}(x)} y_{0,n}, b_n \rangle + \frac{1}{2} \| \Pi_{\mathcal{R}(x)} (b_n) \|_2^2. \end{aligned} $$

Hence, using Cauchy-Schwarz inequality, we obtain for all $x$

$$ |H_n(x) - H_{0,n}(x)| \le \|\Pi_{\mathcal{R}(x)} y_{0,n}\|_2 \|\Pi_{\mathcal{R}(x)} b_n\|_2 + \frac{1}{2} \|\Pi_{\mathcal{R}(x)} (b_n)\|_2^2. $$

Using the facts that $\|\Pi_{\mathcal{R}(x)}b_n\|_2 \le \|b_n\|_2$ and that $\|b_n\|_2 \le \theta\|y_{0,n}\|_2$, we obtain

$$ \|H_n - H_{0,n}\|_{L^\infty(\mathcal{D})} \le \|y_{0,n}\|_2^2 \left( \theta + \frac{1}{2}\theta^2 \right). $$

For the inequality (42) to make sense, we need $4(\theta + \frac{1}{2}\theta^2) \le 1$, which is equivalent to $\theta < \frac{\sqrt{6}}{2} - 1$ and Theorem 3.3 is proven. $\square$

## 5.5 Proof of Theorem 3.4

*Proof.* Assume that $\bar{\mu} = \sum_{n=1}^N \bar{w}_n \delta_{\bar{x}_n}$. Then, we get

$$ \mathbf{y} = \sum_{n=1}^{N} \mathbf{y}_n = \sum_{n=1}^{N} \bar{w}_n \mathbf{E}(\bar{x}_n) \bar{\alpha}. \quad (43) $$

We have

$$ H(x) = H_1(x) + \underbrace{\langle \Pi_{\mathcal{R}(x)} y_1, \Pi_{\mathcal{R}(x)} \sum_{n=2}^{N} y_n \rangle}_{\epsilon_1(x)} + \frac{1}{2} \underbrace{\| \Pi_{\mathcal{R}(x)} \sum_{n=2}^{N} y_n \|_2^2}_{\epsilon_2(x)}, \quad (44) $$

where $H_1(x) = \frac{1}{2} \|\Pi_{\mathcal{R}(x)} y_1\|_2^2$. Under Assumption 2.2, $H_1$ possesses a unique global maximizer at $\bar{x}_1$. Let $\delta > 0$ denote some radius and $B_{\delta/2} = \{x \in \mathbb{R}^D, \|x - \bar{x}_1\|_2 \le \delta/2\}$. Similarly to the previous proof, we now apply Lemma 5.1 with $S = B_{\delta/2}$, $f = H_1$, $\epsilon = \epsilon_1 + \epsilon_2$ and $\varphi(t) = \frac{1}{2}\phi(t)\|y_1\|_2^2$.

Consider a point $x \in B(\delta/2)$. Assumption 2.2 and the hypothesis

$$ \|\bar{x}_1 - \bar{x}_n\|_2 \ge \delta \text{ for all } n \ge 2 $$
allows to derive the following inequalities

$$
\begin{aligned}
|\epsilon_1(\mathbf{x})| &= \left| \left\langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} \mathbf{y}_1, \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} \sum_{n=2}^{N} \mathbf{y}_n \right\rangle \right| \\
&\le \sum_{n=2}^{N} \left| \left\langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} \mathbf{\Pi}_{\mathcal{R}(\bar{\mathbf{x}}_1)} \mathbf{y}_1, \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} \mathbf{\Pi}_{\mathcal{R}(\bar{\mathbf{x}}_n)} \mathbf{y}_n \right\rangle \right| \\
&\le \sum_{n=2}^{N} (1 - \phi(\|\mathbf{x} - \bar{\mathbf{x}}_1\|_2))(1 - \phi(\|\mathbf{x} - \bar{\mathbf{x}}_n\|_2)) \|\mathbf{y}_1\|_2 \|\mathbf{y}_n\|_2 \\
&\le (N-1)c(1 - \phi(\delta)) \|\mathbf{y}_1\|_2^2 \le (N-1)c(1 - \phi(\delta/2)) \|\mathbf{y}_1\|_2^2.
\end{aligned}
$$

Similarly, using Jensen inequality, we obtain

$$
\begin{aligned}
|\epsilon_2(\mathbf{x})| &= \frac{1}{2} \left\| \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} \sum_{n=2}^{N} \mathbf{y}_n \right\|_2^2 \\
&\le \frac{N-1}{2} \sum_{n=2}^{N} (1 - \phi(\delta/2))^2 \|\mathbf{y}_n\|_2^2 \\
&\le \frac{(N-1)^2 c^2}{2} (1 - \phi(\delta/2))^2 \|\mathbf{y}_1\|_2^2.
\end{aligned}
$$

Hence $\|\epsilon\|_{L^\infty(B_{\delta/2})} \le \frac{(N-1)^2 c^2}{2} (1 - \phi(\delta/2))^2 \|\mathbf{y}_1\|_2^2 + (N-1)c(1 - \phi(\delta/2)) \|\mathbf{y}_1\|_2^2$.

Set $Z = (N-1)c$. Lemma 5.1 states that any global maximizer $\hat{\mathbf{x}}$ of $H$ over $B_{\delta/2}$ satisfies

$$
\begin{aligned}
\| \hat{\mathbf{x}} - \bar{\mathbf{x}}_1 \|_2 &\le \phi^{-1} \left( 4 \left( \frac{Z^2}{2} (1 - \phi(\delta/2))^2 + Z(1 - \phi(\delta/2)) \right) \right) \\
&= a \cdot 2^{1/b} \cdot \left[ \frac{Z \cdot (2 + Z + 2(\delta/2a)^b)}{(1 + (\delta/2a)^b)^2 - 2Z^2 - 4Z(1 + (\delta/2a)^b)} \right]^{1/b},
\end{aligned}
$$

where the last line was obtained using a symbolic calculus software. Letting $\tau \ge 5$, the denominator above is positive. Setting $\delta = 2a(\tau Z)^{1/b}$, the above expression becomes

$$
r \geq \frac{a \cdot 2^{1/b}}{\tau} \cdot \left[ \frac{(\tau^2 + \tau)Z + 2\tau}{1 + (\tau^2 - 4\tau - 2)Z + (2\tau - 4)Z} \right].
$$

For $\tau \ge 5$, the above expression is bounded above by $18.4 \cdot \frac{a \cdot 2^{1/b}}{\tau}$, obtained for $\tau = 5$ and $Z \to \infty$. $\square$

## 5.6 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 $\mathbf{b} \sim \mathcal{N}(0, \sigma^2 \mathbf{I}_d)$ is white Gaussian noise of variance $\sigma^2$. Define the following two random processes

$$
\Delta_1(\mathbf{x}) \stackrel{\text{def.}}{=} \langle \Pi_{\mathcal{R}(\mathbf{x})} \mathbf{y}_0, \Pi_{\mathcal{R}(\mathbf{x})} \mathbf{b} \rangle \quad \text{and} \quad \Delta_2(\mathbf{x}) \stackrel{\text{def.}}{=} \frac{1}{2} \| \Pi_{\mathcal{R}(\mathbf{x})} (\mathbf{b}) \|_2^2. \quad (45)
$$
Then under Assumption 2.2

$$ \|\hat{\mathbf{x}} - \bar{\mathbf{x}}\|_2 \le \phi_+^{-1} \left( \frac{2\mathrm{Ampl}(\Delta_1 + \Delta_2)}{\|\mathbf{y}_0\|_2^2} \right), \quad (46) $$

where

$$ \mathrm{Ampl}(f) \stackrel{\mathrm{def.}}{=} \sup_{\mathbf{x} \in \mathcal{D}} f(\mathbf{x}) - \inf_{\mathbf{x} \in \mathcal{D}} f(\mathbf{x}). $$

*Proof.* This proposition derives from point iii) in the previous proof. The inequality (42) can be improved as

$$ \|\hat{\mathbf{x}} - \bar{\mathbf{x}}\|_2 \le \phi_+^{-1} \left( \frac{4 \inf_{c \in \mathbb{R}} \|H - H_0 - c\|_{\infty}}{\|\mathbf{y}_0\|_2^2} \right), \quad (47) $$

since a constant term does not affect the location of a minimizer. We have

$$ \Delta(\mathbf{x}) \stackrel{\mathrm{def.}}{=} H(\mathbf{x}) - H_0(\mathbf{x}) = \langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} y_0, \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} b \rangle + \frac{1}{2} \| \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} (\mathbf{b}) \|_2^2. $$

It is therefore natural to set $c = \frac{1}{2}(\sup_{\mathbf{x} \in \mathbb{R}^D} \Delta(\mathbf{x}) - \inf_{\mathbf{x} \in \mathbb{R}^D} \Delta(\mathbf{x}))$, to minimize in the infinite norm in Problem (47). □

This proposition reveals that the critical quantity to control, to evaluate the localization error is the amplitude of the random process $\Delta_1 + \Delta_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 $\mathrm{Ampl}(\Delta_1 + \Delta_2) \le \mathrm{Ampl}(\Delta_1) + \mathrm{Ampl}(\Delta_2)$. In addition, the random variable $Z_1 = \mathrm{Ampl}(\Delta_1)$ is sub-Gaussian:

$$ \mathbb{P}(|Z_1 - \bar{Z}_1| \ge t) \le 2 \exp(-t^2/(8\sigma^2 \|y_0\|_2^2)) $$

and the random variable $Z_2 = \mathrm{Ampl}(\Delta_2)$ is sub-exponential:

$$ \mathbb{P}(|Z_2 - \bar{Z}_2| \ge t) \le 2 \exp(-Ct/\sigma^2), $$

where $\bar{Z}_1$ and $\bar{Z}_2$ are the expectations of $Z_1$ and $Z_2$ and C is a universal constant.

*Proof.* First notice that

$$ \mathrm{Ampl}(\Delta_1 + \Delta_2) = \sup_{\mathbf{x} \in \mathcal{D}} \Delta_1(\mathbf{x}) + \Delta_2(\mathbf{x}) - \inf_{\mathbf{x} \in \mathcal{D}} \Delta_1(\mathbf{x}) + \Delta_2(\mathbf{x}) $$

$$ \le \sup_{\mathbf{x} \in \mathcal{D}} \Delta_1(\mathbf{x}) + \sup_{\mathbf{x} \in \mathcal{D}} \Delta_2(\mathbf{x}) - \left( \inf_{\mathbf{x} \in \mathcal{D}} \Delta_1(\mathbf{x}) + \inf_{\mathbf{x} \in \mathcal{D}} \Delta_2(\mathbf{x}) \right) $$

$$ = \mathrm{Ampl}(\Delta_1) + \mathrm{Ampl}(\Delta_2). $$

We will treat the two random processes $\Delta_1$ and $\Delta_2$ separately.

i) Consider the function $f_1(\mathbf{b}) \stackrel{\mathrm{def.}}{=} \sup_{\mathbf{x} \in \mathcal{D}} \langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} y_0, \mathbf{b} \rangle$ and define the random variable $V_1^+ = f_1(\mathbf{b})$ with mean $\bar{V}_1^+$. Similarly, define $V_1^- = \inf_{\mathbf{x} \in \mathcal{D}} \langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} y_0, \mathbf{b} \rangle$. In addition, notice that $Z_1 = \mathrm{Ampl}(\Delta_1) \le V_1^+ - V_1^-$. We first show that $f_1$ is Lipschitz continuous. We have

$$ f_1(\mathbf{b} + \boldsymbol{\epsilon}) = \sup_{\mathbf{x} \in \mathcal{D}} \langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} y_0, \mathbf{b} + \boldsymbol{\epsilon} \rangle \le \sup_{\mathbf{x} \in \mathcal{D}} \langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} y_0, \mathbf{b} \rangle + \|y_0\|_2 \|\boldsymbol{\epsilon}\|_2 $$

$$ f_1(\mathbf{b} + \boldsymbol{\epsilon}) = \sup_{\mathbf{x} \in \mathcal{D}} \langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} y_0, \mathbf{b} + \boldsymbol{\epsilon} \rangle \ge \sup_{\mathbf{x} \in \mathcal{D}} \langle \mathbf{\Pi}_{\mathcal{R}(\mathbf{x})} y_0, \mathbf{b} \rangle - \|y_0\|_2 \|\boldsymbol{\epsilon}\|_2 $$
Hence, $|f_1(\mathbf{b}) - f_1(\mathbf{b} + \boldsymbol{\epsilon})| \le \|\mathbf{y}_0\|_2\|\boldsymbol{\epsilon}\|_2$ and $f_1$ is $\|\mathbf{y}_0\|_2$-Lipschitz continuous. Using a Gaussian logarithmic Sobolev inequality [47, Thm 5.6], we obtain that $V_1^+$ is sub-Gaussian with

$$ \mathbb{P}(|V_1^+ - \bar{V}_1^+| \ge t) \le 2 \exp(-t^2/(2\sigma^2 \|\mathbf{y}_0\|^2)). $$

The same result holds for $V_1^-$. Finally, the sum of two dependent sub-Gaussian variables with parameters $\sigma_1$ and $\sigma_2$ is sub-Gaussian with a sub-Gaussian parameter smaller than $\sigma_1 + \sigma_2$, so that

$$ \mathbb{P}(|Z_1 - \bar{Z}_1| \ge t) \le 2 \exp(-t^2/(8\sigma^2 \|\mathbf{y}_0\|^2)). $$

ii) Now, define the random variable $Y_2^+ \stackrel{\text{def.}}{=} \sup_{x \in D} \|\Pi_{\mathcal{R}(x)}\mathbf{b}\|_2$ and $V_2^+ \stackrel{\text{def.}}{=} \frac{1}{2}(Y_2^+)^2$. The function $\mathbf{b} \mapsto \sup_{x \in D} \|\Pi_{\mathcal{R}(x)}\mathbf{b}\|_2$ is 1-Lipschitz continuous. Hence using a Gaussian logarithmic Sobolev inequality again, we obtain that $Y_2^+$ is sub-Gaussian with

$$ \mathbb{P}(|Y_2^+ - \bar{Y}_2^+| \ge t) \le 2 \exp \left( -\frac{t^2}{2\sigma^2} \right). $$

Using [48, Lemma 2.7.6], we conclude that $V_2^+$ is sub-exponential and satisfies

$$ \mathbb{P}(|V_2^+ - \bar{V}_2^+| \ge t) \le 2 \exp(-Ct/\sigma^2), $$

for an absolute constant $C$. We can make a similar proof for the random variable $Y_2^- \stackrel{\text{def.}}{=} \inf_{x \in D} \|\Pi_{\mathcal{R}(x)}\mathbf{b}\|_2$ and conclude as in the previous proof. $\square$

Notice that the averages $\bar{Z}_1$ and $\bar{Z}_2$ 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 $\sigma\|\mathbf{y}_0\|_2$ and the deviation around the mean of $Z_2$ scales as $\sigma^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 $\bar{Z}_1$ and $\bar{Z}_2$. 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 $\sigma\|\mathbf{y}_0\|_2 + \sigma^2$ is negligible.

Unfortunately, the averages $\bar{Z}_1$ and $\bar{Z}_2$ 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 $\bar{Z}_1$ and $\bar{Z}_2$ under additional regularity assumptions.

**Proposition 5.4** (Control of $\bar{Z}_1$). Under the assumptions of Theorem 3.5, we have, for all $\bar{x} \in \mathbb{R}^D$

$$ \bar{Z}_1 \le c \cdot \left(\frac{L}{\beta}\right)^{\frac{\alpha}{\alpha+1}} \cdot \sigma \sqrt{DI} \|\mathbf{y}_0\|_2, \quad (48) $$

where $c$ is an absolute constant.

*Proof.* Here, we wish to control the supremum $\bar{Z}_1$ of the centered Gaussian process $\Delta_1$. A traditional approach to bound it consists in computing Dudley's entropy integral, see e.g., [47, Corollary 13.2].
**Pseudo-metric, covering numbers and Dudley's integral** To this end, we introduce the pseudo-metric:

$$ \mathrm{dist}(\boldsymbol{x}, \boldsymbol{x}') \stackrel{\mathrm{def.}}{=} \sqrt{\mathbb{E}((\Delta_1(\boldsymbol{x}) - \Delta_1(\boldsymbol{x}'))^2)}. \quad (49) $$

Let $\mathcal{B}(c, \delta, \mathrm{dist}) = \{\boldsymbol{x} \in \mathbb{R}^D, \mathrm{dist}(c, \boldsymbol{x}) \le \delta\}$ denote a ball of radius $\delta$ centered at $c$ with respect to $\mathrm{dist}$. Let $S \subseteq D$ denote a set and define the covering number $\Theta(\delta, S)$ as the minimum number of $\delta$-balls needed to cover $S$. We recall the following simple results: for $S \subseteq S'$, we have $\Theta(\delta, S) \le \Theta(\delta, S')$. In addition, if $S$ can be partitioned as $S = S_1 \sqcup S_2$, we have $\Theta(\delta, S) \le \Theta(\delta, S_1) + \Theta(\delta, S_2)$.

Dudley's theorem reads as follows:

$$ \mathbb{E} \left( \sup_{\boldsymbol{x} \in S} \Delta_1(\boldsymbol{x}) \right) \le 24 \int_0^{\eta/4} \sqrt{\log(\Theta(u, S))} du \quad (50) $$

with

$$ \eta \stackrel{\mathrm{def.}}{=} \inf_{\boldsymbol{x} \in S} \sup_{\boldsymbol{x}' \in S} \mathrm{dist}(\boldsymbol{x}, \boldsymbol{x}'). \quad (51) $$

**Bounding the pseudo-metric** Using the fact that $y_0 \in \mathcal{R}(\bar{x})$ and that $y_0 = \Pi_{\mathcal{R}(\bar{x})} y_0$:

$$
\begin{align*}
\mathrm{dist}(\boldsymbol{x}, \boldsymbol{x}')^2 &= \mathbb{E}((\Delta_1(\boldsymbol{x}) - \Delta_1(\boldsymbol{x}'))^2) \\
&= \mathbb{E}((\langle \Pi_{\mathcal{R}(\boldsymbol{x})} y_0, \Pi_{\mathcal{R}(\boldsymbol{x})} b \rangle - \langle \Pi_{\mathcal{R}(\boldsymbol{x}')} y_0, \Pi_{\mathcal{R}(\boldsymbol{x}')} b \rangle)^2) \\
&= \mathbb{E}((\langle \Pi_{\mathcal{R}(\boldsymbol{x})} \Pi_{\mathcal{R}(\bar{x})} y_0, b \rangle - \langle \Pi_{\mathcal{R}(\boldsymbol{x}')} \Pi_{\mathcal{R}(\bar{x})} y_0, b \rangle)^2) \\
&= \mathbb{E}((\langle [\Pi_{\mathcal{R}(\boldsymbol{x})} - \Pi_{\mathcal{R}(\boldsymbol{x}')}] \Pi_{\mathcal{R}(\bar{x})} y_0, b \rangle)^2) \\
&= \mathbb{E}((\langle [\Pi_{\mathcal{R}(\boldsymbol{x})} - \Pi_{\mathcal{R}(\boldsymbol{x}')}] \Pi_{\mathcal{R}(\bar{x})} y_0, \Pi_{\mathcal{R}(\bar{x})} b \rangle)^2) \\
&\le \| [\Pi_{\mathcal{R}(\boldsymbol{x})} - \Pi_{\mathcal{R}(\boldsymbol{x}')}] \Pi_{\mathcal{R}(\bar{x})} \|_{2 \to 2}^2 \| y_0 \|_2^2 \cdot \mathbb{E} (\| \Pi_{\mathcal{R}(\bar{x})} b \|_2^2) \\
&\le \sigma^2 I \| y_0 \|_2^2 \| [\Pi_{\mathcal{R}(\boldsymbol{x})} - \Pi_{\mathcal{R}(\boldsymbol{x}')}] \Pi_{\mathcal{R}(\bar{x})} \|_{2 \to 2}^2.
\end{align*}
$$

In addition:

$$
\begin{align*}
\| (\Pi_{\mathcal{R}(\boldsymbol{x})} - \Pi_{\mathcal{R}(\boldsymbol{x}')} ) \Pi_{\mathcal{R}(\bar{x})} \|_{2 \to 2} &\le \min (1, \| \Pi_{\mathcal{R}(\boldsymbol{x})} \Pi_{\mathcal{R}(\bar{x})} \|_{2 \to 2} + \| \Pi_{\mathcal{R}(\boldsymbol{x}')} \Pi_{\mathcal{R}(\bar{x})} \|_{2 \to 2}) \\
&\le \min \left( 1, \frac{1}{\beta^\alpha \| \boldsymbol{x} - \bar{\boldsymbol{x}} \|_2^\alpha} + \frac{1}{\beta^\alpha \| \boldsymbol{x}' - \bar{\boldsymbol{x}} \|_2^\alpha} \right)
\end{align*}
$$

Letting $\zeta \stackrel{\mathrm{def.}}{=} \sigma \sqrt{I} \|y_0\|_2$ and $\psi(r) = \frac{1}{\beta^\alpha r^\alpha}$, we get:

$$ \mathrm{dist}(\boldsymbol{x}, \boldsymbol{x}') \le \zeta \cdot \min (1, \psi(\|\boldsymbol{x} - \bar{\boldsymbol{x}}\|_2) + \psi(\|\boldsymbol{x}' - \bar{\boldsymbol{x}}\|_2)). \quad (52) $$

In addition, the Lipschitz continuity of $\Pi_{\mathcal{R}(\boldsymbol{x})}$ also implies that:

$$
\begin{align*}
\mathrm{dist}(\boldsymbol{x}, \boldsymbol{x}')^2 &= \mathbb{E}((\Delta_1(\boldsymbol{x}) - \Delta_1(\boldsymbol{x}'))^2) \\
&= \mathbb{E}((\langle [\Pi_{\mathcal{R}(\boldsymbol{x})} - \Pi_{\mathcal{R}(\boldsymbol{x}')}] y_0, \Pi_{\mathcal{R}(\bar{x})} b \rangle)^2) \tag{53}
\end{align*}
$$

$$ \le \mathbb{E} \left( \| \Pi_{\mathcal{R}(\boldsymbol{x})} - \Pi_{\mathcal{R}(\boldsymbol{x}')} \|_{2 \to 2}^2 \| y_0 \|_2^2 \| \Pi_{\mathcal{R}(\bar{x})} b \|_2^2 \right) \quad (54) $$

$$ \le \sigma^2 I \|y_0\|_2^2 L^2 \|\boldsymbol{x} - \boldsymbol{x}'\|_2^2. \quad (55) $$
In what follows, we let

$$ \overline{\mathrm{dist}}(\mathbf{x}, \mathbf{x}') \stackrel{\mathrm{def.}}{=} \zeta \cdot \min (1, L\|\mathbf{x} - \mathbf{x}'\|_2, \psi(\|\mathbf{x} - \bar{\mathbf{x}}\|_2) + \psi(\|\mathbf{x}' - \bar{\mathbf{x}}\|_2)) \quad (56) $$

and $\bar{\Theta}$ denote the corresponding covering number. The inequality $\overline{\mathrm{dist}}(\mathbf{x}, \mathbf{x}') \le \overline{\mathrm{dist}}(\mathbf{x}, \mathbf{x}')$ implies that $\Theta(\delta, S) \le \bar{\Theta}(\delta, S)$ for all $S$ and $\delta$.

**Bounding the integral** We have $\overline{\mathrm{dist}}(\mathbf{x}, \mathbf{x}') \le \zeta$ for all $\mathbf{x}, \mathbf{x}'$. Therefore:

$$ \mathbb{E} \left( \sup_{\mathbf{x} \in \mathcal{D}} \Delta_1(\mathbf{x}) \right) \le \mathbb{E} \left( \sup_{\mathbf{x} \in \mathbb{R}^D} \Delta_1(\mathbf{x}) \right) \le 24 \int_0^{\zeta/4} \sqrt{\log(\bar{\Theta}(u, \mathbb{R}^D))} du. \quad (57) $$

By working over $\mathbb{R}^D$ rather than $\mathcal{D}$, we can assume that $\bar{\mathbf{x}} = 0$ without loss of generality. We have

$$ \psi(\|\mathbf{x} - \bar{\mathbf{x}}\|_2) + \psi(\|\mathbf{x}' - \bar{\mathbf{x}}\|_2) \le \frac{2}{\beta^\alpha \min(\|\mathbf{x}' - \bar{\mathbf{x}}\|_2, \|\mathbf{x} - \bar{\mathbf{x}}\|_2)^\alpha} $$

Hence, if both $\|\mathbf{x}\|_2$ and $\|\mathbf{x}'\|_2$ are larger than $(\frac{2\zeta}{\beta^\alpha u})^{1/\alpha}$, $\overline{\mathrm{dist}}(\mathbf{x}, \mathbf{x}') \le u$. This implies that

$$ \bar{\Theta} \left( u, \left\{ \mathbf{x} \in \mathbb{R}^D, \|\mathbf{x}\|_2 \ge \left( \frac{2\zeta}{\beta^\alpha u} \right)^{1/\alpha} \right\} \right) = 1. $$

It remains to control the covering number of the central Euclidean ball of radius $(\frac{2\zeta}{\beta^\alpha u})^{1/\alpha}$. In this region, we can use the Lipschitz inequality $\overline{\mathrm{dist}}(\mathbf{x}, \mathbf{x}') \le \zeta L \|\mathbf{x} - \mathbf{x}'\|_2$. It implies that is sufficient to cover the central ball with Euclidean balls of radius $\frac{u}{\zeta L}$. Standard results (see e.g., [48, Cor 4.2.13]) yield:

$$ \bar{\Theta} \left( u, \left\{ \mathbf{x} \in \mathbb{R}^D, \|\mathbf{x}\|_2 \le \left( \frac{2\zeta}{\beta^\alpha u} \right)^{1/\alpha} \right\} \right) \le \left( 3 \cdot \left( \frac{2\zeta}{\beta^\alpha u} \right)^{\frac{1}{\alpha}} \cdot \frac{\zeta L}{u} \right)^D. $$

Letting $a = \left(3 \cdot \left(\frac{2}{\beta^\alpha}\right)^{\frac{1}{\alpha}} \cdot L\right)^D$ and $b = D^{\frac{\alpha+1}{\alpha}}$, we obtain:

$$ \bar{\Theta}(u, \mathbb{R}^D) \le a \cdot \left(\frac{u}{\zeta}\right)^{-b} + 1. \quad (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(\bar{\Theta}(u, \mathbb{R}^D))$ only up to the smallest $u$, say $u_s$ such that $\bar{\Theta}(u, \mathbb{R}^D) \ge 2$. Up to this value, $a \cdot (u/\zeta)^{-b}$ is necessarily larger than 1, so that $u_s \le \zeta a^{1/b}$ and we obtain

$$ \begin{aligned} \mathbb{E} \left( \sup_{\mathbf{x} \in \mathbb{R}^D} \Delta_1(\mathbf{x}) \right) &\le 24 \int_0^{u_s} \sqrt{\log(a \cdot u^{-b} + 1)} du \le 24 \int_0^{\zeta a^{1/b}} \sqrt{\log(2a \cdot (u/\zeta)^{-b})} du \\ &= 24\zeta \cdot \sqrt{b} \cdot (2a)^{1/b} \int_0^{2^{-1/b}} \sqrt{\log(v^{-1})} dv \le 24\zeta \cdot \sqrt{b} \cdot (2a)^{1/b} \int_0^1 \sqrt{\log(v^{-1})} dv \\ &= 24\zeta \cdot \sqrt{b} \cdot (2a)^{1/b} \cdot \frac{\sqrt{\pi}}{2}. \end{aligned} $$
Letting $c$ denote an absolute constant, we finally obtain the bound

$$ \mathbb{E} \left( \sup_{x \in \mathbb{R}^D} \Delta_1(\boldsymbol{x}) \right) \leq c \cdot \sqrt{D \frac{\alpha+1}{\alpha}} \cdot \beta^{-\alpha/(\alpha+1)} \cdot L^{\frac{\alpha}{\alpha+1}} \cdot \sigma \sqrt{I} \|\boldsymbol{y}_0\|_2. $$

The term $\sqrt{\frac{\alpha+1}{\alpha}}$ can be absorbed in the constant since $\alpha \ge 1/2$. To conclude, we use the fact that $\bar{Z}_1 \le 2\mathbb{E}(\sup_{x \in \mathbb{R}^D} \Delta_1(\boldsymbol{x}))$. $\square$

**Proposition 5.5 (Control of $\bar{Z}_2$).** Under the assumptions of Theorem 3.5, we have, for all $\bar{x} \in \mathcal{D}$

$$ \begin{aligned} & \mathbb{P} \left( Z_2 \ge c_1 \sigma^2 (D \log(L) + \sqrt{D \log(L)I}) + t \right) \\ & \le 2 \exp \left( -c_2 \min \left( \frac{t^2}{\sigma^4 (D \log(L) + I)}, \frac{t}{\sigma^2} \right) \right) \end{aligned} $$

for some absolute constants $c, c_1, c_2$. This implies

$$ \bar{Z}_2 \le c \cdot \sigma^2 \left( D \log(L) + \sqrt{D \log(L)I} \right). \quad (59) $$

*Proof.* Controlling $\bar{Z}_2$ 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 [49, 50, 51, 47]. We will make use of the following theorem, which has been rewritten using our notation.

**Theorem 5.6 (Theorem 3.1 in [51]).** Let $\mathcal{S} = \{\Pi_{\mathcal{R}(\boldsymbol{x})}, \boldsymbol{x} \in \mathcal{D}\}$ denote the family of projections. Let $d_{2 \to 2}(\mathcal{S}) = \sup_{\mathcal{A} \in \mathcal{S}} \|\mathcal{A}\|_{2 \to 2}$ and $d_F(\mathcal{S}) = \sup_{\mathcal{A} \in \mathcal{S}} \|\mathcal{A}\|_F$. Let us define the following Dudley integral

$$ \gamma = c_1 \int_0^{d_{2 \to 2}(\mathcal{S})} \sqrt{\log \Theta(u, \mathcal{S})} du, $$

where the covering number is evaluated w.r.t. the spectral distance $\| \cdot \|_{2 \to 2}$. Then

$$ \mathbb{P} \left( \sup_{x \in \mathcal{D}} |\Delta_2(x) - \mathbb{E}(\Delta_2(x))| \ge \sigma^2 (c_2 E + t) \right) \le 2 \exp \left( -c_3 \min \left( \frac{t^2}{\sigma^4 V^2}, \frac{t}{\sigma^2 U} \right) \right), $$

with

$$ E = \gamma^2 + \gamma d_F(\mathcal{S}) + d_F(\mathcal{S}) d_{2 \to 2}(\mathcal{S}), \quad U = d_{2 \to 2}^2(\mathcal{S}), \quad V = d_{2 \to 2}(\mathcal{S}) (\gamma + d_F(\mathcal{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

$$ \mathrm{Ampl}(\Delta_2) \le 2 \sup_{x \in \mathcal{D}} |\Delta_2(x) - \mathbb{E}(\Delta_2(x))|. $$

We have $d_F(\mathcal{S}) = \sqrt{I}$ since $\dim(\mathcal{R}(\boldsymbol{x})) = I$ and $d_{2 \to 2}(\mathcal{S}) = 1$ since the matrices $\Pi_{\mathcal{R}(\boldsymbol{x})}$ are projections. The Lipschitz continuity assumption (23) also yields $\Theta(u, \mathcal{S}) \le (c_u^L)^D$.
To evaluate Dudley's integral, we can separate two cases. Let us first assume that $cL > 1$. In that case, we get

$$ 
\begin{aligned}
\gamma &= c_1 \int_{0}^{d_2 \to 2(S)} \sqrt{\log \Theta(u, S)} du \le c_1 \int_{0}^{1} \sqrt{\log \left(\frac{cL}{u}\right)^D} du \\
&= c_1 \sqrt{D} \left( \sqrt{\log(cL)} + \frac{\sqrt{\pi}}{2} cL \cdot \operatorname{erfc}(\sqrt{cL}) \right) \le c_1 \sqrt{D} \left( \sqrt{\log(cL)} + \frac{1}{2\sqrt{\log(cL)}} \right).
\end{aligned}
$$

where we used the inequality $\operatorname{erfc}(z) < \frac{\exp(-z^2)}{\sqrt{\pi z}}$ to obtain the last result. If $L$ is larger than $\exp(1)/L$, the term $\frac{1}{2\sqrt{\log(cL)}}$ is smaller than $\sqrt{\log(cL)}$ and can be discarded, up to changing the multiplicative constant. Therefore, we get

$$ \gamma \le c_1 \sqrt{D} \sqrt{\log(cL)}. $$

Now, let us propose another bound covering the case $cL \le 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 $\Theta(u, S) \ge 2$. The previous inequality yields $(c \frac{L}{u_s})^D \ge 2$ and hence $u_s \le cL/2^{1/D} \le cL$. Therefore, we get that

$$ 
\begin{aligned}
\gamma &= c_1 \int_{0}^{d_2 \to 2(S)} \sqrt{\log \Theta(u, S)} du = c_1 \int_{0}^{u_s} \sqrt{\log \Theta(u, S)} du \\
&\le c_1 \int_{0}^{cL} \sqrt{\log \left(\frac{cL}{u}\right)^D} du = c_1 \sqrt{D} \int_{0}^{cL} \sqrt{\log \frac{cL}{u}} du = cL\sqrt{D}.
\end{aligned}
$$

Overall, if $L$ is larger than a universal constant $c$, we obtain:

$$ 
\begin{aligned}
E &\le c_1 \cdot \left( D \log(L) + \sqrt{D \log(L)} I + \sqrt{I} \right). \\
V &\le c_1 \cdot \left( \sqrt{D \log(L)} + \sqrt{I} \right) \text{ and } U = 1.
\end{aligned}
$$

Recall that for a nonnegative random variable $X$, $E(X) = \int_{0}^{\infty} P(X > t) dt$. Hence, the concentration inequality

$$ P(\text{Ampl}(\Delta_2) \ge \sigma^2(c_2E + t)) \le 2 \exp \left( -c_3 \min \left( \frac{t^2}{\sigma^4 V^2}, \frac{t}{\sigma^2 U} \right) \right), $$

yields

$$ 
\begin{aligned}
\bar{Z}_2 &\stackrel{\text{def.}}{=} E(\text{Ampl}(\Delta_2)) \le c \cdot (\sigma^2 E + \sigma^2 V + \sigma U) \\
&\le c\sigma^2(D \log(L) + \sqrt{D \log(L)}I + \sqrt{I}).
\end{aligned}
$$

□

We now have all the ingredients to prove the main Theorem 3.5. We have for any real $t_1, t_2$

$$ 
\begin{aligned}
\mathbb{P} (Z_1 + Z_2 - \bar{Z}_1 - \bar{Z}_2 \ge t_1 + t_2) &\le \mathbb{P} (([Z_1 - \bar{Z}_1 \le t_1] \cap [Z_2 - \bar{Z}_2 \le t_2])^c) \\
&\le \mathbb{P} ([Z_1 - \bar{Z}_1 \ge t_1] \cup [Z_2 - \bar{Z}_2 \ge t_2]) \le \mathbb{P} (Z_1 - \bar{Z}_1 \ge t_1) + \mathbb{P} (Z_2 - \bar{Z}_2 \ge t_2) \\
&\le 2 \exp \left( -\frac{t_1^2}{8\sigma^2 \|y_0\|_2^2} \right) + 2 \exp \left( -\frac{ct_2}{\sigma^2} \right).
\end{aligned}
$$
where we used Proposition 5.3 to obtain the last inequality. Let $\rho > 0$ denote an arbitrary number and set $t_1 = c\rho\sigma\|\mathbf{y}_0\|_2$, $t_2 = c\sigma^2\rho$. Injecting these values in the above inequality, we obtain

$$ \mathbb{P}(Z_1 + Z_2 - \bar{Z}_1 - \bar{Z}_2 \ge t_1 + t_2) \le \exp(-c\rho^2) + \exp(-c\rho) \quad (60) $$

for some universal constant $c$. In words, the probability that $Z_1 + Z_2$ deviates from $\bar{Z}_1 + \bar{Z}_2$ by more than $\rho(\sigma\|\mathbf{y}_0\|_2 + \sigma^2)$ decays exponentially fast in $\rho$. Now let $\epsilon \in \mathbb{R}_+$ and let $t = \phi(\epsilon) \in [0, 1]$. Using the fact that both $\phi$ and $\phi^{-1}$ are increasing, we obtain

$$ \begin{aligned} \mathbb{P}(\|\hat{\mathbf{x}} - \bar{\mathbf{x}}\|_2 \ge \epsilon) &= \mathbb{P}(\|\hat{\mathbf{x}} - \bar{\mathbf{x}}\|_2 \ge \phi^{-1}(t)) \\ &= \mathbb{P}(\phi(\|\hat{\mathbf{x}} - \bar{\mathbf{x}}\|_2) \ge t) \end{aligned} \quad (61) $$

$$ \begin{aligned} &\le \mathbb{P}\left(\frac{2\text{Ampl}(\Delta_1 + \Delta_2)}{\|\mathbf{y}_0\|_2^2} \ge t\right) \\ &\le \mathbb{P}\left(Z_1 + Z_2 \ge \frac{t\|\mathbf{y}_0\|_2^2}{2}\right). \end{aligned} \quad (62) $$

As previously, set $t_1 = c\rho\sigma\|\mathbf{y}_0\|_2$, $t_2 = c\sigma^2\rho$ and

$$ \frac{t\|\mathbf{y}_0\|_2^2}{2} = \bar{Z}_1 + \bar{Z}_2 + t_1 + t_2. $$

Combining (60) and (62), we obtain

$$ \mathbb{P}(\|\hat{\mathbf{x}} - \bar{\mathbf{x}}\|_2 \ge \epsilon) \le \exp(-c\rho^2) + \exp(-c\rho) $$

for

$$ \epsilon = \phi^{-1}\left(2\frac{\bar{Z}_1 + \bar{Z}_2 + t_1 + t_2}{\|\mathbf{y}_0\|_2^2}\right) $$

It remains to use Propositions 5.4 and 5.5 to obtain

$$ \begin{aligned} &\bar{Z}_1 + \bar{Z}_2 + c\rho(\sigma\|\mathbf{y}_0\|_2 + \sigma^2) \\ &\le c\left(\sigma^2(D\log(L) + \sqrt{DI}\log(L)) + \left(\frac{L}{\beta}\right)^{\frac{\alpha}{\alpha+1}} \cdot \sigma\sqrt{DI}\|\mathbf{y}_0\|_2 + \rho(\sigma\|\mathbf{y}_0\|_2 + \sigma^2)\right) \\ &= c \cdot \left(\sigma\|\mathbf{y}_0\|_2 \left[\left(\frac{L}{\beta}\right)^{\frac{\alpha}{\alpha+1}} \cdot \sqrt{DI} + \rho\right] + \sigma^2 \left[D\log(L) + \sqrt{DI}\log(L) + \rho\right]\right). \end{aligned} $$

## 5.7 Proof of Proposition 3.6

*Proof.* We take the convention of [28]. Define $\text{sinc}(x) \stackrel{\text{def.}}{=} \sin(\pi x)/(\pi x)$ and let $\mathcal{F}$ denote the Fourier transform on $L^2(\mathbb{R})$ defined by

$$ \mathcal{F}(f)(\omega) \stackrel{\text{def.}}{=} \int_{\mathbb{R}} f(x) \exp(-ix\omega) dx. $$

We recall that

$$ \mathcal{F}^{-1}(g)(x) = \frac{1}{2\pi} \int_{\mathbb{R}} g(\omega) \exp(ix\omega) d\omega, $$
and the Parseval identity reads

$$ \langle u, v \rangle_{L^2(\mathbb{R})} = \frac{1}{2\pi} \langle \mathcal{F}(u), \mathcal{F}(v) \rangle_{L^2(\mathbb{R})}. $$

Notice that the normalization by $1/\sqrt{a}$ ensures that $\|e\|_{L^2(\mathbb{R})} = 1$ for all $a > 0$. The Fourier transform of $e$ is given by:

$$ \mathcal{F}(e)(\omega) = \int_{\mathbb{R}} e(x) \exp(-ix\omega) dx = \sqrt{a} \cdot 1_{\Omega_a}(\omega), $$

with $\Omega_a = [-\pi/a, \pi/a]$.

We recall the Shannon-Whittaker interpolation formula. For a given $b \le a$, let

$$ \psi_m(t) \stackrel{\text{def.}}{=} \frac{\text{sinc}(t/b - m)}{\sqrt{b}}. \qquad (63) $$

Then $(\psi_m)_{m \in \mathbb{Z}}$ is an orthogonal basis of $\mathcal{PW}(\Omega_a)$ for the usual scalar product on $L^2(\mathbb{R})$. In addition, for $u \in \mathcal{PW}(\Omega_a)$,

$$ \langle u, \psi_m \rangle_{L^2(\mathbb{R})} = \sqrt{b} \cdot u(bm). $$

Hence for any $u \in \mathcal{PW}(\Omega_a)$, we get the Shannon interpolation formula:

$$ u(x) = \sum_{m \in \mathbb{Z}} \langle u, \psi_m \rangle_{L^2(\mathbb{R})} \psi_m(x) = \sqrt{b} \sum_{m \in \mathbb{Z}} u(mb) \psi_m(x). $$

Therefore, for any $u, v \in \mathcal{PW}(\Omega_a)$, we obtain

$$ \langle u, v \rangle_{L^2(\mathbb{R})} = b \sum_{m \in \mathbb{Z}} u(mb)v(mb). $$

The assumptions imply that $\mathbf{E}(x)$ is the infinite vector in $\ell^2$ defined by $\mathbf{E}(x) = (e(x - z_m))_{m \in \mathbb{Z}}$. Using Parseval identity, we get for all $x, x' \in \mathbb{R}$ that

$$ \begin{aligned} \| \mathbf{\Pi}_{\mathcal{R}(x)} \mathbf{\Pi}_{\mathcal{R}(x')} \|_{2 \to 2} &\le \frac{|\langle \mathbf{E}(x), \mathbf{E}(x') \rangle_{\ell^2}|}{\| \mathbf{E}(x) \|_2 \| \mathbf{E}(x') \|_2} \\ &= \frac{|\int_{\mathbb{R}} e(t-x)e(t-x') dt|}{\| e(\cdot-x) \|_{L^2(\mathbb{R})} \| e(\cdot-x') \|_{L^2(\mathbb{R})}} \\ &= \frac{1}{a} \left| \int_{\mathbb{R}} \text{sinc}((t-x)/a) \text{sinc}((t-x')/a) dt \right| \\ &= \frac{1}{a \cdot 2\pi} \left| \int_{\mathbb{R}} \mathcal{F}(\text{sinc}((\cdot-x)/a))(\omega) \overline{\mathcal{F}(\text{sinc}((\cdot-x')/a))(\omega)} d\omega \right| \\ &= \frac{a}{2\pi} \left| \int_{-\pi/a}^{\pi/a} \exp(-ix\omega) \exp(ix'\omega) d\omega \right| \\ &= \left| \text{sinc} \left( \frac{x-x'}{a} \right) \right| \\ &\le \frac{1}{(|x-x'|/a)}. \end{aligned} $$
Hence, the Assumption of Theorem 3.5 holds with $\alpha = 1$ and $\beta = 1/a$. In addition, we have

$$
\begin{align*}
\| \mathbf{\Pi}_{\mathcal{R}(x)} - \mathbf{\Pi}_{\mathcal{R}(x')} \|_{2 \to 2}^2 &= \left\| \frac{\mathbf{E}(x)}{\|\mathbf{E}(x)\|_2} - \frac{\mathbf{E}(x')}{\|\mathbf{E}(x')\|_2} \right\|_{\ell^2}^2 \\
&= \|e(\cdot - x) - e(\cdot - x')\|_{L^2(\mathbb{R})}^2 \\
&= \frac{1}{2\pi} \|[\exp(-ix\cdot) - \exp(-ix\cdot')]\mathcal{F}(e)\|_{L^2(\mathbb{R})}^2 \\
&= \frac{a}{2\pi} \int_{-\pi/a}^{\pi/a} |\exp(-ix\omega) - \exp(-ix'\omega)|^2 d\omega \\
&= \frac{a}{2\pi} \int_{-\pi/a}^{\pi/a} 2 - 2\cos((x' - x)\omega) d\omega \\
&= \frac{a}{\pi} \left[ \omega - \frac{\sin((x' - x)\omega)}{x' - x} \right]_{-\pi/a}^{\pi/a} \\
&= 2[1 - \mathrm{sinc}((x' - x)/a)]
\end{align*}
$$

Using the Taylor expansion of the sinc, we obtain for all $x \in \mathbb{R}$

$$1 - \mathrm{sinc}(x) \leq \frac{x^2 \pi^2}{6}.$$

This allows to take $\phi(r) = \frac{r^2 \pi^2}{3a^2}$ in Assumption 2.2, $\phi^{-1}(r) = \frac{a\sqrt{3r}}{\pi}$ and

$$\| \mathbf{\Pi}_{\mathcal{R}(x)} - \mathbf{\Pi}_{\mathcal{R}(x')} \|_{2 \to 2}^2 \leq \frac{2(x - x')^2 \pi^2}{6a^2} \text{ implying } L \leq \frac{\pi}{\sqrt{3a}}.$$

The fact that we observe a single source located at $\bar{x}_1 \in \mathcal{D}$ with weight $\bar{w}_1 = 1$ implies that $y_{0,1} = \mathbf{E}(\bar{x}_1)$ and that $\|y_0\|_{\ell^2}^2 = \frac{1}{b}$. Plugging these expressions in $\epsilon$ and ignoring the constants, we get that

$$\epsilon = \phi^{-1} \left( c \left( \sigma \sqrt{b} (1 + \rho) + b \sigma^2 \left( -\log(a) + \sqrt{-\log(a)} + \rho \right) \right) \right). \quad (64)$$

Now, to correctly sample the signal, we can set $b = \tau a$ for some $\tau \le 1$. We obtain using $\sqrt{-\log(t)} \lesssim -\log(t)$ for small $t$:

$$
\begin{align*}
\epsilon &= \phi^{-1} \left( c \left( \sigma \sqrt{\tau a} (1 + \rho) + \tau a \sigma^2 (-\log(a) + \rho) \right) \right) \\
&= ca \sqrt{\left( \sigma \sqrt{\tau a} (1 + \rho) + \tau a \sigma^2 (-\log(a) + \rho) \right)}.
\end{align*}
$$

□

## 5.8 Proof of Theorem 3.7

*Proof.* Let $\mathbf{P}(\mathbf{x}) \stackrel{\text{def.}}{=} \mathbf{E}^*(\mathbf{x})\mathbf{E}(\mathbf{x})$. By definition, we have

$$\hat{\gamma} = \mathbf{P}(\hat{\mathbf{x}})^{-1} \mathbf{E}^*(\hat{\mathbf{x}}) (\mathbf{y}_{0,1} + \mathbf{b}_1) \quad (65)$$

We have $\mathbf{E}(\hat{\mathbf{x}}) = \mathbf{E}(\bar{x}) + \mathbf{\Delta}$ with $\|\mathbf{\Delta}\|_{2 \to 2} \le \sqrt{\sigma_+ L_E} \|\hat{\mathbf{x}} - \bar{x}\|_2$. Hence $\mathbf{P}(\hat{\mathbf{x}}) = \mathbf{P}(\bar{x}) + \mathbf{\Delta}'$ with

$$\|\mathbf{\Delta}'\|_{2 \to 2} \le 2L_E \|\hat{\mathbf{x}} - \bar{x}\|_2 \sqrt{\sigma_+} \|\mathbf{E}(\bar{x})\|_{2 \to 2} + \sigma_+ L_E^2 \|\hat{\mathbf{x}} - \bar{x}\|_2^2. \quad (66)$$
The linear system to recover $\hat{\gamma}$ can be rewritten as

$$ (\mathbf{P}(\bar{x}) + \mathbf{\Delta}')\hat{\gamma} = (\mathbf{E}^*(\bar{x}) + \mathbf{\Delta})(\mathbf{y}_{0,1} + \mathbf{b}_1) = \mathbf{E}^*(\bar{x})\mathbf{y}_{0,1} + \mathbf{\delta} \quad (67) $$

with $\mathbf{\delta} = \mathbf{\Delta}(\mathbf{y}_{0,1} + \mathbf{b}_1) + \mathbf{E}^*(\bar{x})\mathbf{b}_1$. Under Assumption 2.1, the unique solution of $\mathbf{P}(\bar{x})\mathbf{\gamma} = \mathbf{E}^*(\bar{x})\mathbf{y}_{0,1}$ is $\bar{\gamma}$. If $\hat{x}$ is sufficiently close to $\bar{x}$, we have $\|\mathbf{\Delta}'\|_{2\to2} < \sigma_-$ and $\|\mathbf{P}(\bar{x})^{-1}\mathbf{\Delta}'\|_{2\to2} < 1$. We can now use standard results of linear algebra, see e.g. [52, 3.6], to obtain that

$$ \frac{\|\hat{\gamma} - \bar{\gamma}\|_2}{\|\bar{\gamma}\|_2} \le \frac{\|\mathbf{P}(\bar{x})\|_{2\to2} \|\mathbf{P}(\bar{x})^{-1}\|_{2\to2}}{1 - \|\mathbf{P}(\bar{x})^{-1}\mathbf{\Delta}'\|_{2\to2}} \left( \frac{\|\mathbf{\Delta}'\|_{2\to2}}{\|\mathbf{P}(\bar{x})\|_{2\to2}} + \frac{\|\mathbf{\delta}\|_2}{\|\mathbf{E}^*(\bar{x})\mathbf{y}_{0,1}\|_2} \right) \\ \le \frac{\sigma_+}{\sigma_-} \frac{1}{\left(1 - \frac{\|\mathbf{\Delta}'\|_{2\to2}}{\sigma_-}\right)} \left( \frac{\|\mathbf{\Delta}'\|_{2\to2}}{\sigma_-} + \frac{\|\mathbf{\Delta}\|_{2\to2}(\|\mathbf{y}_{0,1}\|_2 + \|\mathbf{b}_1\|_2) + \sqrt{\sigma_+}\|\mathbf{b}_1\|_2}{\sqrt{\sigma_-}\|\mathbf{y}_{0,1}\|_2} \right). $$

Letting $\hat{x} \to \bar{x}$, we see that the noise term contains an error term $\epsilon_1$ defined by

$$ \epsilon_1 \stackrel{\text{def.}}{=} \kappa^{3/2} \frac{\|\mathbf{b}_1\|_2}{\|\mathbf{y}_{0,1}\|_2}, $$

which does not vanish as $\hat{x} \to \bar{x}$ and an additional term $\epsilon_2(\hat{x})$ which does vanish. Using the fact that

$$ \frac{1}{\left(1 - \frac{\|\mathbf{\Delta}'\|_{2\to2}}{\sigma_-}\right)} = 1 + \frac{\|\mathbf{\Delta}'\|_{2\to2}}{\sigma_-} + O(\|\hat{x} - \bar{x}\|_2^2), $$

we obtain

$$ \epsilon_2(\bar{x}) = c\kappa L_E \|\hat{x} - \bar{x}\|_2 \left[ \kappa + \sqrt{\kappa} \left( 1 + \frac{\|\mathbf{b}_1\|_2}{\|\mathbf{y}_{0,1}\|_2} \right) + \kappa^{3/2} \frac{\|\mathbf{b}_1\|_2}{\|\mathbf{y}_{0,1}\|_2} \right] \\ + O(\|\hat{x} - \bar{x}\|_2^2) \\ \le c\kappa^{5/2} L_E \|\hat{x} - \bar{x}\|_2 \left( 1 + \frac{\|\mathbf{b}_1\|_2}{\|\mathbf{y}_{0,1}\|_2} \right) + O(\|\hat{x} - \bar{x}\|_2^2), $$

for some absolute constant $c$. $\square$

## 5.9 Proof of Theorem 3.8

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

## References

[1] Eric Betzig, George H Patterson, Rachid Sougrat, O Wolf Lindwasser, Scott Olenych, Juan S Bonifacino, Michael W Davidson, Jennifer Lippincott-Schwartz, and Harald F Hess. Imaging intracellular fluorescent proteins at nanometer resolution. *Science*, 313(5793):1642–1645, 2006.

[2] William E Moerner and Lothar Kador. Optical detection and spectroscopy of single molecules in a solid. *Physical review letters*, 62(21):2535, 1989.
[3] Daniel Sage, Hagai Kirshner, Thomas Pengo, Nico Stuurman, Junhong Min, Suliana Manley, and Michael Unser. Quantitative evaluation of software packages for single-molecule localization microscopy. *Nature methods*, 12(8):717-724, 2015.

[4] Daniel Sage, Thanh-An Pham, Hazen Babcock, Tomas Lukes, Thomas Pengo, Jerry Chao, Ramraj Velmurugan, Alex Herbert, Anurag Agrawal, Silvia Colabrese, et al. Super-resolution fight club: assessment of 2d and 3d single-molecule localization microscopy software. *Nature methods*, 16(5):387-395, 2019.

[5] Vincent Duval and Gabriel Peyré. Exact support recovery for sparse spikes deconvolution. *Foundations of Computational Mathematics*, 15(5):1315-1355, 2015.

[6] Valentin Debarnot, Paul Escande, and Pierre Weiss. A scalable estimator of sets of integral operators. *Inverse Problems*, 2019.

[7] Valentin Debarnot, Paul Escande, Thomas Mangeat, and Pierre Weiss. Learning low-dimensional models of microscopes. *IEEE Transactions on Computational Imaging*, 2020.

[8] Otmar Scherzer, Markus Grasmair, Harald Grossauer, Markus Haltmeier, and Frank Lenzen. *Variational methods in imaging*. Springer, 2009.

[9] Yohann De Castro and Fabrice Gamboa. Exact reconstruction using beurling minimal extrapolation. *Journal of Mathematical Analysis and applications*, 395(1):336-354, 2012.

[10] Emmanuel J Candès and Carlos Fernandez-Granda. Towards a mathematical theory of super-resolution. *Communications on pure and applied Mathematics*, 67(6):906-956, 2014.

[11] Quentin Denoyelle, Vincent Duval, and Gabriel Peyré. Support recovery for sparse super-resolution of positive measures. *Journal of Fourier Analysis and Applications*, 23(5):1153-1194, 2017.

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

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

[14] Axel Flinth, Frédéric de Gournay, and Pierre Weiss. On the linear convergence rates of exchange and continuous methods for total variation minimization. *Mathematical Programming*, pages 1-37, 2020.

[15] Ali Ahmed, Benjamin Recht, and Justin Romberg. Blind deconvolution using convex programming. *IEEE Transactions on Information Theory*, 60(3):1711-1732, 2013.

[16] Shuyang Ling and Thomas Strohmer. Self-calibration and biconvex compressive sensing. *Inverse Problems*, 31(11):115002, 2015.
[17] Yuejie Chi. Guaranteed blind sparse spikes deconvolution via lifting and convex optimization. *IEEE Journal of Selected Topics in Signal Processing*, 10(4):782–794, 2016.

[18] Peter Jung, Felix Krahmer, and Dominik Stöger. Blind demixing and deconvolution at near-optimal rate. *IEEE Transactions on Information Theory*, 64(2):704–727, 2017.

[19] Ali Ahmed and Laurent Demanet. Leveraging diversity and sparsity in blind deconvolution. *IEEE Transactions on Information Theory*, 64(6):3975–4000, 2018.

[20] Mohamed A Suliman and Wei Dai. Blind two-dimensional super-resolution and its performance guarantee. *arXiv preprint arXiv:1811.02070*, 2018.

[21] Jinchi Chen, Weiguo Gao, Sihan Mao, and Ke Wei. Vectorized hankel lift: A convex approach for blind super-resolution of point sources. *arXiv preprint arXiv:2008.05092*, 2020.

[22] Shun-ichi Amari, Scott C Douglas, Andrzej Cichocki, and Howard H Yang. Multichannel blind deconvolution and equalization using the natural gradient. In *First IEEE Signal Processing Workshop on Signal Processing Advances in Wireless Communications*, pages 101–104. IEEE, 1997.

[23] Filip Sroubek and Jan Flusser. Multichannel blind iterative image restoration. *IEEE Transactions on Image Processing*, 12(9):1094–1106, 2003.

[24] Liming Wang and Yuejie Chi. Blind deconvolution from multiple sparse inputs. *IEEE Signal Processing Letters*, 23(10):1384–1388, 2016.

[25] Yanjun Li and Yoram Bresler. Multichannel sparse blind deconvolution on the sphere. *IEEE Transactions on Information Theory*, 65(11):7415–7436, 2019.

[26] Laixi Shi and Yuejie Chi. Manifold gradient descent solves multi-channel sparse blind deconvolution provably and efficiently. *IEEE Transactions on Information Theory*, 2021.

[27] Yanjun Li, Kiryung Lee, and Yoram Bresler. Identifiability and stability in blind deconvolution under minimal assumptions. *IEEE Transactions on Information Theory*, 63(7):4619–4633, 2017.

[28] Stéphane Mallat. *A wavelet tour of signal processing*. Elsevier, 1999.

[29] Paul Escande and Pierre Weiss. Approximation of integral operators using product-convolution expansions. *Journal of Mathematical Imaging and Vision*, 58(3):333–348, 2017.

[30] Ralf C. Flicker and Francois J. Rigaut. Anisoplanatic deconvolution of adaptive optics images. *J. Opt. Soc. Am. A*, 22(3):504–513, Mar 2005.

[31] Jérémie Bigot, Paul Escande, and Pierre Weiss. Estimation of linear operators from scattered impulse responses. *Applied and Computational Harmonic Analysis*, 47(3):730–758, 2019.
[32] Felix Krahmer and Dominik Stöger. On the convex geometry of blind deconvolution and matrix completion. *arXiv preprint arXiv:1902.11156*, 2019.

[33] Yuxin Chen, Jianqing Fan, Bingyan Wang, and Yuling Yan. Convex and nonconvex optimization are both minimax-optimal for noisy blind deconvolution. *arXiv preprint arXiv:2008.01724*, 2020.

[34] Robert Beinert and Kristian Bredies. Tensor-free proximal methods for lifted bilinear/quadratic inverse problems with applications to phase retrieval. *arXiv preprint arXiv:1907.04875*, 2019.

[35] P-A Absil, Robert Mahony, and Rodolphe Sepulchre. *Optimization algorithms on matrix manifolds*. Princeton University Press, 2009.

[36] Jérôme Bolte, Shoham Sabach, and Marc Teboulle. Proximal alternating linearized minimization for nonconvex and nonsmooth problems. *Mathematical Programming*, 146(1-2):459–494, 2014.

[37] Valerio Cambareri and Laurent Jacques. Through the haze: a non-convex approach to blind gain calibration for linear random sensing models. *Information and Inference: A Journal of the IMA*, 8(2):205–271, 2019.

[38] Yann Traonmilin and Jean-François Aujol. The basins of attraction of the global minimizers of the non-convex sparse spike estimation problem. *Inverse Problems*, 36(4):045003, 2020.

[39] Z. Zhu, Q. Li, G. Tang, and M. B. Wakin. Global optimality in low-rank matrix optimization. *IEEE Transactions on Signal Processing*, 66(13):3614–3628, 2018.

[40] Xiaodong Li, Shuyang Ling, Thomas Strohmer, and Ke Wei. Rapid, robust, and reliable blind deconvolution via nonconvex optimization. *Applied and computational harmonic analysis*, 47(3):893–934, 2019.

[41] Michael Kech and Felix Krahmer. Optimal injectivity conditions for bilinear inverse problems with applications to identifiability of deconvolution problems. *SIAM Journal on Applied Algebra and Geometry*, 1(1):20–37, 2017.

[42] Henrik Eisenmann, Felix Krahmer, Max Pfeffer, and André Uschmajew. Riemannian thresholding methods for row-sparse and low-rank matrix recovery. *arXiv preprint arXiv:2103.02356*, 2021.

[43] Claire Boyer, Antonin Chambolle, Yohann De Castro, Vincent Duval, Frédéric De Gournay, and Pierre Weiss. On representer theorems and convex regularization. *SIAM Journal on Optimization*, 29(2):1260–1281, 2019.

[44] Patrick L Combettes and Jean-Christophe Pesquet. Proximal splitting methods in signal processing. In *Fixed-point algorithms for inverse problems in science and engineering*, pages 185–212. Springer, 2011.

[45] Michael Grant and Stephen Boyd. Cvx: Matlab software for disciplined convex programming, version 2.1, 2014.
[46] GW Stewart. On the continuity of the generalized inverse. *SIAM Journal on Applied Mathematics*, 17(1):33–45, 1969.

[47] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. *Concentration inequalities: A nonasymptotic theory of independence*. Oxford university press, 2013.

[48] Roman Vershynin. *High-dimensional probability: An introduction with applications in data science*, volume 47. Cambridge university press, 2018.

[49] Rafał Latała et al. Estimates of moments and tails of gaussian chaoses. *The Annals of Probability*, 34(6):2315–2331, 2006.

[50] Michel Talagrand. *The generic chaining: upper and lower bounds of stochastic processes*. Springer Science & Business Media, 2006.

[51] Felix Krahmer, Shahar Mendelson, and Holger Rauhut. Suprema of chaos processes and the restricted isometry property. *Communications on Pure and Applied Mathematics*, 67(11):1877–1904, 2014.

[52] Eugene E Tyrtyshnikov. *A brief introduction to numerical analysis*. Springer Science & Business Media, 2012.