Super-resolution via superset selection and pruning

I Introduction

where $A$ is the partial, short and wide Fourier matrix $A_{jk}=e^{2\pi ijk/n}$, $0\leq j<m$, $-n/2\leq k<n/2$, $n$ even, and, say, $e\sim N(0,\sigma^{2}I_{m})$.

When recovery is successful in this scenario of contiguous measurements, we may speak of super-resolution: the spacing between neighboring nonzero components in $x_{0}$ can be much smaller than the Rayleigh limit $n/m$ suggested by Shannon-Nyquist theory. But in contrast to the compressed sensing scenario, where the $m$ values of $j$ are drawn at random from $\{0,\ldots,n-1\}$, super-resolution can be arbitrarily ill-posed. Open questions concern not only recovery bounds, but the very algorithms needed to define good estimators.

In this paper, we discuss a simple algorithm for solving (1) based on

subspace identification as in the matrix pencil method, but without the subsequent eigenvalue computation; and

a removal procedure for tightening the active set, remindful of a step in certain greedy pursuits.

II Noiseless subspace identification

For completeness we start by recalling the classical uniqueness result for (2).

The “superset method” hinges on a special property that the partial Fourier matrix $A$ does not share with arbitrary dictionaries: each column $a_{k}$ is translation-invariant in the sense that any restriction of $a_{k}$ to $s\leq m$ consecutive elements gives rise to the same sequence, up to an overall scalar. In other words, exponentials are eigenfunctions of the translation operator. This structure is important. There is an opportunity cost in ignoring it and treating (1) as a generic compressed sensing problem.

A way to leverage translation invariance is to recognize that it gives access to the subspace spanned by the atoms $a_{k}$ for $k\in T$, such that $y=\sum_{k\in T}(x_{0})_{k}a_{k}$. Algorithmically, one picks a number $1<L<m$ and juxtaposes translated copies of (restrictions of) $y$ into the Hankel matrix $Y=\mbox{Hankel}(y)$, defined as

The range of $Y$ is the subspace we seek.

If $L\geq|T|$, then the rank of $Y$ is $|T|$, and

The lemma suggests a simple recovery procedure in the noiseless case: loop over all the candidate atoms $a_{k}$ for $-n/2\leq k<n/2$ and select those for which the angle

Once the set $T$ is identified, the solution is obtained by solving the determined system

The proofs of lemma 2 and theorem 3 hinge on the fact that $A$ has full spark.

The idea of subspace identification is at the heart of a different method, the matrix pencil, which seeks the rank-reducing numbers $z$ of the pencil

where $\overline{Y}$ is $Y$ with its first row removed, and $\underline{Y}$ is $Y$ with its last row removed. These numbers $z$ are computed as the generalized eigenvalues of the couple $(\underline{Y}^{*}\overline{Y},\underline{Y}^{*}\underline{Y})$. $z$ can also be found via solving the eigenvalues of the matrix $\underline{Y}^{\dagger}\overline{Y}$. When $|T|\leq L\leq m-|T|$, the collection of these generalized eigenvalues includes $e^{2\pi ijk/n}$ for $k\in T$, as well as $m-L-|T|$ zeros. There exist variants that consider a Toeplitz matrix instead of a Hankel matrix, with slightly better numerical stability properties. When $L=|T|$, the matrix pencil method reduces to Prony’s method, a numerically inferior choice that should be avoided in practice if possible.

III Noisy subspace identification

The problem becomes more difficult when the observations are contaminated by noise. In this situation $\mbox{Ran}\,A_{T}^{L}\neq\mbox{Ran}\,Y$, though in low-noise situations we may still be able to recover $T$ from the indices of the smallest angles $\angle(a^{L}_{k},\mbox{Ran}\,Y)$.

Let $y=y_{0}+e$ with $e\sim N(0,\sigma^{2}I_{m})$, and form the corresponding $L\times(m-L)$ matrices $Y$ and $Y_{0}$ as previously. Denote the singular values of $Y^{m-L}_{0}$ by $s_{n,0}$. Then there exists positive $c_{1},C_{1}$ and $c$, such that with probability at least $1-c_{1}m^{-C_{1}}$,

for all indices $k$ in the support set and

Here we sketch the proof of this proposition. We note that $a^{L}_{k}\in\mbox{Ran}\,Y_{0}$ when $k$ is in the true support. Thus

Denote the compact singular value decomposition of $A^{L}_{T}=US^{L}V^{*}$. Recalling that $a^{L}_{k}\in\mbox{Ran}\,Y_{0}$ and a well-known fact that $Y_{0}=A^{L}_{T}D(A^{m-L}_{T})^{*}$ where $D=\operatorname{diag}((x_{0})_{T})$, we can write $a^{L}_{k}=U\alpha=\sum_{i=1}^{|T|}\alpha_{i}u_{i}$. Thus,

Next, since $Y=Y_{0}+E=A^{L}_{T}D(A^{m-L}_{T})^{*}+E$, we have $Y[D(A^{m-L}_{T})^{*}]^{\dagger}=A^{L}_{T}+E[D(A^{m-L}_{T})^{*}]^{\dagger}$ where $A^{\dagger}$ is the pseudo-inverse matrix of $A$. By multiplying both sides by $(\mathcal{P}_{Y^{\perp}}u_{i})^{*}$, we get

Since the vector $\mathcal{P}_{Y^{\perp}}u_{i}$ is orthogonal to $\mbox{Ran}\,Y$, the left hand side is zero. Thus multiplying both sides by $v_{i}$, the $i$-th right singular vector of $A^{L}_{T}$, we have

We can see that $(\mathcal{P}_{Y^{\perp}}u_{i})^{*}A^{L}_{T}v_{i}=(\mathcal{P}_{Y^{\perp}}u_{i})^{*}s^{L}_{i}u_{i}=s^{L}_{i}\left\|\mathcal{P}_{Y^{\perp}}u_{i}\right\|_{2}^{2}$ where $s^{L}_{i}$ is the $i$-th singular value of $A^{L}_{T}$. We therefore obtain

where $s^{m-L}_{|T|}$ is the smallest singular value of $A^{m-L}_{T}$.

Recalling that $a^{L}_{k}=U\alpha$, we have $\alpha_{i}=u^{*}_{i}a^{L}_{k}$. From the SVD of $A^{L}_{T}$, we see that $A^{L}_{T}(A^{L}_{T})^{*}=U(S^{L})^{2}U^{*}$, so that

This identity implies that $\left\|u^{*}_{i}A^{L}_{T}\right\|_{2}=s^{L}_{i}$, and thus, $|\alpha_{i}|\leq s^{L}_{i}$. Combining this result with (III) and (7) yields

Using the matrix Bernstein inequality of one obtains that $\left\|E\right\|\leq\sigma\sqrt{cL\log m}$ with high probability. Finally, writing $Y^{m-L}_{T}$ as $Y^{m-L}_{T}=A^{m-L}_{T}D^{1/2}(D^{1/2})^{*}(A^{m-L}_{T})^{*}$, we have

There are a few unknown quantities involving $\epsilon_{1}$, which can empirically be controlled. The support size $T$ can be estimated by a reasonably large constant, say $m/2$. The dynamic range of the signal can presumably be known if we know in prior the type of underlying signal of interest. The singular value $s_{|T|,0}$ of $Y^{m-L}_{0}$ can be replaced by that of $Y^{m-L}$ via the simple Weyl’s inequality $|s_{i}-s_{i,0}|\leq\|\,\mbox{Hankel}(e)\|$, which can in turn be controlled as $O(\sigma\sqrt{L\log m})$ with high probability.

The subspace identification step now gathers all the values of $k$ such that

The resulting set $\Omega$ of indices is only expected to be a superset of the true support $T$, with high probability.

A second step is now needed to prune $\Omega$ in order to extract $T$. For this purpose, a loop over $k$ is set up where we test the membership of $y$ in $\mbox{Ran}\,A_{\Omega\backslash k}$, the range of $A_{\Omega}$ with the $k$-th column removed. We are now considering a new set of angles where the roles of $y$ and $A$ are reversed: in a noiseless situation, $k\in T$ if and only if

When noise is present, we first filter out the noise off $\Omega$ by projecting $y$ onto the range of $A_{\Omega}$, then estimate $k\in T$ only when the angle is above a certain threshold. It is easier to work directly with projections $\Pi$:

The effect of noise on the left-hand side is as follows.

Let $y=y_{0}+e$ with $e\sim N(0,\sigma^{2}I_{m})$. Let $\Pi_{\Omega}y$ be the projection of $y$ onto $\mbox{Ran}\,A_{\Omega}$, and let $\Delta\Pi=\Pi_{\Omega}-\Pi_{\Omega\backslash k}$. Then there exists $c>0$ such that, with high probability,

Algorithm 1 for the superset method implements the removal step in an iterative fashion, one atom at a time.

IV Experimental Results

In the next simulation, we consider a signal of size $n=1000$ which contains two nearby spikes at locations $$and has magnitudes$1/\sqrt{2}$and$-1/\sqrt{2}$. We empirically investigate the algorithm’s ability to recover the signal from varying measurements$m=\{10,20,...,220\}$and noise levels$log_{10}\sigma=\{-3.5,-3.4,...,-2\}$. For each pair$(m,\sigma)$, we report the frequency of success over$100$random realizations of$e$. The greyscale goes from white (100 successes) to black (100 failures). A trial is declared successful if the recovered$\widehat{x}$satisfies$\left\|\widehat{x}-x_{0}\right\|_{2}/\left\|x_{0}\right\|_{2}<10^{-3}$. The horizontal axis indicates the noise level$\sigma$in log scale, and the vertical axis indicates$\log_{10}(1-\mu)$where$\mu$ is the coherence as earlier.

We note that the coherence is inversely proportional to the amount of measurements $m$ and proportional to the super-resolution factor $n/m$: increasing $m$ (decreasing the super-resolution factor) will reduce the coherence $\mu$. On the vertical axis, smaller values imply higher coherence, or equivalently smaller amount of measurements. As shown in Fig. 2, for reasonably small noise, the algorithm is able to recover the signal exactly even the coherence is nearly $1$.

For reference, we also compare the superset method with the matrix pencil method as set up in . The noise is filtered out by preparing low-rank approximations of $\underline{Y}$ and $\overline{Y}$ where only the singular values above $c\sigma\sqrt{L\log L}$ are kept, for some heuristically optimized constant $c$. Two more signals are considered: (1) a 3-sparse signal consisting of three neighboring spikes, each of magnitude $1/\sqrt{3}$ with alternating signs, and (2) a 4-sparse signal with neighboring spikes of alternating signs and equal magnitude $1/2$. Fig. 2 is a good illustration of the contrasting numerical behaviors of the two methods: the matrix pencil is often the better method in the special case of a signal with 2 spikes, but loses ground to the superset method in various cases of progressively less sparse signals. Understanding the performance of the matrix pencil would require formulating a lower bound on the (typically extremely small) $S$-th eigenvalues of $Y_{0}$ where $S$ is the sparsity of $y_{0}$.

V Conclusion

Empirical evidence is presented for the potential of the superset method as a viable computational method for super-resolution. Further theoretical justifications will be presented elsewhere.

References