Atomic norm denoising with applications to line spectral estimation

Introduction

Extracting the frequencies and relative phases of a superposition of complex exponentials from a small number of noisy time samples is a foundational problem in statistical signal processing. These line spectral estimation problems arise in a variety of applications, including the direction of arrival estimation in radar target identification , sensor array signal processing and imaging systems and also underlies techniques in ultra wideband channel estimation , spectroscopy , molecular dynamics , and power electronics .

While polynomial interpolation using Prony’s technique can estimate the frequency content of a signal exactly from as few as $2k$ samples if there are $k$ frequencies, Prony’s method is inherently unstable due to sensitivity of polynomial root finding. Several methods have been proposed to provide more robust polynomial interpolation (for an extensive bibliography on the subject, see ), and these techniques achieve excellent noise performance in moderate noise. However, the denoising performance is often sensitive to the model order estimated, and theoretical guarantees for these methods are all asymptotic with no finite sample error bounds. Motivated by recent work on atomic norms , we propose a convex relaxation approach to denoise a mixture of complex exponentials, with theoretical guarantees of noise robustness and a better empirical performance than previous subspace based approaches.

Specializing these denoising results to the line spectral estimation, we provide mean-squared-error estimates for denoising line spectra with the atomic norm. The denoising algorithm amounts to soft thresholding the noise corrupted measurements in the atomic norm and we thus refer to the problem as Atomic norm Soft Thresholding (AST). We show, via an appeal to the theory of positive polynomials, that AST can be solved using semidefinite programming (SDP) , and we provide a reasonably fast method for solving this SDP via the Alternating Direction Method of Multipliers (ADMM) . Our ADMM implementation can solve instances with a thousand observations in a few minutes.

While the SDP based AST algorithm can be thought of as solving an infinite dimensional Lasso problem, the computational complexity can be prohibitive for very large instances. To compensate, we show that solving the Lasso problem on an oversampled grid of frequencies approximates the solution of the atomic norm minimization problem to a resolution sufficiently high to guarantee excellent mean-squared error (MSE). The gridded problem reduces to the Lasso, and by leveraging the Fast Fourier Transform (FFT), can be rapidly solved with freely available software such as SpaRSA . A Lasso problem with thousands of observations can be solved in under a second using Matlab on a laptop. The prediction error and the localization accuracy for line spectral estimation both increase as the oversampling factor increases, even if the actual set of frequencies in the line spectral signal are off the Lasso grid.

We compare and contrast our algorithms, AST and Lasso, with classical line spectral algorithms including MUSIC ,and Cadzow’s and Matrix Pencil methods . Our experiments indicate that both AST and the Lasso approximation outperform classical methods in low SNR even when we provide the exact model order to the classical approaches. Moreover, AST has the same complexity as Cadzow’s method, alternating between a least-squares step and an eigenvalue thresholding step. The discretized Lasso-based algorithm has even lower computational complexity, consisting of iterations based upon the FFT and simple linear time soft-thresholding.

where $\conv(\mathcal{A})$ is the convex hull of points in $\mathcal{A}.$ We analyze the denoising performance of an estimate that uses the atomic norm to encourage sparsity in $\mathcal{A}$.

In Section 2, we characterize the performance of the estimate $\hat{x}$ that solves

where $\tau$ is an appropriately chosen regularization parameter. We provide an upper bound on the MSE when the noise statistics are known. Before we state the theorem, we note that the dual norm $\lVert\cdot\rVert_{\mathcal{A}}^{*}$, corresponding to the atomic norm, is given by

This theorem implies that when $\E\lVert w\rVert_{\mathcal{A}}^{*}$ is $o(n)$, the estimate $\hat{x}$ is consistent.

Our lower bound on $\tau$ is in terms of the expected dual norm of the noise process $w$, equal to

That is, the optimal $\tau$ and achievable MSE can be estimated by studying the extremal values of the stochastic process indexed by the atomic set $\mathcal{A}$.

The set $\mathcal{A}$ can be viewed as an infinite dictionary indexed by the continuously varying parameters $f$ and $\phi$. When the number of observations, $n$, is much greater than $k$, $x^{\star}$ is $k$-sparse and thus line spectral estimation in the presence of noise can be thought of as a sparse approximation problem. The regularization parameter for the strongest guarantee in Theorem 1 is given in terms of the expected dual norm of the noise and can be explicitly computed for many noise models. For example, when the noise is Gaussian, we have the following theorem for the MSE:

It is instructive to compare this to the trivial estimator $\hat{x}=y$ which has a per-element MSE of $\sigma^{2}$. In contrast, Theorem 2 guarantees that AST produces a consistent estimate when $k=o(\sqrt{n/\log(n)})$.

The atomic formulation not only offers a way to denoise the line spectral signal, but also provides an efficient frequency localization method. After we obtain the signal estimate $\hat{x}$ by solving (1.1), we also obtain the solution $\hat{z}$ to the dual problem as $\hat{z}=y-\hat{x}$. As we shall see in Corollary 1, the dual solution $\hat{z}$ both certifies the optimality of $\hat{x}$ and reveals the composing atoms of $\hat{x}$. For line spectral estimation, this provides an alternative to polynomial interpolation for localizing the constituent frequencies.

Indeed, when there is no noise, Candés and Fernandez-Granda showed the dual solution recovers these frequencies exactly under mild technical conditions . This frequency localization technique is later extended in to the random undersampling case to yield a compressive sensing scheme that is robust to basis mismatch. When there is noise, numerical simulations show that the atomic norm minimization problem (1.1) gives approximate frequency localization.

A number of Prony-like techniques have been devised that are able to achieve excellent denoising and frequency localization even in the presence of noise. Our experiments in Section 5 demonstrate that our proposed estimation algorithms outperform Matrix Pencil, MUSIC and Cadzow’s methods. Both AST and the discretized Lasso algorithms obtain lower MSE compared to previous approaches, and the discretized algorithm is much faster on large problems.

Abstract Denoising with Atomic Norms

The foundation of our technique consists of extending recent work on atomic norms in linear inverse problems in . In this work, the authors describe how to reconstruct models that can be expressed as sparse linear combinations of atoms from some basic set $\mathcal{A}.$ The set $\mathcal{A}$ can be very general and not assumed to be finite. For example, if the signal is known to be a low rank matrix, $\mathcal{A}$ could be the set of all unit norm rank-$1$ matrices.

We show how to use an atomic norm penalty to denoise a signal known to be a sparse nonnegative combination of atoms from a set $\mathcal{A}$. We compute the mean-squared-error for the estimate we thus obtain and propose an efficient computational method.

The atomic norm $\lVert\cdot\rVert_{\mathcal{A}}$ of $\mathcal{A}$ is the Minkowski functional (or the gauge function) associated with $\conv(\mathcal{A})$ (the convex hull of $\mathcal{A}$):

To set up the atomic norm denoising problem, suppose we observe a signal $y=x^{\star}+w$ and that we know a priori that $x^{\star}$ can be written as a linear combination of a few atoms from $\mathcal{A}$. One way to estimate $x^{\star}$ from these observations would be to search over all short linear combinations from $\mathcal{A}$ to select the one which minimizes $\lVert y-x\rVert_{2}$. However, this could be formidable: even if the set of atoms is a finite collection of vectors, this problem is the NP-hard SPARSEST VECTOR problem .

We now establish some universal properties about the problem (1.1). First, we collect a simple consequence of the optimality conditions in a lemma:

$\hat{x}$ is the solution of (1.1) if and only if (i) $\lVert y-\hat{x}\rVert_{\mathcal{A}}^{*}\leq\tau$, (ii) $\langle y-\hat{x},\hat{x}\rangle=\tau\lVert\hat{x}\rVert_{\mathcal{A}}.$

The supremum in (2.2) is achievable, namely, for any $x$ there is a $z$ that achieves equality. Since $\mathcal{A}$ contains all extremal points of $\{x:\|x\|_{\mathcal{A}}\leq 1\}$, we are guaranteed that the optimal solution will actually lie in the set $\mathcal{A}$:

The dual norm will play a critical role throughout, as our asymptotic error rates will be in terms of the dual atomic norm of noise processes. The dual atomic norm also appears in the dual problem of (1.1)

The dual problem admits a unique solution $\hat{z}$ due to strong concavity of the objective function. The primal solution $\hat{x}$ and the dual solution $\hat{z}$ are specified by the optimality conditions and there is no duality gap: (i) $y=\hat{x}+\hat{z},$ (ii) $\lVert\hat{z}\rVert_{\mathcal{A}}^{*}\leq\tau,$ (iii) $\langle\hat{z},\hat{x}\rangle=\tau\lVert\hat{x}\rVert_{\mathcal{A}}$.

The proofs of Lemma 1 and Lemma 2 are provided in Appendix A. A straightforward corollary of this Lemma is a certificate of the support of the solution to (1.1).

Suppose for some $S\subset\mathcal{A}$, $\hat{z}$ is a solution to the dual problem (2) satisfying

$\langle\hat{z},a\rangle=\tau$ whenever $a\in S,$

$|\langle\hat{z},a\rangle|<\tau$ if $a\not\in S.$

Then, any solution $\hat{x}$ of (1.1) admits a decomposition $\hat{x}=\sum_{a\in S}{c_{a}a}$ with $\lVert\hat{x}\rVert_{\mathcal{A}}=\sum_{a\in S}{c_{a}}.$

Thus the dual solution $\hat{z}$ provides a way to determine a decomposition of $\hat{x}$ into a set of elementary atoms that achieves the atomic norm of $\hat{x}$. In fact, one could evaluate the inner product $\left<\hat{z},a\right>$ and identify the atoms where the absolute value of the inner product is $\tau$. When the SNR is high, we expect that the decomposition identified in this manner should be close to the original decomposition of $x^{\star}$ under certain assumptions.

We are now ready to state a proposition which gives an upper bound on the MSE with the optimal choice of the regularization parameter.

If the regularization parameter $\tau>\lVert w\rVert_{\mathcal{A}}^{*}$, the optimal solution $\hat{x}$ of (1.1) has the MSE

where for (2.7) we have used Lemma 1 and (2.3). The theorem now follows from (2.7) and (2.8) since $\tau>\lVert w\rVert_{\mathcal{A}}^{*}.$ The value of the regularization parameter $\tau$ to ensure the MSE is upper bounded thus, is $\lVert w\rVert_{\mathcal{A}}^{*}.$ ∎

which is the stability result for Lasso reported in assuming no conditions on $\Phi$.

In this section, we provide conditions under which a faster convergence rate can be obtained for AST.

Suppose the set of atoms $\mathcal{A}$ is centrosymmetric and $\lVert w\rVert_{\mathcal{A}}^{*}$ concentrates about its expectation so that $P({\lVert w\rVert_{\mathcal{A}}^{*}}\geq\E{\lVert w\rVert_{\mathcal{A}}^{*}}+t)<\delta(t)$. For $\gamma\in$, define the cone

is strictly positive for some $\gamma>{\E\lVert w\rVert_{\mathcal{A}}^{*}}/{\tau}$. Then

with probability at least $1-\delta(\gamma\tau-\E\lVert w\rVert_{\mathcal{A}}^{*})$.

Since $\lVert w\rVert_{\mathcal{A}}^{*}$ concentrates about its expectation, with probability at least $1-\delta(\gamma\tau-\E\lVert w\rVert_{\mathcal{A}}^{*})$, we have $\lVert w\rVert_{\mathcal{A}}^{*}\leq\gamma\tau$ and hence $\hat{x}-x^{\star}\in C_{\gamma}$. Using (2.6), if $\tau>\lVert w\rVert_{\mathcal{A}}^{*}$,

So, with probability at least $1-\delta(\gamma\tau-\E\lVert w\rVert_{\mathcal{A}}^{*})$:

The main difference between (2.10) and (2.5) is that the MSE is controlled by $\tau^{2}$ rather than $\tau\|x^{*}\|_{\mathcal{A}}$. As we will now see (2.10) provides minimax optimal rates for the examples of sparse vectors and low-rank matrices.

We show in the appendix that in this case $\phi_{\gamma}(x^{\star},\mathcal{A})>\frac{(1-\gamma)}{2\sqrt{k}}$. We also have $\tau_{0}=\E\lVert w\rVert_{\infty}\geq\sigma\sqrt{2\log(n)}.$ Pick $\tau>\gamma^{-1}\tau_{0}$ for some $\gamma<1.$ Then, using our lower bound for $\phi_{\gamma}$ in (2.10), we get a rate of

for the AST estimate with high probability. This bound coincides with the minimax optimal rate derived by Donoho and Johnstone . Note that if we had used (2.5) instead, our MSE would have instead been $O\left(\sqrt{\sigma^{2}k\log n}\|x^{\star}\|_{2}/n\,\right)$, which depends on the norm of the input signal $x^{\star}$.

We prove in the appendix that $\phi_{\gamma}(X^{\star},\mathcal{A})\geq\frac{1-\gamma}{2\sqrt{2r}}$. To obtain an estimate for $\tau$, we note that the spectral norm of the noise matrix satisfies $\|W\|\leq 2\sqrt{n}$ with high probability . Substituting these estimates for $\tau$ and $\phi_{\gamma}$ in (2.10), we get the minimax optimal MSE

2 Expected MSE for Approximated Atomic Norms

We close this section by noting that it may sometimes be easier to solve (1.1) on a different set $\widetilde{\mathcal{A}}$ (say, an $\epsilon$-net of $\mathcal{A}$ instead of $\mathcal{A}$. If for some $M>0,$

holds for every $x$, then Theorem 1 still applies with a constant factor $M$. We will need the following lemma.

${\lVert z\rVert_{\mathcal{A}}^{*}\leq M\lVert z\rVert_{\widetilde{\mathcal{A}}}^{*}}$ for every $z$ iff ${M^{-1}\lVert x\rVert_{\widetilde{\mathcal{A}}}\leq\lVert x\rVert_{\mathcal{A}}}$ for every $z$.

We will show the forward implication – the converse will follow since the dual of the dual norm is again the primal norm. By (2.3), for any $x$, there exists a $z$ with $\lVert z\rVert_{\widetilde{\mathcal{A}}}^{*}\leq 1$ and ${\langle x,z\rangle=\lVert x\rVert_{\widetilde{\mathcal{A}}}}$. So,

Now, we can state the sufficient condition for the following proposition in terms of either the primal or the dual norm:

then under the same conditions as in Theorem 1,

By assumption, $\E\left(\lVert w\rVert_{\mathcal{A}}^{*}\right)\leq\tau$. Now, (2.12) implies $\E\left(\lVert w\rVert_{\widetilde{\mathcal{A}}}^{*}\right)\leq\tau.$ Applying Theorem 1, and using (2.13), we get

Application to Line Spectral Estimation

The infinite set $\mathcal{A}=\{a_{f,\phi}:f\in,\phi\in\}$ forms an appropriate collection of atoms for $x^{\star}$, since $x^{\star}$ in (1.2) can be written as a sparse nonnegative combination of atoms in $\mathcal{A}.$ In fact, $x^{\star}=\sum_{l=1}^{k}c_{l}^{\star}a_{f_{l}^{\star},0}=\sum_{l=1}^{k}|c_{l}^{\star}|a_{f_{l}^{\star},\phi_{l}},$ where $c_{l}^{\star}=|c_{l}^{\star}|e^{i2\pi\phi_{l}}.$

The corresponding dual norm takes an intuitive form:

In other words, $\|v\|_{\mathcal{A}}^{*}$ is the maximum absolute value attained on the unit circle by the polynomial $\zeta\mapsto\sum_{l=0}^{n-1}v_{l}\zeta^{l}$. Thus, in what follows, we will frequently refer to the dual polynomial as the polynomial whose coefficients are given by the dual optimal solution of the AST problem.

In this section, we present a semidefinite characterization of the atomic norm associated with the line spectral atomic set $\mathcal{A}=\{a_{f,\phi}|f\in,\phi\in\}$. This characterization allows us to rewrite (1.1) as an equivalent semidefinite programming problem.

Let $T^{*}$ denote the adjoint of the map $T$. Then we have the following succinct characterization

[35, Theorem 4.24] For any given causal trigonometric polynomial $V(f)=\sum_{l=0}^{n-1}v_{l}e^{-2\pi ilf}$, $|V(f)|\leq\tau$ if and only if there exists complex Hermitian matrix $Q$ such that

Here, ${e}_{1}$ is the first canonical basis vector with a one at the first component and zeros elsewhere and $v^{*}$ denotes the Hermitian adjoint (conjugate transpose) of $v$.

Using Lemma 4, we rewrite the atomic norm $\|x\|_{\mathcal{A}}=\sup_{\|v\|_{\mathcal{A}}^{*}\leq 1}\left<x,v\right>$ as the following semidefinite program:

The dual problem of (3.3) (after a trivial rescaling) is then equal to the atomic norm of $x$:

Therefore, the atomic denoising problem (1.1) for the set of trigonometric atoms is equivalent to

The semidefinite program (3.4) can be solved by off-the-shelf solvers such as SeDuMi and SDPT3 . However, these solvers tend to be slow for large problems. For the interested reader, we provide a reasonably efficient algorithm based upon the Alternating Direction Method of Multipliers (ADMM) in Appendix

2 Choosing the regularization parameter

The choice of the regularization parameter is dictated by the noise model and we show the optimal choice for white gaussian noise samples in our analysis. As noted in Theorem 1, the optimal choice of the regularization parameter depends on the dual norm of the noise. A simple lower bound on the expected dual norm occurs when we consider the maximum value of $n$ uniformly spaced points in the unit circle. Using the result of , the lower bound whenever $n\geq 5$ is

Using a theorem of Bernstein and standard results on the extreme value statistics of Gaussian distribution, we can also obtain a non-asymptotic upper bound on the expected dual norm of noise for $n>3$:

(See Appendix D for a derivation of both the lower and upper bound). If we set the regularization parameter $\tau$ equal to an upper bound on the expected dual atomic norm, i.e.,

an application of Theorem 1 yields the asymptotic result in Theorem 2.

3 Determining the frequencies

As shown in Corollary 1, the dual solution can be used to identify the frequencies of the primal solution. For line spectra, a frequency $f\in$ is in the support of the solution $\hat{x}$ of (1.1) if and only if

That is, $f$ is in the support of $\hat{x}$ if and only if it is a point of maximum modulus for the dual polynomial. Thus, the support may be determined by finding frequencies $f$ where the dual polynomial attains magnitude $\tau$.

Figure 1 shows the dual polynomial for (1.1) with $n=64$ samples and $k=6$ randomly chosen frequencies. The regularization parameter $\tau$ is chosen as described in Section 3.2.

A recent result by Candes and Fernandez-Granda establishes that in the noiseless case, the frequencies localized by the dual polynomial are exact provided the minimum separation between the frequencies is at least $4/n$ where $n$ is the number of samples in the line spectral signal. Under similar separation condition, numerical simulations suggest that (1.1) achieves approximate frequency location in the noisy case.

4 Discretization and Lasso

When the number of samples is larger than a few hundred, the running time of our ADMM method is dominated by the cost of computing eigenvalues and is usually expensive . For very large problems, we now propose using Lasso as an alternative to the semidefinite program (3.4). To proceed, pick a uniform grid of $N$ frequencies and form $\mathcal{A}_{N}=\left\{a_{m/N,\phi}~\middle|~0\leq m\leq N-1\right\}\subset\mathcal{A}$ and solve (1.1) on this grid. i.e., we solve the problem

Because of the relatively simple structure of the atomic set, the optimal solution $\hat{x}$ for (3.6) can be made arbitrarily close to (3.4) by picking $N$ a constant factor larger than $n$. In fact, we show that the atomic norms on $\mathcal{A}$ and $\mathcal{A}_{N}$ are equivalent (See Appendix C):

Using Proposition 3 and (3.5), we conclude

Due to the efficiency of the FFT, the discretized approach has a much lower algorithmic complexity than either Cadzow’s alternating projections method or the ADMM method described in Appendix E, which each require computing an eigenvalue decomposition at each iteration. Indeed, fast solvers for (3.7) converge to an $\epsilon$ optimal solution in no more than $1/\sqrt{\epsilon}$ iterations. Each iteration requires a multiplication by $\Phi$ and a simple “shrinkage” step. Multiplication by $\Phi$ or $\Phi^{*}$ requires $O(N\log N)$ time and the shrinkage operation can be performed in time $O(N)$.

As we discuss below, this fast form of basis pursuit has been proposed by several authors. However, analyzing this method with tools from compressed sensing has proven daunting because the matrix $\Phi$ is nowhere near a restricted isometry. Indeed, as $N$ tends to infinity, the columns become more and more coherent. However, common sense says that a larger grid should give better performance, for both denoising and frequency localization! Indeed, by appealing to the atomic norm framework, we are able to show exactly this point: the larger one makes $N$, the closer one approximates the desired atomic norm soft thresholding problem. Moreover, we do not have to choose $N$ to be too large in order to achieve nearly the same performance as the AST.

Related Work

The classical methods of line spectral estimation, often called linear prediction methods, are built upon the seminal interpolation method of Prony . In the noiseless case, with as little as $n=2k$ measurements, Prony’s technique can identify the frequencies exactly, no matter how close the frequencies are. However, Prony’s technique is known to be sensitive to noise due to instability of polynomial rooting . Following Prony, several methods have been employed to robustify polynomial rooting method including the Matrix Pencil algorithm , which recasts the polynomial rooting as a generalized eigenvalue problem and cleverly uses extra observations to guard against noise. The MUSIC and ESPRIT algorithms exploit the low rank structure of the autocorrelation matrix.

Cadzow proposed a heuristic that improves over MUSIC by exploiting the Toeplitz structure of the matric of moments by alternately projecting between the linear space of Toeplitz matrices and the space of rank $k$ matrices where $k$ is the desired model order. Cadzow’s technique is very similar to a popular technique in time series literature called Singular Spectrum Analysis , which uses autocorrelation matrix instead of the matrix of moments for projection. Both these techniques may be viewed as instances of structured low rank approximation which exploit additional structure beyond low rank structure used in subspace based methods such as MUSIC and ESPRIT. Cadzow’s method has been identified as a fruitful preprocessing step for linear prediction methods . A survey of classical linear prediction methods can be found in and an extensive list of references is given in .

Most, if not all of the linear prediction methods need to estimate the model order by employing some heuristic and the performance of the algorithm is sensitive to the model order. In contrast, our algorithms AST and the Lasso based method, only need a rough estimate of the noise variance. In our experiments, we provide the true model order to Matrix Pencil, MUSIC and Cadzow methods, while we use the estimate of noise variance for AST and Lasso methods, and still compare favorably to the classical line spectral methods.

In contrast to linear prediction methods, a number of authors have suggested using compressive sensing and viewing the frequency estimation as a sparse approximation problem. For instance, notes that the Lasso based method has better empirical localization performance than the popular MUSIC algorithm. However, the theoretical analysis of this phenomenon is complicated because of the need to replace the continuous frequency space by an oversampled frequency grid. Compressive sensing based results (see, for instance, ) need to carefully control the incoherence of their linear maps to apply off-the-shelf tools from compressed sensing. It is important to note that the performance of our algorithm improves as the grid size increases. But this seems to contradict conventional wisdom in compressed sensing because our design matrix $\Phi$ becomes more and more coherent, and limits how fine we can grid for the theoretical guarantees to hold.

We circumvent the problems in the conventional compresssive sensing analysis by directly working in the continuous parameter space and hence step away from such notions as coherence, focussing on the geometry of the atomic set as the critical feature. By showing that the continuous approach is the limiting case of the Lasso based methods using the convergence of the corresponding atomic norms, we justify denoising line spectral signals using Lasso on a large grid. Since the original submission of this manuscript, Candès and Fernandez-Granda showed that our SDP formulation exactly recovers the correct frequencies in the noiseless case.

Experiments

AST needs an estimate of the noise variance $\sigma^{2}$ to pick the regularization parameter according to (3.5). In many situations, this variance is not known to us a priori. However, we can construct a reasonable estimate for $\sigma$ when the phases are uniformly random. It is known that the autocorrelation matrix of a line spectral signal (see, for example Chapter 4 in ) can be written as a sum of a low rank matrix and $\sigma^{2}I$ if we assume that the phases are uniformly random. Since the empirical autocorrelation matrix concentrates around the true expectation, we can estimate the noise variance by averaging a few smallest eigenvalues of the empirical autocorrelation matrix. In the following experiments, we form the empirical autocorrelation matrix using the MATLAB routine corrmtx using a prediction order $m=n/3$ and averaging the lower $25\%$ of the eigenvalues. We used this estimate in equation (3.5) to determine the regularization parameter for both our AST and Lasso experiments.

We implemented Cadzow’s method as described by the pseudocode in , the Matrix Pencil as described in and MUSIC using the MATLAB routine rootmusic. All these algorithms need an estimate of the number of sinusoids. Rather than implementing a heuristic to estimate $k$, we fed the true $k$ to our solvers. This provides a huge advantage to these algorithms. Neither AST or the Lasso based algorithm are provided the true value of $k$, and the noise variance $\sigma^{2}$ required in the regularization parameter is estimated from $y$.

In Figure 2, we show MSE vs SNR plots for a subset of experiments when $n=128$ time samples are taken to take a closer look at the differences. It can be seen from these plots that the performance difference between classical algorithms such as MUSIC and Cadzow with respect to the convex optimization based AST and Lasso is most pronounced at lower sparsity levels. When the noise dominates the signal (SNR $\leq 0$ dB), all the algorithms are comparable. However, AST and Lasso outperform the other algorithms in almost every regime.

We note that the denoising performance of Lasso improves with increased grid size as shown in the MSE vs SNR plot in Figure 3(a). The figure shows that the performance improvement for larger grid sizes is greater at high SNRs. This is because when the noise is small, the discretization error is more dominant and finer gridding helps to reduce this error. Figures 3(a) and (b) also indicate that the benefits of increasing discretization levels are diminishing with the grid sizes, at a higher rate in the low SNR regime, suggesting a tradeoff among grid size, accuracy, and computational complexity.

Finally, in Figure 3(b), we provide numerical evidence supporting the assertion that frequency localization improves with increasing grid size. Lasso identifies more frequencies than the true ones due to basis mismatch. However, these frequencies cluster around the true ones, and more importantly, finer discretization improves clustering, suggesting over-discretization coupled with clustering and peak detection as a means for frequency localization for Lasso. This observation does not contradict the results of where the authors look at the full Fourier basis ($N=n$) and the noise-free case. This is the situation where discretization effect is most prominent. We instead look at the scenario where $N\gg n$.

From the performance profile in Figure 4(a), we see that AST is the best performing algorithm, with Lasso coming in second. Cadzow does not perform as well as AST, even though it is fed the true number of sinusoids. When Cadzow is fed an incorrect $k$, even off by $1$, the performance degrades drastically, and never provides adequate mean-squared error. Figure 4(b) shows that the denoising performance improves with grid size.

Conclusion and Future Work

The Atomic norm formulation of line spectral estimation provides several advantages over prior approaches. By performing the analysis in the continuous domain we were able to derive simple closed form rates using fairly straightforward techniques.We only grid the unit circle at the very end of our analysis and determine the loss incurred from discretization. This approach allowed us to circumvent some of the more complicated theoretical arguments that arise when using concepts from compressed sensing or random matrix theory.

This work provides several interesting possible future directions, both in line spectral estimation and in signal processing in general. We conclude with a short outline of some of the possibilities.

Determining checkable conditions on the cones in Section 2.1 for the atomic norm problem is a major open problem. Our experiments suggest that when the frequencies are spread out, AST performs much better with a slightly larger regularization parameter. This observation was also made in the model-based compressed sensing literature . Moreover, Candès and Fernandez Granda also needed a spread assumption to prove their theories.

This evidence together suggests the fast rate developed in Section 2.1 may be active for signals with well separated frequenies. Determining concrete conditions on the signal $x^{\star}$ that ensure this fast rate require techniques for estimating the parameter $\phi$ in (2.9). Such an investigation should be accompanied by a determination of the minimax rates for line spectral estimation. Such minimax rates would shed further light on the rates achievable for line spectral estimation.

Our work also naturally extends to moment problems where the atomic measures are supported on the unit disk in the complex plane. These problems arise naturally in controls and systems theory and include model order reduction, system identification, and control design. Applying the standard program developed in Section 2 provides a new look at these classic operator theory problems in control theory. It would be of significant importance to develop specialized atomic-norm denoising algorithms for control theoretic problems. Such an approach could yield novel statistical bounds for estimation of rational functions and $\mathcal{H}_{\infty}$-norm approximations.

Our abstract denoising results in Section 2 apply to any atomic models and it is worth investigating their applicability for other models in statistical signal processing. For instance, it might be possible to pose a scheme for denoising a signal corrupted by multipath reflections. Here, the atoms might be all time and frequency shifted versions of some known signal. It remains to be seen what new insights in statistical signal processing can be gleaned from our unified approach to denoising.

Acknowledgements

The authors would like to thank Vivek Goyal, Parikshit Shah, and Joel Tropp for many helpful conversations and suggestions on improving this manuscript. This work was supported in part by NSF Award CCF-1139953 and ONR Award N00014-11-1-0723.

References

Appendix A Optimality Conditions

The function $f(x)=\tfrac{1}{2}\lVert y-x\rVert_{2}^{2}+\tau\lVert x\rVert_{\mathcal{A}}$ is minimized at $\hat{x}$, if for all $\alpha\in(0,1)$ and all $x$,

Since $\lVert\cdot\rVert_{\mathcal{A}}$ is convex, we have

for all $x$ and for all $\alpha\in(0,1)$. Thus, by letting $\alpha\to 0$ in (A.1), we note that $\hat{x}$ minimizes $f(x)$ only if, for all $x$,

However if (A.2) holds, then, for all $x$

implying $f(x)\geq f(\hat{x}).$ Thus, (A.2) is necessary and sufficient for $\hat{x}$ to minimize $f(x)$.

The condition (A.2) simply says that $\tau^{-1}\left(y-\hat{x}\right)$ is in the subgradient of $\lVert\cdot\rVert_{\mathcal{A}}$ at $\hat{x}$ or equivalently that $0\in\partial f(\hat{x})$.

But by definition of the dual atomic norm,

where $I_{A}(\cdot)$ is the convex indicator function. Using this in (A.3), we find that $\hat{x}$ is a minimizer if and only if $\lVert y-\hat{x}\rVert_{\mathcal{A}}^{*}\leq\tau$ and $\langle y-\hat{x},\hat{x}\rangle\geq\tau\lVert\hat{x}\rVert_{\mathcal{A}}$. This proves the theorem. ∎

A.0.2 Proof of Lemma 2

We can rewrite the primal problem (1.1) as a constrained optimization problem:

Now, we can introduce the Lagrangian function

where the first infimum follows by completing the squares and the second infimum follows from (A.4). Thus the dual problem of maximizing $g(z)$ can be written as in (2).

The solution to the dual problem is the unique projection $\hat{z}$ of $y$ on to the closed convex set $C=\{z:\lVert z\rVert_{\mathcal{A}}^{*}\leq\tau\}$. By projection theorem for closed convex sets, $\hat{z}$ is a projection of $y$ onto $C$ if and only if $\hat{z}\in C$ and $\langle z-\hat{z},y-\hat{z}\rangle\leq 0$ for all $z\in C$, or equivalently if $\langle\hat{z},y-\hat{z}\rangle\geq\sup_{z}\,\langle z,y-\hat{z}\rangle=\tau\lVert y-\hat{z}\rVert_{\mathcal{A}}.$ These conditions are satisfied for $\hat{z}=y-\hat{x}$ where $\hat{x}$ minimizes $f(x)$ by Lemma 1. Now the proof follows by the substitution $\hat{z}=y-\hat{x}$ in the previous lemma. The absence of duality gap can be obtained by noting that the primal objective function at $\hat{x},$

Appendix B Fast Rate Calculations

Let $z\in C_{\gamma}(x^{\star},\mathcal{A}).$ For some $\alpha>0$ we have,

In the above inequality, set $z=z_{T}+z_{T^{c}}$ where $z_{T}$ are the components on the support of $T$ and $z_{T^{c}}$ are the components on the complement of $T$. Since $x^{\star}+z_{T}$ and $z_{T^{c}}$ have disjoint supports, we have,

that is, $z$ satisfies the null space property with a constant of $\tfrac{1+\gamma}{1-\gamma}.$ Thus,

Now we can turn to the case of low rank matrices.

and let $\mathcal{P}_{T_{0}}$, $\mathcal{P}_{T}$, and $\mathcal{P}_{T^{\perp}}$ be projection operators that respectively map onto the subspaces $T_{0}$, $T$, and the orthogonal complement of $T$. Now, if $Z\in C_{\gamma}(X^{\star},\mathcal{A})$, then for some $\alpha>0$, we have

Since $\lVert\mathcal{P}_{T_{0}}(Z)\rVert_{*}\leq\lVert\mathcal{P}_{T}(Z)\rVert_{*}$, we have

Putting these computations together gives the estimate

That is, we have $\phi_{\gamma}(X^{\star},\mathcal{A})\geq\frac{1-\gamma}{2\sqrt{2r}}$ as desired. ∎

Appendix C Approximation of the Dual Atomic Norm

This section furnishes the proof that the atomic norms induced by $\mathcal{A}$ and $\mathcal{A}_{N}$ are equivalent. Note that the dual atomic norm of $w$ is given by

i.e., the maximum modulus of the polynomial $W_{n}$ defined by

Treating $W_{n}$ as a function of $f$, with a slight abuse of notation, define

We show that we can approximate the maximum modulus by evaluating $W_{n}$ in a uniform grid of $N$ points on the unit circle. To show that as $N$ becomes large, the approximation is close to the true value, we bound the derivative of $W_{n}$ using Bernstein’s inequality for polynomials.

Let $p_{n}$ be any polynomial of degree $n$ with complex coefficients. Then,

Note that for any $f_{1},f_{2}\in,$ we have

where the last inequality follows by Bernstein’s theorem. Letting $s$ take any of the $N$ values $0,\tfrac{1}{N}\ldots,\tfrac{N-1}{N}$, we see that,

Since the maximum on the grid is a lower bound for maximum modulus of $W_{n}$, we have

Appendix D Dual Atomic Norm Bounds

This section derives non-asymptotic upper and lower bounds on the expected dual norm of gaussian noise vectors, which are asymptotically tight upto $\log\log$ factors. Recall that the dual atomic norm of $w$ is given by $\sqrt{n}\sup_{f\in}|W_{f}|$ where

Thus, the $n$ samples $\left\{W_{m/n}\right\}_{m=0}^{n-1}$ are uncorrelated and thus independent because of their joint gaussianity. This gives a simple non-asymptotic lower bound using the known result for maximum value of $n$ independent gaussian random variables whenever $n>5$:

We will show that the lower bound is asymptotically tight neglecting $\log\log$ terms. Since the dual norm induced by $\mathcal{A}_{N}$ approximates the dual norm induced by $\mathcal{A}$, (See C), it is sufficient to compute an upper bound for $\lVert w\rVert_{\mathcal{A}_{N}}^{*}.$ Note that $|W_{f}|^{2}$ has a chi-square distribution since $W_{f}$ is a Gaussian process. We establish a simple lemma about the maximum of chi-square distributed random variables.

Let $x_{1},\ldots,x_{N}$ be complex gaussians with unit variance. Then,

Let $x_{1},\ldots,x_{N}$ be complex Gaussians with unit variance: $\E[|x_{i}|^{2}]=1$. Note that $2|x_{i}|^{2}$ is a chi-squared random variable with two degrees of freedom. Using Jensen’s inequality, also observe that

Now let $z_{1},\ldots,z_{n}$ be chi-squared random variables with $2$ degrees of freedom. Then we have

Setting $\delta=2\log(N)$ gives $\E\left[\max_{1\leq i\leq N}z_{i}\right]\leq 2\log{N}+2$. Plugging this estimate into (D.1) gives $\E\left[\max_{1\leq i\leq N}|x_{i}|\right]\leq\sqrt{\log{N}+1}$. ∎

Plugging in $N=4\pi n\log(n)$ and using (C.1) and (C.4) establishes a tight upper bound.

Appendix E Alternating Direction Method of Multipliers for AST

A thorough survey of the ADMM algorithm is given in . We only present the details essential to the implementation of atomic norm soft thresholding. To put our problem in an appropriate form for ADMM, rewrite (3.4) as

and dualize the equality constraint via an Augmented Lagrangian:

The updates with respect to $t$, $x$, and $u$ can be computed in closed form:

Here $W$ is the diagonal matrix with entries

The $Z$ update is simply the projection onto the positive definite cone

Projecting a matrix $Q$ onto the positive definite cone is accomplished by forming an eigenvalue decomposition of $Q$ and setting all negative eigenvalues to zero.

To summarize, the update for $(t,u,x)$ requires averaging the diagonals of a matrix (which is equivalent to projecting a matrix onto the space of Toeplitz matrices), and hence operations that are $O(n)$. The update for $Z$ requires projecting onto the positive definite cone and requires $O(n^{3})$ operations. The update for $\Lambda$ is simply addition of symmetric matrices.

Note that the dual solution $\hat{z}$ can be obtained as $\hat{z}=y-\hat{x}$ from the primal solution $\hat{x}$ obtained from ADMM by using Lemma 2.