Towards a Mathematical Theory of Super-Resolution

Introduction

Super-resolution is a word used in different contexts mainly to design techniques for enhancing the resolution of a sensing system. Interest in such techniques comes from the fact that there usually is a physical limit on the highest possible resolution a sensing system can achieve. To be concrete, the spatial resolution of an imaging device may be measured by how closely lines can be resolved. For an optical system, it is well known that resolution is fundamentally limited by diffraction. In microscopy, this is called the Abbe diffraction limit and is a fundamental obstacle to observing sub-wavelength structures. This is the reason why resolving sub-wavelength features is a crucial challenge in fields such as astronomy , medical imaging , and microscopy . In electronic imaging, limitations stem from the lens and the size of the sensors, e. g. pixel size. Here, there is an inflexible limit to the effective resolution of a whole system due to photon shot noise which degrades image quality when pixels are made smaller. Some other fields where it is desirable to extrapolate fine scale details from low-resolution data—or resolve sub-pixel details—include spectroscopy , radar , non-optical medical imaging and geophysics . For a survey of super-resolution techniques in imaging, see and the references therein.

This paper is about super-resolution, which loosely speaking is the process whereby the fine scale structure of an object is retrieved from coarse scale information only.We discuss here computational super-resolution methods as opposed to instrumental techniques such as interferometry. A useful mathematical model may be of the following form: start with an object $x(t_{1},t_{2})$ of interest, a function of two spatial variables, and the point-spread function $h(t_{1},t_{2})$ of an optical instrument. This instrument acts as a filter in the sense that we may observe samples from the convolution product

In the frequency domain, this equation becomes

where $\hat{x}$ is the Fourier transform of $x$, and $\hat{h}$ is the modulation transfer function or simply transfer function of the instrument. Now, common optical instruments act as low-pass filters in the sense that their transfer function $\hat{h}$ vanishes for all values of $\omega$ obeying $|\omega|\geq\Omega$ in which $\Omega$ is a frequency cut-off; that is,

In microscopy with coherent illumination, the bandwidth $\Omega$ is given by $\Omega=2\pi\text{NA}/\lambda$ where NA is the numerical aperture and $\lambda$ is the wavelength of the illumination light. For reference, the transfer function in this case is simply the indicator function of a disk and the point-spread function has spherical symmetry and a radial profile proportional to the ratio between the Bessel function of the first order and the radius. In microscopy with incoherent light, the transfer function is the Airy function and is proportional to the square of the coherent point-spread function. Regardless, the frequency cut-off induces a physical resolution limit which is roughly inversely proportional to $\Omega$ (in microscopy, the Rayleigh resolution distance is defined to be $0.61\times 2\pi/\Omega$).

The daunting and ill-posed super-resolution problem then consists in recovering the fine-scale or, equivalently, the high-frequency features of $x$ even though they have been killed by the measurement process. This is schematically represented in Figure 1, which shows a highly resolved signal together with a low-resolution of the same signal obtained by convolution with a point-spread function. Super-resolution aims at recovering the fine scale structure on the left from coarse scale features on the right. Viewed in the frequency domain, super-resolution is of course the problem of extrapolating the high-end and missing part of the spectrum from the low-end part, as seen in Figure 2. For reference, this is very different from a typical compressed sensing problem in which we wish to interpolate—and not extrapolate—the spectrum.

2 Models and methods

For concreteness, consider a continuous-time model in which the signal of interest is a weighted superposition of spikes

where $\{t_{j}\}$ are locations in $$and$\delta_{\tau}$is a Dirac measure at$\tau$. The amplitudes$a_{j}$may be complex valued. Expressed differently, the signal$x$is an atomic measure on the unit interval putting complex mass at time points$t_{1}$,$t_{2}$,$t_{3}$and so on. The information we have available about$x$is a sample of the lower end of its spectrum in the form of the lowest$2f_{c}+1$Fourier series coefficients ($f_{c}$ is an integer):

For simplicity, we shall use matrix notations to relate the data $y$ and the object $x$ and will write (1.3) as $y=\mathcal{F}_{n}\,x$ where $\mathcal{F}_{n}$ is the linear map collecting the lowest $n=2f_{c}+1$ frequency coefficients. It is important to bear in mind that we have chosen this model mainly for ease of exposition. Our techniques can be adapted to settings where the measurements are modeled differently, e. g. by sampling the convolution of the signal with different low-pass kernels. The important element is that just as before, the frequency cut-off induces a resolution limit inversely proportional to $f_{c}$; below we set $\lambda_{c}=1/f_{c}$ for convenience.

Let $T=\{t_{j}\}$ be the support of $x$. If the minimum distance obeys

then $x$ is the unique solution to (1.4). This holds with the proviso that $f_{c}\geq 128$. If $x$ is known to be real-valued, then the minimum gap can be lowered to $1.87\,\lambda_{c}$.

We find this result particularly unexpected. The total-variation norm makes no real assumption about the structure of the signal. Yet, not knowing that there are any spikes, let alone how many there are, total-variation minimization locates the position of those spikes with infinite precision! Even if we knew that (1.4) returned a spike train, there is no reason to expect that the locations of the spikes would be infinitely accurate from coarse scale information only. In fact, one would probably expect the fitted locations to deviate at least a little from the truth. This is not what happens.

An interesting aspect of this theorem is that it cannot really be tested numerically. Indeed, one would need a numerical solver with infinite precision to check that total-variation minimization truly puts the spikes at exactly the right locations. Although this is of course not available, Section 4 shows how to solve the minimum total-variation problem (1.4) by using ideas from semidefinite programming and allowing recovery of the support with very high precision. Also, we demonstrate through numerical simulations in Section 5 that Theorem 1.2 is fairly tight in the sense that a necessary condition is a separation of at least $\lambda_{c}=1/f_{c}$.

Viewed differently, one can ask: how many spikes can be recovered from $n=2f_{c}+1$ low-frequency samples? The answer given by Theorem 1.2 is simple. At least $n/4$ provided we have the minimum separation discussed above. A classical argument shows that any method whatsoever would at least need two samples per unknown spike so that the number of spikes cannot exceed half the number of samples, i. e. $n/2$. This is another way of showing that the theorem is reasonably tight.

3 Super-resolution in higher dimensions

Our results extend to higher dimensions and reveal the same dependence between the minimum separation and the measurement resolution as in one dimension. For concreteness, we discuss the $2$-dimensional setting and emphasize that the situation in $d$ dimensions is similar. Here, we have a measure

as before but in which the $t_{j}\in^{2}$. We are given information about $x$ in the form of low-frequency samples of the form

This again introduces a physical resolution of about $\lambda_{c}=1/f_{c}$. In this context, we may think of our problem as imaging point sources in the 2D plane—such as idealized stars in the sky—with an optical device with resolution about $\lambda_{c}$—such as a diffraction limited telescope. Our next result states that it is possible to locate the point sources without any error whatsoever if they are separated by a distance of $2.38\,\lambda_{c}$ simply by minimizing the total variation.

Then if $x$ is real valued, it is the unique minimum total-variation solution among all real objects obeying the data constraints (1.7). Hence, the recovery is exact. For complex measures, the same statement holds but with a slightly different constant.

Whereas we have tried to optimize the constant in one dimensionThis is despite the fact that the authors have a proof—not presented here—of a version of Theorem 1.2 with a minimum separation at least equal to $1.85\,\lambda_{c}$., we have not really attempted to do so here in order to keep the proof reasonably short and simple. Hence, this theorem is subject to improvement.

4 Discrete super-resolution

the connection with the previous sections is obvious since $x$ might be interpreted as samples of a discrete signal on a grid $\{t/N\}$ with $t=0,1,\ldots,N-1$. In fact, the continuous-time setting is the limit of infinite resolution in which $N$ tends to infinity while the number of samples remains constant ($f_{c}$ fixed). Instead, we can choose to study the regime in which the ratio between the actual resolution of the signal $1/N$ and the resolution of the data defined as $1/f_{c}$ is constant. This gives the corollary below.

Let $T\subset\{0,1,\ldots,N-1\}$ be the support of $\{x_{t}\}_{t=0}^{N-1}$ obeying

in which $F_{n}$ is the partial Fourier matrix in (1.9) is exact.

5 The super-resolution factor

In the discrete framework, we wish to resolve a signal on a fine grid with spacing $1/N$. However, we only observe the lowest $n=2f_{c}+1$ Fourier coefficients so that in principle, one can only hope to recover the signal on a coarser grid with spacing only $1/n$ as shown in Figure 4. Hence, the factor $N/n$, or equivalently, the ratio between the spacings in the coarse and fine grids, can be interpreted as a super-resolution factor (SRF). Below, we set

The reason for introducing the SRF is that with inexact data, we obviously cannot hope for infinite resolution. Indeed, noise will ultimately limit the resolution one can ever hope to achieve and, therefore, the question of interest is to study the accuracy one might expect from a practical super-resolution procedure as a function of both the noise level and the SRF.

6 Stability

for some $\delta\geq 0$, where $F_{n}$ is as before. Letting $P_{n}$ be the orthogonal projection of a signal onto the first $n$ Fourier modes, $P_{n}=\frac{1}{N}F_{n}^{*}F_{n}$, we can view (1.13) as an input noise model since with $w=F_{n}z$, we have

since the high-frequency part of $z$ is filtered out by the measurement process. Finally, with $s={N}^{-1}F_{n}^{*}y$, (1.13) is equivalent to

We propose studying the relaxed version of the noiseless problem (1.11)

We show that this recovers $x$ with a precision inversely proportional to $\delta$ and to the square of the super-resolution factor.

Assume that $x$ obeys the separation condition (1.6). Then with the noise model (1.14), the solution $\hat{x}$ to (1.15) obeys

This theorem, which shows the simple dependence upon the super-resolution factor and the noise level, is proved in Section 3. Clearly, plugging in $\delta=0$ in (1.16) gives Corollary 1.4.

Versions of Theorem 1.5 hold in the continuous setting as well, where the locations of the spikes are not assumed to lie on a given fine grid but can take on a continuum of values. The arguments are more involved than those needed to establish (1.16) and we leave a detailed study to a future paper.

7 Sparsity and stability

Researchers in the field know that super-resolution under sparsity constraints alone is hopelessly ill posed. In fact, without a minimum distance condition, the support of sparse signals can be very clustered, and clustered signals can be nearly completely annihilated by the low-pass sensing mechanism. The extreme ill-posedness can be understood by means of the seminal work of Slepian on discrete prolate spheroidal sequences. This is surveyed in Section 3.2 but we give here a concrete example to drive this point home.

To keep things simple, we consider the ‘analog version’ of (1.14) in which we observe

we essentially have $\mathcal{P}_{W}=P_{n}$ where the equality is true in the limit where $N\rightarrow\infty$ (technically, $\mathcal{P}_{W}$ is the convolution with the sinc kernel while $P_{n}$ uses the Dirichlet kernel). Set a mild level of super-resolution to fix ideas,

Now the work of Slepian shows that there is a $k$-sparse signal supported on $[0,\ldots,k-1]$ obeying

Even knowing the support ahead of time, how are we going to recover such signals from noisy measurements? Even for a very mild super-resolution factor of just $\text{SRF}=1.05$ (we only seek to extend the spectrum by 5%), (1.17) becomes

which implies that there exists a unit-norm signal with at most $256$ consecutive nonzero entries such that $\left|\left|\mathcal{P}_{W}x\right|\right|_{2}\leq 1.2\times 10^{-15}$. Of course, as the super-resolution factor increases, the ill-posedness gets worse. For large values of SRF, there is $x$ obeying (1.17) with

8 Comparison with related work

Finally, we would like to mention an alternative approach to the super-resolution of pointwise events from coarse scale data. Leveraging ideas related to error correction codes and spectral estimation, shows that it is possible to recover trains of Dirac distributions from low-pass measurements at their rate of innovation (in essence, the density of spikes per unit of time). This problem, however, is extraordinarily ill posed without a minimum separation assumption as explained in Sections 1.7 and 3.2. Moreover, the proposed reconstruction algorithm in needs to know the number of events ahead of time, and relies on polynomial root finding. As a result, it is highly unstable in the presence of noise as discussed in , and in the presence of approximate sparsity. Algebraic techniques have also been applied to the location of singularities in the reconstruction of piecewise polynomial functions from a finite number of Fourier coefficients (see and references therein). The theoretical analysis of these methods proves their accuracy up to a certain limit related to the number of measurements. Corollary 1.6 takes a different approach, guaranteeing perfect localization if there is a minimum separation between the singularities.

9 Connections to sparse recovery literature

Theorem 1.2 and Corollary 1.4 can be interpreted in the framework of sparse signal recovery. For instance, by swapping time and frequency, Corollary 1.4 asserts that one can recover a sparse superposition of tones with arbitrary frequencies from $n$ time samples of the form

where the frequencies are of the form $\omega_{j}=j/N$. Since the spacing between consecutive frequencies is not $1/n$ but $1/N$, we may have a massively oversampled discrete Fourier transform, where the oversampling ratio is equal to the super-resolution factor. In this context, a sufficient condition for perfectly super-resolving these tones is a minimum separation of $4/n$. In addition, Theorem 1.2 extends this to continuum dictionaries where tones $\omega_{j}$ can take on arbitrary real values.

The matrix with normalized columns $f_{j}=\{e^{-i2\pi t\omega_{j}}/\sqrt{n}\}_{t=0}^{n-1}$ does not obey the restricted isometry property since a submatrix composed of a very small number of contiguous columns is already very close to singular, see and Section 3.2 for related claims. For example, with $N=512$ and a modest SRF equal to 4, the smallest singular value of submatrices formed by eight consecutive columns is $3.32\;10^{-5}$.

If $n<N/2$, i.e. $\text{SRF}>2$, this says that $|T|$ must be zero. In other words, to recover one spike, we would need at least half of the Fourier samples.

Other guarantees based on the coherence of the dictionary yield similar results. A popular condition requires that

where $M$ is the coherence of the system defined as $\max_{i\neq j}|\langle f_{i},f_{j}\rangle|$. When $N=1024$ and $\text{SRF}=4$, $M\approx 0.9003$ so that this becomes $|T|\leq 1.055$, and we can only hope to recover one spike.

The condition $\text{WERC}\left(T\right)<1$ guarantees exact recovery. Considering three spikes and using Taylor expansions to bound the sine function, the minimum distance needed to ensure that $\text{WERC}\left(T\right)<1$ may be lower bounded by $24\text{SRF}^{3}/\pi^{3}-2\text{SRF}$. This is achieved by considering three spikes at $\omega\in\{0,\pm\Delta\}$, where $\Delta=(k+1/2)/n$ for some integer $k$; we omit the details. If $N=20,000$ and the number of measurements is $1,000$, this allows for the recovery of at most $3$ spikes, whereas Corollary 1.4 implies that it is possible to reconstruct at least $n/4=250$. Furthermore, the cubic dependence on the super-resolution factor means that if we fix the number of measurements and let $N\rightarrow\infty$, which is equivalent to the continuous setting of Theorem 1.2, the separation needed becomes infinite and we cannot guarantee the recovery of even two spikes.

10 Extensions

Standard Fourier analysis gives that the $k$th Fourier coefficient of this measure is given by

If $T=\{t_{j}\}$ obeys (1.6), $x$ is determined exactly from $y$ by solving (1.22).

Extensions to non-periodic functions, other types of discontinuities and smoothness assumptions are straightforward.

11 Organization of the paper

The remainder of the paper is organized as follows. We prove our main noiseless result in Section 2. There, we introduce our techniques which involve the construction of an interpolating low-frequency polynomial. Section 3 proves our stability result and argues that sparsity constraints cannot be sufficient to guarantee stable super-resolution. Section 4 shows that (1.4) can be cast as a finite semidefinite program. Numerical simulations providing a lower bound for the minimum distance that guarantees exact recovery are presented in Section 5. We conclude the paper with a short discussion in Section 6.

Noiseless Recovery

This result follows from elementary measure theory and is included in Section A of the Appendix for completeness. Constructing a bounded low-frequency polynomial interpolating the sign pattern of certain signals becomes increasingly difficult if the minimum distance separating the spikes is too small. This is illustrated in Figure 5, where we show that if spikes are very near, it would become in general impossible to find an interpolating low-frequency polynomial obeying (2.2).

2 Proof of Theorem 1.2

Theorem 1.2 is a direct consequence of the proposition below, which establishes the existence of a valid dual polynomial provided the elements in the support are sufficiently spaced.

The remainder of this section proves this proposition. Our method consists in interpolating $v$ on $T$ with a low-frequency kernel and correcting the interpolation to ensure that the derivative of the dual polynomial is zero on $T$. The kernel we employ is

and $K(0)=1$. If $f_{c}$ is even, $K(t)$ is the square of the Fejér kernel which is a trigonometric polynomial with frequencies obeying $|k|\leq f_{c}/2$. As a consequence, $K$ is of the form (2.1). The careful reader might remark that the choice of the interpolation kernel seems somewhat arbitrary. In fact, one could also use the Fejér kernel or any other power of the Fejér kernel using almost identical proof techniques. We have found that the second power nicely balances the trade-off between localization in time and in frequency, and thus yields a good constant.

To construct the dual polynomial, we interpolate $v$ with both $K$ and its derivative $K^{\prime}$,

whereas in order to obey $|q(t)|<1$ for $t\in T^{c}$, we impose $q^{\prime}(t_{k})=0$,

As we will see, this implies that the magnitude of $q$ reaches a local maximum at those points, which in turn can be used to show that (2.2) holds.

The proof of Proposition 2.1 consists of three lemmas, which are the object of the following section. The first one establishes that if the support is spread out, it is possible to interpolate any sign pattern exactly.

Under the hypotheses of Proposition 2.1, there exist coefficient vectors $\alpha$ and $\beta$ obeying

such that (2.5)–(2.6) hold. Further, if $v_{1}=1$,

To complete the proof, Lemmas 2.3 and 2.4 show that $|q\left(t\right)|<1$.

Fix $\tau\in T$. Under the hypotheses of Proposition 2.1, $|q(t)|<1$ for $\left|t-\tau\right|\in(0,0.1649\,\lambda_{c}]$.

Fix $\tau\in T$. Then under the hypotheses of Proposition 2.1, $|q(t)|<1$ for $\left|t-\tau\right|\in[0.1649\,\lambda_{c},\Delta/2]$. This can be extended as follows: letting $\tau_{+}$ be the closest spike to the right, i. e. $\tau_{+}=\min\{t\in T:t>\tau\}$. Then $\left|q(t)\right|<1$ for all $t$ obeying $0<t-\tau\leq(\tau_{+}-\tau)/2$, and likewise for the left side.

Finally, we record a useful lemma to derive stability results.

If $\Delta\left(T\right)\geq 2.5\,\lambda_{c}$, then for any $\tau\in T$,

Further, for $\min_{\tau\in T}\,\left|t-\tau\right|>0.1649\,\lambda_{c}$, $\left|q\left(t\right)\right|$ is upper bounded by the right-hand side above evaluated at $0.1649\,\lambda_{c}$.

Section 2.5 describes how the proof can be adapted to obtain a slightly smaller bound on the minimum distance for real-valued signals.

3 Proofs of Lemmas

The proofs of the three lemmas above make repeated use of the fact that the interpolation kernel and its derivatives decay rapidly away from the origin. The intermediate result below proved in Section B of the Appendix quantifies this.

where $H_{0}^{\infty}=1$, $H_{1}^{\infty}=4$, $H_{2}^{\infty}=18$, $H_{3}^{\infty}=77$,

This lemma is used to control quantities of the form $\sum_{t_{i}\in T\setminus\{\tau\}}\left|K\left(t-t_{i}\right)\right|$ ($\tau\in T$) as shown below.

Suppose $0\in T$. Then for all $t\in[0,\Delta/2]$,

Proof We consider the sum over positive $t_{i}\in T$ first and denote by $t_{+}$ the positive element in $T$ closest to . We have

Let us assume $t_{+}<2\Delta_{\text{min}}$ (if $t_{+}>2\Delta_{\text{min}}$ the argument is very similar). Note that the assumption that $f_{c}\geq 128$ implies $21\Delta_{\text{min}}<0.33<\sqrt{2}/\pi$. By Lemma 2.6 and the minimum separation condition, this means that the second term in the right-hand side is at most

which can be upper bounded using the fact that

the first inequality holds because $t<\Delta_{\text{min}}$ and the last because the Riemann zeta function is equal to $\pi^{4}/90$ at 4. Also,

To verify the claim about the monotonicity w.r.t. $\Delta$, observe that both terms

Now set $t^{\prime}>t$. Then by Lemma 2.6,

Finally, a last fact we shall use is that $K\left(0\right)=1$ is the global maximum of $K$ and $\left|K^{\prime\prime}\left(0\right)\right|=\left|-\pi^{2}f_{c}\left(f_{c}+4\right)/3\right|$ the global maximum of $\left|K^{\prime\prime}\right|$.

where $j$ and $k$ range from $1$ to $\left|T\right|$. With this, (2.5) and (2.6) become

A standard linear algebra result asserts that this system is invertible if and only if $D_{2}$ and its Schur complement $D_{0}-D_{1}D_{2}^{-1}D_{1}$ are both invertible. To prove that this is the case we can use the fact that a symmetric matrix $M$ is invertible if

where $\|A\|_{\infty}$ is the usual infinity norm of a matrix defined as $\|A\|_{\infty}=\max_{\|x\|_{\infty}=1}\|Ax\|_{\infty}=\max_{i}\sum_{j}|a_{ij}|$. This follows from $M^{-1}=(I-H)^{-1}=\sum_{k\geq 0}H^{k}$, $H=I-M$, where the series is convergent since $\left|\left|H\right|\right|_{\infty}<1$. In particular,

We also make use of the inequalities below, which follow from Lemma 2.7

Note that $D_{2}$ is symmetric because the second derivative of the interpolation kernel is symmetric. The bound (2.16) and the identity $K^{\prime\prime}\left(0\right)=-\pi^{2}f_{c}\left(f_{c}+4\right)/3$ give

which implies the invertibility of $D_{2}$. The bound (2.13) then gives

Combining this with (2.14) and (2.15) yields

Note that the Schur complement of $D_{2}$ is symmetric because $D_{0}$ and $D_{2}$ are both symmetric whereas $D_{1}^{T}=-D_{1}$ since the derivative of the interpolation kernel is odd. This shows that the Schur complement of $D_{2}$ is invertible and, therefore, the coefficient vectors $\alpha$ and $\beta$ are well defined.

There just remains to bound the interpolation coefficients, which can be expressed as

where $C$ is the Schur complement. The relationships (2.13) and (2.18) immediately give a bound on the magnitude of the entries of $\alpha$

Similarly, (2.15), (2.17) and (2.18) allow to bound the entries of $\beta$:

Finally, with $v_{1}=1$, we can use (2.18) to show that $\alpha_{1}$ is almost equal to 1. Indeed,

$\left|\gamma_{1}\right|\leq\left|\left|\text{\em I}-C^{-1}\right|\right|_{\infty}$, and

3.2 Proof of Lemma 2.3

We assume without loss of generality that $\tau=0$ and $q(0)=1$. By symmetry, it suffices to show the claim for $t\in(0,0.1649\,\lambda_{c}]$. Since $q^{\prime}(0)=0$, local strict concavity would imply that $\left|q\left(t\right)\right|<1$ near the origin. We begin by showing that the second derivative of $\left|q\right|$ is strictly negative in the interval $\left(0,0.1649\,\lambda_{c}\right)$. This derivative is equal to

where $q_{R}$ is the real part of $q$ and $q_{I}$ the imaginary part. As a result, it is sufficient to show that

as long as $\left|q\left(t\right)\right|$ is bounded away from zero. In order to bound the different terms in (2.19), we use the series expansions of the interpolation kernel and its derivatives around the origin to obtain the inequalities, which hold for all $t\in[-1/2,1/2]$,

The lower bounds are decreasing in $t$, while the upper bounds are increasing in $t$, so we can evaluate them at $0.1649\,\lambda_{c}$ to establish that for all $t\in[0,0.1649\,\lambda_{c}]$,

We combine this with Lemmas 2.7 and 2.2 to control the different terms in (2.19) and begin with $q_{R}\left(t\right)$. Here,

These bounds allow us to conclude that $|q|^{\prime\prime}$ is negative on $[0,0.1649\lambda_{c}]$ since

3.3 Proof of Lemma 2.4

As before, we assume without loss of generality that $\tau=0$ and $q(0)=1$. We use Lemma 2.7 again to bound the absolute value of the dual polynomial on $[0.1649\lambda_{c},\Delta/2]$ and write

Note that we are assuming adversarial sign patterns and as a result we are unable to exploit cancellations in the coefficient vectors $\alpha$ and $\beta$. To control $|K(t)|$ and $|K^{\prime}(t)|$ between $0.1649\,\lambda_{c}$ and $0.7559\,\lambda_{c}$, we use series expansions around the origin which give

This derivative is strictly negative between $0.1649\,\lambda_{c}$ and $0.7559\,\lambda_{c}$, which implies that $L_{1}\left(t\right)$ is decreasing in this interval. Put

By Lemma 2.7, this function is increasing. With (2.28), this gives the crude bound

Table 2 shows that taking $\left\{t_{1},t_{2}\right\}=\left\{0.1649\,\lambda_{c},0.4269\,\lambda_{c}\right\}$ and then $\left\{t_{1},t_{2}\right\}=\left\{0.4269\,\lambda_{c},0.7559\,\lambda_{c}\right\}$ proves that $\left|q(t)\right|<1$ on $[0.1649\lambda_{c},0.7559\lambda_{c}]$. For $0.7559\lambda_{c}\leq t\leq\Delta/2$, we apply Lemma 2.6 and obtain

here, the second step follows from the monotonicity of $B_{0}$ and $B_{1}$. Finally, for $\Delta/2\leq t\leq t_{+}/2$, this last inequality applies as well. This completes the proof.

4 Proof of Lemma 2.5

Replacing $\Delta=1.98\,\lambda_{c}$ by $\Delta=2.5\,\lambda_{c}$ and going through exactly the same calculations as in Sections 2.3.1 and 2.3.2 yields that for any $t$ obeying $0\leq\left|t-\tau\right|\leq 0.1649\,\lambda_{c}$,

At a distance of $0.1649\,\lambda_{c}$, the right-hand side is equal to $0.9909$. The calculations in Section 2.3.3 with $\Delta=2.5\,\lambda_{c}$ imply that the magnitude of $q(t)$ at locations at least $0.1649\,\lambda_{c}$ away from an element of $T$ is bounded by $0.9843$. This concludes the proof.

5 Improvement for real-valued signals

Stability

This section proves Theorem 1.5 and we begin by establishing a strong form of the null-space property. In the remainder of the paper $P_{T}$ is the orthogonal projector onto the linear space of vectors supported on $T$, namely, $(P_{T}x)_{i}=x_{i}$ if $i\in T$ and is zero otherwise.

Under the assumptions of Theorem 1.5, any vector $h$ such that $F_{n}h=0$ obeys

for some numerical constant $\rho$ obeying $0<\rho<1$. This constant is of the form $1-\rho=\alpha/\text{\em SRF}^{2}$ for some positive $\alpha>0$. If $\text{SRF}\geq 3.03$, we can take $\alpha=0.0883$.

Proof Let $P_{T}h_{t}=\left|P_{T}h_{t}\right|e^{i\phi_{t}}$ be the polar decomposition of $P_{T}h$, and consider the low-frequency polynomial $q(t)$ in Proposition 2.1 interpolating $v_{t}=e^{-i\phi_{t}}$. We shall abuse notations and set $q=\{q_{t}\}_{t=0}^{N-1}$ where $q_{t}=q(t/N)$. For $t\notin T$, $|q(t/N)|=|q_{t}|\leq\rho<1$. By construction $q=P_{n}q$, and thus $\langle q,h\rangle=\langle q,P_{n}h\rangle=0$. Also,

For the numerical constant, we use Lemma 2.5 which says that if $0\in T$ and $1/N\leq 0.1649\lambda_{c}$, which is about the same as $1/\text{SRF}\approx 2f_{c}/N\leq 2\times 0.1649$ or $\text{SRF}>3.03$, we have

This applies directly to any other $t$ such that $\min_{\tau\in T}\left|t-\tau\right|=1$. Also, for all $t$ at distance at least 2 from $T$, Lemma 2.5 implies that $|q(t/N)|\leq\rho$. This completes the proof.

The proof is a fairly simple consequence of Lemma 3.1. Set $h=\hat{x}-x$ and decompose the error into its low- and high-pass components

The high-frequency part is in the null space of $P_{n}$ and (3.1) gives

where the last inequality follows from (3.2). Hence,

where the last inequality follows from (3.3).

Since from Lemma 3.1, we have $1-\rho=\alpha/\text{SRF}^{2}$ for some numerical constant $\alpha$, the upper bound is of the form $4\alpha^{-1}\,\text{SRF}^{2}\,\delta$. For $\Delta(T)\geq 2.5\lambda_{c}$, we have $\alpha^{-1}\approx 11.235$.

2 Sparsity is not enough

For $\text{SRF}=16$ this is true of a subspace of dimension 36, two thirds of the total dimension. Such signals can be completely canceled out by perturbations of norm $5.02\;10^{-8}$, so that even at signal-to-noise ratios (SNR) of more than 145 dB, recovery is impossible by any method whatsoever.

Interestingly, the sharp transition shown in Figure 7 between the first singular values almost equal to one and the others, which rapidly decay to zero, can be characterized asymptotically by using the work of Slepian on prolate spheroidal sequences . Introduce the operator $\mathcal{T}_{k}$, which sets the value of an infinite sequence to zero on the complement of an interval $T$ of length $k$. With the notation of Section 1.7, the eigenvectors of the operator $\mathcal{P}_{W}\mathcal{T}_{k}$ are the discrete prolate spheroidal sequences $\{s_{j}\}_{j=1}^{k}$ introduced in ,

Set $v_{j}=\mathcal{T}_{k}s_{j}/\sqrt{\lambda_{j}}$, then by (3.4), it is not hard to see that

Therefore, for a fixed value of $\text{SRF}=1/2W$, and $k\geq 20$, the small eigenvalues are equal to zero for all practical purposes. In particular, for $\text{SRF}=4$ and $\text{SRF}=1.05$ we obtain (1.17) and (1.19) in Section 1.7 respectively. Additionally, a Taylor series expansion of $\gamma$ for large values of SRF yields (1.20).

Since $\|\mathcal{P}_{W}v_{j}\|_{L_{2}}=\sqrt{\lambda_{j}}$, the bound on $\lambda_{j}$ for $j$ near $k$ directly implies that some sparse signals are essentially zeroed out, even for small super-resolution factors. However, Figure 7 suggests an even stronger statement: as the super-resolution factor increases not only some, but most signals supported on $T$ seem to be almost completely suppressed by the low pass filtering. Slepian provides an asymptotic characterization for this phenomenon. Indeed, just about the first $2kW$ eigenvalues of $\mathcal{P}_{W}\mathcal{T}_{k}$ cluster near one, whereas the rest decay abruptly towards zero. To be concrete, for any $\epsilon>0$ and $j\geq 2kW\left(1+\epsilon\right)$, there exist positive constants $C_{0}$, and $\gamma_{0}$ (depending on $\epsilon$ and $W$) such that

This holds for all $k\geq k_{0}$, where $k_{0}$ is some fixed integer. This implies that for any interval $T$ of length $k$, there exists a subspace of signals supported on $T$ with dimension asymptotically equal to $\left(1-1/\text{SRF}\right)k$, which is obliterated by the measurement process. This has two interesting consequences. First, even if the super-resolution factor is just barely above one, asymptotically there will always exist an irretrievable vector supported on $T$. Second, if the super-resolution factor is two or more, most of the information encoded in clustered sparse signals is lost. Consider for instance a random sparse vector $x$ supported on $T$ with i.i.d. entries. Its projection onto a fixed subspace of dimension about $\left(1-1/\text{SRF}\right)k$ (corresponding to the negligible eigenvalues) contains most of the energy of the signal with high probability. However, this component is practically destroyed by low-pass filtering. Hence, super-resolving almost any tightly clustered sparse signal in the presence of noise is hopeless. This justifies the need for a minimum separation between nonzero components.

Minimization via semidefinite programming

At first sight, finding the solution to the total-variation norm problem (1.4) might seem quite challenging, as it requires solving an optimization problem over an infinite dimensional space. It is of course possible to approximate the solution by discretizing the support of the signal, but this could lead to an increase in complexity if the discretization step is reduced to improve precision. Another possibility is to try approximating the solution by estimating the support of the signal in an iterative fashion . Here, we take a different route and show that (1.4) can be cast as a semidefinite program with just $\left(n+1\right)^{2}/2$ variables, and that highly accurate solutions can be found rather easily. This formulation is similar to that in which concerns a related infinite dimensional convex program. Our exposition is less formal here than in the rest of the paper.

the constraint says that the trigonometric polynomial $(\mathcal{F}_{n}^{\ast}\,c)(t)=\sum_{|k|\leq f_{c}}c_{k}e^{i2\pi kt}$ has a modulus uniformly bounded by $1$ over the interval $$. The interior of the feasible set contains the origin and is consequently non empty, so that strong duality holds by a generalized Slater condition . The cost function involves a finite vector of dimension$n$, but the problem is still infinite dimensional due to the constraints. A corollary to Theorem 4.24 in allows to express this constraint as the intersection between the cone of positive semidefinite matrices$\{X:X\succeq 0\}$ and an affine hyperplane.

Returning to (4.1), the polynomial $e^{i2\pi f_{c}t}\,(\mathcal{F}_{n}^{\ast}\,c)(t)$ is causal and has the same magnitude as $(\mathcal{F}_{n}^{\ast}\,c)(t)$. Hence, the dual problem is equivalent to

The careful reader will observe that we have just shown how to compute the optimal value of (1.4), but not how we could obtain a solution. To find a primal solution, we abuse notation by letting $c$ be the solution to (4.3) and consider the trigonometric polynomial

which implies that the trigonometric polynomial $\mathcal{F}_{n}^{\ast}c$ is exactly equal to the sign of $\hat{x}$ when $\hat{x}$ is not vanishing. This is illustrated in Figure 8. Thus, to recover the support of the solution to the primal problem, we must simply locate the roots of $p_{2n-2}$ on the unit circle, for instance by computing the eigenvalues of its companion matrix . As shown in Table 5, this scheme allows to recover the support with very high precision. Having obtained the estimate for the support $\hat{T}$, the amplitudes of the signal can be reconstructed by solving the system of equations $\sum_{t\in\hat{T}}e^{-i2\pi kt}a_{t}=y_{k}$, $\left|k\right|\leq f_{c}$, using the method of least squares. There is a unique solution as we have at most $n-1$ columns which are linearly independent since one can add columns to form a Vandermonde system.The set of roots contains the support of a primal optimal solution; if it is a strict superset, then some amplitudes will vanish.Figure 9 illustrates the accuracy of this procedure; a Matlab script reproducing this example is available at http://www-stat.stanford.edu/~candes/superres_sdp.m.

which shows that $c$ is a solution to the dual (4.3) that does not carry any information about the support of $x$. Fortunately, this situation is in practice highly unusual. In fact, it does not occur as long as

and we use interior point methods as in SDPT3 to solve (4.3). (Our simulations use CVX which in turn calls SDPT3.) This phenomenon is explained below. At the moment, we would like to remark that Condition (4.5) is sufficient for the primal problem (1.4) to have a unique solution, and holds except in very special cases. To illustrate this, suppose $y$ is a random vector, not a measurement vector corresponding to a sparse signal. In this case, we typically observe dual solutions as shown in Figure 10 (non-vanishing polynomials with at most $n-1$ roots). To be sure, we have solved 400 instances of (4.3) with different values of $f_{c}$ and random data $y$. In every single case, condition (4.5) held so that we could construct a primal feasible solution $x$ with a duality gap below $10^{-8}$, see Figure 11. In all instances, the support of $x$ was constructed by determining roots of $p_{2n-2}(z)$ at a distance at most $10^{-4}$ from the unit circle.

Interior point methods approach solutions from the interior of the feasible set by solving a sequence of optimization problems in which an extra term, a scaled barrier function, is added to the cost function . To be more precise, in our case (4.3) would become

where $t$ is a positive parameter that is gradually reduced towards zero in order to approach a solution to (4.3). Let $\lambda_{k}$, $1\leq k\leq n$, denote the eigenvalues of $Q-cc^{\ast}$. By Schur’s formula (Theorem 1.1 in ) we have

To conclude, we have merely presented an informal discussion of a semidefinite programming approach to the minimum-total variation problem (1.4). It is beyond the scope of this paper to rigorously justify this approach—for example, one would need to argue that the root finding procedure can be made stable, at least under the conditions of our main theorem—and we leave this to future work along with extensions to noisy data.

Numerical experiments

For a super-resolution factor $\text{SRF}=N/n$, we work with a partial DFT matrix $F_{n}$ with frequencies up to $f_{c}=\lfloor n/2\rfloor$. Fix a candidate minimum distance $\Delta$.

Using a greedy algorithm, construct an adversarial support with elements separated by at least $\Delta$ by sequentially adding elements to the support. Each new element is chosen to minimize the condition number formed by the columns corresponding to the selected elements.

Take the signal $x$ to be the singular vector corresponding to the smallest singular value of $F_{n}$ restricted to $T$.

This construction of an adversarial signal was found to be better adapted to the structure of our measurement matrix than other methods proposed in the literature such as . We used this scheme and a simple binary search to determine a lower bound for the minimum distance that guarantees exact recovery for $N=4096$, super-resolution factors of 8, 16, 32 and 64 and support sizes equal to 2, 5, 10, 20 and 50. The simulations were carried out in Matlab, using CVX to solve the optimization problem. Figure 12 shows the results, which suggest that on the discrete grid we need at least a minimum distance equal to twice the super-resolution factor in order to guarantee reconstruction of the signal (red curve). Translated to the continuous setting, in which the signal would be supported on a grid with spacing $1/N$, this implies that $\Delta\gtrsim\lambda_{c}$ is a necessary condition for exact recovery.

Discussion

In this paper, we have developed the beginning of a mathematical theory of super-resolution. In particular, we have shown that we can super-resolve ‘events’ such as spikes, discontinuity points, and so on with infinite precision from just a few low-frequency samples by solving convenient convex programs. This holds in any dimension provided that the distance between events is proportional to $1/f_{c}=\lambda_{c}$, where $f_{c}$ is the highest observed frequency; for instance, in one dimension, a sufficient condition is that the distance between events is at least $2\lambda_{c}$. Furthermore, we have proved that when such condition holds, stable recovery is possible whereas super-resolution—by any method whatsoever—is in general completely hopeless whenever events are at a distance smaller than about $\lambda_{c}/2$.

2 Extensions

We have focused in this paper on the super-resolution of point sources, and by extension of discontinuity points in the function value, or in the derivative and so on. Clearly, there are many other models one could consider as well. For instance, we can imagine collecting low-frequency Fourier coefficients of a function

where $\{\varphi_{j}(t)\}$ are basis functions. Again, $f$ may have lots of high-frequency content but we are only able to observe the low-end of the spectrum. An interesting research question is this: suppose the coefficient sequence $x$ is sparse, then under what conditions is it possible to super-resolve $f$ and extrapolate its spectrum accurately? In a different direction, it would be interesting to extend our stability results to other noise models and error metrics. We leave this to further research.

Acknowledgements

E. C. is partially supported by NSF via grant CCF-0963835 and the 2006 Waterman Award, by AFOSR under grant FA9550-09-1-0643 and by ONR under grant N00014-09-1-0258. C. F. is supported by a Caja Madrid Fellowship and was previously supported by a La Caixa Fellowship. E. C. would like to thank Mikhail Kolobov for fruitful discussions, and Mark Davenport, Thomas Strohmer and Vladislav Voroninski for useful comments about an early version of the paper. C. F. would like to thank Armin Eftekhari for a remark on Lemma 2.6.

References

Appendix A Background on the recovery of complex measures

For further details, we refer the reader to .

Proof The proof is a variation on the well-known argument for finite signals, and we note that a proof for continuous-time signals, similar to that below, can be found in . Let $\hat{x}$ be a solution to (1.4) and set $\hat{x}=x+h$. Consider the Lebesgue decomposition of $h$ relative to $\left|x\right|$,

The existence of $q$ suffices to establish a valuable inequality between the total-variation norms of $h_{T}$ and $h_{T^{c}}$. Begin with

with a strict inequality if $h\neq 0$. Assuming $h\neq 0$, we have

This is a contradiction and thus $h=0$. In other words, $x$ is the unique minimizer.

Appendix B Proof of Lemma 2.6

The first inequality in the lemma holds due to two lower bounds on the sine function:

The proof for these expressions, which we omit, is based on concavity of the sine function and on a Taylor expansion around the origin. Put $f=f_{c}/2+1$ for short. Some simple calculations give $K^{\prime}(0)=0$ and for $t\neq 0$,

Further calculations show that the value of the second derivative of $K$ at the origin is $-\pi^{2}f_{c}\left(f_{c}+4\right)/3$, and for $t\neq 0$,

It is also possible to check that the third derivative of $K$ is zero at the origin, and for $t\neq 0$,

Appendix C Proof of Theorem 1.3

Theorem 1.3 follows from Proposition C.1 below, which guarantees the existence of a dual certificate. In this section, we write $\Delta=\Delta(T)\geq\Delta_{\text{min}}=2.38\,\lambda_{c}$. Unless specified otherwise, $|r-r^{\prime}|$ is the $\infty$ distance.

The proof is similar to that of Proposition 2.1 in that we shall construct the dual polynomial $q$ by interpolation with a low-pass, yet rapidly decaying two-dimensional kernel. Here, we consider

obtained by tensorizing the square of the Fejer kernel (2.3). (For reference, if we had data in which $y(k)$ is observed if $\|k\|_{2}\leq f_{c}$, we would probably use a radial kernel.) Just as before, we have fixed $K$ somewhat arbitrarily, and it would probably be possible to optimize this choice to improve the constant factor in the expression for the minimum distance. We interpolate the sign pattern using $K^{\text{2D}}$ and its partial derivatives, denoted by $K^{\text{2D}}_{\left(1,0\right)}$ and $K^{\text{2D}}_{\left(0,1\right)}$ respectively, as follows:

and we fit the coefficients so that for all $t_{j}\in T$,

The first intermediate result shows that the dual polynomial is well defined, and also controls the magnitude of the interpolation coefficients.

Under the hypotheses of Proposition C.1, there are vectors $\alpha$, $\beta_{1}$ and $\beta_{2}$ obeying (C.2) and

where $\beta=(\beta_{1},\beta_{2})$. Further, if $v_{1}=1$,

Proposition C.1 is now a consequence of the two lemmas below which control the size of $q$ near a point $r_{0}\in T$. Without loss of generality, we can take $r_{0}=0$.

Assume $0\in T$. Then under the hypotheses of Proposition C.1, $\left|q\left(r\right)\right|<1$ for all $0<|r|\leq 0.2447\,\lambda_{c}$.

Assume $0\in T$. Then under the conditions of Proposition C.1, $\left|q\left(r\right)\right|<1$ for all $r$ obeying $0.2447\,\lambda_{c}\leq\left|r\right|\leq\Delta/2$. This also holds for all $r$ that are closer to $0\in T$ (in the $\infty$ distance) than to any other element in $T$.

To express the interpolation constraints in matrix form, define

We split this sum into different regions corresponding to whether $|x_{j}|$ or $|y_{j}|\leq\Delta/2$ and to $\min(|x_{j}|,|y_{j}|)\geq\Delta/2$. First,

This holds because the $x_{j}$’s must be at least $\Delta$ apart, $B_{0}$ is nonincreasing and the absolute value of $K^{\text{2D}}$ is bounded by one. The region $\{r_{j}\neq 0,\,|x_{j}|<\Delta/2\}$ yields the same bound. Now observe that Lemma C.5 below combined with Lemma 2.6 gives

To bound this expression, we apply the exact same technique as for (2.11) in Section 2.3, starting at $j=0$ and setting $j_{0}=20$. This gives

In turn, the same upper-bounding technique yields

where we have used the fact that $\|K^{\prime}\|_{\infty}\leq 2.08\left(f_{c}+2\right)$, which follows from combining Lemma 2.6 with (2.21). Likewise,

since $\|K^{\prime\prime}\|_{\infty}=\pi^{2}f_{c}\left(f_{c}+4\right)/3$, as $\left|K^{\prime\prime}\right|$ reaches its global maximum at the origin.

We use these estimates to show that the system (C.5) is invertible and to show that the coefficient sequences are bounded. To ease notation, set

Note that $S_{1}$ is a Schur complement of $D_{\left(0,2\right)}$ and that a standard linear algebra identity gives

Applying (2.13) from Section 2.3.1, we obtain

which together with $K^{\text{2D}}_{\left(2,0\right)}\left(0\right)=-\pi^{2}f_{c}\left(f_{c}+4\right)/3$ and (C.8) imply

Another application of (2.13) then yields

Next, (C.7), (C.8) and (C.10) allow to bound $S_{2}$,

which combined with (C.6), (C.7), (C.10) and (C.11) implies

The results above allow us to derive bounds on the coefficient vectors by applying (2.13) one last time, establishing

where the last lower bound holds if $v_{1}=1$. The derivation for $\left|\left|\beta_{2}\right|\right|_{\infty}$ is identical and we omit it.

C.2 Proof of Lemma C.3

Since $v$ is real valued, $\alpha$, $\beta$ and $q$ are all real valued. For $|r|\leq 0.2447\,\lambda_{c}$, we show that the Hessian matrix of $q$,

is negative definite. In what follows, it will also be useful to establish bounds on the kernel and its derivatives near the origin. Using (2.20)–(2.24), we obtain

These bounds are all monotone in $x$ and $y$ so we can evaluate them at $x=0.2447\,\lambda_{c}$ and $y=0.2447\,\lambda_{c}$ to show that for any $|r|\leq 0.2447\,\lambda_{c}$,

Similarly, the contribution from the bands where either $|r_{j,1}|$ or $|r_{j,2}|\leq\Delta/2$ obeys

Since $\operatorname{Tr}(H)=q_{\left(2,0\right)}+q_{\left(0,2\right)}<0$ and $\det(H)=|q_{\left(2,0\right)}||q_{\left(0,2\right)}|-|q_{\left(1,1\right)}|^{2}>0$, both eigenvalues of $H$ are strictly negative.

We have shown that $q$ decreases along any segment originating at . To complete the proof, we must establish that $q>-1$ in the square. Similar computations show

C.3 Proof of Lemma C.4

For $0.2447\,\lambda_{c}\leq\left|r\right|\leq\Delta/2$,

Using the series expansion around the origin of $K$ and $K^{\prime}$ (2.29), we obtain that for $t_{1}\leq|r|\leq t_{2}$,

The same bound holds for $K^{\text{2D}}_{\left(0,1\right)}$. Now set

where $\alpha^{\infty}$ and $\beta^{\infty}$ are the upper bounds from Lemma C.2. The quantities reported in Table 7 imply that setting $\left\{t_{1},t_{2}\right\}$ to $\left\{0.1649\,\lambda_{c},0.27\,\lambda_{c}\right\}$, $\left\{0.27\,\lambda_{c},0.36\,\lambda_{c}\right\}$, $\left\{0.36\,\lambda_{c},0.56\,\lambda_{c}\right\}$ and $\left\{0.56\,\lambda_{c},0.84\,\lambda_{c}\right\}$ yields $\left|q\right|<0.9958$, $\left|q\right|<0.9929$, $\left|q\right|<0.9617$ and $\left|q\right|<0.9841$ respectively in the corresponding intervals. Finally, for $\left|r\right|$ between $0.84\,\lambda_{c}$ and $\Delta/2$, applying Lemma (2.6) yields $W\left(r\right)\leq 0.5619$, $Z_{\left(0,0\right)}\left(0.84\,\lambda_{c}\right)\leq 0.3646$ and $Z_{\left(0,1\right)}\left(0.84\,\lambda_{c}\right)\leq 0.6502\,f_{c}$ , so that $\left|q\right|\leq 0.9850$. These last bounds also apply to any location beyond $\Delta/2$ closer to than to any other element of $T$ because of the monotonicity of $B_{0}$ and $B_{1}$. This concludes the proof.