{'fraction_non_alphanumeric': 0.039626115941353346, 'fraction_numerical': 0.01442856809864572, 'mean_word_length': 4.423767383059419, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Adaptive Superresolution in Deconvolution of Sparse Peaks

The aim of this paper is to investigate superresolution in deconvolution driven by sparsity priors. The observed signal is a convolution of an original signal with a continuous kernel.With the prior knowledge that the original signal can be considered as a sparse combination of Dirac delta peaks, we seek to estimate the positions and amplitudes of these peaks by solving a finite dimensional convex problem on a computational grid. Because, the support of the original signal may or may not be on this grid, by studying the discrete deconvolution of sparse peaks using L1-norm sparsity prior, we confirm recent observations that canonically the discrete reconstructions will result in multiple peaks at grid points adjacent to the location of the true peak. Owning to the complexity of this problem, we analyse carefully the de-convolution of single peaks on a grid and gain a strong insight about the dependence of the reconstructed magnitudes on the exact peak location. This in turn allows us to infer further information on recovering the location of the exact peaks i.e. to perform super-resolution. We analyze in detail the possible cases that can appear and based on our theoretical findings, we propose an self-driven adaptive grid approach that allows to perform superresolution in one-dimensional and multi-dimensional spaces. With the view that the current study can provide a further step in the development of more robust algorithms for the detection of single molecules in fluorescence microscopy or identification of characteristic frequencies in spectral analysis, we demonstrate how the proposed approach can recover sparse signals using simulated clusters of point sources (peaks) of low-resolution in one and two-dimensional spaces.


A. Deconvolution of sparse peaks on discrete grids
In a wide range of imaging applications, the signal of interest comprises a sequence of sparse peaks (or point sources) for instance in fluoresce microscopy [45], [55], [72], astronomy [52], ultrasound or Doppler technology [78], [79], [4], [6], medical imaging [76], [14] and computational neuroscience [33]. In these applications, one has frequently to solve the problem that the signal of interest cannot be observed directly, but has to be inferred from other quantities, often low spatial resolution observations which mathematically can be described as the convolution of the original signal with a smooth kernel.
In this article, we study the superresolution problem, known as sparse peak deconvolution [29], where one seeks to estimate the positions and amplitudes of the underlying sparse peaks from a set of blurred observations. As idealized data we consider the convolution of a measure µ on Ω ⊆ R d with a known symmetric and smooth (with infinite support) kernel G which attains its maximum at 0, i.e.
f (x) = (G * µ)(x) := Ω G(x − y) dµ(y) . (I.1) Here we consider the convolution operator from M(Ω) to L 2 (Ω), which is well-defined by the Fourier convolution theorem (cf. [19]) and we are interested in the reconstruction of sparse peaks when the corresponding original signal is of the form where L is the total number of peaks and δ ξ l = δ(x−ξ l ) denotes a concentrated measure (expressed through the Dirac-delta function δ) at location ξ l : Ω → R d with amplitude γ l .
In order to obtain a sparse reconstruction it is nowadays standard to employ the well-established ℓ 1 -norm minimization approaches (also known as Basic Pursuit or LASSO) [24], [67], [20] which, in addition to sparse promoting solutions, allow the linearization of the original problem, the direct application of fast convex optimization solvers (e.g. [38]) and do not require application of Fourier transform [16]. So, instead of solving a continuous deconvolution problem [21], [15], [13], [16], the aim is to reconstruct µ via a discrete set of concentrated measures, i.e to look for a discrete July  Since the support of signal µ (I.2) may or may not be on the computational grid {x k } k=1,...,N , three chief questions for the ℓ 1 -norm estimates arise: • How the error between the original signal µ and discrete signal µ N is quantified based on the discretization?
• What are the expected patterns of the discrete estimation µ N on an arbitrary grid?
• Can the locations and amplitudes of the original signal µ be approximated with the help of the solution µ N ?
In this article, we investigate and answer these questions with the help of convex optimization theory and standard numerical analysis. We anticipate that the understanding of the effects of sparsity promoting solvers when the computational grid and the support of the original sparse signal do not coincide will allow the development of more robust algorithms required in applications such as fluorescence microscopy [34], [70], [69], [40], [37], [72], [43], [61].

B. Related works
Sparsity prior driven deconvolution approaches in continuous domains have been studied in several works including [21], [15], [13], [27], [36]. The signal (sum of Dirac functions) to be recovered is not a finite-dimensional vector (as in (I.3)) but a Radon measure and the minimization problem is formulated with the help of the total variation (TV) term and this problem is referred to as Beurling-Lasso (BLASSO) [21], [27], [65]. An extensive theoretical analysis of BLASSO in the case of one dimensional Fourier measurements was provided in [16]. Particularly, it was shown July 14, 2020 DRAFT that if the spikes are separated enough, then the exact recovery is possible (when the fraction of the measurement noise and regularization parameter tends to zero). Robustness to noise under this separation condition was studied by [3], [27], while the effect of the positivity constraint was analyzed in [8], [56], [22], [31], [30]. The observation sampling and the exact support recovery was theoretically studied in [25]. In [65], the BLASSO problem was analyzed for measures in higher dimensional spaces revealing that the kernel and arrangement of the original peaks affect the stability in the estimates.
BLASSO is a convex but infinite dimensional optimization problem. As shown in [16], [9], [75], [7], [36], solvers exist for ideal lowpass filters (i.e. Dirichlet type of kernels) when the observations are transformed into the Fourier domain and consider a finite number of frequencies in one dimensional signal spaces. Particularly, in theses cases, the primal (BLASSO) problem is expressed via its (Fenchel-Rockafellar) associated finite dimensional dual problem, which for the numerical computations is encoded as a finite semi-definite program (SDP) [16], [18]. The core of these approaches rely on the duality between peak locations and the existence of an interpolating trigonometric polynomial (often referred to as dual certificate) in the measurement (dual) space (which is bounded by 1 in magnitude at locations indicating the underlying peaks [17]). However, apart from the one dimensional spaces (line and torus), there is not a canonical extension or exact SPD formulation in higher dimensional spaces. We refer for instance to [18] (and the references therein) for relaxed SDP versions in higher dimensions.
In arbitrary spaces and for general kernels, one has to approximate the BLASSO problem by first introducing discrete grids and then solve a finite dimensional minimization problem (e.g. LASSO or basic pursuit). Several authors have proposed approximation or non-convex optimization steps to be included in the standard LASSO to recover the exact locations and amplitudes. In particular, in [32], [29], the continuous basis-pursuit which involves a first order Taylor approximation of the kernel in the fidelity term of the minimization problem accompanied by the ℓ 1 -norm regularization term have been utilized to improve the accuracy in the peak localization. However, multiple peaks around the original peak is a common result (as also theoretically justified in [29]). More robust iterative approaches using the Frank-Wolfe algorithm (also known as conditional gradient method) has been proposed in [13], [10]. These include an alternation between two steps. In the first step, the computational support is renewed by generating a new peak location using the conditional gradient method and then a non-convex step follows where only the locations and amplitudes are computed July 14, 2020 DRAFT while the number of peak locations stays fixed.
Even though the current paper is focusing on superresolution using convex optimization methods and especially the ℓ 1 -norm regularization, we would like to mention that there is also a vast literature on spectral superresolution algorithms that rely on Prony's concept (for a general review see [49], [73]), for example MUSIC [71], ESPRIT [66]) or pencil method [44]. These methods perform well in noiseless setting and do not require a minimum separation condition to fully recover positive and negative peaks; however they rely strongly on the signal, noise and measurement modelling and their extension to higher dimension is not trivial see e.g. [53], [63], [50], [2], [51], [23].

C. Contributions
In this work, the aim is to find new connections between the super-resolution algorithms which impose sparsity assumptions on the signal to be recovered [72] and theoretical studies (e.g. [56], [5], [29]) which have been developed rather separately so far. To that end, we first explain how convex optimization techniques [12], [41], [11] and, more precisely, ℓ 1 -norm sparsity constraints affect the solution of such inverse problems as the deconvolution of sparse peaks (or point sources) on discrete grids (or meshes) when the convolution kernel is smooth (admissible) [8] and then we propose an adaptive super-resolution scheme. In particular, our contributions are two-fold and are summarized as follows:

1) Theoretical:
• With the help of the first order optimality condition of the ℓ 1 -norm minimization problem, we show that the numerical solution consists of one or multiple peaks at grid points (or nodes) adjacent to the location of the actual peak. Our conclusions are inline with recent results presented in [29] for one dimensional spaces; but, the methodology employed, as well as the form in which the problems in question are expressed, are different. Previous approaches study the properties of the ℓ 1 -norm numerical solution by introducing the extended computational support notion [29] or by deriving dual certificates that fulfils particular properties [64] which both were used nicely to obtain asymptotic properties of the signal support. Here, we take a step forward to characterize the values of the numerical solution on its support also in dependence of the exact peak locations. We use the optimality condition of the finite dimensional ℓ 1 -norm minimization problem to investigate the conditions under which a single or multiple peaks are recovered in one dimensional spaces. Then, we define an optimality curve, directly related to July 14, 2020 DRAFT the optimality condition of the problem (respectively the dual certificate), whose shape allows us not only to justify the patterns of the expected numerical solutions on fixed computational grids both in one and higher dimensional spaces, but also infer further information on the location of the exact peaks.
• We show that we can explicitly approximate the locations and amplitudes of the exact peaks based on a set of linear equations derived from the associated normal equations of the ℓ 1 -norm problem.
• We derive an a-posterior error between the original signal (I.2) and its discrete version (I.3) by employing the Bregman distance [42]. We show that the recovery error depends on the relative distance between the computational grid points and the locations of the original peaks. The proposed adaptive algorithm shares some similarities with other superresolution algorithms e.g. [55], [81]; however, our approach embeds an automatic adaptation scheme [77] since it restricts and refines the grid in an unsupervised manner only in the neighborhoods where there is indication that a peak exists. This additionally allows to solve iteratively a small to medium size linear problem using convex optimization techniques [12], instead of a big size problem as in [55].
Overall, our analysis provides theoretical insights on the effects of gridding (a.k.a discetization of the parameter space), and it can help the algorithmic development in the direction of avoiding heuristic post-processing steps by using information about the convolution kernel properties and the formulation of the minimization problem rather than resorting in unreliable approximations as often happens in application papers (e.g. in [74]).

II. DECONVOLUTION OF SPARSE PEAKS BY CONVEX OPTIMIZATION
In the following we discuss the theoretical basis of sparse peak deconvolution using convex optimization approaches and put it in perspective with classical discretization issues in numerical analysis.

A. Sparsity over the continuum and its discretization
Let us start by formulating the problem over the continuum, following [13], [27], which is the underlying ideal sparse peak deconvolution to which we expect minimizers of (I.5) to converge to.
For a Radon measure µ on Ω we denote its total variation by The convex variational problem solved for sparse peak deconvolution in a continuum setting is then given by Now (I.5) can be interpreted as a discretization on a given grid, it can indeed be rephrased as It is straight-forward to show that the functionals J N Γ-converge to J ∞ , but one can also ask for more quantitative error estimates, which we shall discuss below.
By standard arguments we can verify the following result for the discretized problem (cf. [14], [13], [27] for analogous results on the continuum problem (II.2)): Proposition II.1. For λ ≥ 0 there exists a solution of (I.5). If λ > A * f ∞ then the unique minimizer is given by c = 0. If λ < A * f ∞ , then each solution is different from zero.
Proof. Convexity, coercivity, and nonnegativity immediately imply the existence of a minimizer.
July 14, 2020 DRAFT Now let λ > A * f ∞ , then with inequality only for c = 0. Hence, c = 0 is the unique minimizer. In the case λ < A * f ∞ we choose c = ǫA * f with ǫ > 0 sufficiently small to verify that there exists a c yielding a functional value lower than 1 2 f 2 .

B. Optimality conditions
As a next step we state the optimality conditions for (II.2) and the discrete version (I.5). Those are important for error estimates and further analysis in this paper.
Let us start with the sub-differential of the total variation norm, which is given by (cf. [13]) Here µ = µ + − µ − denotes the standard Jordan decomposition of the signed measure µ. Since the quadratic part of the functional J ∞ is differentiable and G is continuous, i.e. the convolution maps into the pre-dual of the space of Radon measures, we obtain the optimality condition where H = G * G. On the other hand, the optimality condition of the discrete problem (I.5) is where p ∈ R N is contained in the sub-differential of c 1 . The right hand side of the previous The optimality 1 The convolution of signal µ N with a kernel, e.g. Gaussian G, is G * µ N = N k=1 c k G(x − x k ) and in matrix form this can be expressed as [Ac]j = N k=1 c k G(xj − x k ). Moreover, G * G = H and since function G is symmetric and the convolution is July 14, 2020 DRAFT condition can be written as for k = j, p j ∈ sign(c j ) when c j = 0 and |p j | < 1 when c j = 0.
In order to highlight the connection with the continuum formulation, we rewrite the optimality solely for the measure µ N and deduce that We see that the main difference to the optimality condition in the continuum is that the first equality only holds on the grid points x k and not in the whole domain Ω. Note that due to the continuity of G and H one will expect (at least for sufficiently small grid size) that if |H * µ N −G * f | λ − 1 is strictly less than zero in a set of neighbouring grid points, then it remain less than zero in the area bounded by these points (further details are given in section III-B2). Hence, the main violation of the continuum optimality condition considered for µ N will appear close to grid points where it equals zero, usually corresponding to non-zero coefficients c k . This yields a first idea for using an adaptive computational grid. As we shall see below this can be further improved and backed up by a-posteriori error estimation.

C. A-Posteriori error estimate
In order to derive suitable error estimates for non-smooth convex variational problems such as problem (I.5), it is now a standard approach to use the Bregman distance as proposed in [60] (we refer to [42] for an overview). The Bregman distance for the total variation distance is given by for a subgradient q ∈ ∂ µ T V . Given a subgradientq ∈ ∂ μ T V , we will denote by the symmetric Bregman distance.
July 14, 2020 DRAFT The key idea here is to use the difference in the optimality conditions and take a duality product with the difference of the measures. For this sake, we use the following notation It is straightforward to see that q N ∈ ∂ µ N (x) T V and hence from (II.5), (II.9) and (II.13) we (II.14) Now the announced duality product with µ − µ N implies an a-posterior error estimate of the form Thus, we observe that only regions with r N = 0 contribute to the error between µ and µ N .
Moreover, via the optimality condition (II.9) of problem (I.5), we have that r N (x k ) = 0 for any grid point x k and thus r N , µ N = 0. So, we can write The previous expression shows that by reducing the supremum norm of r N is crucial for reducing the global error. This can be achieved by introducing finer computational grids. In Figure 1 we can observe that the value of r N ∞ decreases with respect to the size of the computational grid.
However, we remark that by introducing fixed fine grids (as proposed for instance in [81], [57], [55]) the computational performance and stability can be affected. To overcome computational limitations, we later propose to introduce a progressive grid refinement.

III. NUMERICAL SOLUTIONS IN THE CASES OF SINGLE PEAKS
Here, as started in [39] we where ξ ∈ R d and d ≥ 1 2 . We can easily interpret Proposition (II.1) in this case as ..,N and A * f ∞ > λ in order to obtain nonzero solutions. Hence, we need λ < γ max j {H(x j − ξ)} and since H attains its maximum at zero, a simple sufficient condition is given by which is also necessary in the case of ξ coinciding with a grid point. We will thus assume condition (III.2) throughout the whole section without further notice.

A. Exact recovery
The simplest case to start with, which can directly be treated in arbitrary dimensions, is that ξ coincides with one of the grid points. In this case we obviously expect perfect reconstruction, which is confirmed by the following result: Then there exists a one-sparse solution 2 Our analysis is based on a single positive peak which is often the case in image processing applications. We note that the theorems/conclusions presented in this section are valid also for a negative peak. July 14, 2020 DRAFT Proof. Without restriction of generality assume that γ > 0. In order to prove the assertion, we have to check whether the optimality condition of (I.5) holds under the assumptions mentioned above.
From our prior computations (II.8), the optimality condition (II.8) reduces to We have to differentiate between the cases where j = K and j = K.
For j = K the optimality condition (III.3) is due to the fact that H attains its maximum at zero. Hence in both cases the optimality condition is fulfilled and we obtain the assertion.
Therefore, the reconstruction of the support of a delta peak is exact if the position of the peak coincides with a grid point and the regularization parameter is small enough.

B. Recoveries for off-the-grid peaks
Let us consider the more frequent case where µ = γδ ξ is located among a set of grid points Here, with the help of the optimality conditions (II.9), we define a so-called optimality curve p(x) given by where f = G * µ, H = G * G (smooth and symmetric) and µ N = x k ∈N c k δ x k is the nontrivial numerical solution obtained from the minimization problem (I.5). The optimality curve p(x) can be regarded as analogous to the T V dual certificate [65] for the ℓ 1 -norm minimization problem that will allow us to understand the expected patterns of the numerical solutions around ξ.
To ease our analysis, we rewrite p(x) as (III.6) July 14, 2020 DRAFT We note that p(x j ) < 0 when c j = 0 and p(x j ) = 0 when c j > 0.

1) Single spatial dimension:
In one dimensional spaces, we consider that signal µ = γδ ξ is located between two grid points, i.e. ξ ∈ (x K , x K+1 ). For the following consideration, the interval length h will be defined as By employing function p(x) (III.6) in the vicinity of peak ξ, we can prove the following Theorem.
Let H ∈ C 3 (R) be nonnegative with a unique maximum at zero and let h be sufficiently small.
there exists a solution of (I.5), which can be written nonzero and of the same sign as γ.
The proof of Theorem III.2 is given in Appendix A. Figure 2 illustrates the assertion of Theorem III.2. Note that due to the symmetry of H, the analogous claim holds for ξ in the other half of the interval. Figure 3 depicts the optimality curve p(x) for a positive peak when H(x) is a Gaussian Fig. 2: If ξ is in the blue interval, the reconstructed solution µ N consists of only one peak. In the case that ξ is located in the red interval, then one recovered peak is not sufficient.
kernel. The curve is downward concave in the area around ξ which implies that there are at most two points on x-axis where p(x) = 0. From these points, at least one is the grid point with the nonzero coefficient of µ N . We can observe that the number of the recovered peaks depends on the distance between the location of the exact peak (denoted by red x) and the neighboring grid points (given fixed λ = 0.1λ max ). 2) Higher spatial dimension: In higher spatial dimensions, the topological structure is more complicated which makes a rigorous proof by analogous arguments impossible. However, we can at least make some formal arguments and computational experiments concerning the optimality curve p(x) (III.6). First of all we expect that for γ > 0, p(x) is concave around ξ and nonzero entries are only found in the convex hull of ξ on the grid, i.e. the largest convex hull C ξ that can be formed of a set of grid points surrounding the peak ξ such that no other grid point is contained The following observations can be made: • Given any location x far from ξ, we have that p(x) < 0, since the positive term H(x − ξ) is • Since the kernel H(x) is smooth, and considering that its width is greater than the resolution of the selected grid (which is often the case for low resolution images), then p(x) < 0 in the area bounded by grid points where c j = 0.
• Given {x j | c j > 0} forms a small neighborhood of ξ, consisting of N ≥ 1 points, we can make a local Taylor expansion similar to the one-dimensional case. First of all we have that p(x j ) = 0 for all such j. By summing those with respect to j we get July 14, 2020 DRAFT Using the lowest order approximation for small arguments we find i.e. to first order Using this approximation the Hessian can be computed to leading order as which is negative definite due to our assumptions on H. Hence, p is concave in a neighbourhood of ξ, which implies that its level sets are convex. The points x j with p(x j ) = 0 are on the Thus, from those arguments we see that the active grid points (c j > 0) are to be expected in the convex hull of ξ on the grid. This can be made rigorous under the assumption that the local grid size around ξ is small and there are no active grid points at large distance from ξ, which is confirmed in all our numerical experiments. Figure 4 illustrates this behaviour by showing the shape of the function p and its relationship to the nonzero coefficients of the reconstructed signal µ N . In this figure, the small blue dots depict the computational grid, the big blue circles show the grid points with nonzero entries (i.e. estimated peaks). For the computations, the regularization parameter was set λ = 0.1λ max and H(x) was Gaussian with standard deviation σ = 1.5 √ 2h (where h was the grid resolution). Based on the previous analysis, we can see in Figure 5 that the numerical solution depends on λ and the properties of kernel H. As expected, the number of active grid points increases as λ decreases which is effectively a property of the finite-dimensional where all the c k have same sign s ∈ {+1, −1}.
July 14, 2020 DRAFT The associated normal equations are we can write the associated normal equations as For a small neighbourhood around ξ (and h → 0 a bound for the grid size) we can perform a Taylor-expansion and obtain where we have used ∇H(0) = 0. We observe that all equations have the same leading order term, which yields up to order two In order to access higher-order terms we can exploit the fact that set N of nonzero coefficients has more than one grid point and thus we can estimate differences of equation (IV.3) for pairs of grid points x i , x j ∈ N . This yields we can uniquely solve for the location ξ and obtain a second order approximation in h (considering the contribution fromf negligible).

A. Examples of peak recoveries
To demonstrate the previous theoretical results, we present some examples in one and two dimensional spaces.

1) 1D spaces:
In the following examples, we consider a signal with three peaks with amplitudes γ 1 , γ 2 and γ 3 at positions ξ 1 , ξ 2 and ξ 3 where ξ l ∈ (0, 1) for l = 1, 2, 3. The signal is given by Moreover, we choose a Gaussian convolution kernel G with standard deviation σ = 0.03. The continuous convolved data can be expressed analytically as For the estimation of the numerical solution µ N the domain [0, 1] is discretized and the ℓ 1 -norm minimization problem (I.5) is solved with λ = 0.01 A * f ∞ on a uniform grid of size N.
We first consider a grid that includes ξ 1 , ξ 2 and ξ 3 . In Figure 6, we observe that the exact recovery is feasible (which is in accordance with preposition (III.1)).
Now we can consider the case of Theorem III.2 where the three peaks of µ are located between the grid points. Figure 7   We note that in the case where the numerical solution yields to a single nonzero coefficient c k at x k , thenγ l = c k andξ l = x k . Table I summarizes the values of the amplitudes and locations of the original and estimated peaks for the two different computational grids of Figure 7 and 8. Based on these results, the reconstructions in a fine grid are slightly more accurate than ones obtained using a coarse grid, which is inline with the a-posteriori error analysis presented in section II-C.

Based on equation (IV.5), the peak location is approximated aŝ
2) 2D spaces: In higher dimensions, the exact signal is µ = L l=1 γ l δ ξ l where ξ l ∈ Ω ⊂ R d (d > 1) and the estimated µ N solution has nonzero values clustered in grid points around the locations of the actual peaks ξ l . If a cluster of grid points with nonzero coefficients around peak ξ l is denoted by N l = {x k } 1:N l , we approximate the peak amplitude according tô    Fig. 7 and Fig 8. The first and second column show the locations ξ and amplitudes γ of the underlying peaks, then there are the estimated locations and amplitudes for the cases where a coarse grid (N=16 points) and a fine grid (N=51 points) were used.
The location is approximated similarly as in (IV.5). Particularly, if there are at least two grid points x i and x j ∈ N l , we havê where  Figure 9). The middle image of Figure 9 shows the numerical result obtained solving the ℓ 1 -norm minimization problem. The intense blue circles illustrate the locations where nonzero entries appeared. We can observe that there are four distinctive clusters of grid points with nonzero coefficients.
Therefore, four peaks are approximated, one for each cluster using (IV.9). The peak approximations are shown in the right image of Figure 9. Table II

V. RECOVERIES IN THE CASE OF MULTIPLE PEAK SIGNALS
A coarse computational grid, even though reduces the computational cost, imposes some limitations to detect and separate neighbouring peaks. For example, there is always a possibility that there are more than one positive peak in an interval between two grid points (see one dimensional example of Figure 10.A), or two or more original peaks may be located in adjacent intervals (e.g. Figure 10.B). Then the numerical solution µ N of the ℓ 1 -norm minimization problem may not be accurate enough.
Following similar analysis as in Theorem (III.2), we can easily show that there exists a ℓ 1 -norm as depicted in Figure 10.A even though the original peaks are two. Additionally, when the original peaks are distributed between two intervals then we can expect up to three reconstructed peaks as in Figure 10.B and Figure 10.C. A natural way to improve the estimates is by refining the grid. Figure 11 illustrates how by performing local refinements on the grid (and fitting the input data with a solution in the updated grid), we can achieve a separation of the underlying peaks.
A. Adaptive super-resolution for sparse signal 1) Overview: In the following context, to ease our analysis and to proceed with the domain refinement in higher dimensions instead of using grid/points we use the mesh/nodes notion as in the finite element methods. Hence, the computational domain can be described by a mesh consisting of a set of nodes (equivalent to grid points) and elements (e.g. line segments in one dimension or July 14, 2020 DRAFT Fig. 11: By restricting and refining iteratively the computational grid around the nonzero coefficients of the numerical solution µ N , we can determine small disjoint intervals where the two original peaks (in red) belong to and thus separate them. Fig.(A) shows the observations, original peaks (in red) and the locations of the used computational grid (small blue dots). A more analytical description of the approach is given in Appendix B. We need to mention that the proposed scheme is applicable for peaks of similar sign or when the positive and negative peaks satisfy the separation criteria (which is based on the kernel's width and noise type/level) as studied in some cases for example in [29], [64].

VI. RESULTS AND DISCUSSION
To demonstrate how the proposed superresolution approch can be used, we reconstruct super-

A. Simulated data
The simulations were carried out in a two dimensional square domain Ω = [0, 1] 2 . The aim was to approximate the locations and amplitudes of an original signal µ = L l=1 γ l δ ξ l from low resolution images W ∈ R M ×M where As a convolution kernel, we use the one from study [55], given by Also, we considered a low additive measurement noise ε = scε ∈ R M ×M whereε ∈ R M ×M was sampled from a Gaussian distribution with zero mean and variance one. The scaling parameter sc July 14, 2020 DRAFT was estimated based on the level of the signal-to-noise ratio (SNR), SNR = 10 log 10 In the following simulations, we used SNR = 40 dB.

B. Details about the adaptive superresolution approach
We estimated the locations and the amplitudes of the underlying peaks by employing the proposed scheme of section V. The initial estimation (by solving the ℓ 1 -norm minimization) was performed in a uniform mesh of N × N nodes. Then, the mesh was updated automatically around the nonzero entries of vector c. In practice, to avoid small numerical inaccuracies, the new domain was defined by keeping the nodes with absolute values of the estimated peaks greater than a small threshold (i.e. 0.5% of the maximum |c| of vector c). In the current implementations, the ℓ 1 -norm minimization problem was solved using the hierarchical adaptive lasso (HAL) [58]. Other algorithms e.g. [10], [48] could be used as well. Here, we used HAL to reduce the amplitude shrinkage of the estimated nonzero coefficient given a regularization parameter λ. In the following examples, λ = 0.1λ max .
Moreover, the incorporation of a Bregman iteration [80] could be considered in the future for cases with relative high measurement noise.
The updates of the computational support terminated when the distance between the existing nodes and the additional nodes became small. In the following examples, we used as a criterion for adding a new node, the minimum distance of this candidate node from the existing nodes, h min = 0.25h M (approximately 0.125 of the Gaussian kernel's standard deviation). This choice was made to enable a computational efficiency (i.e. a reasonable number of iterations) and to allow a good approximation of the peaks using small clusters of the nonzero coefficients and (IV.9).

C. Comparison metrics
In tests with only few peaks, we used: For dense distributions of peaks, we employed the earth mover's distance (EMD) (Wasserstein metric) as a measure of dissimilarity between the original and the estimated peaks locations [68], [62].

D. Examples
Three different examples are presented to demonstrate the different stages of the proposed scheme.
In the first example, we show step by step the estimation of the locations and amplitudes of five positive peaks (see Figure 12). In the second example, presented in Figure 13, we use the proposed apprach to recover both positive and negative peaks. In the last example in Figure 14, we illustrate the potential of the algorithm to deal with denser peak distributions.
In the example of Figure 12  The reconstruction results (amplitudes and locations) are summarized in Table III  To show that the proposed scheme can be used to recover both positive and negative peaks, Figure 13 illustrates the estimation of two positive and two negative peaks. The original peaks are marked with + for the positive peaks and x for the negative peaks and the estimated ones with circles and squares respectively in the right-hand side image of Figure 13. Table IV Figure 12, we present the progressive steps of the approach to separate the underlying peaks. Also, the image on the lowermost left side shows the grid points corresponding to non-zero coefficients in the clusters (in grey) of the 9 th iteration of the approach indicate that, in a low noise regime, the proposed superresolution adaptive scheme can recover as many peaks as the exact number of them in most of the cases when their distance does not violate an underlying separation condition e.g. minimum distance ∆ min ≥ Cσ for Gaussian kernels [64].
Our numerical simulations showed that, when the noise level was SNR=40dB, the correct number of peaks could be recovered if ∆ min /σ > 0.2 given h M . We remark that even though the effect of  the measurement noise (either Gaussian or Poisson) and the measurement/observation sampling can affect the recovery of a multi-dimensional signal, these mathematical questions admit of different July 14, 2020 DRAFT analysis than the currently addressed questions and will be considered in a follow-up study.

E. Discussion and future prospects
Questions regarding the deconvolution of sparse peaks present great mathematical difficulties with some of them investigated in studies such as in [26], [27], [28], [22], [65]. However, even though very important, these studies often do not accommodate easily accessible solutions to software developers and engineers working on superresolution applications. The current work aimed to shed light into some of these theoretical findings and put them into perspective with practical solutions in superresolution algorithms.
In the theoretical part of this study, we explained why clusters of nonzero peaks appear around the locations of the original peaks when we solve the ℓ 1 -norm minimization problem on a discrete grid and in which parts of the grid these nonzero peaks are more likely to appear. Moreover, we showed how the locations of the underlying peaks are connected with the numerical solution using the optimality condition of the ℓ 1 -norm minimization problem. One important remark here is that the distribution of the nonzero coefficients of the numerical solution depend on the properties of the convolution kernel H = G * G (III.6) and not on kernel G (I.1). Therefore thinking in more general terms, for an inverse problem with forward operator R, it is only important that R * R is a convolution to anticipate a numerical solution following similar pattern as in the current problem.
This could be true for example in classical tomography e.g. filtered back-projection [59].
In general, we envision that similar superresolution schemes can be performed in a wide variety of inverse problems in the fields of geophysics, astronomy and spectroscopy [47], [35], [46] because many of these applications share the same characteristics and properties as this deconvolution problem. However, some new aspects need to be investigated. For example, in neuroimaging, the EEG source imaging problem, even though it shares seemingly similarities with the current deconvolution, problem, is a severely ill-posed problem where the forward operator has a singularity and its computational version has a matrix with a large null space [54]. Therefore, special design of the prior model (e.g. weighting) is required whereas the expected pattern of the numerical reconstructions has to be studied carefully.
In the application part, the main two novelties of the proposed approach were a) the automatic adaptation of the computational domain (using elements) and b) the approximation of the underlying peaks using a numerical approximation of the ℓ 1 norm optimality condition that stemmed from the findings of our theoretical analysis. As a natural next step though, within the microscopy field July 14, 2020 DRAFT we expect to compare the proposed approach with other state-of-the art algorithms [72] which lie either on the variational or spectral framework, e.g. Alternating Descent Conditional Gradient Method [10] or MUSICAL [1] respectively.
Moreover, we are considering extension that could improve the algorithmic performance, for example, the incorporation of a non-convex step as in [10] to possibly speed up the convergence.
In that step possibly, the optimality curve (III.6) could guide the update of the computational domain simultaneously in multiple locations. Furthermore, the idea of employing the ensemble learning or committee method [58] (which allows to estimate a weighted solution in each discretization level by solving multiple deconvolution problems in a parallel fashion), could help to reduce possible bias introduced due to the regularization or high measurement noise.

VII. CONCLUSIONS
The current work bridges the gap between theoretical studies and implementations of algorithms that impose sparsity constraints on the signal to be recovered. First, we studied theoretically the deconvolution of single peaks using the ℓ 1 -norm and we confirmed recent observations that a discrete reconstruction yields to multiple peaks at grid points adjacent to the location of the actual peak. We showed that by using these adjacent peaks and the first order optimality condition of this convex problem, we can obtain a set of linear equations to approximate the location of the actual peak. We quantified the errors between the continuous (TV) problem (that allows exact peak recoveries) and the finite ℓ 1 -norm minimization problem, which designated that the accuracy of the numerical estimates depends on the discretization that can be improved by applying finer gridding.
Second, using the previous theoretical finding we proposed an iterative scheme in which automated local refinement on the computational grid was performed to identify the areas where the true peaks were located. Then, with the help of the equations from the optimality condition, the peak locations and amplitudes were estimated. Finally, low resolution images, obtained using simulated focal sources convoluted with a smooth kernel, were used to show that our approach can increase the spatial resolution by allowing the separation and localization of these focal sources.

A. Proof of Theorem III.2
Proof.
Let us first consider that the reconstructed signal is µ N = αδ x K with a > 0. When j = K, p K = 1 July 14, 2020 DRAFT in the optimality condition (II.8) since p K ∈ ∂c K and a = c K > 0. Thus, the optimality condition (II.8) reduces to Now let us consider the case where j = K + 1, then (II.8) becomes This is equivalent to Thus, we obtain which is true since we have λ < γH(x K − ξ) and H(x K − ξ) < H(0). Now assume that h is sufficiently small in the latter case and make the Ansatz µ N = c K δ x K µ N = + c K+1 δ x K+1 . Without restriction of generality we consider γ > 0 hence we look for c K > 0 and c K+1 > 0, the other sign is analogous. We extend the vector c by c j = 0 for j / ∈ {K, K + 1} and verify that it is a minimizer of J in (I.5) by constructing an appropriate subgradient in the optimality condition (II.7).
In particular, we have The conditions p K = p K+1 = 1 lead to the following 2 × 2 system Steps 1-2 are repeated until the computational support is not longer updated (the number of nodes and elements stays fixed). This happens when the distance between the nodes becomes small. As a minimum distance (between two nodes) we can use a limit for peak separation presented for some convolution kernels in [64], [65]). Alternatively, prior information about the expected size of the underlying peaks (e.g. in microscopy the sizes of the molecules) can be considered. Then, we recover as many peaks as the number of the disjoint clusters Ωl, forl = 1, . . . ,L (wherê L is the total number of disjoint clusters estimated in the last iteration). The amplitude denoted byγl in cluster Ωl follows from equation (IV.8) using the coefficients entries of that cluster. For the approximation of the the peak location in cluster Ωl, we first check the number of nonzero July 14, 2020 DRAFT coefficients denoted by Nl. If Nl > d (where d is the dimensionality of the problem), the peak location,ξ l , is estimated by solving a linear system formed using the expression (IV.9). Now, if 2 ≤ Nl ≤ d then the peak location can be approximated with the help of linear basis functions, φ k .
If the approximated location is expressed asξl = Nl k=1 φ k x k then by inserting the previous linear representation forξl in expression (IV.9), we can obtain an approximation for the peak location. If Nl = 1, thenξl equals to the value of the nonzero node (following from (IV.3)).
ACKNOWLEDGEMENT AK was supported by the Academy of Finland Postdoctoral Researcher program (No 316542).