Near Minimax Line Spectral Estimation
Gongguo Tang, Badri Narayan Bhaskar, Benjamin Recht
Introduction
Spectrum estimation is one of the fundamental problems in statistical signal processing. Despite of hundreds of years of research on this subject, there still remain several fundamental open questions in this area. This paper addresses a central one of these problems: how well can we determine the locations and magnitudes of spectral lines from noisy temporal samples? In this paper, we establish lower bounds on how well we can recover such signals and demonstrate that these worst case bounds can be nearly saturated by solving a convex programming problem. Moreover, we prove that the estimator approximately localizes the frequencies of the true spectral lines.
Suppose the line spectral signal is given by (1.1) and we observe noisy consecutive samples where is i.i.d. complex Gaussian with variance . If the frequencies in satisfy a minimum separation condition
with the distance metric on the torus, then we can determine an estimator satisfying
with high probability by solving a semidefinite programming problem.
Note that if we exactly knew the frequencies , the best rate of estimation we could achieve would be . Our upper bound is merely a logarithmic factor larger than this rate. On the other hand, we will demonstrate via minimax theory that a logarithmic factor is unavoidable when the support is unknown. Hence, our estimator is nearly minimax optimal.
It is instructive to compare our stability rate to the optimal rate achievable for estimating a sparse signal from a finite, discrete dictionary . In the case that there are incoherent dictionary elements, no method can estimate a -sparse signal from measurements corrupted by Gaussian noise at a rate less than . In our problem, there are an infinite number of candidate dictionary elements and it is surprising that we can still achieve such a fast rate of convergence with our highly coherent dictionary. We emphasize that none of the standard techniques from sparse approximation can be immediately generalized to our case. Not only is our dictionary infinite, but also it does not satisfy the usual assumptions such as restricted eigenvalue conditions or coherence conditions that are used to derive stability results in sparse approximation. Nonetheless, in terms of mean-square error performance, our results match those obtained when the frequencies are restricted to lie on a discrete grid.
In the absence of noise, polynomial interpolation can exactly recover a line spectral signal of arbitrary frequencies with as few as equispaced measurements. In the light of our minimum frequency separation requirement (1.2), why should one favor convex techniques for line spectral estimation? Our stability result coupled with minimax optimality establish that no method can perform better than convex methods when the frequencies are well-separated. And, while polynomial interpolation and subspace methods do not impose any resolution limiting assumptions on the constituent frequencies, these methods are empirically highly sensitive to noise. To the best of our knowledge, there is no result similar to Theorem 1 that provides finite sample guarantees about the noise robustness of polynomial interpolation techniques.
Additionally, little is known about how well spectral lines can be localized from noisy observations. The frequencies estimated by any method will never exactly coincide with the true frequencies in the signal in the presence of noise. However, we can characterize the localization performance of our convex programming approach, and summarize this performance in Theorem 2.
Let be the solution to the same semidefinite programming (SDP) problem as referenced in Theorem 1 and . Let and form the decomposition of into coefficients and frequencies, as revealed by the SDP. Then, there exist fixed numerical constants and such that with high probability
.
If for any frequency , the corresponding amplitude , then with high probability there exists a corresponding frequency in the recovered signal such that,
Part (i) of Theorem 2 shows that the estimated amplitudes corresponding to frequencies far from the support are small. In practice, we note that we rarely find any spurious frequencies in the far region, suggesting that our bound (i) is conservative. Parts (ii) and (iii) of the theorem show that in a neighborhood of each true frequency, the recovered signal has amplitude close to the true signal. Part (iv) shows that the larger a particular coefficient is, the better our method is able to estimate the corresponding frequency. In particular, note that if , then . In all four parts, note that the localization error goes to zero as the number of samples grows.
We proceed as follows. In Section 2, we begin by contextualizing our result in the canon of line spectral estimation. We emphasize the advantages and shortcomings of prior art, and describe the methods on which our analysis is built upon. We then in Section 3 describe the semidefinite programming approach to line spectral estimation, originally introduced in , and explain how it relates to other recent spectrum estimation algorithms. We present minimax lower-bounds for line spectral estimation in Section 4. We then provide the proofs of our main results in Section 5. Finally, in Section 6, we empirically demonstrate that the semidefinite programming approach outperforms MUSIC and Cadzow’s technique in terms of the localization metrics defined by parts (i), (ii) and (iii) of Theorem 2.
Prior Art in Line Spectral Estimation
To date, line spectral analysis may be broadly classified into two camps. Subspace methods build upon polynomial interpolation and exploit certain low rank structure in the spectrum estimation problem for denoising. Research on subspace approaches has yielded several standard algorithms that are widely deployed and shown to achieve Cramér-Rao bound asymptotically . However, the sensitivity to noise and model order is not well understood, and there are few guarantees of how these algorithms perform given a limited number of noisy measurements. For a review of many of these classical approaches, see for example .
More recently, approaches based on convex optimization have gained favor and have been demonstrated to perform well on a variety of spectrum estimation tasks . These convex programming methods restrict the frequencies to lie on a finite grid of points and view line spectral signals as a sparse combination of single frequencies. While these methods are reported to have significantly better localization properties than subspace methods (see for example, ) and admit fast and robust algorithms, they have two significant drawbacks. First, while finer gridding may lead to better performance, very fine grids are often numerically unstable. Furthermore, traditional compressed sensing theory does not adequately characterize the performance of fine gridding in these algorithms as the dictionary becomes highly coherent.
Some very recent work bridges the gap between the performant discretized algorithms and continuous subspace approaches by developing a new theory of convex relaxations for infinite continuous dictionary of frequencies. Our work in applies the atomic norm framework proposed by Chandrasekaran et al to the line spectral estimation problem. There, we established stability results on the denoising error and demonstrated empirically that our algorithm compared favorably with both the classical and recent convex approaches which assume the frequencies are on an oversampled DFT grid. Our prior results made no assumption about the separation between frequencies. When the frequencies are well separated, the current work demonstrates that much faster convergence rates are achieved.
Our work is closely related to recent results established by Candès and Fernandez-Granda on exact recovery using convex methods and their recent work on exploiting the robustness of their dual polynomial construction to show super-resolution properties of convex methods. The total variation norm formulation used in is equivalent to the atomic norm specialized to the line spectral estimation problem.
Robustness bounds were established in both our earlier work and in the work of Candès and Fernandez-Granda . In , a slow convergence rate was established with no assumptions about the separation of frequencies in the true signal. In , the authors provide guarantees on the energy of error in the frequency domain in the case that the frequencies are well separated. The noise is assumed to be adversarial with a small spectral energy. In contrast, our paper shows near minimax denoising error under Gaussian noise. It is also not clear that there is a computable formulation for the optimization problem analyzed in . While the guarantees the authors derive in are not comparable with our results, several of their mathematical constructions are used in our proofs here.
Additional recent work derives conditions for approximate support recovery under the Gaussian noise model using the Beurling-Lasso . There, the authors show that there is a true frequency in the neighborhood of every estimated frequency with large enough amplitude. We note that the Beurling-Lasso is equivalent to the atomic norm algorithm that we analyze in this paper. A more recent paper by Fernandez-Granda improves this result by giving conditions on recoverability in terms of the true signal instead of the estimated signal and prove a theorem similar to Theorem 2, but use a worst case bound on the noise samples. Here, we improve these recent results in our proof of Theorem 2, providing tighter guarantees under the Gaussian noise model.
Frequency Localization using Atomic Norms
Denote by the dimensional vector composed of equispaced Nyquist samples for .
for , where is i.i.d. circularly symmetric complex Gaussian noise.
where is the phase of the th component. So, the target signal may be viewed as a sparse non-negative combination of elements from the atomic set given by
For a general atomic set , the atomic norm of a vector is defined as the gauge function associated with the convex hull of atoms:
In this paper, we analyze the performance of the atomic norm soft thresholding (AST) estimate:
where the atomic norm corresponds to the atomic set in (3.3), and is a suitably chosen regularization parameter. The corresponding dual problem is interesting because it gives a way of localizing the frequencies in an atomic norm achieving decomposition of . The dual problem of AST is given by the following semi-infinite program:
to denote the decomposition of given by the dual polynomial .
We show in that a good choice of for obtaining accelerated convergence rates is
for some . We shall use this choice of regularization parameter throughout this paper.
As shown in Section III.A of our prior work, problem (3.5) is equivalent to the semidefinite programming problem
where denotes a Hermitian Toeplitz matrix with as its first row and is the first component of . Similarly, the dual semi-infinite program (3.1) is equivalent to the dual semidefinite program of (3.9).
What is the best rate we can expect?
Using results about minimax achievable rates for linear models , we can deduce that the convergence rate stated in (1.3) is near optimal. Define the set of well separated frequencies as
The expected minimax denoising error for a line spectral signal with frequencies from is defined as the lowest expected denoising error rate for any estimate for the worst case signal with support . Note that we can lower bound by restricting the set of candidate frequencies to smaller set. To that end, suppose we restrict the signal to have frequencies only drawn from an equispaced grid on the torus . Note that any set of frequencies from are pairwise separated by at least . If we denote by a partial DFT matrix with (unnormalized) columns corresponding to frequencies from , we can write for some with . Thus,
Here, the first inequality is the restriction of . The second inequality follows because we project out all components of that do not lie in the span of . Such projections can only reduce the Euclidean norm. The third inequality uses the fact that the minimum singular value of is since . Now we may directly apply the lower bound for estimation error for linear models derived by Candés and Davenport. Namely, Theorem 1 of states that
for some constant that is independent of , , and .
This theorem and Theorem 1 certify that AST is nearly minimax optimal for spectral estimation of well separated frequencies.
Proofs of Main Theorems
We describe the preliminaries and notations, and restate some recent results we used before sketching the proof of Theorems 1 and 2.
The sample may be regarded as the th trigonometric moment of the discrete measure given by (3.1):
for . Thus, the problem of extracting the frequencies and amplitudes from noisy observations may be regarded as the inverse problem of estimating a measure from noisy trigonometric moments.
We can write the vector of observations in terms of an atomic decomposition
or equivalently in terms of a corresponding representing measure given by (3.1) satisfying
There is a one-one correspondence between atomic decompositions and representing measures. Note that there are infinite atomic decompositions of and also infinite corresponding representing measures. However, since every collection of atoms is linearly independent, forms a full spark frame and therefore the problem of finding the sparsest decomposition of is well-posed if there is a decomposition which is at least sparse.
2 Dual Certificate and Exact Recovery
whenever .
The authors of not only explicitly constructed such a certificate characterized by the dual polynomial , but also showed that their construction satisfies some stability conditions, which is crucial for showing that denoising using the atomic norm provides stable recovery in the presence of noise.
For each , interpolates the sign vector so that
In each neighborhood corresponding to defined by , the polynomial behaves like a quadratic and there exist constants so that
When , there is a numerical constant such that
We use results in and (reproduced in Appendix D for convenience) and borrow several ideas from the proofs in , with nontrivial modifications to establish the error rate of atomic norm regularization.
3 Proof of Theorem 1
Let be the representing measure for the solution of (3.5) with minimum total variation norm, that is,
Using a Taylor series approximation in each of the near regions , we first show that the denoising error (or in general any integral of a trigonometric polynomial against the difference measure) can be controlled in terms of an integral in the far region and the zeroth, first, and second moments of the difference measure in the near regions. The precise result is presented in the following lemma, whose proof is given in Appendix A.
Then for any th order trigonometric polynomial , we have
Applying Lemma 1 to the error function, we get
As a consequence of our choice of in (3.8), we can show that with high probability. In fact, we have
The second inequality follows from the optimality conditions for (3.5). It is shown in Appendix C of that the penultimate inequality holds with high probability.
Therefore, to complete the proof, it suffices to show that the other terms on the right hand side of (5.3) are . While there is no exact frequency recovery in the presence of noise, we can hope to get the frequencies approximately right. Hence, we expect that the integral in the far region can be well controlled and the local integrals of the difference measure in the near regions are also small due to cancellations. Next, we utilize the properties of the dual polynomial in Theorems 4 and another polynomial given in Theorem 5 in Appendix B to show that the zeroth and first moments of may be controlled in terms of the other two quantities in (5.3) to upper bound the error rate. The following lemma is similar to Lemmas 2.2 and 2.3 in , but we have made several modifications to adapt it to our signal and noise model. For completeness, we provide the proof in Appendix C.
There exists numeric constants and such that
Let . If is large enough, then there exists a numerical constant such that, with high probability
Putting together Lemmas 1, 2 and 3, we finally prove our main theorem:
The first three inequalities come from successive applications of Lemmas 1, 2 and 3 respectively. The fourth inequality follows from (5.4) and the fifth by our choice of according to Eq. (3.8). This completes the proof of Theorem 1.
4 Proof of Theorem 2
The first two statements in Theorem 2 are direct consequences of Lemma 3. For (iii.), we follow and use the dual polynomial constructed in Lemma 2.2 of which satisfies
We note that . Then, by applying triangle inequality several times,
We upper bound the first term using Lemma 5 in Appendix B which yields
The other terms can be controlled using the properties of :
Using Lemma 3, both of the above are upper bounded by . Now, by combining these upper bounds, we finally have
This shows part (iii) of the theorem. Part (iv) can be obtained by combining parts (ii) and (iii).
Experiments
In , we demonstrated with extensive experiments that AST outperforms classical subspace algorithms in terms of mean squared estimation error. In the experiments here, we focus on frequency localization and compare the performance of AST, MUSIC and Cadzow’s method under various choices of number of frequencies, number of samples and signal to noise ratios (SNRs).
AST needs an estimate of the noise variance to pick the regularization parameter according to (3.8). In our experiments, we do not provide our algorithm with the true noise variance. Instead, we can construct an estimate for with the following heuristic. We formed the empirical autocorrelation matrix using the MATLAB routine corrmtx using a prediction order and averaging the lower of the eigenvalues. We then use this estimate in equation (3.8) to determine the regularization parameter. See for more details.
We implemented AST using the Alternating Direction Method of Multipliers (ADMM, see for example, , or for the specific details). We used the stopping criteria described in and set for all experiments. We use the dual solution to determine the support of the optimal solution . Once the frequencies are extracted, we ran the least squares problem where to obtain debiased estimates of the amplitudes.
We implemented Cadzow’s method as described by the pseudocode in , and MUSIC using the MATLAB routine rootmusic. These algorithms need an estimate of the number of sinusoids. Rather than implementing a heuristic to estimate , we fed the true to our solvers. This provides a significant advantage to these algorithms. On the contrary, AST is not provided the true value of , and the noise variance required in the regularization parameter is estimated from .
Let and denote the amplitudes and frequencies estimated by any of the algorithms - AST, MUSIC or Cadzow. We use the following error metrics to characterize the frequency localization of various algorithms:
Sum of the absolute value of amplitudes in the far region ,
The weighted frequency localization error,
Error in approximation of amplitudes in the near region,
These are precisely the quantities that we prove tend to zero in Theorem 2.
To summarize the results, we first provide performance profiles to summarize the behavior of the various algorithms across all of the parameter settings. Performance profiles provide a good visual indicator of the relative performance of many algorithms under a variety of experimental conditions . Let be the set of experiments and let be the value of the error measure of experiment using the algorithm . Then the ordinate of the graph at specifies the fraction of experiments where the ratio of the performance of the algorithm to the minimum error across all algorithms for the given experiment is less than , i.e.,
The performance profiles in Figure 1 show that AST is the best performing algorithm for all the three metrics. AST in fact outperforms MUSIC and Cadzow by a substantial margin for metrics and .
In Figure 2, we display how the error metrics vary with increasing SNR for AST, MUSIC and Cadzow. We restrict these plots to the experiments with samples. These plots demonstrate that AST localizes frequencies substantially better than MUSIC and Cadzow even for low signal to noise ratios as there is very little energy in the far region of the frequencies () and has the smallest weighted mean square frequency deviation (). Although we have plotted the average value in these plots, we observed spikes in the plots for Cadzow’s algorithm as the average is dominated by the worst performing instances. These large errors are due to the numerical instability of polynomial root finding.
Conclusion and Future Work
In this paper, we demonstrated stability of atomic norm regularization by analysis of specific properties of the atomic set of moments and the associated dual space of trigonometric polynomials. The key to our analysis is the existence and properties of various trigonometric polynomials associated with signals with well separated frequencies.
Though we have made significant progress at understanding the theoretical limits of line-spectral estimation and superresolution, our bounds could still be improved. For instance, it remains open as to whether the logarithmic term in Theorem 1 can be improved to . Deriving such an upper bound or improving our minimax lower bound would provide an interesting direction for future work.
Additionally, it is not clear if our localization bounds in Theorem 2 have the optimal dependence on the number of sinusoids . For instance, we expect that the condition on signal amplitudes for approximate support recovery should not depend on , by comparison with similar guarantees that have been established for Lasso . We additionally conjecture that for a large enough regularization parameter, there will be no spurious recovered frequencies in the solution. That is, there should be no non-zero coefficients in the “far region” in Theorem 2. Future work should investigate whether better guarantees on frequency localization are possible.
References
Appendix A Proof of Lemma 1
We first split the domain of integration into the near and far regions.
by using Hölder’s inequality for the last inequality. Using Taylor’s theorem, we may expand the integrand around as
where for the last inequality we have used a theorem of Bernstein for trigonometric polynomials (see, for example ):
Substituting back into (A.1) yields the desired result.
Appendix B Some useful lemmas
In addition to Theorem 4, we recall another result in where the authors show the existence of a trigonometric polynomial that is linear in each which is also an essential ingredient in our proof.
For every there exists a numerical constant such that
For , there exists a numerical constant such that
We will also need the following straightforward consequence of the constructions of the polynomials in Theorem 4, Theorem 5, and Section 5.4.
There exists a numerical constant such that the constructed in Theorem 4, in Theorem 5, and in Section 5.4 satisfy respectively
We will give a detailed proof of (B.3), and list the necessary modifications for proving (B.4) and (B.5). The dual polynomial constructed in is of the form
where is the squared Fejér kernel (recall that )
for some numerical constants and . Using (B.6) and triangle inequality, we bound as follows:
To continue, note that where is the Fejér kernel, since is the squared Fejér kernel. We can write
where . Now, by using Parseval’s identity, we obtain
for some numerical constant when .
Now let us turn our attention to . Since , we have
We have already established that and we will now show that . Differentiating the expression for in (B.9), we get
Therefore, by applying Parseval’s identity again, we get
for some constant . Combining (B.12) and (B.10) with (B.8) gives the desired result in (B.3).
The dual polynomial is also of the form (B.6) with coefficient vectors and , which satisfy [2, Proof of Lemma 2.7]
Combining the above two bounds with (B.7), (B.12) and (B.10) gives the desired result in (B.4).
The last polynomial also has the form (B.6) with coefficient vectors and . According to [22, Proof of Lemma 2.2], these coefficients satisfy
which yields (B.5) following the same argument leading to (B.3).
Using Lemma 4, we can derive the estimates we need in the following lemma.
Let be the difference measure. Then, there exists numerical constant such that
Here we use Parseval’s identity in the second to last step and Hölder’s inequality in the last inequality. Then, the result follows by using Lemma 4 and (5.4). ∎
We also need the following consequence of the optimality condition of AST from [7, Lemma 2]:
Appendix C Proof of Lemma 2
Set and let be the dual polynomial promised by Theorem 4 for this . Then, we have
We use a similar argument for bounding but this time use the dual polynomial guaranteed by Theorem 5. Again, start with the polar form
Set in Theorem 5 to obtain
For the first inequality, we have used (B.1) and triangle inequality, and for the last inequality, we have used (B.14) and (B.2). Equations (C.1) and (C.2) complete the proof.
Appendix D Proof of Lemma 3
Denote by the projection of the difference measure on the support set of so that is supported on . Then, setting the polynomial in Theorem 4 that interpolates the sign of , we have
where for the first inequality we used triangle inequality and for the last inequality we used (B.13). The integration over is can be bounded using Hölder’s inequality
Now, we appeal to Proposition 1 and obtain
with high probability, where for the penultimate inequality we used our choice of and with high probability, a fact shown in Appendix C of . Substituting (D) in (D.2), we get
As a consequence of (D.1) and (D.5), we get,
whence the result follows for large enough