Super-Resolution from Noisy Data

Introduction

It is often of great interest to study the fine details of a signal at a scale beyond the resolution provided by the available measurements. In a general sense, super-resolution techniques seek to recover high-resolution information from coarse-scale data. There is a gigantic literature on this subject as researchers try to find ways of breaking the diffraction limit—a fundamental limit on the possible resolution—imposed by most imaging systems. Examples of applications include conventional optical imaging , astronomy , medical imaging , and microscopy . In electronic imaging, photon shot noise limits the pixel size, making super-resolution techniques necessary to recover sub-pixel details . Among other fields demanding and developing super-resolution techniques, one could cite spectroscopy , radar , non-optical medical imaging and geophysics .

In many of these applications, the signal we wish to super-resolve is a superposition of point sources; depending upon the situation, these may be celestial bodies in astronomy , molecules in microscopy , or line spectra in speech analysis . A large part of the literature on super-resolution revolves around the problem of distinguishing two blurred point sources that are close together, but there has been much less analysis on the conditions under which it is possible to super-resolve the location of a large number of point sources with high precision. This question is of crucial importance, for instance, in fluorescence microscopy. Techniques such as photoactivated localization microscopy (PALM) or stochastic optical reconstruction microscopy (STORM) are based on the use of probes that switch randomly between a fluorescent and a non-fluorescent state. To super-resolve a certain object, multiple frames are gathered and combined. Each frame consists of a superposition of blurred light sources that correspond to the active probes and are mostly well separated.

In the companion article , the authors studied the problem of recovering superpositions of point sources in a noiseless setting, where one has perfect low-frequency information. In contrast, the present paper considers a setting where the data are contaminated with noise, a situation which is unavoidable in practical applications. In a nutshell, proves that with noiseless data, one can recover a superposition of point sources exactly, namely, with arbitrary high accuracy, by solving a simple convex program. This phenomenon holds as long as the spacing between the sources is on the order of the resolution limit. With noisy data now, it is of course no longer possible to achieve infinite precision. In fact, suppose the noise level and sensing resolution are fixed. Then one expects that it will become increasingly harder to recover the fine details of the signal as the scale of these features become finer. The goal of this paper is to make this vague statement mathematically precise; we shall characterize the estimation error as a function of the noise level and of the resolution we seek to achieve, showing that it is in fact possible to super-resolve point sources from noisy data with high precision via convex optimization.

To formalize matters, we have observations about an object $x$ of the form

where $t$ is a continuous parameter (time, space, and so on) belonging to the $d$-dimensional cube $^{d}$. Above, $z$ is a noise term which can either be stochastic or deterministic, and $Q_{\text{lo}}$ is a bandlimiting operator with a frequency cut-off equal to $f_{\text{lo}}=1/\lambda_{\text{lo}}$. Here, $\lambda_{\text{lo}}$ is a positive parameter representing the finest scale at which $x$ is observed. To make this more precise, we take $Q_{\text{lo}}$ to be a low-pass filter of width $\lambda_{\text{lo}}$ as illustrated at the top of Figure 1; that is,

such that in the frequency domain the convolution equation becomes

Here and henceforth we denote the usual Fourier transform of a measure or function $g$, provided that it exists, by $\hat{g}(f)=\int e^{-i2\pi\langle f,t\rangle}g(\text{d}t)$. The spectrum of the low-pass kernel $\widehat{K}_{\text{lo}}(f)$ vanishes outside of the cell $[-f_{\text{lo}},f_{\text{lo}}]^{d}$.

Our goal is to resolve the signal $x$ at a finer scale $\lambda_{\text{hi}}\ll\lambda_{\text{lo}}$. In other words, we would like to obtain a high-resolution estimate $x_{\text{est}}$ such that $Q_{\text{hi}}\,x_{\text{est}}\approx Q_{\text{hi}}\,x$, where $Q_{\text{hi}}$ is a bandlimiting operator with cut-off frequency $f_{\text{hi}}:=1/\lambda_{\text{hi}}>f_{\text{lo}}$. This is illustrated at the bottom of Figure 1, which shows the convolution between $K_{\text{hi}}$ and $x$. A different way to pose the problem is as follows: we have noisy data about the spectrum of an object of interest in the low-pass band $\left[-f_{\text{lo}},f_{\text{lo}}\right]$, and would like to estimate the spectrum in the possibly much wider band $\left[-f_{\text{hi}},f_{\text{hi}}\right]$. We introduce the super-resolution factor (SRF) as:

in words, we wish to double the resolution if the SRF is equal to two, to quadruple it if the SRF equals four, and so on. Given the notorious ill-posedness of spectral extrapolation, a natural question is how small the error at scale $\lambda_{\text{hi}}$ between the estimated and the true signal $K_{\text{hi}}\ast(x_{\text{est}}-x)$ can be? In particular, how does it scale with both the noise level and the SRF? This paper addresses this important question.

2 Models and methods

As mentioned earlier, we are interested in superpositions of point sources modeled as

The measurement error $z$ is otherwise arbitrary and can be adversarial. For concreteness, we set $K_{\text{lo}}$ to be the periodic Dirichlet kernel

3 Main result

Our objective is to approximate the signal up until a certain resolution determined by the width of the smoothing kernel $\lambda_{\text{hi}}>\lambda_{\text{lo}}$ used to compute the error. To fix ideas, we set

to be the Fejér kernel with cut-off frequency $f_{\text{hi}}=1/\lambda_{\text{hi}}$. Figure 2 shows this kernel together with its spectrum.

As explained in Section 3.2 of , no matter what method is used to achieve super-resolution, it is necessary to introduce a condition about the support of the signal, which prevents the sources from being too clustered together. Otherwise, the problem is easily shown to be hopelessly ill-posed by leveraging Slepian’s work on prolate spheroidal sequences . In this paper, we use the notion of minimum separation.

Our model (1.3) asserts that we can achieve a low-resolution error obeying

but that we cannot do better as well. The main question is: how does this degrade when we substitute the low-resolution with the high-resolution kernel?

Assume that the support $T$ of $x$ obeys the separation condition

where $C_{0}$ is a positive numerical constant.

Thus, minimizing the total-variation norm subject to data constraints yields a stable approximation of any superposition of Dirac measures obeying the minimum-separation condition. When $z=0$, setting $\delta=0$ and letting SRF$\rightarrow\infty$, this recovers the result in which shows that $x_{\text{est}}=x$, i.e. we achieve infinite precision. What is interesting here is the quadratic dependence of the estimation error in the super-resolution factor.

We have chosen to analyze problem (1.5) and a perturbation with bounded $L_{1}$ norm for simplicity, but our techniques can be adapted to other recovery schemes and noise models. For instance, suppose we observe noisy samples of the spectrum

where $\epsilon_{k}$ is an iid sequence of complex-valued $\mathcal{N}(0,\sigma^{2})$ variables (this means that the real and imaginary parts are independent $\mathcal{N}(0,\sigma^{2})$ variables). This is equivalent to a line-spectra estimation problem with additive Gaussian white noise, as we explain below. In order to super-resolve the signal under this model, we propose the following convex program

which can be implemented using off-the-shelf software as discussed in Section 3. A corollary to our main theorem establishes that with high probability solving this problem allows to super-resolve the signal despite the added perturbation with an error that scales with the square of the super-resolution factor and is proportional to the noise level.

Fix $\gamma>0$. Under the stochastic noise model (1.8), the solution to problem (1.9) with $\delta=\left(1+\gamma\right)\sigma\sqrt{4f_{\emph{\text{lo}}}+2}$ obeys

with probability at least $1-e^{-2f_{\emph{\text{lo}}}\gamma^{2}}$.

This result is proved in Section C of the appendix.

4 Extensions

Other high-resolution kernels. We work with the high-resolution Fejér kernel but our results hold for any symmetric kernel that obeys the properties (1.11) and (1.12) below, since our proof only uses these simple estimates. The first reads

This is to make sure that (2.6) holds. (For the Fejér kernel, we can take $f$ to be quadratic in $\left[0,1/2-\lambda_{\text{hi}}\right]$ and constant in $\left[1/2-\lambda_{\text{hi}},1/2\right]$.)

Higher dimensions. Our techniques can be applied to establish robustness guarantees for the recovery of point sources in higher dimensions. The only parts of the proof of Theorem 1.2 that do not generalize directly are Lemmas 2.4, 2.5 and 2.7. However, the methods used to prove these lemmas can be extended without much difficulty to multiple dimensions as described in Section D of the Appendix.

Spectral line estimation. Swapping time and frequency, Theorem 1.2 can be immediately applied to the estimation of spectral lines in which we observe

where $\alpha$ is a vector of complex-valued amplitudes and $z$ is a noise term. Here, our work implies that a non-parametric method based on convex optimization is capable of approximating the spectrum of a multitone signal with arbitrary frequencies, as long as these frequencies are sufficiently far apart, and furthermore that the reconstruction is stable. In this setting, the smoothed error quantifies the quality of the approximation windowed at a certain spectral resolution.

5 Related work

Earlier work on the super-resolution problem in the presence of noise studied under which conditions recovery is not hopelessly ill-posed, establishing that sparsity is not sufficient even for signals supported on a grid . More recently, studies the local stability of the problem in a continuous domain. These works, however, do not provide any tractable algorithms to perform recovery.

Since at least the work of Prony , parametric methods based on polynomial rooting have been a popular approach to the super-resolution of trains of spikes and, equivalently, of line spectra. These techniques are typically based on the eigendecomposition of a sample covariance matrix of the data . The theoretical analysis available for these methods is based on an asymptotic characterization of the sample covariance matrices under Gaussian noise , which unfortunately does not allow to obtain explicit guarantees on the recovery error beyond very simple cases involving one or two spikes. Other works extend these results to explore the trade-off between resolution and signal-to-noise ratio for the detection of two closely-spaced line spectra or light sources . A recent reference , which focuses mainly on the related problem of imaging point scatterers, analyzes the performance of a parametric method in the case of signals sampled randomly from a discrete grid under the assumption that the sample covariance matrix is close enough to the true one. In general, parametric techniques require prior knowledge of the model order and rely heavily on the assumption that the noise is white or at least has known spectrum (see Chapter 4 of ). An alternative approach that overcomes the latter drawback is to perform nonlinear least squares estimation of the model parameters . Unfortunately, the resulting optimization problem has an extremely multimodal cost function, which makes it very sensitive to initialization .

Proof of Theorem 1.2

It is useful to first introduce various objects we shall need in the course of the proof. We let $T=\{t_{j}\}$ be the support of $x$ and define the disjoint subsets

here, $\lambda\in\left\{\lambda_{\text{lo}},\lambda_{\text{hi}}\right\}$, and $j$ ranges from 1 to $\left|T\right|$. We write the union of the sets $S_{\text{near}}^{\lambda}\left(j\right)$ as

for any measure $\mu$ and $\lambda\in\left\{\lambda_{\text{lo}},\lambda_{\text{hi}}\right\}$. Finally, we reserve the symbol $C$ to denote a numerical constant whose value may change at each occurrence.

Set $h=x-x_{\text{est}}$. The error obeys

and has bounded total-variation norm since $\left|\left|h\right|\right|_{\text{TV}}\leq\left|\left|x\right|\right|_{\text{TV}}+\left|\left|x_{\text{est}}\right|\right|_{\text{TV}}\leq 2\left|\left|x\right|\right|_{\text{TV}}$. Our aim is to bound the $L_{1}$ norm of the smoothed error $e:=K_{\text{hi}}*h$,

We begin with a lemma bounding the total-variation norm of $h$ ‘away’ from $T$.

Under the conditions of Theorem 1.2, there exist positive constants $C_{a}$ and $C_{b}$ such that

This lemma is proved in Section 2.1 and relies on the existence of a low-frequency dual polynomial constructed in to guarantee exact recovery in the noiseless setting.

To develop a bound about $\|e\|_{L_{1}}$, we begin by applying the triangle inequality to obtain

By a corollary of the Radon-Nykodim Theorem (see Theorem 6.12 in ), it is possible to perform the polar decomposition $P_{S_{\text{far}}^{\lambda_{\text{hi}}}}\left(h\right)\left(\text{d}\tau\right)=e^{i2\pi\theta\left(\tau\right)}\left|P_{S_{\text{far}}^{\lambda_{\text{hi}}}}\left(h\right)\right|\left(\text{d}\tau\right)$ such that $\theta\left(\tau\right)$ is a real function and $\left|P_{S_{\text{far}}^{\lambda_{\text{hi}}}}\left(h\right)\right|$ is a positive measure. Then

where we have applied Fubini’s theorem and (1.11) (note that the total-variation norm of $\left|P_{S_{\text{far}}^{\lambda_{\text{hi}}}}\left(h\right)\right|$ is bounded by $2\left|\left|x\right|\right|_{\text{TV}}<\infty$).

In order to control the second term in the right-hand side of (2.1), we use a first-order approximation of the super-resolution kernel provided by the Taylor series expansion of $\psi\left(\tau\right)=K_{\text{hi}}\left(t-\tau\right)$ around $t_{j}$: for any $\tau$ such that $\left|\tau-t_{j}\right|\leq 0.16\lambda_{\text{hi}}$, we have

Applying this together with the triangle inequality, and setting $t_{j}=0$ without loss of generality, give

(To be clear, we do not lose generality by setting $t_{j}=0$ since the analysis is invariant by translation; in particular by a translation placing $t_{j}$ at the origin. To keep things as simple as possible, we shall make a frequent use of this argument.) We then combine Fubini’s theorem with (1.11) to obtain

Some simple calculations show that (1.11) and (1.12) imply

for a positive constant $C_{4}$. This together with Fubini’s theorem yield

for any $u$. In order to make use of these bounds, it is necessary to control the local action of the measure $h$ on a constant and a linear function. The following two lemmas are proved in Sections 2.2 and 2.3.

Take $T$ as in Theorem 1.2 and any measure $h$ obeying $\left|\left|Q_{\emph{\text{lo}}}h\right|\right|_{L_{1}}\leq 2\delta$. Then

Take $T$ as in Theorem 1.2 and any measure $h$ obeying $\left|\left|Q_{\emph{\text{lo}}}h\right|\right|_{L_{1}}\leq 2\delta$. Then

We may now conclude the proof of our main theorem. Indeed, the inequalities (2.2), (2.3), (2.4), (2.5) and (2.7) together with $I_{S_{\text{near}}^{\lambda_{\text{hi}}}}\left(h\right)\leq I_{S_{\text{near}}^{\lambda_{\text{lo}}}}\left(h\right)$ imply

where the second inequality follows from Lemma 2.1.

The proof relies on the existence of a certain low-frequency polynomial, characterized in the following lemma which recalls results from Proposition 2.1 and Lemma 2.5 in .

Invoking a corollary of the Radon-Nykodim Theorem (see Theorem 6.12 in ), it is possible to perform a polar decomposition of $P_{T}h$,

Next, since $q$ interpolates $e^{-i\phi\left(t\right)}$ on $T$,

Applying (2.10) in Lemma 2.4 and Hölder’s inequality, we obtain

Set $t_{j}=0$ without loss of generality. The triangle inequality and (2.9) in Lemma 2.4 yield

Combining (2.12), (2.13) and (2.14) gives

Observe that we can substitute $\lambda_{\text{lo}}$ with $\lambda_{\text{hi}}$ in (2.12) and (2.14) and obtain

This follows from using (2.9) instead of (2.10) to bound the magnitude of $q$ on $S_{\text{far}}^{\lambda_{\text{hi}}}$.

These inequalities can be interpreted as a generalization of the strong null-space property used to obtain stability guarantees for super-resolution on a discrete grid (see Lemma 3.1 in ). Combined with the fact that $\hat{x}$ has minimal total-variation norm among all feasible points, they yield

2 Proof of Lemma 2.2

The proof of this lemma relies upon the low-frequency polynomial from Lemma 2.4 and the fact that $q(t)$ is close to the chosen sign pattern when $t$ is near any element of the support. The following intermediate result is proved in Section A of the Appendix.

There is a polynomial $q$ satisfying the properties from Lemma 2.4 and, additionally,

where $\theta_{j}\in[0,2\pi)$. We set $v_{j}=e^{i\theta_{j}}$ in Lemma 2.4 and apply the triangular inequality to obtain

for all $t_{j}\in T$. By another application of the triangle inequality and (2.11)

To bound the remaining term in (2.15), we apply Lemma 2.5 with $t_{j}=0$ (this is no loss of generality),

It follows from this, (2.15) and (2.16) that

3 Proof of Lemma 2.3

Proof Note that in the interval $\left[0.16\lambda_{\text{hi}},0.16\lambda_{\text{lo}}\right]$, $t/0.16\,\lambda_{\text{hi}}\geq 1$, whence

We now turn our attention to the proof of Lemma 2.3. By the triangle inequality,

The second term is bounded via Lemma 2.6. For the first, we use an argument very similar to the proof of Lemma 2.2. Here, we exploit the existence of a low-frequency polynomial that is almost linear in the vicinity of the elements of $T$. The result below is proved in Section B of the Appendix.

where $\theta_{j}\in[0,2\pi)$, $t_{j}\in T$, and set $v_{j}=e^{i\theta_{j}}$ in Lemma 2.7. Again, suppose $t_{j}=0$. Then

The inequality (2.18) and Hölder’s inequality allow to bound the first term in the right-hand side of (2.20),

Another application of the triangular inequality yields

We employ Hölder’s inequality, (2.11), (2.18) and (2.19) to bound each of the terms in the right-hand side. First,

Combining (2.17) with these estimates gives

Numerical implementation

In this section we discuss briefly how to solve problem (1.9) by semidefinite programming. The dual problem of (1.9) takes the form

where $F_{\text{lo}}$ denotes the linear operator that maps a function to its first $n:=2f_{\text{lo}}+1$ Fourier coefficients as in (1.8) so that $Q_{\text{lo}}=F_{\text{lo}}^{*}F_{\text{lo}}$. The dual can be recast as the semidefinite program (SDP)

where $X$ is an $n\times n$ Hermitian matrix, leveraging a corollary to Theorem 4.24 in (see also ). In most cases, this allows to solve the primal problem with high accuracy. The following lemma suggests how to obtain a primal solution from a dual solution.

By Hölder’s inequality and the constraint on $F_{\text{lo}}^{\ast}u_{\text{est}}$, $\left|\left|x_{\text{est}}\right|\right|_{\text{TV}}\geq\left\langle x_{\text{est}},F_{\text{lo}}^{\ast}u_{\text{est}}\right\rangle$ so that equality holds. This is only possible if $F_{\text{lo}}^{\ast}u_{\text{est}}$ equals the sign of $x_{\text{est}}$ at every point where $x_{\text{est}}$ is nonzero.

This result implies that it is usually possible to determine the support of the primal solution by locating those points where the polynomial $q(t)=(F_{\text{lo}}^{\ast}u_{\text{est}})(t)$ has modulus equal to one. Once the support is estimated accurately, a solution to the primal problem can be found by solving a discrete problem. Figure 3 shows the result of applying this scheme to a simple example. We omit further details and defer the analysis of this approach to future work.

Discussion

In this work we introduce a theoretical framework that provides non-asymptotic stability guarantees for tractable super-resolution of multiple point sources in a continuous domain. More precisely, we show that it is possible to extrapolate the spectrum of a superposition of point sources by convex programming and that the extrapolation error scales quadratically with the super-resolution factor. This is a worst case analysis since the noise has bounded norm but is otherwise arbitrary. Natural extensions would include stability studies using other error metrics and noise models. For instance, an analysis tailored to a stochastic model might be able to sharpen Corollary 1.3 and be more precise in its findings. In a different direction, our techniques may be directly applicable to related problems. An example concerns the use of the total-variation norm for denoising line spectra . Here, it would be interesting to see whether our methods allow to prove better denoising performance under a minimum-separation condition. Another example concerns the recovery of sparse signals from a random subset of their low-pass Fourier coefficients . Here, it is likely that our work would yield stability guarantees from noisy low-frequency data.

E. C. is partially supported by AFOSR under grant FA9550-09-1-0643, by ONR under grant N00014-09-1-0258 and by a gift from the Broadcom Foundation. C. F. is supported by a Fundación Caja Madrid Fellowship. We thank Carlos Sing-Long for useful feedback about an earlier version of the manuscript.

References

Appendix A Proof of Lemma 2.5

We use the construction described in Section 2 of . In more detail,

Without loss of generality we consider $t_{j}=0$ and bound $q\left(t\right)-v_{j}$ in the interval $\left[0,0.16\lambda_{\text{lo}}\right]$. To ease notation, we define $w(t)=q\left(t\right)-v_{j}=w_{R}(t)+i\,w_{I}(t)$, where $w_{R}$ is the real part of $w$ and $w_{I}$ the imaginary part. Leveraging different results from Section 2 in (in particular the equations in (2.25) and Lemmas 2.2 and 2.7), we have

The same bound holds for $w_{I}$. Since $w_{R}(0)$, $w_{R}^{\prime}(0)$, $w_{I}(0)$ and $w_{I}^{\prime}(0)$ are all equal to zero, this implies $\left|w_{R}(t)\right|\leq C^{\prime}f_{\text{lo}}^{2}t^{2}$ and $\left|w_{I}(t)\right|\leq C^{\prime}f_{\text{lo}}^{2}t^{2}$ in the interval of interest, which allows the conclusion

Appendix B Proof of Lemma 2.7

The proof is similar to that of Lemma 2.4 (see Section 2 of ), where a low-frequency kernel and its derivative are used to interpolate an arbitrary sign pattern on a support satisfying the minimum-distance condition. More precisely, we set

In order to satisfy (2.18) and (2.19), we constrain $q_{1}$ as follows: for each $t_{j}\in T$,

Intuitively, this forces $q_{1}$ to approximate the linear function $v_{j}\left(t-t_{j}\right)$ around $t_{j}$. These constraints can be expressed in matrix form,

and $j$ and $k$ range from $1$ to $\left|T\right|$. It is shown in Section 2.3.1 of that under the minimum-separation condition this system is invertible, so that $\alpha$ and $\beta$ are well defined. These coefficient vectors can consequently be expressed as

where $S$ is the Schur complement. Inequality (B.2) implies

where $\kappa=\left|G^{\left(2\right)}(0)\right|=\pi^{2}f_{\text{lo}}(f_{\text{lo}}+4)/3$.

Let $\|M\|_{\infty}$ denote the usual infinity norm of a matrix $M$ defined as $\|M\|_{\infty}=\max_{\|x\|_{\infty}=1}\|Mx\|_{\infty}=\max_{i}\sum_{j}|a_{ij}|$. Then, if $\left|\left|I-M\right|\right|_{\infty}<1$, the series $M^{-1}=\left(I-\left(I-M\right)\right)^{-1}=\sum_{k\geq 0}\left(I-M\right)^{k}$ is convergent and we have

This, together with (B.4), (B.5) and (B.6) implies

for a certain positive constant $C_{\kappa}$. Note that due to the numeric upper bounds on the constants in (B.2) $C_{\kappa}$ is indeed a positive constant as long as $f_{\text{lo}}\geq 1$. Finally, we obtain a bound on the magnitude of the entries of $\alpha$

where $C_{\alpha}=C_{\kappa}C_{1}/\left(1-C_{0}\right)$, and on the entries of $\beta$

for a positive constant $C_{\beta}=C_{\kappa}$. Combining these inequalities with (B.3) and the fact that the absolute values of $G(t)$ and $G^{\left(1\right)}(t)$ are bounded by one and $7f_{\text{lo}}$ respectively (see the proof of Lemma C.5 in ), we have that for any $t$

where $t_{i}$ denotes the element in $T$ nearest to $t$ (note that all other elements are at least $\Delta/2$ away). Thus, (2.19) holds.

The proof is completed by the following lemma, which proves (2.18).

Proof We assume without loss of generality that $t_{j}=0$. By symmetry, it suffices to show the claim for $t\in(0,0.16\,\lambda_{\text{lo}}]$. To ease notation, we define $w(t)=v_{j}t-q_{1}(t)=w_{R}(t)+i\,w_{I}(t)$, where $w_{R}$ is the real part of $w$ and $w_{I}$ the imaginary part. Leveraging (B.7), (B.8) and (B.2) together with the fact that $G^{\left(2\right)}(t)$ and $G^{\left(3\right)}(t)$ are bounded by $4f_{\text{lo}}^{2}$ and $6f_{\text{lo}}^{3}$ respectively if $\left|t\right|\leq 0.16\lambda_{\text{lo}}$ (see the proof of Lemma 2.3 in ), we obtain

The same bound applies to $w_{I}$. Since $w_{R}(0)$, $w_{R}^{\prime}(0)$, $w_{I}(0)$ and $w_{I}^{\prime}(0)$ are all equal to zero, this implies $\left|w_{R}(t)\right|\leq Cf_{\text{lo}}t^{2}$—and similarly for $\left|w_{I}(t)\right|$—in the interval of interest. Whence, $\left|w(t)\right|\leq Cf_{\text{lo}}t^{2}$.

Appendix C Proof of Corollary 1.3

The proof of Theorem 1.2 relies on two identities

To prove the corollary, we show that (C.1) and (C.2) hold. Due to the fact that $\left|\left|\epsilon\right|\right|_{2}^{2}$ follows a $\chi^{2}$-distribution with $4f_{\text{lo}}+2$ degrees of freedom, we have

for any positive $\gamma$ by a concentration inequality (see [24, Section 4]). By Parseval, this implies that with high probability $\left|\left|Q_{\text{lo}}x-y\right|\right|_{\mathcal{L}_{2}}=\left|\left|\epsilon\right|\right|_{2}\leq\delta$. As a result, $x_{\text{est}}$ is feasible, which implies (C.1) and furthermore

since by the Cauchy-Schwarz inequality $\left|\left|f\right|\right|_{L_{1}}\leq\left|\left|f\right|\right|_{\mathcal{L}_{2}}$ for any function $f$ with bounded $L_{2}$ norm supported on the unit interval. Thus, (C.2) also holds and the proof is complete.

Appendix D Extension to multiple dimensions

The extension of the proof hinges on establishing versions of Lemmas 2.4, 2.5 and 2.7 for multiple dimensions. These lemmas construct bounded low-frequency polynomials which interpolate a sign pattern on a well-separated set of points $S$ and have bounded second derivatives in a neighborhood of $S$. In the multidimensional case, we need the directional derivative of the polynomials to be bounded in any direction, which can be ensured by bounding the eigenvalues of their Hessian matrix evaluated on the support of the signal. To construct such polynomials one can proceed in a way similar to the proof of Lemmas 2.4 and 2.7, namely, by using a low-frequency kernel constructed by tensorizing several squared Fejér kernels to interpolate the sign pattern, while constraining the first-order derivatives to either vanish or have a fixed value. As in the one-dimensional case, one can set up a system of equations and prove that it is well conditioned using the rapid decay of the interpolation kernel away from the origin. Finally, one can verify that the construction satisfies the required conditions by exploiting the fact that the interpolation kernel and its derivatives are locally quadratic and rapidly decaying. This is spelled out in the proof of Proposition C.1 in to prove a version of Lemma 2.4 in two dimensions. In order to clarify further how to adapt our techniques to a multidimensional setting we provide below a sketch of the proof of the analog of Lemma 2.1 in two dimensions. In particular, this illustrates how the increase in dimension does not change the exponent of the SRF in our recovery guarantees.

The proof relies on the existence of a low-frequency polynomial

As in one dimension, we perform a polar decomposition of $P_{T}h$,

and work with $v_{j}=e^{-i\phi(t_{j})}$. The rest of the proof is almost identical to the 1D case. Since $q$ is low frequency,

Next, since $q$ interpolates $e^{-i\phi\left(t\right)}$ on $T$,

Applying (D.3) and Hölder’s inequality, we obtain

Setting $t_{j}=(0,0)$ without loss of generality, the triangle inequality and (D.2) yield

By the same argument as in the 1D case, the fact that $\hat{x}$ has minimal total-variation norm is now sufficient to establish