MUSIC for Single-Snapshot Spectral Estimation: Stability and Super-resolution

Introduction

The field of Compressive Sensing (CS) has provided us with a new technology of reconstructing a signal from a small number of linear measurements. With a new exceptions, signals considered in the compressive sensing community are assumed to be sparse under a discrete, finite-dimensional dictionary.

However, signals arising in applications such as radar , sonar and remote sensing are represented by few parameters on a continuous domain. These signals are usually not sparse under any discrete dictionary but can be approximately sparsely represented by indicator functions on a discrete domain. An approximation error, called gridding error or basis mismatch exists, manifesting the gap between the continuous world and the discrete world. This issue is well illustrated by the spectral estimation problem as follows.

Suppose a signal $y(t)$ consists of linear combinations of $s$ time-harmonic components from the set

where $\varepsilon(t)$ is the external noise.

be the imaging vector of size $M+1$ at the frequency $\omega$ and define

The single-snapshot formulation of spectral estimation takes the form

One can attempt to linearize (3) by expanding the matrix $\Phi^{M}$ via setting up a grid

where $N$ is some large integer, and writing the spectral estimation problem in the form a linear inversion problem

Discretizing $[0,1)$ as in (4) amounts to rounding frequencies on the continuum to the nearest grid points in $\mathcal{G}$, giving rise to a gridding error which is roughly proportional to the grid spacing. On the other hand, as $N$ increases, correlation among adjacent columns of $A$ also increases dramatically .

A key unit of frequency separation is the Rayleigh Length, roughly the minimum resolvable separation of two objects with equal intensities in classical resolution theory . Mathematically, the Rayleigh Length (RL) is the distance between the center and the first zero of the Dirichlet kernel

The ratio $F=N/M$ between RL and the grid spacing is called the refinement factor in and super-resolution factor in . The higher $F$ is, the more coherent the measurement matrix $A$ becomes.

In this paper, to circumvent the gridding problem, we reformulate the spectral estimation problem (3) in the form of multiple measurement vectors that is suitable for the application of the MUltiple Signal Classification (MUSIC) algorithm , widely used in signal processing and array imaging .

Most state-of-the-art spectral estimation methods ( and references therein) assume many snapshots of array measurement as well as statistical assumptions on measurement noise. In contrast, we pursue below a deterministic approach to spectral estimation with a single snapshot of array measurement in common with .

Fixing a positive integer $1\leq L<M$, we form the Hankel matrix

Since its first appearance in Prony’s method the Hankel data matrix (6) plays an important role in modern methods such as the state space method and the matrix pencil method .

It is straightforward to verify that ${\rm Hankel}(y)$ with $y=\Phi^{M}x$ admits the Vandermonde decomposition

Here we use a special property of Fourier measurements: a time translation corresponds to a frequency phase modulation.

Let $H^{\varepsilon}={\rm Hankel}(y^{\varepsilon})$ and $E={\rm Hankel}(\varepsilon)$. The multiple measurement vector formulation of spectral estimation takes the form

More specifically, let the Singular Value Decomposition (SVD) of $H$ be written as

with the singular values $\sigma_{1}\geq\sigma_{2}\geq\sigma_{3}\geq\cdots\sigma_{s}>0.$ The signal and noise spaces are exactly the column spaces of $U_{1}$ and $U_{2}$ respectively.

The following fact is the basis for noiseless MUSIC (See Appendix A for proof).

Suppose $\omega_{k}\neq\omega_{l}\ \forall k\neq l$. If

$\omega\in\mathcal{S}\Longleftrightarrow R(\omega)=0\Longleftrightarrow J(\omega)=\infty$.

Condition (9) says that the number of measurement data $(M+1)\geq 2s$ suffices to guarantee exact reconstruction by the MUSIC algorithm.

For the noisy data matrix $H^{\varepsilon}$ let the SVD be written as

with the singular values $\sigma^{\varepsilon}_{1}\geq\sigma^{\varepsilon}_{2}\geq\sigma^{\varepsilon}_{3}\geq\cdots.$ The noise-space correlation function and imaging function become

respectively with $\mathcal{P}^{\varepsilon}_{2}=U^{\varepsilon}_{2}(U^{\varepsilon}_{2})^{\star}$.

Figure 1 shows a noise-space correlation function and an imaging function in the noise-free case. True frequencies are exactly located where the noise-space correlation function vanishes and the imaging function peaks.

The MUSIC algorithm as formulated above requires the number of frequencies $s$ as an input. There are some techniques for evaluating how many objects are present in the event that such information is not available. When $\sigma_{s}\gg 2\|E\|_{2}$, $s$ can be easily estimated based on the singular value distribution of $H^{\varepsilon}$ due to Weyl’s theorem .

$|\sigma^{\varepsilon}_{j}-\sigma_{j}|\leq\|E\|_{2},\ j=1,2,\ldots$

As a result, $\sigma^{\varepsilon}_{j}\leq\|E\|_{2},\ \forall j\geq s+1$ and $\sigma^{\varepsilon}_{s}\geq\sigma_{s}-\|E\|_{2}$. Hence $\sigma^{\varepsilon}_{s}\gg\sigma^{\varepsilon}_{s+1}$, creating a gap between $\sigma_{s}^{\varepsilon}$ and $\{\sigma^{\varepsilon}_{j}:j\geq s+1\}$. An example is shown in Figure 3.

For simplicity, we denote $\Phi^{M}=\Phi^{0\rightarrow M}.$

2 Contribution of the present work

The main contribution of the paper is a stability analysis for the MUSIC algorithm with respect to general support set $\mathcal{S}$ and external noise.

In the MUSIC algorithm frequency candidates are identified at the $s$ smallest local minima of the noise-space correlation which measures how much an imaging vector is correlated with the noise space. In noise-free case, the noise-space correlation function $R(\omega)$ vanishes exactly on $\mathcal{S}$. For the noisy case we prove

To make the bounds (10) explicit and more meaningful, we prove the discrete Ingham inequalities (Corollary 1) which implies

Furthermore, we prove that for every $\omega_{j}\in\mathcal{S}$, there exists a local minimizer $\hat{\omega}_{j}$ of $R^{\varepsilon}$ such that $\hat{\omega}_{j}\to\omega_{j}$ as noise decreases to .

To relax the restriction on the minimum separation between adjacent frequencies, condition (13) suggests that $L$ should be about $M/2$ and then the resolving power of the present form of MUSIC is as good as $2/M=$ 2 RL.

By the results of , the spectral norm of the random Hankel matrix $E$ constructed from a zero mean, independently and identically distributed (i.i.d.) sequence of a finite variance is on the order of $\sqrt{M\log M}$ for $M\gg 1$ while $\sigma_{s}$ is on the order of $M$ (with $L\approx M/2$). In this case the factor $\alpha$ in (10) is almost always positive for sufficiently large $M$ regardless of the variance of noise and

Also the super-resolution effect of MUSIC is studied. When the minimum separation between frequencies drops below 1 RL, we show that the noise level that MUSIC can tolerate obeys a power law with respect to the minimum separation with an exponent smaller than an estimate established by Donoho.

Our analysis can be easily extended to other settings where the MUSIC algorithm can be applied, such as the estimation of Directions of Arrivals (DOA) and inverse scattering .

3 Comparison with other works

Other closely related work includes and where Vandermonde decomposition of the Hankel matrix (6) are used to design different algorithms.

In Demanet et al. proposed an approach to spectral estimation with a selection step of the support set followed by a pruning step. In the selection step, any $\omega$ satisfying $\sin\measuredangle(\phi^{L}(\omega),{\rm Range}H^{\varepsilon})$ $\leq\eta$ for some judicious choice of $\eta>0$ is kept as a frequency candidate based on their estimate

In Chen and Chi exploited the low-rank property of the Hankel matrix $H$ and applied the matrix completion technique to recover a spectrally sparse signal from its partial time-domain samples. The focus of is on the completion and denoising of of data from the partial noisy samples while MUSIC is designed for frequency recovery. A combination of and our work constitutes a new framework for single-snapshot spectral estimation with compressive noisy measurements which is to be discussed in Section 6.

As for frequency recovery, recent progresses center around greedy algorithms and Total Variation (TV) minimization.

The challenge of applying greedy algorithms to (5) while $N\gg M$ lies in the high coherence and ill conditioning of the sensing matrix $A$. In order to mitigate this effect, we exploited the coherence pattern of $A$ and introduced the techniques of Band exclusion and Local Optimization (BLO) to enhance standard compressive sensing algorithms. The performance guarantee in assumes $q\geq 3$ RL and ensures reconstruction of $\mathcal{S}$ to the accuracy of $1$ RL.

In , Candès and Fernandez-Granda proposed TV minimization and showed that, under the assumption of $q\geq 4$ RL, the TV minimizer yields an $L_{1}$ reconstruction error linearly proportional to noise with a magnification factor proportional to $F^{2}$ where $F$ is the refinement/super-resolution factor. Inspired by this approach, Tang et al. developed an atomic norm (equivalent to the TV norm in 1D) minimization for the completion of $y$ from its partial samples and showed exact reconstruction in the noise-free case. Like , a main emphasis in is on sparse measurements. Unfortunately, the effect of noise is not considered in . For numerical implementation a SemiDefinite Programming (SDP) on the dual problem is solved in where numerical efficiency and stable retrieval of primal solutions may become a problem.

Historically, Prony was the first to address the problem of spectral estimation . Unfortunately, Prony’s method is numerically unstable and numerous modifications were attempted to improve its numerical behavior. Approximate Prony Method (APM) proposed by Beylkin and Monzón in is a major breakthrough for function approximation by exponential sums. Specifically, Beylking and Monzón considered the following problem: given $2N+1$ values of function $f(t)$ on a uniform grid on $$and a target accuracy$\varepsilon>0$, they find the minimal number$s$of complex weights$w_{j}$and complex nodes$\gamma_{j}$ such that

Many interesting examples were provided in . For instance, the Bessel function $J_{0}(100\pi t)$ in $$is approximated by exponential sums of$28$complex nodes with accuracy$\varepsilon=10^{-10}$ by APM.

In comparison the spectral estimation problem (1) is the identification of $\{\gamma_{j}\}$ from noisy data, instead of approximation of the signal. For spectral estimation with noisy data, APM’s stability may be questionable. The numerical examples of spectral estimation by APM in all have low NSR $=\mathcal{O}(10^{-\delta})$ where $\delta\geq 4$. In contrast our simulations in Section 5 are performed with NSR as large as $0.5$. Furthermore, the super-resolution effect of MUSIC is quantitatively documented in Section 4 while it has not been reported in literature whether APM has the capability of localizing closely spaced frequencies.

In terms of discrete Ingham inequalities, Theorem 2 is the first result in which both the gap condition and the upper/lower bounds are explicitly given. In comparison semi-discrete Ingham inequalities in [32, 39, Lemma 3.1] give the correct scaling with respect to the size of the Vandermonde matrix $\Phi^{L}$ without an explicit estimate for the constants (cf. (17) and (20) in Section 2). In other words the previous Ingham inequalities affirm only that the matrix $\Phi^{L}$ has a finite condition number under certain gap condition of $\mathcal{S}$ without an explicit estimate on the magnitude of the condition number.

Detailed numerical comparisons of the MUSIC algorithm with Band-excluded Locally Optimized Orthogonal Matching Pursuit (BLOOMP) of , SemiDefinite Programming (SDP) of and Matched filtering using prolates enhanced by the Band-excluded and Locally Optimized technique are presented in Section 5.

Since the SVD step is its primary computational cost, MUSIC has low computational complexity compared to other existing methods. As we will also see, MUSIC is also among the most accurate algorithms. Finally, MUSIC is the only algorithm that can resolve frequencies with complex amplitudes closely spaced below 1 RL. Indeed, the resolution of MUSIC can be arbitrarily small for sufficiently small noise.

The paper is organized as follows. We estimate nonzero singular values of rectangular Vandermonde matrices with nodes on the unit disk in Section 2. Perturbation theory for MUSIC is presented in in Section 3 and super-resolution effect of MUSIC is studied in Section 4. Numerical experiments are provided in Section 5. We finally conclude and discuss extensions of our current work in Section 6.

Vandermonde matrices with nodes on the unit circle

Performance of the MUSIC algorithm in the presence of noise is crucially dependent on $\sigma_{1}$ and $\sigma_{s}$, the maximum and minimum nonzero singular values of the noiseless Hankel data matrix. To pave a way for the stability analysis, we discuss singular values of the rectangular Vandermonde matrix $\Phi^{L}$ in this section.

then the system of complex exponentials $\{e^{-2\pi i\omega_{j}t},\ t\in[-T/2,T/2],\ j=1,\ldots,s\}$ form a Riesz basis of its span in $L^{2}[-T/2,T/2]$, i.e.,

Ingham inequalities can be considered as a generalization of the Parseval’s identity for non-harmonic Fourier series. The gap condition is necessary for a positive lower bound in (15) but the upper bound always holds.

We prove a discrete version of Ingham inequalities.

Suppose $\mathcal{S}$ satisfies the gap condition

Proof of Theorem 2 is provided in Appendix B.

The difference between the bounds of the discrete and the continuous Ingham inequalities is $\mathcal{O}(1/L)$ which is negligible when $L$ is large. The upper bound in (17) holds even when the gap condition (16) is violated; however, (16) is necessary for the positivity of the lower bound.

Some form of discrete Ingham inequalities are developed in for the analysis of the control/observation properties of numerical schemes of the 1-d wave equation. The main result therein is that when time integrals in (15) are replaced by discrete sums on a discrete mesh, discrete Ingham inequalities converge to the continuous one as the mesh becomes infinitely fine. Their asymptotic analysis, however, do not provide the non-asymptotic results stated in Theorem 2.

Perturbation of noise-space correlation

In this section we use tools in classical matrix perturbation theory and develop a perturbation estimate on the noise-space correlation function, the key ingredient of the MUSIC algorithm. Our main results are presented in Theorem 3, Corollary 1 and Theorem 4 and proofs are provided in Appendix C.1 and C.2.

Suppose $L\geq s$, $M-L+1\geq s$ and $\|E\|_{2}<\sigma_{s}$. Then

In particular, for $\omega_{j}\in\mathcal{S}$, $R(\omega_{j})=0$ and

Suppose noise vector $\varepsilon$ contains i.i.d. random variables of variance $\sigma^{2}$, $\|E\|_{2}\leq\|E\|_{F}=\mathcal{O}(\sigma)$. For fixed $M$ and $\mathcal{S}$,

Theorem 3 holds for all signal models with any support set $\mathcal{S}$. In view of the Vandermonde decomposition (7) we next derive explicit bounds for the perturbation of noise-space correlation by combining Theorem 2 and 3.

Let $L$ and $M-L$ be even integers. Suppose $\mathcal{S}$ satisfies the following gap condition

Corollary 1 is stated in the case that both $L$ and $M-L$ are even integers. However, the results hold in other cases with a slightly different $\alpha_{1}$ based on (17) and (20).

As noted in Section 1.2, for i.i.d. noise $\|E\|_{2}$ grows like $\sqrt{M\log M}$ which is much smaller than $\sqrt{L(M-L)}$ with $L\sim M$ as $M\to\infty$. As a consequence,

How close are the MUSIC estimates, namely the $s$ lowest local minimizers of $R^{\varepsilon}(\omega)$, to the true frequencies which are the zeros of $R(\omega)$? While we can not at the moment answer this question directly, the following asymptotic result says that every true frequency has a near-by strict local minimizer of $R^{\varepsilon}$ in its vicinity converging to it. Denote $[\phi^{L}(\omega)]^{\prime}={d\phi^{L}(\omega)}/{d\omega}$.

When $\|E\|_{2}$ is sufficiently small (details in (63) and (64)), there exists a strict local minimizer $\hat{\omega}_{j}$ of $R^{\varepsilon}(\omega)$ near $\omega_{j}$ such that

where $Q(\omega)=R^{2}(\omega)$, $\eta(L)=2\pi\sqrt{1^{2}+2^{2}+\ldots+L^{2}}/\sqrt{L+1}$ and

For i.i.d. random variables of variance $\sigma^{2}$, $\|E\|_{2}=\mathcal{O}(\sigma)$ for fixed $M$. Hence

cf. (59), assumption (27) is a generic condition that says the true frequencies are not degenerate minimizers of $R^{2}$. Figure 2 shows $\|\mathcal{P}_{2}[\phi^{L}(\omega)]^{\prime}\|_{2}/\|[\phi^{L}(\omega)]^{\prime}\|_{2}$ for various $M$ and suggests that for $q\geq$4 RL,

for some constants $C_{1},C_{2}>0$ independent of $M$.

Under the assumption (23), the right hand side of (28) scales like $\sigma\sqrt{M\log M}$ with $L=M/2$ for i.i.d. random noise of variance $\sigma^{2}$. Under the assumption (29),

As shown in the proof of Theorem 4 (Step 1),

Super-resolution effect of MUSIC

Super-resolution refers to the capability of resolving frequencies separated below 1 RL. The super-resolution effect of MUSIC has been numerically demonstrated in various applications , but theoretical guarantees are lacking. In this section we aim to analyze the super-resolution effect in light of the preceding results and that of .

With $2s$ noiseless data MUSIC guarantees to exactly recover $s$ distinct frequencies. With noisy data, the resulting perturbation of the noise-space correlation function depends crucially on $\sigma_{1}$ and $\sigma_{s}$ (Theorem 3). As

the larger $\sigma_{\rm min}(\Phi^{L}),x_{\text{min}},\sigma_{\rm min}(\Phi^{M-L})$ and the smaller $\sigma_{\rm max}(\Phi^{L}),x_{\text{max}},\sigma_{\rm max}(\Phi^{M-L})$ are, the less sensitive MUSIC is to noise. It follows that a close-to-unity dynamic range $x_{\text{max}}/x_{\text{min}}$ and good conditioning of $\Phi^{L}$ and $\Phi^{M-L}$ are conducive to the stability of MUSIC.

In particular, the denominator $(\sigma_{s}-\|E\|_{2})^{2}$ on the right hand side of (21) indicates that the amount of noise that can be tolerated by MUSIC is approximately $\sigma_{s}$, which by (31) is at least $x_{\text{min}}\sigma_{\rm min}(\Phi^{L})\sigma_{\rm min}(\Phi^{M-L})$.

Hence, to understand the super-resolution effect of MUSIC it is essential to estimate the smallest nonzero singular value of $\Phi^{L}$ with closely spaced frequencies. Under certain weakened gap conditions proposed in , we provide an explicit upper bound on $\sigma_{\rm max}(\Phi^{L})$ and discuss the possible implication of the bounds on $\sigma_{\rm min}(\Phi^{L})$ in on the super-resolution effect.

Suppose $\mathcal{S}$ is an $s$-periodic sequence and satisfies the following weakened gap condition

By choosing $R$ and $\rho$ such that $R\rho=2/L$, in which case there are at most $R$ frequencies in any interval of 4 RL (1 RL = 1/(2L) when L = M/2), we can maintain $B(R\rho,L)$ roughly independent of $L\gg 1$ and obtain the (asymptotic) upper bound ${17\sqrt{2}\over 4\pi}R$ for $\sigma_{\rm max}^{2}(\Phi^{L})/L$. It is noteworthy that the minimum separation between two consecutive frequencies significantly affect the least nonzero singular value but not the largest one.

Presently we can not prove an explicit lower bound for $\sigma_{\rm min}(\Phi^{L})$ when two or more frequencies are spaced below 2 RL (1 RL = 1/(2L) when L = M/2).

In other words, $R_{*}$ is the size of the largest cluster whose members are separated from each other by less than $4$ RL. Define

where $\alpha(L,R_{*})$ is some positive constant depending on $L$ and $R_{*}$. No algorithm is proposed for support reconstruction in .

The definition (36) is the continuous analog of

Hence by making the identification $\Delta=\min_{i\neq j}d(\omega_{i},\omega_{j})=q$, it is plausible that, even with discrete data and without the lattice substrate, the bound

Our preceding analysis in Sections 2 and 3 is for the case $R_{*}=1$.

As commented above, for a support set with Rayleigh index $R_{*}$, the amount of noise that can be tolerated by MUSIC is approximately $\sigma_{s}$ which, according to (39), decays at worst like $x_{\text{min}}q^{4R_{*}+2}$ as $q\rightarrow 0$. Our numerical experiments in Section 5.4 show that the noise level that MUSIC can tolerate obeys $q^{e(R_{*})}$ for $R_{*}=2,3,4,5$ and $e(2)=3.6691$, $e(3)=6.0565$, $e(4)=8.3861$ and $e(5)=11.2392$, suggesting that $e(R_{*})\approx 2.504R_{*}-1.4262$.

Numerical experiments

A systematic numerical simulation is performed on MUSIC, BLOOMP, SDP and Matched Filtering using prolates in this section, showing that MUSIC combines the advantages of strong stability and low computation complexity for the detection of well-separated frequencies and furthermore only MUSIC yields an exact reconstruction in the noise-free case regardless of the distribution of true frequencies and processes the capability of resolving closely spaced frequencies.

We compare the performances of various algorithms on the spectral estimation problem (1) with $M=100$ and i.i.d. Gaussian noise, i.e. $\varepsilon\sim N(0,\sigma^{2}I)+iN(0,\sigma^{2}I)$. Define the

We test and compare the following algorithms.

The MUSIC algorithm: As suggested by (24) in Corollary 1 we set $M=2L$.

When frequencies are separated by at least 4 RL, we choose $r=$ 1 RL.

Matched Filtering (MF) using prolates: In , Eftekhari and Wakin use matched filtering windowed by the Discrete Prolate Spheroidal (Slepian) Sequence for the same problem while frequencies are extracted by band-excluded and locally optimized thresholding proposed in . In its current form, MF using prolates can not deal with complex-valued amplitudes so it is tested with real-valued amplitudes only.

Reconstruction error is measured by Hausdorff distance between the exact ($\mathcal{S}$) and the recovered ($\hat{\mathcal{S}}$) sets of frequencies:

2 Noise-free case

In the noise-free case only MUSIC processes a theory of exact reconstruction regardless of the distribution of true frequencies. In theory, BLOOMP requires a separation of 3 RL for approximate support recovery while SDP requires a separation of 4 RL for exact recovery. In this test we use the four algorithms to recover $15$ real-valued amplitudes separated by $1$ RL. Figure 4 shows that MUSIC achieves the accuracy of about $0.004$ RL while BLOOMP, SDP and MF using prolates essentially fail, which implies that certain separation condition is necessary for BLOOMP, SDP and MF using prolates.

3 Detection of well-separated frequencies

Figure 5 shows reconstructions of $15$ real-valued frequencies separated by $4$ RL. By extracting $15$ largest local maxima of the imaging function $J^{\varepsilon}(\omega)$, MUSIC yields a reconstruction distance of $0.06$ RL. As predicted by the theory in , every recovered object of BLOOMP is within 1 RL distance from a true one. Indeed, in this simulation BLOOMP achieves the best accuracy of $0.05$ RL among tested algorithms. The primal solution of SDP is usually not $s$-sparse and the recovered frequencies tend to cluster around the true ones which degrades the accuracy. The Hausdorff distance between the recovered spikes with the $s$ strongest amplitudes and the true frequencies is 3.94 RL in this simulation. The BET technique can be used to enhance the accuracy of reconstruction and achieve the accuracy of 0.13 RL. Similarly the BLO technique introduced in can be applied to improve the result of Matched filtering windowed by the DPSS sequence (the blue curve in Figure 5(d)) and achieve the accuracy of $0.10$ RL.

Figure 6 shows the average errors of 100 trials by SDP with HT, BET-enhanced SDP, BLOOMP and MUSIC for complex-valued objects separated between 4 RL and 5 RL (Fig. 6(a)(b)) or separated between 2 RL and 3 RL (Fig. 6(c)(d)) versus NSR when dynamic range = 1 (Fig. 6(a)(c)) and when dynamic range = 10 (Fig. 6(b)(d)). In this simulation $[0,1)$ are fully occupied by frequencies satisfying the separation condition and amplitudes $x$ are complex-valued with random phases. Refinement factor $F$ in BLOOMP is adaptive according to the rule: $F=\max(5,\min(1/{\rm NSR},20))$. Figure 6 shows that BLOOMP is the stablest algorithm while frequencies are separated above 4 RL and MUSIC becomes the stablest one while frequencies are separated between 2 RL and 3 RL. Simply extracting $s$ largest amplitudes from the SDP solution (black curve) is not a good idea and the BET technique (green curve) can mitigate the problem with SDP. The average running time in Figure 6 shows that MUSIC takes about 0.33s for one experiment and is the most efficient one among all methods being tested. SDP needs about 20.5s for one experiment and is computationally most expensive. Running time of BLOOMP is dependent on sparsity $s$ and refinement factor $F$. The running time of BLOOMP in Fig. 6(c)(d) is more than the time in Fig. 6(a)(b) as $s\in$ in Fig. 6(c)(d) and $s\in$ in Fig. 6(a)(b).

4 Super-resolution of MUSIC

Theory in Section 4 implies that MUSIC has super-resolution effect and moreover the noise level that MUSIC can handle follows a power law with respect to the minimum separation of the frequencies. We numerically investigate the $2,3,4,5$-point resolution of MUSIC here as numerical verification.

In Figure 7, we consider support set $\mathcal{S}$ containing two, three, four and five equally spaced frequencies. We run MUSIC algorithm on reconstructions of randomly phased complex objects supported on $\mathcal{S}$ with varied separation $q$ and varied NSR for 100 trials and record the average of $d(\mathcal{S},\hat{\mathcal{S}})/q$. Figure 7 (a)-(d) displays the color plot of the logarithm to the base 2 of average $d(\mathcal{S},\hat{\mathcal{S}})/q$ with respect to NSR (y-axis) and $q$ (x-axis) in the unit of RL. Frequency localization is considered successful if $d(\mathcal{S},\hat{\mathcal{S}})/q<1/2$.

A phase transition occurs in (a)-(d), manifesting MUSIC’s capability of resolving two, three, four and five closely spaced complex-valued objects if NSR is below certain level. Theory in Section 4 indicates the noise level that MUSIC can handle scales at worst like $q^{4R_{*}+2}$ where $R_{*}=2$ in Figure 7(a), $R_{*}=3$ in Figure 7(b), $R_{*}=4$ in Figure 7(c) and $R_{*}=5$ in Figure 7(d). The borderline between successful recovery and failure defined by $d(\mathcal{S},\hat{\mathcal{S}})/q=1/2$ are marked out in black in Figure 7 (a)-(d). The phase transition curves for $R_{*}=2,3,4,5$ are shown in (e) in the ordinary scale and in (f) in log-log scale. It appears that the transition curves can be fitted to a constant times $q^{e(R_{*})}$ with $e(2)=3.6691$, $e(3)=6.0565$, $e(4)=8.3861$ and $e(5)=11.2392$, suggesting that a much smaller exponent $e(R_{*})\approx 2.504R_{*}-1.4262$ than $4R_{*}+2$.

Conclusion and extension

We have provided a stability analysis of the MUSIC algorithm for single-snapshot spectral estimation off the grid. We have proved that perturbation of the noise-space correlation by external noise is roughly proportional to the spectral norm of the noise Hankel matrix with a magnification factor given in terms of maximum and minimum nonzero singular values of the Hankel matrix constructed from the noiseless measurements. Under the assumption of frequency separation roughly $\geq$ 2 RL, the magnification factor is explicitly estimated by means of a new version of discrete Ingham inequalities.

A systematic numerical study has shown that the MUSIC algorithm enjoys strong stability and low computation complexity for the reconstruction of well-separated frequencies. MUSIC is the only algorithm that can recover arbitrarily closely spaced frequencies as long as the noise is sufficiently small. And we have numerically documented the super-resolution effect of MUSIC in terms of the relationship among the minimum separation, the Rayleigh index (the size of largest cluster) and the noise. The results conform to the optimal bound conjectured (and partially proved) by Donoho .

Finally we discuss a possible extension of the present work. We became aware of the reference after completing the first draft of this work (arXiv:1404.1484). In , Chen and Chi used the matrix completion technique to obtain a stable approximation of $\{y(k),k=0,\ldots,M\}$ from its partial noisy samples. This can be used as the preprocessing denoising step before invoking the single-snapshot MUSIC. Together and the present work constitute a framework for single-snapshot spectral estimation with compressive noisy measurements.

For the noisy compressive data $\mathcal{P}_{\Lambda}y^{\varepsilon}$ satisfying $\|\mathcal{P}_{\Lambda}(y^{\varepsilon}-y)\|_{2}\leq\delta$, proposes the following denoising strategy of Hankel matrix completion

where $\|\cdot\|_{*}$ denotes the nuclear norm.

The total procedure of compressive spectral estimation is given in the following table.

The following estimate on the difference $\hat{R}(\omega)-R(\omega)$ is obtained by combining [7, Theorem 2] and Theorem 3.

Let $R(\omega)$ and $\hat{R}(\omega)$, respectively, be the noise-space correlation functions for the noiseless $y$ and denoised data $\hat{y}$.

Let $\Lambda$ of size $m$ be uniformly sampled at random from $\{0,\ldots,M\}$. Suppose $\|\mathcal{P}_{\Lambda}(y^{\varepsilon}-y)\|_{2}\leq\delta$. Then there exists a universal constant $C>0$ such that

with probability exceeding $1-(M+1)^{-2}$ provided that

Theorem 6 implies that compressive spectral estimation is stable with matrix completion and MUSIC whenever $\Phi^{L},\Phi^{M-L}$ are well-conditioned and the sample size is sufficiently large.

Furthermore with the discrete Ingham inequalities we can give explicit estimate for the right hand side of (43) and (44). In particular, with $L\approx M/2$ and well-separated ($>$ 2 RL) frequencies, $\mu$ and $\gamma$ scale like a constant and $m=\mathcal{O}(s\log^{3}M)$ suffices for any sufficiently small $\delta$. In this case, the right hand side of (43) is

showing enhanced stability as $M/m\to 1$ where $M/m$ is the compression ratio, $x_{\text{max}}/x_{\text{min}}$ the object’s peak-to-trough ratio and $\delta/x_{\text{min}}$ the noise-to-object ratio.

Acknowledgement

Wenjing Liao would like to thank Armin Eftekhari for providing their codes and helpful discussions at SAMSI.

Appendix A Proof of Theorem 1

In the noise-free case ${\rm Range}(H)$ and ${\rm Range}(\Phi^{L})$ coincide if the matrix $X(\Phi^{M-L})^{T}$ has full row rank, i.e., $\text{Rank\,}(\Phi^{M-L})=s$, which is guaranteed on the condition that $M-L+1\geq s$ and the frequencies in $\mathcal{S}$ are pairwise distinct.

If $\text{Rank\,}(\Phi^{M-L})=s,$ then ${\rm Range}(H)={\rm Range}(\Phi^{L})$.

$\text{Rank\,}(\Phi^{L})=s$ if $L+1\geq s$ and $\omega_{k}\neq\omega_{l},\ \forall k\neq l$.

If $L+1\geq s$, then $s\times s$ square submatrix $\Psi$ of $\Phi^{L}$ is a square Vandermonde matrix whose determinant is given by

Clearly, $\det{(\Psi)}\neq 0$ if and only if $\omega_{i}\neq\omega_{j},i\neq j$. Hence $\text{Rank\,}(\Psi)=s$ which implies $\text{Rank\,}(\Phi^{L})=s$. ∎

Similarly, if $L+1\geq s+1$, the extended matrix $\Phi^{L}_{\omega}=[\Phi^{L}\ \phi^{L}(\omega)]$ has full column rank for any $\omega\notin\mathcal{S}$. As a consequence, $\omega\in\mathcal{S}$ if and only if $\phi^{L}(\omega)$ belongs to ${\rm Range}(\Phi^{L})$.

Appendix B Proof of Theorem 2

The proof of Theorem 2 combines techniques used in and . We take

Graphs of $g(t)$, $|G(\omega)|/L$ and the real part of $G(\omega)/L$ are shown in Figure 8. Function $G$ has the following properties.

$G(-\omega)=e^{-2\pi iL\omega}G(\omega)$ and $|G(-\omega)|=|G(\omega)|$.

$L(\frac{2}{\pi}-\frac{1}{L})\leq G(0)\leq L(\frac{2}{\pi}+\frac{1}{L}).$

$|G(\omega)|\leq\frac{2}{\pi}\frac{L}{|1-4L^{2}\omega^{2}|}+\frac{8}{\pi L}$ for $\omega\in[0,1/2]$.

According to the Poisson summation formula,

In (46) the difference between the discrete and the continuous case lies in

which is bounded above by $8/(\pi L^{2})$, and is therefore negligible when $L$ is sufficiently large.

We start with the following lemma, which paves the way for the proof of Theorem 2.

Suppose objects in $\mathcal{S}$ satisfy the gap condition

It follows from the triangle inequality that

where $\displaystyle\sum_{j=1}^{s}\sum_{l\neq j}|G(\omega_{j}-\omega_{l})\overline{{{\rm c}}_{j}}{{\rm c}}_{l}|$ can be estimated through Property 4 in Lemma 3.

where $\lfloor{s\over 2}\rfloor$ denotes the nearest integer smaller than or equal to $s/2$. Kernel $g$ in (48) is crucial for the convergence of the series in (49).

The equation above along with Property 3 in Lemma 3 yields

The gap condition (47) is derived from the positivity condition of the lower bound, i.e.,

Given that $\mathcal{S}=\{\omega_{1},\ldots,\omega_{s}\}\subset[0,1)$ and frequencies are separated above $1/L$, there are no more than $L$ frequencies in $\mathcal{S}$, i.e., $s<L$.

The lower bound in Theorem 4 follows from Lemma 4 as

The gap condition (16) in Theorem 4 is derived from the positivity condition of the lower bound, i.e.,

We prove the upper bound in Theorem 4 in two cases: $L$ is even or $L$ is odd.

First we substitute $L$ with $2L$ in (48) and obtain

Let $D^{\frac{L}{2}}={\rm diag}(e^{-2\pi i\omega_{1}\frac{L}{2}},e^{-2\pi i\omega_{2}\frac{L}{2}},\ldots,e^{-2\pi i\omega_{s}\frac{L}{2}})$ and $D^{-\frac{L}{2}}=(D^{\frac{L}{2}})^{-1}$. On the one hand,

Eq. (51) above follows from Case 1 as $L+1$ is an even integer. In summary, when $L$ is an odd integer,

Appendix C Proof of Theorems in Section 3

Let $\mathcal{H}=HH^{\star}$, $\mathcal{H}^{\varepsilon}=H^{\varepsilon}{H^{\varepsilon}}^{\star}$ and $\mathcal{E}=HE^{\star}+EH^{\star}+EE^{\star}$. Then

where $\Sigma_{1}={\rm diag}(\sigma_{1},\ldots,\sigma_{s}),\Sigma^{\varepsilon}_{1}={\rm diag}(\sigma^{\varepsilon}_{1},\ldots,\sigma^{\varepsilon}_{s}),$ and $\Sigma^{\varepsilon}_{2}={\rm diag}(\sigma^{\varepsilon}_{s+1},\sigma^{\varepsilon}_{s+2},\ldots).$

On the one hand, the (2,1) entries on both sides of (55) are equal such that

Based on (53), we obtain $U_{2}^{\star}\mathcal{H}=\mathbf{0}$ and therefore

On the other hand, the (1,2) entries on both sides of (55) are equal such that

Based on (53), we obtain $U_{1}^{\star}\mathcal{H}=\Sigma_{1}\Sigma_{1}^{\star}U_{1}^{\star}$ and therefore

Eq. (58) together with (56) and (57) imply

According to Proposition 1, $\sigma^{\varepsilon}_{s}\geq\sigma_{s}-\|E\|_{2}$ and $\sigma^{\varepsilon}_{s+1}\leq\|E\|_{2}$. Then

In particular, while $\omega$ is restricted on $\mathcal{S}$, a sharper upper bound in (22) is derived as follows:

Since $M-L+1\geq s$ and true frequencies are pairwise distinct, $X(\Phi^{M-L})^{T}$ has full row rank. Denote $Y^{\varepsilon}=U^{\varepsilon}_{1}\Sigma^{\varepsilon}_{1}{V^{\varepsilon}_{1}}^{\star}+U^{\varepsilon}_{2}\Sigma^{\varepsilon}_{2}{V^{\varepsilon}_{2}}^{\star}$ where $\Sigma_{1}^{\varepsilon}=\text{diag}(\sigma^{\varepsilon}_{1},\ldots,\sigma^{\varepsilon}_{s})$ and $\Sigma_{2}^{\varepsilon}=\text{diag}(\sigma^{\varepsilon}_{s+1},\sigma^{\varepsilon}_{s+2},\ldots)$. Multiplying ${U^{\varepsilon}_{2}}^{\star}$ on the left and the pseudo-inverse of $X(\Phi^{M-L})^{T}$ on the right of

C.2 Proof of Theorem 4

Both $Q(\omega)$ and $Q^{\varepsilon}(\omega)$ are smooth functions and

Let $D(\omega)=Q^{\varepsilon}(\omega)-Q(\omega)$ and then

First, we derive an upper bound of $|D^{\prime}(\omega)|$ and $|D^{\prime\prime}(\omega)|$ in terms of $\alpha,L$ and $\|E\|_{2}$.

Since $[\phi^{L}(\omega)]^{\prime}=-2\pi i[0\ e^{-2\pi i\omega}\ 2e^{-2\pi i2\omega}\ 3e^{-2\pi i3\omega}\ \ldots\ Le^{-2\pi iL\omega}]^{T}$, we have $\|[\phi^{L}(\omega)]^{\prime}\|_{2}=\eta(L)\sqrt{L+1}$ where $\eta(L)=2\pi\sqrt{1^{2}+2^{2}+\ldots+L^{2}}/\sqrt{L+1}$ and then

Applying the same technique to $D^{\prime\prime}(\omega)$ yields

where $\zeta(L)=(2\pi)^{2}\sqrt{1^{4}+2^{4}+\ldots+L^{4}}/\sqrt{L+1}$.

Next we prove that for each $\omega_{j}\in\mathcal{S}$, there exists a strict local minimizer $\hat{\omega}_{j}$ of $R^{\varepsilon}$ near $\omega_{j}$ satisfying (28).

Since $\omega_{j}$ is a strict local minimizer of $Q(\omega)$ and

We break the following argument into several steps:

Let $Q^{\prime\prime}(\omega_{j})=2m_{j}$. $m_{j}>0$ due to assumption (62). Thanks to the smoothness of $Q^{\prime\prime}$, there exists $\delta_{j}>0$ such that

Consider $Q^{\prime}(\omega_{j}-\delta_{j})$ and $Q^{\prime}(\omega_{j}+\delta_{j})$. There exist $\kappa_{1},\kappa_{2}\in(0,\delta_{j})$ such that

From Step 1, $Q^{\prime\prime}(\omega_{j}-\kappa_{1})>m_{j}$ and $Q^{\prime\prime}(\omega_{j}+\kappa_{1})>m_{j}$, so $Q^{\prime}(\omega_{j}-\delta_{j})<-m_{j}\delta_{j}$ and $Q^{\prime}(\omega_{j}+\delta_{j})>m_{j}\delta_{j}$. According to (60), when noise is sufficiently small such that

$[{Q^{\varepsilon}}(\omega)]^{\prime}$ is a smooth function so there exists $\hat{\omega}_{j}\in(\omega_{j}-\delta_{j},\omega_{j}+\delta_{j})$ such that $[{Q^{\varepsilon}}(\hat{\omega}_{j})]^{\prime}=0$ by intermediate value theorem.

From Step 2 and Step 3 we have obtained an open interval containing $\omega_{j}$: $(\omega_{j}-\delta_{j},\omega_{j}+\delta_{j})$ such that

Also there exists $\hat{\omega}_{j}\in(\omega_{j}-\delta_{j},\omega_{j}+\delta_{j})$ such that

so $\hat{\omega}_{j}$ is a strict local minimizer of $Q^{\varepsilon}(\omega)$ and $R^{\varepsilon}(\omega)$.

through Taylor expansion of $Q^{\prime}(\omega)$ at $\omega=\omega_{j}$. Since $Q^{\prime}(\omega_{j})=0$,

Proof of Theorem 4 ends here. Next we provide the argument for the validation of Remark 9 and Remark 11.

Suppose noise vector $\varepsilon$ contains i.i.d. random variables of variance $\sigma^{2}$. For fixed $M$, $\|E\|_{2}\leq\|E\|_{F}^{2}=\mathcal{O}(\sigma)$. Since $\min_{\xi\in(\omega_{j},\hat{\omega}_{j})}|Q^{\prime\prime}(\xi)|>m_{j}$, we have $|\hat{\omega}_{j}-\omega_{j}|\rightarrow 0$ as $\sigma\rightarrow 0$ and

The asymptotic rate of $\hat{\omega}_{j}\rightarrow\omega_{j}$ as $M=2L\rightarrow\infty$ in the case of $q\geq$ 4 RL is discussed in Remark 11. Here we show that condition (63) and (64) hold as $M=2L\rightarrow\infty$ under assumption (29).

The left hand side (l.h.s.) of (63) scales like $M\sqrt{M\log M}$. On the right hand side (r.h.s.),

Therefore the r.h.s. of (63) grows faster than the l.h.s. and (63) holds as $M=2L\rightarrow\infty$.

Next we show that (64) is also valid as $M\rightarrow\infty$. First

In Step 1, for all $\omega\in(\omega_{j}-\delta_{j},\omega_{j}+\delta_{j})$, $|Q^{\prime\prime}(\omega)-Q^{\prime\prime}(\omega_{j})|$ can be as large as $m_{j}=\mathcal{O}(M^{2})$. Meanwhile $Q^{\prime\prime}(\omega)=Q^{\prime\prime}(\omega_{j})+(\omega-\omega_{j})Q^{\prime\prime\prime}(\kappa)$ for some $\kappa\in(\omega_{j},\omega)$ and then $|Q^{\prime\prime}(\omega)-Q^{\prime\prime}(\omega_{j})|=|\omega-\omega_{j}||Q^{\prime\prime\prime}(\kappa)|$. Since $|Q^{\prime\prime}(\omega)-Q^{\prime\prime}(\omega_{j})|$ can be as large as $\mathcal{O}(M^{2})$ and $Q^{\prime\prime\prime}(\kappa)\leq\mathcal{O}(M^{3})$, $|\omega-\omega_{j}|$ can be as large as $\mathcal{O}(1/M)$. In other words, $\delta_{j}=\mathcal{O}(1/M)$ in Step 1.

In (64), the l.h.s. scales like $\sqrt{M\log M}$ and the r.h.s. $=m_{j}\delta_{j}/2=\mathcal{O}(M^{2}/M)=\mathcal{O}(M)$ as $M=2L\rightarrow\infty$. As a result, (64) holds while $M=2L\rightarrow\infty$ under assumption (29).

Appendix D Proof of Theorem 5

As $|z_{1}+\ldots+z_{R}|^{2}\leq R(|z_{1}|^{2}+\ldots+|z_{R}|^{2})$,

References