---
abstract: 'The mathematical theory of super-resolution developed recently by Candès and Fernandes-Granda states that a continuous, sparse frequency spectrum can be recovered with infinite precision via a (convex) atomic norm technique given a set of uniform time-space samples. This theory was then extended to the cases of partial/compressive samples and/or multiple measurement vectors via atomic norm minimization (ANM), known as off-grid/continuous compressed sensing (CCS). However, a major problem of existing atomic norm methods is that the frequencies can be recovered only if they are sufficiently separated, prohibiting commonly known high resolution. In this paper, a novel (nonconvex) sparse metric is proposed that promotes sparsity to a greater extent than the atomic norm. Using this metric an optimization problem is formulated and a locally convergent iterative algorithm is implemented. The algorithm iteratively carries out ANM with a sound reweighting strategy which enhances sparsity and resolution, and is termed as reweighted atomic-norm minimization (RAM). Extensive numerical simulations are carried out to demonstrate the advantageous performance of RAM with application to direction of arrival (DOA) estimation.'
author:
- 'Zai Yang, [*Member, IEEE*]{}, and Lihua Xie, [*Fellow, IEEE*]{} [^1]'
title: Enhancing Sparsity and Resolution via Reweighted Atomic Norm Minimization
---

Continuous compressed sensing (CCS), DOA estimation, frequency estimation, gridless sparse method, high resolution, reweighted atomic norm minimization (RAM).

Introduction
============

Compressed sensing (CS) [@candes2006compressive; @donoho2006compressed] refers to a technique of reconstructing a high dimensional signal from far fewer samples and has brought significant impact on signal processing and information theory in the past decade. In conventional wisdom, the signal of interest needs to be sparse under a finite discrete dictionary for successful reconstruction, which limits its applications, for example, to array processing, radar and sonar, where the dictionary is typically specified by one or more continuous parameters. In this paper, we are concerned about a compressed sensing problem with a continuous dictionary which arises in line spectral estimation and array processing [@krim1996two; @stoica2005spectral]. In particular, we are interested in recovering $L$ discrete sinusoidal signals which compose the data matrix $\m{Y}^o\in\bC^{N\times L}$ with its $(j,t)$th element (corrupted by noise in practice) where $i=\sqrt{-1}$, $f_k\in\bT\triangleq\left[0,1\right]$, $s_{kt}\in\bC$ and $\mbra{N}=\lbra{1,2,\dots,N}$. This means that each column of $\m{Y}^o$ is superimposed by $K$ discrete sinusoids with frequencies $\lbra{f_k}$ and amplitudes $\lbra{s_{kt}}$. To recover $\m{Y}^o$ (and the frequencies in many applications), however, we are only given partial/compressive samples on its rows indexed by $\m{\Omega} \subset \mbra{N}$ (of size $M < N$), denoted by $\m{Y}_{\m{\Omega}}^o$. This problem is referred to as off-grid or continuous compressed sensing (CCS) according to [@tang2012compressed; @yang2014continuous] differing from the existing CS framework in the sense that every frequency $f_k$ can take any continuous value in $\bT$ rather than constrained on a finite discrete grid.

The CCS problem in the case of $L=1$ (a.k.a. the single-measurement-vector (SMV) case) is usually known as line spectral estimation in which frequency recovery though is of main interest. The use of compressive data can lead to efficient sampling and/or energy saving. It also can be caused by data missing due to adversary environmental effects. The multiple-measurement-vector (MMV) case with $L>1$ is common in array processing where one estimates directions of a few narrowband sources using outputs of an antenna array. Readers are referred to [@krim1996two; @stoica2005spectral] for derivation of the model in . Therein $\m{Y}^o$ consists of outputs of a virtual $N$-element uniform linear array (ULA), in which adjacent antennas are spaced by half a wavelength, over $L$ time snapshots. In particular, each column of $\m{Y}^o$ corresponds to one snapshot of the ULA and each row consists of outputs of a single antenna. The fact that we have only access to $\m{Y}_{\m{\Omega}}^o$ means that we actually use a sparse linear array (SLA) that is obtained by retaining the antennas of the ULA indexed by $\m{\Omega}$. Therefore, the index set $\m{\Omega}$ refers to geometry of the SLA and a smaller $M$ means use of fewer antennas (note that SLAs are common in practice for obtaining a large aperture from a limited number of antennas, see, e.g., [@linebarger1993difference] and the references therein). Each frequency component corresponds to one source. The value of $f_k$ uniquely determines the direction of source $k$, and vice versa. Consequently, the problem of direction of arrival (DOA) estimation using a SLA $\m{\Omega}$ is exactly the frequency estimation problem in CCS given the measurement matrix $\m{Y}_{\m{\Omega}}^o$.

Due to its connections to line spectral estimation and DOA estimation, studies of the CCS problem have a long history while frequency estimation has been mainly focused on. Well known conventional methods include periodogram (or beamforming), Capon’s beamforming and subspace methods like MUSIC (see the review in [@stoica2005spectral]). Periodogram suffers from the so-called leakage problem and the Fourier resolution limit of $\frac{1}{N}$ even in the full data case when $M=N$ [@stoica2005spectral]. It therefore has difficulties in resolving two closely spaced frequencies. The situation becomes even worse in the compressive data case. Capon’s beamforming and MUSIC are *high resolution* methods in the sense that they can break the aforementioned resolution limit. Since they are covariance-based methods sufficient snapshots are required to estimate the data covariance. Moreover, they are sensitive to source correlations. With the development of sparse signal representation and later the CS concept, sparse methods have been popular in the last decade which exploit the prior knowledge that the number of frequency components $K$ is small [@malioutov2005sparse; @hyder2010direction; @stoica2011new; @stoica2011spice; @wei2012doa; @fannjiang2012coherence; @hu2013doa; @duarte2013spectral; @liu2013sparsity; @hu2012compressed; @yang2012robustly; @yang2013off; @hu2013fast; @austin2013dynamic; @tan2014joint]. In these methods, however, the frequency domain $\bT$ has to be gridded/discretized into a finite set, resulting in the grid mismatch problem that limits the estimation accuracy as well as brings challenges to the theoretical performance analysis [@chi2011sensitivity; @yang2013off]. Though modified, off-grid estimation methods [@hu2012compressed; @yang2012robustly; @yang2013off; @hu2013fast; @austin2013dynamic; @tan2014joint] have been implemented to alleviate these drawbacks, overall they are still based on gridding of the frequency domain.

A mathematical theory of super-resolution was recently introduced by Candès and Fernandes-Granda [@candes2013towards]. They studied frequency estimation in the SMV, full data case and proposed a *gridless* convex optimization method based on the atomic norm (or the total variation norm) [@chandrasekaran2012convex]. In addition, they proved that the frequencies can be recovered with infinite precision in the absence of noise once they are mutually separated by at least $\frac{4}{N}$. This theory was then extended to the cases of compressive data and MMVs by Tang *et al.* [@tang2012compressed] and the authors [@yang2014continuous; @yang2014exact], showing that the signal and the frequencies can be exactly recovered with high probability via atomic norm minimization (ANM) provided $M\geq O\sbra{K\ln K\ln N}$ and the same frequency separation condition holds. Other related papers include [@bhaskar2013atomic; @candes2013super; @tang2015near; @yang2014discretization; @yang2015gridless; @azais2014spike; @condat2015cadzow; @chen2014robust; @tan2014direction; @mishra2014off; @lu2015distributed]. While the atomic norm techniques completely eliminate grid mismatches of earlier grid-based sparse methods, a major problem is that the frequencies have to be sufficiently separated for successful recovery, prohibiting *high resolution*.[^2]

In this paper, we propose a high resolution gridless sparse method for signal and frequency recovery in CCS. Our method is motivated by the formulations and properties of the atomic $\ell_0$ norm and the atomic norm in [@yang2014continuous; @yang2014exact]. In particular, the atomic $\ell_0$ norm directly exploits sparsity and has no resolution limit but is NP hard to compute. To the contrary, as a convex relaxation the atomic norm can be efficiently computed but suffers from a resolution limit as mentioned above. We propose a novel sparse metric and theoretically show that the new metric fills the gap between the atomic $\ell_0$ norm and the atomic norm. It approaches the former under appropriate parameter setting and breaks the resolution limit. Using this sparse metric we formulate a nonconvex optimization problem for signal and frequency recovery. A locally convergent iterative algorithm is presented to solve the problem. Some further analysis shows that the algorithm iteratively carries out ANM with a sound reweighting strategy that determines preference of frequency selection based on the latest estimate and enhances sparsity and resolution. The resulting algorithm is termed as reweighted atomic-norm minimization (RAM). Extensive numerical simulations are carried out to demonstrate the performance of RAM with application to DOA estimation compared to existing art.

We note that the idea of reweighted optimization for enhancing sparsity is not new. For example, reweighted $\ell_1$ algorithms have been introduced for discrete CS [@lobo2007portfolio; @candes2008enhancing; @wipf2010iterative; @stoica2014weighted], and reweighted trace minimization for low rank matrix recovery (LRMR) [@fazel2003log; @mohan2012iterative]. However, it is unclear how to implement a reweighting strategy in the continuous dictionary setting until this paper. Furthermore, besides sparsity we show that the proposed reweighted algorithm enhances resolution that is of great importance in CCS.

Notations used in this paper are as follows. $\bR$ and $\bC$ denote the sets of real and complex numbers respectively. $\bT$ denotes the unit circle $\left[0,1\right]$ by identifying the beginning and the ending points. Boldface letters are reserved for vectors and matrices. For an integer $N$, $[N]\triangleq\lbra{1,\cdots,N}$. $\abs{\cdot}$ denotes cardinality of a set, amplitude of a scalar, or determinant of a squared matrix. $\onen{\cdot}$, $\twon{\cdot}$ and $\frobn{\cdot}$ denote the $\ell_1$, $\ell_2$ and Frobenius norms respectively. $\m{A}^T$ and $\m{A}^H$ are the matrix transpose and conjugate transpose of $\m{A}$ respectively. $x_j$ is the $j$th entry of a vector $\m{x}$. Unless otherwise stated, $\m{x}_{\m{\Omega}}$ and $\m{A}_{\m{\Omega}}$ respectively reserve the entries of $\m{x}$ and the rows of $\m{A}$ indexed by a set $\m{\Omega}$. For a vector $\m{x}$, $\diag\sbra{\m{x}}$ is a diagonal matrix with $\m{x}$ being its diagonal. $\rank\sbra{\cdot}$ denotes the rank and $\tr\sbra{\cdot}$ the trace. $\m{A}\geq\m{0}$ means that $\m{A}$ is positive semidefinite (PSD).

The rest of the paper is organized as follows. Section \[sec:preliminary\] revisits preliminary gridless sparse methods that motivate this paper. Section \[sec:novelmetric\] presents a novel sparse metric for signal and frequency recovery. Section \[sec:RAM\] introduces the RAM algorithm. Section \[sec:implementation\] presents some algorithm implementation strategies for accuracy and speed considerations. Section \[sec:simulation\] provides extensive numerical simulations to demonstrate the performance of RAM. Section \[sec:conclusion\] concludes this paper.

Preliminary Gridless Sparse Methods by Exploiting Sparsity {#sec:preliminary}
==========================================================

Unless otherwise stated, we assume in this paper that the observed data $\m{Y}_{\m{\Omega}}^o$ is contaminated by noise whose Frobenius norm is bounded by $\eta\geq0$. It is clear that $\eta=0$ refers to the noiseless case. The CCS problem is solved by exploiting sparsity in the sense that the number of frequency components $K$ is small. In particular, we seek a sparse candidate $\m{Y}$ that is composed of a few frequency components and is meanwhile consistent with the observed data by imposing that $\m{Y}\in\cS$, where Therefore, we first define a sparse metric of $\m{Y}$ and then optimize the metric over $\cS$ for its solution. The frequencies $\lbra{f_k}$ are estimated using the frequency components composing $\m{Y}$.

A direct sparse metric is the smallest number of frequency components composing $\m{Y}$, known as the atomic $\ell_0$ norm and denoted by $\norm{\m{Y}}_{\cA,0}$ [@tang2012compressed; @yang2014continuous; @yang2014exact]: where $\m{a}\sbra{f}=\mbra{1,e^{i2\pi f},\dots,e^{i2\pi\sbra{N-1}f}}^T\in\bC^N$ denotes a discrete sinusoid with frequency $f\in\bT$ and $\m{s}_k\in\bC^{L\times 1}$ is the coefficient vector of the $k$th sinusoid. Following from [@tang2012compressed; @yang2014continuous; @yang2014exact], $\norm{\m{Y}}_{\cA,0}$ can be characterized as the following rank minimization problem: Throughout this paper we use the following identity whenever $\m{R}\in \bC^{N\times N}$ is positive semidefinite: The first constraint in (\[formu:atom0norm\]) imposes that $\m{Y}$ lies in the range space of a (Hermitian) Toeplitz matrix where $u_j$ is the $j$th entry of $\m{u}\in\bC^N$. The frequencies composing $\m{Y}$ are encoded in $T\sbra{\m{u}}$. Once an optimizer of $\m{u}$, say $\m{u}^*$, is obtained the frequencies can be retrieved from $T\sbra{\m{u}^*}$ using the Vandermonde decomposition lemma (see, e.g., [@stoica2005spectral]), which states that any PSD Toeplitz matrix $T\sbra{\m{u}^*}$ can be decomposed as where the order $K^*=\rank\sbra{T\sbra{\m{u}^*}}$ and $p_k^*>0$ (note that this decomposition is unique if $K^*<N$ and a computational method can be found in [@yang2015gridless Appendix A]). Therefore, by (\[formu:atom0norm\]) the CCS problem is reformulated as a LRMR problem in which the matrix $T\sbra{\m{u}}$ is Toeplitz and PSD and its range space contains $\m{Y}$.

The atomic $\ell_0$ norm exploits sparsity to the greatest extent possible; however, it is nonconvex and NP-hard to compute according to the rank minimization formulation and it thus encourages computationally feasible alternatives. In this spirit, the atomic ($\ell_1$) norm, denoted by $\atomn{\m{Y}}$, is introduced as a convex relaxation of $\norm{\m{Y}}_{\cA,0}$ [@tang2012compressed; @yang2014continuous; @yang2014exact]: which is a continuous counterpart of the $\ell_{2,1}$ norm utilized for joint sparse recovery in discrete CS (see, e.g., [@malioutov2005sparse; @van2010theoretical]). $\atomn{\m{Y}}$ is a norm and has the following semidefinite formulation [@tang2012compressed; @yang2014continuous; @yang2014exact]: From the perspective of LRMR, (\[formu:AN\_SDP\]) attempts to recover the low rank matrix $T\sbra{\m{u}}$ by relaxing the pseudo rank norm in (\[formu:atom0norm\]) to the nuclear norm (or the trace norm for a PSD matrix). Again, the frequencies are encoded in $T\sbra{\m{u}}$ and can be obtained using the Vandermonde decomposition once the optimization problem is solved within a polynomial time. The atomic norm is computationally advantageous compared to the atomic $\ell_0$ norm while it suffers from a resolution limit due to the relaxation which is not shared by the latter [@candes2013towards; @tang2012compressed; @yang2014exact].

Enhancing Sparsity and Resolution via A Novel Sparse Metric {#sec:novelmetric}
===========================================================

Inspired by the link between CCS and LRMR demonstrated above, we propose the following sparse metric of $\m{Y}$: where $\epsilon>0$ is a regularization parameter that avoids the first term being $-\infty$ when $T\sbra{\m{u}}$ is rank deficient. Note that the log-det heuristic $\log\abs{\cdot}$ has been widely used as a smooth surrogate for the rank of a PSD matrix (see, e.g., [@david1994algorithms; @fazel2003log; @mohan2012iterative]). Also, a similar logarithmic penalty has been adopted to approximate the pseudo $\ell_0$ norm for discrete sparse recovery [@rao1999affine; @chartrand2008iteratively; @candes2008enhancing]. From the perspective of LRMR, the atomic $\ell_0$ norm minimizes the number of nonzero eigenvalues of $T\sbra{\m{u}}$ while the atomic norm minimizes the sum of the eigenvalues. In contrast, the new metric $\cM^{\epsilon}\sbra{\m{Y}}$ puts penalty on $\sum_{k=1}^N\ln\abs{\lambda_k+\epsilon}$, where $\lbra{\lambda_k}_{k=1}^N$ denotes the eigenvalues. We plot the function $h(\lambda)=\ln\abs{\lambda+\epsilon}$ with different $\epsilon$’s in Fig. \[Fig:sparsity\] together with the constant function (except at $\lambda=0$) and the identity function corresponding to the $\ell_0$ and $\ell_1$ norms respectively, where $h(\lambda)$ is translated and scaled for better illustration without altering its sparsity-enhancing property. Intuitively, $h(\lambda)$ gets close to the $\ell_1$ norm for large $\epsilon$ while it approaches the $\ell_0$ norm as $\epsilon\rightarrow0$. Therefore, we expect that the new metric $\cM^{\epsilon}\sbra{\m{Y}}$ bridges $\atomn{\m{Y}}$ and $\norm{\m{Y}}_{\cA,0}$ when $\epsilon$ varies from $+\infty$ to $0$. Formally, we have the following results.

![Illustration of the sparsity-promoting property of $\cM^{\epsilon}\sbra{\cdot}$ with respect to $\epsilon$. The plotted curves include the $\ell_0$ and $\ell_1$ norms corresponding to $\norm{\cdot}_{\cA,0}$ and $\norm{\cdot}_{\cA}$ respectively, and $\ln\abs{\lambda+\epsilon}$ corresponding to $\cM^{\epsilon}\sbra{\cdot}$ with $\epsilon=10, 1, 0.1, 10^{-3}$ and $10^{-6}$. $\ln\abs{\lambda+\epsilon}$ is translated and scaled such that it equals 0 and 1 at $\lambda=0$ and 1 respectively for better illustration.[]{data-label="Fig:sparsity"}](Fig_sparsity.pdf){width="3in"}

Let $\epsilon\rightarrow+\infty$. Then, i.e., they are equivalent infinitesimals. \[thm:epsilontoinf\]

See Appendix \[sec:append:epsilontoinf\].

Let $r=\norm{\m{Y}}_{\cA,0}$ and $\epsilon\rightarrow0$. Then, we have the following results:

1.  If $r\leq N-1$, then i.e., they are equivalent infinities. Otherwise, $\cM^{\epsilon}\sbra{\m{Y}}$ is a positive constant depending only on $\m{Y}$;

2.  Let $\m{u}_{\epsilon}^*$ be the (global) optimizer of $\m{u}$ to the optimization problem in (\[formu:nonconvexrelax\]). Then, the smallest $N-r$ eigenvalues of $T\sbra{\m{u}_{\epsilon}^*}$ are either zero or approach zero as fast as $\epsilon$;

3.  For any cluster point of $\m{u}_{\epsilon}^*$ at $\epsilon=0$, denoted by $\m{u}_0^*$, there exists an atomic decomposition $\m{Y}=\sum_{k=1}^r\m{a}\sbra{f_k}\m{s}_k$ such that $T\sbra{\m{u}_0^*}=\sum_{k=1}^r\twon{\m{s}_k}^2\m{a}\sbra{f_k}\m{a}\sbra{f_k}^H$.[^3]

\[thm:epsilontozero\]

See Appendix \[sec:append:epsilontozero\].

Note that the term $\ln\frac{1}{\epsilon}$ in that becomes unbounded as $\epsilon\rightarrow0$ is not problematic in the optimization problem to introduce, since the objective function $\cM^{\epsilon}\sbra{\m{Y}}$ can be re-scaled by $\sbra{\ln\frac{1}{\epsilon}}^{-1}$ for any $\epsilon>0$ without altering the optimizer. By similar arguments we see that $\epsilon^{-\frac{1}{2}}$ in is not problematic as well.

Theorem \[thm:epsilontoinf\] shows that the new metric $\cM^{\epsilon}\sbra{\m{Y}}$ plays the same role as $\atomn{\m{Y}}$ as $\epsilon\rightarrow+\infty$, while Theorem \[thm:epsilontozero\] states that it approaches $\norm{\m{Y}}_{\cA,0}$ as $\epsilon\rightarrow0$. Consequently, it bridges $\atomn{\m{Y}}$ and $\norm{\m{Y}}_{\cA,0}$ and is expected to enhance sparsity and resolution compared to $\atomn{\m{Y}}$. Moreover, Theorem \[thm:epsilontozero\] characterizes the properties of the optimizer $\m{u}_{\epsilon}^*$ as $\epsilon\rightarrow0$ including the convergent speed of the smallest $N-\norm{\m{Y}}_{\cA,0}$ eigenvalues and the limiting form of $T\sbra{\m{u}^*}$ via the Vandermonde decomposition. In fact, we always observe via simulations that the smallest $N-\norm{\m{Y}}_{\cA,0}$ eigenvalues of $T\sbra{\m{u}^*}$ become zero once $\epsilon$ is appropriately small.

In DOA estimation, a difficult scenario is when the source signals $\lbra{\m{s}_k}_{k=1}^K$ are highly or even completely correlated (the latter case is usually called coherent). For example, covariance-based methods like Capon’s beamforming and MUSIC cannot produce satisfactory results since a faithful covariance estimate is unavailable. In contrast, the proposed sparse metric is robust to source correlations by Theorem \[thm:epsilontozero\] in which we have not made any assumption for the sources.

According to Theorem \[thm:epsilontozero\], the solution $T\sbra{\m{u}^*}$ can be interpreted as the data covariance of $\m{Y}$ after removing correlations among the sources. \[rem:datacov\]

Reweighted Atomic-Norm Minimization (RAM) {#sec:RAM}
=========================================

A Locally Convergent Iterative Algorithm
----------------------------------------

Using the proposed sparse metric $\cM^{\epsilon}\sbra{\m{Y}}$ we solve the following optimization problem for signal and frequency recovery: or equivalently, Note that $\ln\abs{T\sbra{\m{u}}+\epsilon\m{I}}$ is a concave function of $\m{u}$ since $\ln\abs{\m{R}}$ is a concave function of $\m{R}$ on the positive semidefinite cone [@boyd2004convex]. It follows that the problem in is nonconvex and no efficient algorithms can guarantee to obtain the global optimizer. A popular locally convergent approach to minimization of such a concave $+$ convex function is the majorization-maximization (MM) algorithm (see, e.g., [@fazel2003log]). Let $\m{u}_j$ denote the $j$th iterate of the optimization variable $\m{u}$. Then, at the $\sbra{j+1}$th iteration we replace $\ln\abs{T\sbra{\m{u}}+\epsilon\m{I}}$ by its tangent plane at the current value $\m{u}=\m{u}_j$: where $c_j$ is a constant independent of $\m{u}$. As a result, the optimization problem at the $\sbra{j+1}$th iteration becomes Since $\ln\abs{T\sbra{\m{u}}+\epsilon\m{I}}$ is strictly concave in $\m{u}$, at each iteration its value decreases by an amount greater than the decrease of its tangent plane. It follows that by iteratively solving the objective function in (\[formu:problem\]) monotonically decreases and converges to a local minimum.

Interpretation as RAM  {#sec:WAN_SDP}
----------------------

To interpret the optimization problem in (\[formu:problem\_j\]), let us define a weighted continuous dictionary w.r.t. the original continuous dictionary $\lbra{\m{a}\sbra{f}:\; f\in\bT}$, where $w\sbra{f}\geq 0$ is a weighting function. For $\m{Y}\in\bC^{N\times L}$, we define its weighted atomic norm w.r.t. $\cA^w$ as its atomic norm induced by $\cA^w$: According to the definition above, $w\sbra{f}$ specifies preference of the atoms $\lbra{\m{a}\sbra{f}}$. To be specific, an atom $\m{a}\sbra{f_0}$, $f_0\in\bT$, is more likely selected if $w\sbra{f_0}$ is larger. Moreover, the atomic norm is a special case of the weighted atomic norm with a constant weighting function (i.e., without any preference). Similar to the atomic norm, the proposed weighted atomic norm also admits a semidefinite formulation when assigned an appropriate weighting function, which is stated in the following theorem.

Suppose that $w\sbra{f}= \frac{1}{\sqrt{\m{a}\sbra{f}^H\m{W}\m{a}\sbra{f}}}$ with $\m{W}\in\bC^{N\times N}$. Then, \[thm:weightAN\]

See Appendix \[sec:Append\_weightAN\].

Let $\m{W}_j=\frac{1}{N}\sbra{T\sbra{\m{u}_j}+\epsilon\m{I}}^{-1}$ and $w_j\sbra{f} = \frac{1}{\sqrt{\m{a}\sbra{f}^H \m{W}_j \m{a}\sbra{f}}}$. It follows from Theorem \[thm:weightAN\] that the optimization problem in (\[formu:problem\_j\]) can be exactly written as the following weighted atomic norm minimization problem: As a result, the proposed iterative algorithm can be interpreted as reweighted atomic-norm minimization (RAM), where the weighting function is updated based on the latest solution of $\m{u}$. If we let $w_0(f)$ be a constant function or equivalently, $\m{u}_0=\m{0}$, such that no preference of the atoms is specified at the first iteration, then the first iteration coincides with ANM. From the second iteration on, the preference is defined by the weighting function $w_j\sbra{f}$ given above. Note that $w_j^2(f)$ is in fact the power spectrum of Capon’s beamforming provided that $T\sbra{\m{u}_j}$ is interpreted as the noiseless data covariance following from Remark \[rem:datacov\] and $\epsilon$ as the noise variance. Therefore, the reweighting strategy makes the frequencies around those produced by the current iteration preferable at the next iteration and thus enhances sparsity. At the same time, the preference results in finer details of the frequency spectrum in those areas and therefore enhances resolution. Empirical evidences will be provided in Section \[sec:simulation\].

Computationally Efficient Implementations {#sec:implementation}
=========================================

Optimization Using Standard SDP Solver
--------------------------------------

At each iteration of RAM, we need to solve the SDP in (\[formu:problem\_j\]) as follows: where $\m{W}=\sbra{T\sbra{\m{u}_j}+\epsilon\m{I}}^{-1}$. Its dual problem is given as follows by a standard Lagrangian analysis (see Appendix \[sec:append\_duality\]): where $\Re$ takes the real part of the argument and $Z_{n,j}$ denotes the $(n,j)$th entry of $\m{Z}$. We empirically find that the dual problem (\[formu:duality\]) can be solved more efficiently than the primal problem (\[formu:weightedSDP\]) with a standard SDP solver SDPT3 [@toh1999sdpt3]. Note that the optimizer to (\[formu:weightedSDP\]) is given for free via duality when we solve (\[formu:duality\]). As a result, the reweighted algorithm can be iteratively implemented.

A First-order Algorithm via ADMM
--------------------------------

A reasonably fast approach for ANM is based on ADMM [@boyd2011distributed; @bhaskar2013atomic; @yang2015gridless], which is a first-order algorithm and guarantees global optimality. To derive the ADMM algorithm, we reformulate the SDP in (\[formu:weightedSDP\]) as follows: which is very similar to the SDP solved in [@yang2015gridless]. Then we can write the augmented Lagrangian function and iteratively update $\sbra{\m{u},\m{X},\m{Y}}$, $\m{\cQ}$ and the Lagrangian multiplier in closed-form expressions until convergence. We omit the details since all the formulae and derivations are similar to those in [@yang2015gridless], to which interested readers are referred. We mention that an eigen-decomposition of a matrix of order $N+L$ (the order of $\m{\cQ}$) is required at each iteration. Note that the ADMM converges slowly to an extremely accurate solution while moderate accuracy is typically sufficient in practical applications [@boyd2011distributed].

Dimension Reduction for Large $L$
---------------------------------

The number of measurement vectors $L$ can be large, possibly with $L\gg M$, in DOA estimation, which increases considerably the computational workload. We provide the following result to reduce this number from $L$ to $\rank\sbra{\m{Y}_{\m{\Omega}}^o}\leq \min(L,M)$.

Let $r=\rank\sbra{\m{Y}_{\m{\Omega}}^o}\leq \min(L,M)$. Find a unitary matrix $\m{Q}\in\bC^{L\times L}$ (for example, by QR decomposition) satisfying that $\m{Y}_{\m{\Omega}}^{o}\m{Q}= \begin{bmatrix}\m{Y}_{\m{\Omega}}^{o}\m{Q}_1 & \m{0}\end{bmatrix}$, where $\m{Q}_1\in\bC^{L\times r}$. If we make the substitutions $\m{Y}_{\m{\Omega}}^{o}\rightarrow\m{Y}_{\m{\Omega}}^{o}\m{Q}_1$ and $\m{Y}\in\bC^{N\times L}\rightarrow\m{Z}\in\bC^{N\times r}$ in (\[formu:problem\_j\]) and denote by $\sbra{\m{Z}^*,\m{u}^*}$ the optimizer to the resulting optimization problem, then the optimizer to (\[formu:problem\_j\]) is given by $\sbra{\m{Z}^*\m{Q}_1^H,\m{u}^*}$. The same result holds for the nonconvex optimization problem in (\[formu:criterion3\]). \[prop:dimreduce\]

See Appendix \[sec:Append\_dimreduce\].

If only the frequencies are of interest, e.g., in DOA estimation, we can replace $\m{Y}_{\m{\Omega}}^{o}\m{Q}_1$ by $\sbra{\m{Y}_{\m{\Omega}}^{o}\m{Y}_{\m{\Omega}}^{oH}}^{\frac{1}{2}}$ (in fact, any matrix $\widetilde{\m{Y}}$ satisfying that $\widetilde{\m{Y}}\widetilde{\m{Y}}^H=\m{Y}_{\m{\Omega}}^{o}\m{Y}_{\m{\Omega}}^{oH}$) to obtain the same solution of $\m{u}$. \[rem:sameusolution\]

With a similar proof, the dimension reduction technique in Proposition \[prop:dimreduce\] can be extended to a more general linear model with observations expressed by $\m{\Phi}\m{Y}^o$ + noise, where $\m{\Phi}$ denotes a sensing matrix. Also, it can be applied to conventional discrete dictionary models, in which, for example, the $\ell_{2,p}$ norm, $0\leq p\leq1$, needs to be optimized. Note that the dimension reduction approach introduced here is different from that in [@malioutov2005sparse]. In particular, the approach in [@malioutov2005sparse] requires the knowledge of the model order $K$, which is unknown in practical scenarios, and gives an optimization problem which is an approximation of the original one. In contrast, our approach does not need $K$ but produces a dimension-reduced, equivalent problem.

By Proposition \[prop:dimreduce\], when $L>M$ we can reduce the order of the PSD matrix in (\[formu:weightedSDP\]) \[and (\[formu:duality\]), (\[formu:weightedSDP2\])\] from $N+L$ to $N+r\leq N+M$. Therefore, the resulting problem dimension depends only on $M$ and $N$. Both the SDPT3 and ADMM implementations of RAM above can be reasonably fast when $M$ and $N$ are small though they may not possess good scalability, especially for SDPT3. An application at hand is DOA estimation in which the array size $M$ and aperture $N$ are typically small (on the order of 10) though $L$ can be a few hundred or even greater. Note also that the dimension reduction technique takes $O\sbra{M^2L}$ flops in DOA estimation according to Remark \[rem:sameusolution\]. Extensive numerical simulations will be provided in Section \[sec:simulation\] to demonstrate usefulness of our method.

Remarks on Algorithm Implementation
-----------------------------------

According to the discussions in Section \[sec:WAN\_SDP\], we can always start with $\m{u}_0=\m{0}$ and the first iteration coincides with ANM. When $L$ is large, we can also implement a weighting function in the first iteration inspired by Capon’s beamforming for faster convergence (see an example in Section \[sec:simulation\_DOA\]). Moreover, we gradually decrease $\epsilon$ during the algorithm and define the weighting function using the latest solution for avoiding local minima (note that the first iteration with $\m{u}_0=\m{0}$ essentially corresponds to $\epsilon=+\infty$ following from Theorem \[thm:epsilontoinf\]). In fact, this is like an aggressive continuation strategy in which we attempt to solve the nonconvex optimization problem in (\[formu:criterion3\]) at decreasing values of $\epsilon$. The convergence of the reweighted algorithm is retained if we fix $\epsilon$ when it is sufficiently small. In the ADMM implementation, we can further accelerate the algorithm by adopting loose convergence criteria in the first few iterations of RAM. Finally, note that we need to trade off the algorithm performance for computational time by keeping the number of iterations of RAM being small.

Numerical Simulations {#sec:simulation}
=====================

Implementation Details of RAM
-----------------------------

In our implementation of RAM, we first scale the measurements and the noise such that $\frobn{\m{Y}_{\m{\Omega}}}^2=M$ (the noise energy becomes $\eta'^2=\frac{M\eta^2}{\frobn{\m{Y}}^2}$) and compensate the recovery afterwards. We start with $\m{u}_0=\m{0}$ and $\epsilon=1$ as default. We halve $\epsilon$ when beginning a new iteration until $\epsilon=\frac{1}{2^{10}}$ or $\epsilon<\frac{\eta'^2}{10}$. When $\eta=0$ we terminate RAM if the relative change (in the Frobenius norm) of the solution $\m{Y}^*$ at two consecutive iterations is less than $10^{-6}$ or the maximum number of iterations, set to 20, is reached. All simulations are carried out in Matlab v.8.1.0 on a PC with a Windows 7 system and a 3.4 GHz CPU.

An Illustrative Example
-----------------------

We provide a simple example in this subsection to illustrate the iterative process of RAM. In particular, we consider a sparse frequency spectrum consisting of $K=5$ spikes located at 0.1, 0.108, 0.125, 0.2 and 0.5. We randomly generate the complex amplitudes and randomly select $M=30$ samples among $N=64$ consecutive time-space measurements, with $L=1$. Then, we run the RAM algorithm to reconstruct the frequency spectrum from the samples. Note that the first three frequencies are mutually separated by only about $\frac{0.51}{N}$ and $\frac{1.09}{N}$. Implemented with SDPT3 the RAM algorithm converges in four iterations. We plot the simulation results in Fig. \[Fig:illusexample\], where the first subfigure presents variation of the eigenvalues of $T\sbra{\m{u}^*}$ during the iterations, the second row presents the recovered spectra of the first three iterations, and the last row plots the weighting functions used in the first three iterations. Note that the first iteration, which exploits a constant weighting function and coincides with ANM, can detect a rough area where the first three spikes are located but cannot accurately determine their locations and number. In the second iteration, a weighting function is implemented based on the previous estimate to provide preference of the frequencies around those produced in the first iteration. As a result, the third spike is identified while the first two are still not. Following from the same reweighting process, all the frequencies are correctly determined in the next iteration and the algorithm converges after that. It is worth noting that $T\sbra{\m{u}^*}$ becomes rank-5 and the remaining eigenvalues become zero (within numerical precision) after three iterations, where $\epsilon=0.25$. Finally, we report that the relative mean squared error (MSE) of signal recovery improves during the iterations from $1.28\times10^{-4}$ to $2.01\times 10^{-7}$, $3.16\times10^{-19}$ and $3.01\times10^{-21}$. Each iteration takes about 1.7s.

![image](Fig_illusexample.pdf){width="6in"}

Sparsity-Separation Phase Transition
------------------------------------

In this subsection, we study the success rate of RAM in signal and frequency recovery compared to ANM. In particular, we fix $N=64$ and $M=30$ with the sampling index set $\m{\Omega}$ being generated uniformly at random. We vary the duo $\sbra{K,\Delta_f}$ and for each combination we randomly generate $K$ frequencies such that they are mutually separated by at least $\Delta_f$. We randomly generate the amplitudes $\lbra{s_{kt}}$ independently and identically from a standard complex normal distribution. After obtaining the noiseless samples, we carry out signal reconstruction and frequency recovery using ANM and RAM, both implemented by SDPT3. The recovery is called successful if both the relative MSE of signal recovery and the MSE of frequency recovery are less than $10^{-12}$. For each combination $\sbra{K,\Delta_f}$, the success rate is measured over 20 Monte Carlo runs.

We plot the success rates of ANM and RAM with $L=1$ in Fig. \[Fig:phasetrans\_1\], where it is shown that successful recovery can be obtained with more ease in the case of a smaller $K$ and a larger frequency separation $\Delta_f$, leading to a phase transition in the sparsity-separation domain. By comparing the two images, we see that RAM significantly enlarges the success phase and enhances sparsity and resolution. It is worth noting that the phase transitions of both ANM and RAM are not sharp. One reason is that, a set of *well separated* frequencies can be possibly generated at a small value of $\Delta_f$ while we only control that the frequencies are separated by *at least* $\Delta_f$. It is also observed that RAM tends to converge in less iterations with a smaller $K$ and a larger $\Delta_f$.

![Sparsity-separation phase transition of ANM (left) and RAM (right) with $L=1$, $N=64$ and $M=30$. The grayscale images present the success rates, where white and black indicate complete success and complete failure, respectively.[]{data-label="Fig:phasetrans_1"}](Fig_phasetrans_ANM_1.pdf "fig:"){width="1.65in"} ![Sparsity-separation phase transition of ANM (left) and RAM (right) with $L=1$, $N=64$ and $M=30$. The grayscale images present the success rates, where white and black indicate complete success and complete failure, respectively.[]{data-label="Fig:phasetrans_1"}](Fig_phasetrans_RWANM_1.pdf "fig:"){width="1.63in"}

We also consider the MMV case with $L=5$. The success rates of ANM and RAM are presented in Fig. \[Fig:phasetrans\_5\]. Again, remarkable improvement is obtained by the proposed RAM compared to ANM. In fact, we did not find a single failure in our simulation whenever $K\leq20$ and $\Delta_f\geq\frac{0.3}{N}$. By comparing the results in Figs. \[Fig:phasetrans\_1\] and \[Fig:phasetrans\_5\], it can be observed that improved signal and frequency recovery performance can be obtained by increasing $L$, as reported in [@yang2014continuous; @yang2014exact].

![Sparsity-separation phase transition of ANM (left) and RAM (right) with $L=5$, $N=64$ and $M=30$. The grayscale images present the success rates, where white and black indicate complete success and complete failure, respectively.[]{data-label="Fig:phasetrans_5"}](Fig_phasetrans_ANM_5.pdf "fig:"){width="1.65in"} ![Sparsity-separation phase transition of ANM (left) and RAM (right) with $L=5$, $N=64$ and $M=30$. The grayscale images present the success rates, where white and black indicate complete success and complete failure, respectively.[]{data-label="Fig:phasetrans_5"}](Fig_phasetrans_RWANM_5.pdf "fig:"){width="1.63in"}

Application to DOA Estimation {#sec:simulation_DOA}
-----------------------------

We apply the proposed RAM method to DOA estimation. In particular, we consider a 10-element SLA that is obtained from a virtual 20-element ULA, in which adjacent antennas are spaced by half a wavelength, by retaining the antennas indexed by $\m{\Omega}=\lbra{1,2,5,6,8,12,15,17,19,20}$. Hence, we have that $N=20$ and $M=10$. Consider that $K=4$ narrowband sources impinge on the array from directions corresponding to frequencies $0.1$, $0.11$, $0.2$ and $0.5$, and powers $10$, $10$, $3$ and $1$, respectively. Therefore, it is challenging to separate the first two sources which are separated by only $\frac{0.2}{N}$. We consider both cases of uncorrelated and correlated sources. In the latter case, sources 1 and 3 are set to be coherent (completely correlated). Assume that $L=200$ data snapshots are collected which are corrupted by i.i.d. Gaussian noise of unit variance. In our simulation, ANM and RAM are implemented using both SDPT3 and ADMM and based on the proposed dimension reduction technique. A nontrivial weighting function is implemented in the first iteration of RAM with $\m{W}=\m{\Gamma}_{\m{\Omega}}^T \sbra{\frac{1}{L}\m{Y}_{\m{\Omega}}^o\m{Y}_{\m{\Omega}}^{oH}+\epsilon\m{I}}^{-1} \m{\Gamma}_{\m{\Omega}}$, where $\m{\Gamma}_{\m{\Omega}}\in\lbra{0,1}^{M\times N}$ has 1 in the $j$th row only at the $\Omega_j$th position. The weighting function corresponds to Capon’s beamforming with $\frac{1}{L}\m{Y}_{\m{\Omega}}^o\m{Y}_{\m{\Omega}}^{oH}$ being the sample covariance and $\epsilon$ a regularization parameter. We terminate RAM within maximally 10 iterations. We consider MUSIC and ANM for comparison. Assume that the noise variance $\sigma^2=1$ is given for ANM and RAM and the source number $K$ is provided for MUSIC. We set $\eta^2=\sbra{ML+2\sqrt{ML}}\sigma^2$ (mean + twice standard deviation) to upper bound the noise energy with high probability in ANM and RAM.

![Frequency spectra of MUSIC (top), ANM (middle) and RAM (bottom) with uncorrelated (left) and correlated (right) sources in 100 Monte Carlo runs. Sources 1 and 3 are coherent in the case of correlated sources. The area around the first two sources are zoomed in in each subfigure. Only results of the first 20 runs are presented for MUSIC for clearer illustration.[]{data-label="Fig:noisy_noiselevel"}](Fig_noisy_uncor_K_MUSIC.pdf "fig:"){width="1.65in"} ![Frequency spectra of MUSIC (top), ANM (middle) and RAM (bottom) with uncorrelated (left) and correlated (right) sources in 100 Monte Carlo runs. Sources 1 and 3 are coherent in the case of correlated sources. The area around the first two sources are zoomed in in each subfigure. Only results of the first 20 runs are presented for MUSIC for clearer illustration.[]{data-label="Fig:noisy_noiselevel"}](Fig_noisy_cor_K_MUSIC.pdf "fig:"){width="1.63in"} ![Frequency spectra of MUSIC (top), ANM (middle) and RAM (bottom) with uncorrelated (left) and correlated (right) sources in 100 Monte Carlo runs. Sources 1 and 3 are coherent in the case of correlated sources. The area around the first two sources are zoomed in in each subfigure. Only results of the first 20 runs are presented for MUSIC for clearer illustration.[]{data-label="Fig:noisy_noiselevel"}](Fig_noisy_uncor_noiselevel_ANM.pdf "fig:"){width="1.65in"} ![Frequency spectra of MUSIC (top), ANM (middle) and RAM (bottom) with uncorrelated (left) and correlated (right) sources in 100 Monte Carlo runs. Sources 1 and 3 are coherent in the case of correlated sources. The area around the first two sources are zoomed in in each subfigure. Only results of the first 20 runs are presented for MUSIC for clearer illustration.[]{data-label="Fig:noisy_noiselevel"}](Fig_noisy_cor_noiselevel_ANM.pdf "fig:"){width="1.63in"} ![Frequency spectra of MUSIC (top), ANM (middle) and RAM (bottom) with uncorrelated (left) and correlated (right) sources in 100 Monte Carlo runs. Sources 1 and 3 are coherent in the case of correlated sources. The area around the first two sources are zoomed in in each subfigure. Only results of the first 20 runs are presented for MUSIC for clearer illustration.[]{data-label="Fig:noisy_noiselevel"}](Fig_noisy_uncor_noiselevel_RAM.pdf "fig:"){width="1.65in"} ![Frequency spectra of MUSIC (top), ANM (middle) and RAM (bottom) with uncorrelated (left) and correlated (right) sources in 100 Monte Carlo runs. Sources 1 and 3 are coherent in the case of correlated sources. The area around the first two sources are zoomed in in each subfigure. Only results of the first 20 runs are presented for MUSIC for clearer illustration.[]{data-label="Fig:noisy_noiselevel"}](Fig_noisy_cor_noiselevel_RAM.pdf "fig:"){width="1.63in"}

Our simulation results of 100 Monte Carlo runs are presented in Fig. \[Fig:noisy\_noiselevel\] (only the first 20 runs are presented for MUSIC for better illustration). In the absence of source correlations, MUSIC has satisfactory performance in most scenarios. However, its power spectrum exhibits only a single peak around the first two sources (i.e., the two sources cannot be separated) in at least 3 out of the first 20 runs (indicated by the arrows). Moreover, MUSIC is sensitive to source correlations and cannot detect source 1 when it is coherent with source 3. ANM cannot separate the first two sources in the uncorrelated source case and always produces many spurious sources. In contrast, the proposed RAM always correctly detects 4 sources near the true locations, demonstrating its capabilities in enhancing sparsity and resolution. It is interesting to note that the first two sources tend to merge together in RAM. This is reasonable since in the case of heavy noise it is even possible that we can only detect a single source around the first two and 3 sources in total. Note also that the results of ANM and RAM presented in Fig. \[Fig:noisy\_noiselevel\] are produced by the ADMM implementations, while those by SDPT3 are very similar and omitted. In computational time, the SDPT3 versions of ANM and RAM take $0.87$s and $7.31$s on average, respectively, while these numbers are decreased to $0.20$s and $2.65$s for the ADMM ones.

Conclusion {#sec:conclusion}
==========

In this paper, we studied the signal and frequency recovery problem in CCS. Motivated by its connection to the topic of LRMR, we proposed reweighted atomic-norm minimization (RAM) for enhancing sparsity and resolution compared to currently prominent atomic norm minimization (ANM) and validated its advantageous performance via numerical simulations. As a byproduct, we have established a framework for applying LRMR techniques to CCS. In future studies, we may try other methods for matrix rank minimization in the literature and propose more computationally efficient algorithms for CCS. While LRMR represents a 2D counterpart of sparse recovery in discrete CS, this work sheds light on connections between discrete CS, continuous CS and LRMR.

Proof of Theorem \[thm:epsilontoinf\] {#sec:append:epsilontoinf}
-------------------------------------

Note that Let Then, according to (\[formu:AN\_SDP\]) we have that Consider the value of the objective function in (\[formu:resid\]) at $\m{u}=\epsilon^{\frac{1}{2}}\m{u}^*$. It holds that On the other hand, we denote by $\m{u}_{\epsilon}^*$ the optimizer to the optimization problem in (\[formu:resid\]). We first argue that $\tr\sbra{T\sbra{\m{u}_{\epsilon}^*}}\epsilon^{-1}=o\sbra{1}$. Otherwise, by (\[formu:resid\]) $\cM^{\epsilon}\sbra{\m{Y}}-N\ln\epsilon\geq \ln\abs{\epsilon^{-1}T\sbra{\m{u}_{\epsilon}^*} +\m{I} }$ is not an infinitesimal, contradicting (\[formu:residleq\]). Therefore, Combining (\[formu:residleq\]) and the second equality in (\[formu:residgeq\]) yields that $\tr\sbra{T\sbra{\m{u}_{\epsilon}^*}}\epsilon^{-1}=O\sbra{\epsilon^{-\frac{1}{2}}}$. Then, the last equality in (\[formu:residgeq\]) gives that The conclusion is finally drawn by combining (\[formu:residleq\]) and (\[formu:residgeq2\]).

Proof of Theorem \[thm:epsilontozero\] {#sec:append:epsilontozero}
--------------------------------------

Our proof is given in four steps. Let $\lbra{\lambda_{\epsilon,k}}_{k=1}^N$ be the eigenvalues of $T\sbra{\m{u}_{\epsilon}^*}$ that are sorted descendingly. In *Step 1*, we attempt to show that there exists a constant $c>0$ such that $\lambda_{\epsilon,r}\geq c$ holds uniformly for $\epsilon\in(0,1]$. Let $T\sbra{\m{u}_{\epsilon}^*}=\sum_{k=1}^N\lambda_{\epsilon,k}\m{q}_{\epsilon,k}\m{q}_{\epsilon,k}^H= \m{Q}\diag\sbra{\lambda_{\epsilon,1},\dots,\lambda_{\epsilon,N}}\m{Q}^H$ be the eigen-decomposition, where $\m{q}_{\epsilon,k}$ is the $k$th column of $\m{Q}$ and $\m{Q}\m{Q}^H=\m{I}$. Then, where $\overline p_{\epsilon,k} \triangleq\twon{\m{q}_{\epsilon,k}^H\m{Y}}^2$. According to the optimality of $\m{u}_{\epsilon}^*$, the right hand side of the equation above obtains its minimum at $\lambda_{\epsilon,k}$. Since its derivative at $\lambda_{\epsilon,k}$ equals $\frac{1}{\lambda_{\epsilon,k}+\epsilon} -\frac{\overline{p}_{\epsilon,k}}{\lambda_{\epsilon,k}^2}$, we have that Therefore, provided that $\epsilon\leq 1$. It follows that $\m{u}_{\epsilon}^*$ and $\lbra{\lambda_{\epsilon,k}}$ are bounded.

On the other hand, let $T\sbra{\m{u}_{\epsilon}^*}=\sum_{k=1}^{N} p_{\epsilon,k}\m{a}\sbra{f_{\epsilon,k}}\m{a}\sbra{f_{\epsilon,k}}^H=\m{A}\m{P}\m{A}^H$ be any Vandermonde decomposition, where $\lbra{p_{\epsilon,k}}_{k=1}^N$ are sorted descendingly (note that, if $r_{\epsilon}^*=\rank\sbra{T\sbra{\m{u}_{\epsilon}^*}}<N$, then this decomposition is unique and only the first $r_{\epsilon}^*$ elements in $\lbra{p_{\epsilon,k}}_{k=1}^N$ are nonzero). Following from the fact that $\m{Y}$ lies in the range space of $\m{A}$, we have that $\m{Y}=\m{A}\m{S}$ for some $\m{S}$. Let $\m{s}_{\epsilon,k}$ be its $k$th row of $\m{S}$. Then we have According to the optimality of $\m{u}_{\epsilon}^*$, the right hand side of the equation above obtains its minimum at $p_{\epsilon,k}$. As a result, its derivative at $p_{\epsilon,k}$ equals 0, i.e., and so $p_{\epsilon,k}>\twon{\m{s}_{\epsilon,k}}^2$ since provided that $\epsilon>0$. Since $\tr\sbra{T\sbra{\m{u}_{\epsilon}^*}}=N\sum_{k=1}^{N}p_{\epsilon,k}$ and that $\m{u}_{\epsilon}^*$ is bounded as shown previously, $\lbra{p_{\epsilon,k}}$ and $\lbra{\m{s}_{\epsilon,k}}$ are bounded.

We now prove that $\lambda_{\epsilon,r}\geq c$ for some constant $c$. Otherwise, for any $c_j=\frac{1}{j}$, $j=1,2,\dots$, there exists $\epsilon_j\in(0,1]$ such that $\lambda_{\epsilon_j,r}< c_j=\frac{1}{j}$. Since the sequence $\lbra{\sbra{\m{u}_{\epsilon_j}^*, \lambda_{\epsilon_j,k}, \m{q}_{\epsilon_j,k}, p_{\epsilon_j,k}, f_{\epsilon_j,k}, \m{s}_{\epsilon_j,k}}}_{j=1}^{\infty}$ is bounded, there must exist a convergent subsequence. Without loss of generality, we assume that the sequence is convergent itself and denote by $\sbra{\m{u}^*, \lambda_{k}, \m{q}_{k}, p_{k}, f_{k}, \m{s}_{k}}$ the limit point. Since $\lambda_r=\lim_{j\rightarrow\infty}\lambda_{\epsilon_j,r}=0$, we have $\lambda_{k}=0$ for all $k=r+1,\dots,N$. It follows that and $\rank\sbra{T\sbra{\m{u}^*}}\leq r-1$. As a result, at most $r-1$ $f_k$’s are retained in the decomposition $T\sbra{\m{u}^*}=\sum_{k=1}^{N} p_{k}\m{a}\sbra{f_{k}}\m{a}\sbra{f_{k}}^H$ if we remove repetitive $f_k$’s and those with $p_k=0$. Note that $\m{s}_k=\m{0}$ if $p_k=0$ since we have shown that $p_{\epsilon,k}\geq \twon{\m{s}_{\epsilon,k}}^2$. Then, by a similar operation we can reduce the order of the decomposition $\m{Y}=\sum_{k=1}^N \m{a}\sbra{f_k}\m{s}_k$ to maximally $r-1$, i.e., we obtain an atomic decomposition of order at most $r-1$, which contradicts the fact that $\norm{\m{Y}}_{\cA,0}=r$ and leads to the conclusion.

In *Step 2* we prove the first part of the theorem. According to (\[formu:equ1\]) and the bound $\lambda_{\epsilon,r}\geq c$ shown in *Step 1*, we have that On the other hand, we consider an atomic decomposition of $\m{Y}$ of order $r$, $\m{Y}=\sum_{k=1}^r\m{a}\sbra{f_k}\m{s}_k$. Let $T\sbra{\m{u}}=\sum_{k=1}^r\twon{\m{s}_k}^2\m{a}\sbra{f_k}\m{a}\sbra{f_k}^H$, and $\lbra{\lambda_k}_{k=1}^r$ be the $r$ nonzero eigenvalues of $T\sbra{\m{u}}$. Note that $\lbra{\lambda_k}_{k=1}^r$ are constants independent of $\epsilon$. Then, provided $\epsilon\leq 1$ we have that Combining (\[formu:Mgeqside\]) and (\[formu:Mleqside\]), it yields that $\cM^{\epsilon}\sbra{\m{Y}}\sim (N-r)\ln\epsilon$ as $\epsilon\rightarrow0$.

In *Step 3* we prove the second part of the theorem. Based on (\[formu:Mgeqside\]) and (\[formu:Mleqside\]) we have that where $c_1$ and $c_2$ are constants independent of $\epsilon$. Therefore, it must hold that It follows that i.e., $\lambda_{\epsilon,k}=O\sbra{\epsilon}$, $k=r+1,\dots,N$.

Finally, we show the last part of the theorem. For any cluster point $\m{u}_0^*$ of $\m{u}_{\epsilon}^*$ at $\epsilon=0$, there exists a sequence $\lbra{\m{u}_{\epsilon_j}^*}_{j=1}^{\infty}$ converging to $\m{u}_0^*$, where $\epsilon_j\rightarrow0$ as $j\rightarrow\infty$. It must hold that $\rank\sbra{T\sbra{\m{u}_0^*}}=r$ since the smallest $N-r$ eigenvalues of $T\sbra{\m{u}_{\epsilon_j}^*}$ approach 0. Moreover, the eigen-decomposition of $T\sbra{\m{u}_{\epsilon_j}^*}$ converge to that of $T\sbra{\m{u}_0^*}$, where again we denote their eigenvalues by $\lbra{\lambda_{\epsilon_j,k}}$ and $\lbra{\lambda_{k}}$ respectively and use the other notations similarly. Then, according to (\[formu:lambdapbar\]) we have that $\overline{p}_{\epsilon_j,k}=\frac{\lambda_{\epsilon_j,k}^2}{\lambda_{\epsilon_j,k}+\epsilon}\rightarrow \lambda_k=\overline{p}_{k}$, as $j\rightarrow\infty$. Therefore, On the other hand, we similarly write the Vandermonde decomposition of $T\sbra{\m{u}_{\epsilon}^*}$ and let $T\sbra{\m{u}_0^*}=\sum_{k=1}^r p_k\m{a}\sbra{f_k}\m{a}\sbra{f_k}^H$, with $\m{Y}=\sum_{k=1}^r\m{a}\sbra{f_k}\m{s}_k$. It is easy to show that $p_{k}\geq\twon{\m{s}_{k}}^2$ based on the inequality $p_{\epsilon,k}>\twon{\m{s}_{\epsilon,k}}^2$, though $p_{\epsilon,k}$ and $\m{s}_{\epsilon,k}$ do not necessarily converge to $p_{k}$ and $\m{s}_{k}$ (consider the case where an accumulation point of $\lbra{f_{\epsilon_j,1},\dots,f_{\epsilon_j,N}}_{j=1}^{\infty}$ contains identical elements). Then, where the equality holds iff $p_{k}=\twon{\m{s}_{k}}^2$. So, we complete the proof.

Proof of Theorem \[thm:weightAN\] {#sec:Append_weightAN}
---------------------------------

The conclusion is a direct result of the following equalities: where the first equality applies the Vandermonde decomposition, and the second follows the equality (see [@yang2014exact]) given the expression of $\m{R}$ in (\[formu:weightAN\_derivation\]). Note that in the first equality in we did not specify the order of the Vandermonde decomposition of $T(\m{u})$. This means that the proof holds true for any possible decomposition (whenever $T(\m{u})$ is invertible or not).

Lagrangian Analysis of the Dual Problem (\[formu:duality\]) {#sec:append_duality}
-----------------------------------------------------------

Let $\m{\Lambda}=\begin{bmatrix}\m{U}&\m{V}^H \\ \m{V} & \m{Z}\end{bmatrix}\geq\m{0}$. The Lagrangian function of (\[formu:weightedSDP\]) is given as follows: Minimizing $\cL$ with respect to $\sbra{\m{u},\m{X},\m{Y}}$ gives the dual objective which equals $-\lambda^{-1}\frobn{\m{V}_{\m{\Omega}}}^2 -\lambda\eta^2 -2\Re\tr\sbra{\m{Y}_{\m{\Omega}}^{oH} \m{V}_{\m{\Omega}}}$, if or $-\infty$, otherwise, where $T^*\sbra{\cdot}$ denotes the adjoint operator of $T\sbra{\cdot}$. Therefore, we obtain the dual problem in (\[formu:duality\]) by noting that

Proof of Proposition \[prop:dimreduce\] {#sec:Append_dimreduce}
---------------------------------------

Regarding (\[formu:problem\_j\]) and (\[formu:problem\]) we consider the following optimization problem: where $\m{C}\geq\m{0}$ is fixed. We replace the optimization variable $\m{Y}$ by $\m{Z}=\m{Y}\m{Q}$. Since $\m{Q}$ is a unitary matrix, the problem becomes and equivalently, where $\m{Z}_j=\m{Y}\m{Q}_j$, $j=1,2$. Denote the optimizer by $\sbra{\m{Z}_1^*,\m{Z}_2^*}$. It is obvious that $\m{Z}_{2}^*=\m{0}$ since $\tr\sbra{\m{Z}_2^H\m{C}\m{Z}_2}\geq0$. Then the problem becomes which is a dimension reduced version of (\[formu:probY\]) and has the same optimal function value. Moreover, given $\m{Z}^*=\begin{bmatrix}\m{Z}_1^*&\m{0}\end{bmatrix}$ we have the optimizer to (\[formu:probY\]) $\m{Y}^*=\m{Z}^*\m{Q}^H= \m{Z}_1^* \m{Q}_1^H$. Now the conclusion can be easily drawn.

[10]{} \[1\][\#1]{} url@samestyle \[2\][\#2]{} \[2\][[l@\#1=l@\#1\#2]{}]{}

Z. Yang and L. Xie, “Achieving high resolution for super-resolution via reweighted atomic norm minimization,” in *IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, 2015, pp. 3646–3650.

E. Cand[è]{}s, “Compressive sampling,” in *Proceedings of the International Congress of Mathematicians*, vol. 3, 2006, pp. 1433–1452.

D. Donoho, “[Compressed sensing]{},” *IEEE Transactions on Information Theory*, vol. 52, no. 4, pp. 1289–1306, 2006.

H. Krim and M. Viberg, “Two decades of array signal processing research: The parametric approach,” *IEEE Signal Processing Magazine*, vol. 13, no. 4, pp. 67–94, 1996.

P. Stoica and R. L. Moses, *Spectral analysis of signals*.1em plus 0.5em minus 0.4emPearson/Prentice Hall Upper Saddle River, NJ, 2005.

G. Tang, B. N. Bhaskar, P. Shah, and B. Recht, “Compressed sensing off the grid,” *IEEE Transactions on Information Theory*, vol. 59, no. 11, pp. 7465–7490, 2013.

Z. Yang and L. Xie, “Continuous compressed sensing with a single or multiple measurement vectors,” in *IEEE Workshop on Statistical Signal Processing (SSP)*, 2014, pp. 308–311.

D. A. Linebarger, I. H. Sudborough, and I. G. Tollis, “Difference bases and sparse sensor arrays,” *IEEE Transactions on Information Theory*, vol. 39, no. 2, pp. 716–721, 1993.

D. Malioutov, M. Cetin, and A. Willsky, “A sparse signal reconstruction perspective for source localization with sensor arrays,” *IEEE Transactions on Signal Processing*, vol. 53, no. 8, pp. 3010–3022, 2005.

M. Hyder and K. Mahata, “Direction-of-arrival estimation using a mixed $\ell_{2,0}$ norm approximation,” *IEEE Transactions on Signal Processing*, vol. 58, no. 9, pp. 4646–4655, 2010.

P. Stoica, P. Babu, and J. Li, “New method of sparse parameter estimation in separable models and its use for spectral analysis of irregularly sampled data,” *IEEE Transactions on Signal Processing*, vol. 59, no. 1, pp. 35–47, 2011.

P. Stoica, P. Babu, and J. Li, “[SPICE]{}: A sparse covariance-based estimation method for array processing,” *IEEE Transactions on Signal Processing*, vol. 59, no. 2, pp. 629–638, 2011.

X. Wei, Y. Yuan, and Q. Ling, “[DOA]{} estimation using a greedy block coordinate descent algorithm,” *IEEE Transactions on Signal Processing*, vol. 60, no. 12, pp. 6382–6394, 2012.

A. Fannjiang and W. Liao, “Coherence pattern-guided compressive sensing with unresolved grids,” *SIAM Journal on Imaging Sciences*, vol. 5, no. 1, pp. 179–202, 2012.

N. Hu, Z. Ye, X. Xu, and M. Bao, “[DOA]{} estimation for sparse array via sparse signal reconstruction,” *IEEE Transactions on Aerospace and Electronic Systems*, vol. 49, no. 2, pp. 760–773, 2013.

M. F. Duarte and R. G. Baraniuk, “Spectral compressive sensing,” *Applied and Computational Harmonic Analysis*, vol. 35, no. 1, pp. 111–129, 2013.

Z. Liu, Z. Huang, and Y. Zhou, “Sparsity-inducing direction finding for narrowband and wideband signals based on array covariance vectors,” *IEEE Transactions on Wireless Communications*, vol. 12, no. 8, pp. 3896–3907, 2013.

L. Hu, Z. Shi, J. Zhou, and Q. Fu, “Compressed sensing of complex sinusoids: An approach based on dictionary refinement,” *IEEE Transactions on Signal Processing*, vol. 60, no. 7, pp. 3809–3822, 2012.

Z. Yang, C. Zhang, and L. Xie, “Robustly stable signal recovery in compressed sensing with structured matrix perturbation,” *IEEE Transactions on Signal Processing*, vol. 60, no. 9, pp. 4658–4671, 2012.

Z. Yang, L. Xie, and C. Zhang, “Off-grid direction of arrival estimation using sparse [Bayesian]{} inference,” *IEEE Transactions on Signal Processing*, vol. 61, no. 1, pp. 38–43, 2013.

L. Hu, J. Zhou, Z. Shi, and Q. Fu, “A fast and accurate reconstruction algorithm for compressed sensing of complex sinusoids,” *IEEE Transactions on Signal Processing*, vol. 61, no. 22, pp. 5744–5754, 2013.

C. Austin, J. Ash, and R. Moses, “Dynamic dictionary algorithms for model order and parameter estimation,” *IEEE Transactions on Signal Processing*, vol. 61, no. 20, pp. 5117–5130, 2013.

Z. Tan, P. Yang, and A. Nehorai, “Joint sparse recovery method for compressed sensing with structured dictionary mismatch,” *IEEE Transactions on Signal Processing*, vol. 62, no. 19, pp. 4997–5008, 2014.

Y. Chi, L. Scharf, A. Pezeshki, and A. Calderbank, “Sensitivity to basis mismatch in compressed sensing,” *IEEE Transactions on Signal Processing*, vol. 59, no. 5, pp. 2182–2195, 2011.

E. J. Cand[è]{}s and C. Fernandez-Granda, “Towards a mathematical theory of super-resolution,” *Communications on Pure and Applied Mathematics*, vol. 67, no. 6, pp. 906–956, 2014.

V. Chandrasekaran, B. Recht, P. A. Parrilo, and A. S. Willsky, “The convex geometry of linear inverse problems,” *Foundations of Computational Mathematics*, vol. 12, no. 6, pp. 805–849, 2012.

Z. Yang and L. Xie, “Exact joint sparse frequency recovery via optimization methods,” 2014. \[Online\]. Available: <http://arxiv.org/abs/1405.6585>

B. N. Bhaskar, G. Tang, and B. Recht, “Atomic norm denoising with applications to line spectral estimation,” *IEEE Transactions on Signal Processing*, vol. 61, no. 23, pp. 5987–5999, 2013.

E. J. Cand[è]{}s and C. Fernandez-Granda, “Super-resolution from noisy data,” *Journal of Fourier Analysis and Applications*, vol. 19, no. 6, pp. 1229–1254, 2013.

G. Tang, B. N. Bhaskar, and B. Recht, “Near minimax line spectral estimation,” *IEEE Transactions on Information Theory*, vol. 61, no. 1, pp. 499–512, 2015.

Z. Yang, L. Xie, and C. Zhang, “A discretization-free sparse and parametric approach for linear array signal processing,” *IEEE Transactions on Signal Processing*, vol. 62, no. 19, pp. 4959–4973, 2014.

Z. Yang and L. Xie, “On gridless sparse methods for line spectral estimation from complete and incomplete data,” *IEEE Transactions on Signal Processing*, vol. 63, no. 12, pp. 3139–3153, 2015.

J.-M. Azais, Y. De Castro, and F. Gamboa, “Spike detection from inaccurate samplings,” *Applied and Computational Harmonic Analysis*, vol. 38, no. 2, pp. 177–195, 2014.

L. Condat and A. Hirabayashi, “[Cadzow denoising upgraded: A new projection method for the recovery of Dirac pulses from noisy linear measurements]{},” *Sampling Theory in Signal and Image Processing*, vol. 14, no. 1, pp. 17–47, 2015.

Y. Chen and Y. Chi, “Robust spectral compressed sensing via structured matrix completion,” *IEEE Transactions on Information Theory*, vol. 60, no. 10, pp. 6576–6601, 2014.

Z. Tan, Y. C. Eldar, and A. Nehorai, “Direction of arrival estimation using co-prime arrays: A super resolution viewpoint,” *IEEE Transactions on Signal Processing*, vol. 62, no. 21, pp. 5565–5576, 2014.

K. V. Mishra, M. Cho, A. Kruger, and W. Xu, “Off-the-grid spectral compressed sensing with prior information,” in *IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, 2014, pp. 1010–1014.

Z. Lu, R. Ying, S. Jiang, P. Liu, and W. Yu, “Distributed compressed sensing off the grid,” *IEEE Signal Processing Letters*, vol. 22, no. 1, pp. 105–109, 2015.

M. S. Lobo, M. Fazel, and S. Boyd, “Portfolio optimization with linear and fixed transaction costs,” *Annals of Operations Research*, vol. 152, no. 1, pp. 341–365, 2007.

E. J. Candes, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighted $\ell_1$ minimization,” *Journal of Fourier Analysis and Applications*, vol. 14, no. 5-6, pp. 877–905, 2008.

D. Wipf and S. Nagarajan, “Iterative reweighted $\ell_1$ and $\ell_2$ methods for finding sparse solutions,” *IEEE Journal of Selected Topics in Signal Processing*, vol. 4, no. 2, pp. 317–329, 2010.

P. Stoica, D. Zachariah, and J. Li, “[Weighted SPICE: A unifying approach for hyperparameter-free sparse estimation]{},” *Digital Signal Processing*, vol. 33, pp. 1–12, 2014.

M. Fazel, H. Hindi, and S. P. Boyd, “[Log-det heuristic for matrix rank minimization with applications to Hankel and Euclidean distance matrices]{},” in *American Control Conference*, vol. 3, 2003, pp. 2156–2162.

K. Mohan and M. Fazel, “Iterative reweighted algorithms for matrix rank minimization,” *The Journal of Machine Learning Research*, vol. 13, no. 1, pp. 3441–3473, 2012.

E. Van Den Berg and M. Friedlander, “Theoretical and empirical results for recovery from multiple measurements,” *IEEE Transactions on Information Theory*, vol. 56, no. 5, pp. 2516–2527, 2010.

J. David, “Algorithms for analysis and design of robust controllers,” Ph.D. dissertation, Kat. Univ. Leuven, 1994.

B. D. Rao and K. Kreutz-Delgado, “An affine scaling methodology for best basis selection,” *IEEE Transactions on Signal Processing*, vol. 47, no. 1, pp. 187–200, 1999.

R. Chartrand and W. Yin, “Iteratively reweighted algorithms for compressive sensing,” in *IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, 2008, pp. 3869–3872.

S. P. Boyd and L. Vandenberghe, *Convex optimization*.1em plus 0.5em minus 0.4emCambridge University Press, London, 2004.

K.-C. Toh, M. J. Todd, and R. H. T[ü]{}t[ü]{}nc[ü]{}, “[SDPT3–a MATLAB software package for semidefinite programming, version 1.3]{},” *Optimization Methods and Software*, vol. 11, no. 1-4, pp. 545–581, 1999.

S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” *Foundations and Trends in Machine Learning*, vol. 3, no. 1, pp. 1–122, 2011.

[^1]: Manuscript August 2014; revised March and August 2015; accepted October 2015. Parts of this paper were presented at the 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Brisbane, Australia, April 2015 [@yang2015achieving].

    The authors are with the School of Electrical and Electronic Engineering, Nanyang Technological University, 639798, Singapore (e-mail: {yangzai, elhxie}@ntu.edu.sg).

[^2]: The frequency separation $\frac{4}{N}$ is sufficient but not necessary. Empirical studies in [@tang2012compressed] suggest that this value is about $\frac{1}{N}$ in the SMV case, while [@yang2014exact] shows that it also depends on other factors like $K$, $M$ and $L$.

[^3]: $\m{u}^*$ is called a cluster point of a vector-valued function $\m{u}(x)$ at $x=x_0$ if there exists a sequence $\lbra{x_n}_{n=1}^{+\infty}$, $\lim_{n\rightarrow+\infty}x_n=x_0$, satisfying that $\lim_{n\rightarrow+\infty}\m{u}(x_n)\rightarrow \m{u}^*$.
