Robust Spectral Compressed Sensing via Structured Matrix Completion

I Introduction

A large class of practical applications features high-dimensional signals that can be modeled or approximated by a superposition of spikes in the spectral (resp. time) domain, and involves estimation of the signal from its time (resp. frequency) domain samples. Examples include acceleration of medical imaging , target localization in radar and sonar systems , inverse scattering in seismic imaging , fluorescence microscopy , channel estimation in wireless communications , analog-to-digital conversion , etc. The data acquisition devices, however, are often limited by hardware and physical constraints, precluding sampling with the desired resolution. It is thus of paramount interest to reduce sensing complexity while retaining recovery accuracy.

In this paper, we investigate the spectral compressed sensing problem, which aims to recover a spectrally sparse signal from a small number of randomly observed time domain samples. The signal of interest $x\left(\boldsymbol{t}\right)$ with ambient dimension $n$ is assumed to be a weighted sum of multi-dimensional complex sinusoids at $r$ distinct frequencies $\{\boldsymbol{f}_{i}\in[0,1)^{K}:1\leq i\leq r\}$, where the underlying frequencies can assume any continuous values on the unit interval.

Spectral compressed sensing is closely related to the problem of harmonic retrieval, which seeks to extract the underlying frequencies of a signal from a collection of its time domain samples. Conventional methods for harmonic retrieval include Prony’s method , ESPRIT , the matrix pencil method , the Tufts and Kumaresan approach , the finite rate of innovation approach , etc. These methods routinely exploit the shift invariance of the harmonic structure, namely, a consecutive segment of time domain samples lies in the same subspace irrespective of the starting point of the segment. However, one weakness of these techniques if that they require prior knowledge of the model order, that is, the number of underlying frequency spikes of the signal or at least an estimate of it. Besides, these techniques heavily rely on the knowledge of the noise spectra, and are often sensitive against noise and outliers .

Another line of work is concerned with Compressed Sensing (CS) over a discrete domain, which suggests that it is possible to recover a signal even when the number of samples is far below its ambient dimension, provided that the signal enjoys a sparse representation in the transform domain. In particular, tractable algorithms based on convex surrogates become popular due to their computational efficiency and robustness against noise and outliers . Furthermore, they do not require prior information on the model order. Nevertheless, the success of CS relies on sparse representation or approximation of the signal of interest in a finite discrete dictionary, while the true parameters in many applications are actually specified in a continuous dictionary. The basis mismatch between the true frequencies and the discretized grid results in loss of sparsity due to spectral leakage along the Dirichlet kernel, and hence degeneration in the performance of conventional CS paradigms.

The performance of EMaC depends on an incoherence condition that depends only on the frequency locations regardless of the amplitudes of their respective coefficients. The incoherence measure is characterized by the reciprocal of the smallest singular value of some Gram matrix, which is defined by sampling the Dirichlet kernel at the wrap-around differences of all frequency pairs. The signal of interest is said to obey the incoherence condition if the Gram matrix is well conditioned, which arises over a broad class of spectrally sparse signals including but not restricted to signals with well-separated frequencies. We demonstrate that, under this incoherence condition, EMaC enables exact recovery from $\mathcal{O}(r\log^{4}n)$ random samplesThe standard notation $f(n)=\mathcal{O}\left(g(n)\right)$ means that there exists a constant $c>0$ such that $f(n)\leq cg(n)$; $f(n)=\Theta\left(g(n)\right)$ indicates that there are numerical constants $c_{1},c_{2}>0$ such that $c_{1}g(n)\leq f(n)\leq c_{2}g(n)$., and is stable against bounded noise. Moreover, EMaC admits perfect signal recovery from $\mathcal{O}(r^{2}\log^{3}n)$ random samples even when a constant proportion of the samples are corrupted with arbitrary magnitudes. Finally, numerical experiments validate our theoretical findings, and demonstrate the applicability of EMaC in super resolution.

Along the way, we provide theoretical guarantees for low-rank matrix completion of Hankel matrices and Toeplitz matrices, which is of great importance in control, natural language processing, and computer vision. To the best of our knowledge, our results provide the first theoretical guarantees for Hankel matrix completion that are close to the information theoretic limit.

I-B Connection and Comparison to Prior Work

The $K$-fold Hankel structure, which plays a central role in the EMaC algorithm, roots from the traditional spectral estimation technique named Matrix Enhancement Matrix Pencil (MEMP) for multi-dimensional harmonic retrieval. The conventional MEMP algorithm assumes fully observed equi-spaced time domain samples for estimation, and require prior knowledge on the model order. Cadzow’s denoising method also exploits the low-rank structure of the matrix pencil form for denoising line spectrum, but the method is non-convex and lacks performance guarantees.

More recently, Candï¿œs and Fernandez-Granda proposed a total-variation norm minimization algorithm to super-resolve a sparse signal from frequency samples at the low end of the spectrum. This algorithm allows accurate super-resolution when the point sources are sufficiently separated, and is stable against noise . Inspired by this approach, Tang et. al. then developed an atomic norm minimization algorithm for line spectral estimation from $\mathcal{O}(r\log r\log n)$ random time domain samples, which enables exact recovery when the frequencies are separated by at least $4/n$ with random amplitude phases. Similar performance guarantees are later established in for multi-dimensional frequencies. However, these results are established under a random signal model, i.e. the complex signs of the frequency spikes are assumed to be i.i.d. drawn from a uniform distribution. The robustness of the method against noise and outliers is not established either. In contrast, our approach yields deterministic conditions for multi-dimensional frequency models that guarantee perfect recovery with noiseless samples and are provably robust against noise and sparse corruptions. We will provide detailed comparison with the approach of Tang et. al. after we formally present our results. Numerical comparison will also be provided in Section V-C for the line spectrum model.

Our algorithm is inspired by recent advances of Matrix Completion (MC) , which aims at recovering a low-rank matrix from partial entries. It has been shown that exact recovery is possible via nuclear norm minimization, as soon as the number of observed entries exceeds the order of the information theoretic limit. This line of algorithms is also robust against noise and outliers , and allows exact recovery even in the presence of a constant portion of adversarially corrupted entries , which have found numerous applications in collaborative filtering , medical imaging , etc. Nevertheless, the theoretical guarantees of these algorithms do not apply to the more structured observation models associated with the proposed multi-fold Hankel structure. Consequently, direct application of existing MC results delivers pessimistic sample complexity, which far exceeds the degrees of freedom underlying the signal.

Preliminary results of this work have been presented in , where an additional strong incoherence condition was introduced that bore a similar role as the traditional strong incoherence parameter in MC but lacked physical interpretations. This paper removes this condition and further improves the sample complexity.

I-C Organization

The rest of the paper is organized as follows. The signal and sampling models are described in Section II. By restricting our attention to two-dimensional (2-D) frequency models, we present the enhanced matrix form and the associated structured matrix completion algorithms. The extension to multi-dimensional frequency models is discussed in Section III-C. The main theoretical guarantees are summarized in Section III, based on the incoherence condition introduced in Section III-A. We then discuss the extension to low-rank Hankel and Toeplitz matrix completion in Section IV. Section V presents the numerical validation of our algorithms. The proofs of Theorems 1 and 3 are based on duality analysis followed by a golfing scheme, which are supplied in Section VI and Section VII, respectively. Section VIII concludes the paper with a short summary of our findings as well as a discussion of potential extensions and improvements. Finally, the proofs of auxiliary lemmas supporting our results are deferred to the appendices.

II Model and Algorithm

Assume that the signal of interest $x\left(\boldsymbol{t}\right)$ can be modeled as a weighted sum of $K$-dimensional complex sinusoids at $r$ distinct frequencies $\boldsymbol{f}_{i}\in[0,1)^{K}$, $1\leq i\leq r$, i.e.

It is assumed throughout that the frequencies $\boldsymbol{f}_{i}$’s are normalized with respect to the Nyquist frequency of $x(\boldsymbol{t})$ and the time domain measurements are sampled at integer values. We denote by $d_{i}$’s the complex amplitudes of the associated coefficients, and $\left\langle\cdot,\cdot\right\rangle$ represents the inner product. For concreteness, our discussion is mainly devoted to a 2-D frequency model when $K=2$. This subsumes line spectral estimation as a special case, and indicates how to address multi-dimensional models. The algorithms for higher dimensional scenarios closely parallel the 2-D case, which will be briefly discussed in Section III-C.

Consider a data matrix $\boldsymbol{X}=[X_{k,l}]_{0\leq k<n_{1},0\leq l<n_{2}}$ of ambient dimension $n:=n_{1}n_{2}$, which is obtained by sampling the signal (1) on a uniform grid. From (1) each entry $X_{k,l}$ can be expressed as

where for any $i$ ($1\leq i\leq r$) we define

for some frequency pairs $\left\{\boldsymbol{f}_{i}=\left(f_{1i},f_{2i}\right)\mid 1\leq i\leq r\right\}$. We can then express $\boldsymbol{X}$ in a matrix form as follows

The above form (3) is sometimes referred to as the Vandemonde decomposition of $\boldsymbol{X}$.

Suppose that there exists a location set $\Omega$ of size $m$ such that the $X_{k,l}$ is observed if and only if $\left(k,l\right)\in\Omega$. It is assumed that $\Omega$ is sampled uniformly at random. Define $\mathcal{P}_{\Omega}(\boldsymbol{X})$ as the orthogonal projection of $\boldsymbol{X}$ onto the subspace of matrices that vanish outside $\Omega$. We aim at recovering $\boldsymbol{X}$ from $\mathcal{P}_{\Omega}(\boldsymbol{X})$.

II-B Matrix Enhancement

One might naturally attempt recovery by applying the low-rank MC algorithms , arguing that when $r$ is small, perfect recovery of $\boldsymbol{X}$ is possible from partial measurements since $\boldsymbol{X}$ is low rank if $r\ll\min\{n_{1},n_{2}\}$. Specifically, this corresponds to the following algorithm:

where $\left\|\boldsymbol{M}\right\|_{*}$ denotes the nuclear norm (or sum of all singular values) of a matrix $\boldsymbol{M}=[M_{k,l}]$. This is a convex relaxation paradigm with respect to rank minimization. However, naive MC algorithms require at least the order of $r\max\left(n_{1},n_{2}\right)\log\left(n_{1}n_{2}\right)$ samples in order to allow perfect recovery, which far exceeds the degrees of freedom (which is $\Theta\left(r\right)$) in our problem. What is worse, since the number $r$ of spectral spikes can be as large as $n_{1}n_{2}$, $\boldsymbol{X}$ might become full-rank once $r>\min\left(n_{1},n_{2}\right)$. This motivates us to seek other forms that better capture the harmonic structure.

In this paper, we adopt one effective enhanced form of $\boldsymbol{X}$ based on the following two-fold Hankel structure. The enhanced matrix $\boldsymbol{X}_{\text{e}}$ with respect to $\boldsymbol{X}$ is defined as a $k_{1}\times\left(n_{1}-k_{1}+1\right)$ block Hankel matrix

where $1\leq k_{2}\leq n_{2}$ is another pencil parameter. This enhanced form allows us to express each block asNote that the $l$th ($0\leq l<n_{1}$) row $\boldsymbol{X}_{l*}$ of $\boldsymbol{X}$ can be expressed as $\boldsymbol{X}_{l*}=\left[y_{1}^{l},\cdots,y_{r}^{l}\right]\boldsymbol{D}\boldsymbol{Z}^{\top}=\left[y_{1}^{l}d_{1},\cdots,y_{r}^{l}d_{r}\right]\boldsymbol{Z}^{\top},$ and hence we only need to find the Vandemonde decomposition for $\boldsymbol{X}_{0}$ and then replace $d_{i}$ by $y_{i}^{l}d_{i}$.

where $\boldsymbol{Z}_{\text{L}}$, $\boldsymbol{Z}_{\text{R}}$ and $\boldsymbol{Y}_{\text{d}}$ are defined respectively as

Substituting (10) into (8) yields the following:

where $\boldsymbol{E}_{\text{L}}$ and $\boldsymbol{E}_{\text{R}}$ span the column and row space of $\boldsymbol{X}_{\text{e}}$, respectively. This immediately implies that $\boldsymbol{X}_{\text{e}}$ is low-rank, i.e.

This form is inspired by the traditional matrix pencil approach proposed in to estimate harmonic frequencies if all entries of $\boldsymbol{X}$ are available. Thus, one can extract all underlying frequencies of $\boldsymbol{X}$ using methods proposed in , as long as $\boldsymbol{X}$ can be faithfully recovered.

II-C The EMaC Algorithm in the Absence of Noise

We then attempt recovery through the following Enhancement Matrix Completion (EMaC) algorithm:

where $\boldsymbol{M}_{\text{e}}$ denotes the enhanced form of $\boldsymbol{M}$. In other words, EMaC minimizes the nuclear norm of the enhanced form over all matrices compatible with the samples. This convex program can be rewritten into a semidefinite program (SDP)

which can be solved using off-the-shelf solvers in a tractable manner (see, e.g., ). It is worth mentioning that EMaC has a similar computational complexity as the atomic norm minimization method when restricted to the 1-D frequency model.

Careful readers will remark that the performance of EMaC must depend on the choices of the pencil parameters $k_{1}$ and $k_{2}$. In fact, if we define a quantity

II-D The Noisy-EMaC Algorithm with Bounded Noise

In practice, measurements are often contaminated by a certain amount of noise. To make our model and algorithm more practically applicable, we replace our measurements by $\boldsymbol{X}^{\text{o}}=[{X}_{k,l}^{\text{o}}]_{0\leq k<n_{1},0\leq l<n_{2}}$ through the following noisy model

where ${X}_{k,l}^{\text{o}}$ is the observed $(k,l)$-th entry, and $\boldsymbol{N}=[N_{k,l}]_{0\leq k<n_{1},0\leq l<n_{2}}$ denotes some unknown noise. We assume that the noise magnitude is bounded by a known amount $\left\|\mathcal{P}_{\Omega}\left(\boldsymbol{N}\right)\right\|_{\text{F}}\leq\delta$, where $\left\|\cdot\right\|_{\text{F}}$ denotes the Frobenius norm. In order to adapt our algorithm to such noisy measurements, one wishes that small perturbation in the measurements should result in small variation in the estimate. Our algorithm is then modified as follows

That said, the algorithm searches for a candidate with minimum nuclear norm among all signals close to the measurements.

II-E The Robust-EMaC Algorithm with Sparse Outliers

An outlier is a data sample that can deviate arbitrarily from the true data point. Practical data samples one collects may contain a certain portion of outliers due to abnormal behavior of data acquisition devices such as amplifier saturation, sensor failures, and malicious attacks. A desired recovery algorithm should be able to automatically prune all outliers even when they corrupt up to a constant portion of all data samples.

Specifically, suppose that our measurements $\boldsymbol{X}^{\text{o}}$ are given by

where ${X}_{k,l}^{\text{o}}$ is the observed $(k,l)$-th entry, and $\boldsymbol{S}=[S_{k,l}]_{0\leq k<n_{1},0\leq l<n_{2}}$ denotes the outliers, which is assumed to be a sparse matrix supported on some location set $\Omega^{\text{dirty}}\subseteq\Omega$. The sampling model is formally described as follows.

Suppose that $\Omega$ is obtained by sampling $m$ entries uniformly at random, and define $\rho:=\frac{m}{n_{1}n_{2}}$.

Conditioning on $(k,l)\in\Omega$, the events $\left\{(k,l)\in\Omega^{\text{dirty}}\right\}$ are independent with conditional probability

for some small constant corruption fraction $0<\tau<1$.

Define $\Omega^{\text{clean}}:=\Omega\backslash\Omega^{\text{dirty}}$ as the location set of uncorrupted measurements.

EMaC is then modified as follows to accommodate sparse outliers:

II-F Notations

Before continuing, we introduce a few notations that will be used throughout. Let the singular value decomposition (SVD) of $\boldsymbol{X}_{\text{e}}$ be $\boldsymbol{X}_{\text{e}}=\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{V}^{*}$. Denote by

the tangent space with respect to $\boldsymbol{X}_{\text{e}}$, and $T^{\perp}$ the orthogonal complement of $T$. Denote by $\mathcal{P}_{U}$ (resp. $\mathcal{P}_{V}$, $\mathcal{P}_{T}$) the orthogonal projection onto the column (resp. row, tangent) space of $\boldsymbol{X}_{\text{e}}$, i.e. for any $\boldsymbol{M}$,

We let $\mathcal{P}_{T^{\perp}}=\mathcal{I}-\mathcal{P}_{T}$ be the orthogonal complement of $\mathcal{P}_{T}$, where $\mathcal{I}$ denotes the identity operator.

On the other hand, we denote by $\Omega_{\text{e}}(k,l)$ the set of locations of the enhanced matrix $\boldsymbol{X}_{\text{e}}$ containing copies of $X_{k,l}$. Due to the Hankel or multi-fold Hankel structures, one can easily verify the following: each location set $\Omega_{\text{e}}(k,l)$ contains at most one index in any given row of the enhanced form, and at most one index in any given column. For each $\left(k,l\right)\in[n_{1}]\times[n_{2}]$, we use $\boldsymbol{A}_{(k,l)}$ to denote a basis matrix that extracts the average of all entries in $\Omega_{\text{e}}\left(k,l\right)$. Specifically,

III Main Results

This section delivers the following encouraging news: under mild incoherence conditions, EMaC enables faithful signal recovery from a minimal number of time-domain samples, even when the samples are contaminated by bounded noise or a constant portion of arbitrary outliers.

In general, matrix completion from a few entries is hopeless unless the underlying structure is sufficiently uncorrelated with the observation basis. This inspires us to introduce certain incoherence measures. To this end, we define the 2-D Dirichlet kernel as

where $\boldsymbol{f}=(f_{1},f_{2})\in[0,1)^{2}$. Fig. 1 (a) illustrates the amplitude of $\mathcal{D}(k_{1},k_{2},\boldsymbol{f})$ when $k_{1}=k_{2}=6$. The value of $|\mathcal{D}(k_{1},k_{2},\boldsymbol{f})|$ decays inverse proportionally with respect to the frequency $\boldsymbol{f}$. Set $\boldsymbol{G}_{\text{L}}$ and $\boldsymbol{G}_{\text{R}}$ to be two $r\times r$ Gram matrices such that their entries are specified respectively by

where the difference $\boldsymbol{f}_{i}-\boldsymbol{f}_{l}$ is understood as the wrap-around distance in the interval $[-1/2,1/2)^{2}$. Simple manipulation reveals that

Our incoherence measure is then defined as follows.

A matrix $\boldsymbol{X}$ is said to obey the incoherence property with parameter $\mu_{1}$ if

The incoherence measure $\mu_{1}$ only depends on the locations of the frequency spikes, irrespective of the amplitudes of their respective coefficients. The signal is said to satisfy the incoherence condition if $\mu_{1}$ scales as a small constant, which occurs when $\boldsymbol{G}_{\text{L}}$ and $\boldsymbol{G}_{\text{R}}$ are both well-conditioned. Our incoherence condition naturally requires certain separation among all frequency pairs, as when two frequency spikes are closely located, $\mu_{1}$ gets undesirably large. As shown in [43, Theorem 2], a separation of about $2/n$ for line spectrum is sufficient to guarantee the incoherence condition to hold. However, it is worth emphasizing that such strict separation is not necessary as required in , and thereby our incoherence condition is applicable to a broader class of spectrally sparse signals.

To give the reader a flavor of the incoherence condition, we list two examples below. For ease of presentation, we assume below 2-D frequency models with $n_{1}=n_{2}$. Note, however, that the asymmetric cases and general $K$-dimensional frequency models can be analyzed in the same manner.

Random frequency locations: suppose that the $r$ frequencies are generated uniformly at random, then the minimum pairwise separation can be crudely bounded by $\Theta\left(\frac{1}{r^{2}\log n_{1}}\right)$. If $n_{1}\gg r^{2.5}\log n_{1}$, then a crude bound reveals that $\forall i_{1}\neq i_{2},$

holds with high probability, indicating that the off-diagonal entries of $\boldsymbol{G}_{\text{L}}$ and $\boldsymbol{G}_{\text{R}}$ are much smaller than $1/r$ in magnitude. Simple manipulation then allows us to conclude that $\sigma_{\min}\left(\boldsymbol{G}_{\text{L}}\right)$ and $\sigma_{\min}\left(\boldsymbol{G}_{\text{R}}\right)$ are bounded below by positive constants. Fig. 1 (b) shows the minimum eigenvalue of $\boldsymbol{G}_{\text{L}}$ for different $k=k_{1}=k_{2}=6,36,72$ when the spikes are randomly generated and the number of spikes is given as the sparsity level. The minimum eigenvalue of $\boldsymbol{G}_{\text{L}}$ gets closer to one as $k$ grows, confirming our argument.

Small perturbation off the grid: suppose that all frequencies are within a distance at most $\frac{1}{n_{1}r^{1/4}}$ from some grid points $\left(\frac{l_{1}}{k_{1}},\frac{l_{2}}{k_{2}}\right)$ $\left(0\leq l_{1}<k_{1},0\leq l_{2}<k_{2}\right)$. One can verify that $\forall i_{1}\neq i_{2}$,

and hence the magnitude of all off-diagonal entries of $\boldsymbol{G}_{\text{L}}$ and $\boldsymbol{G}_{\text{R}}$ are no larger than $1/(4r)$. This immediately suggests that $\sigma_{\min}\left(\boldsymbol{G}_{\text{L}}\right)$ and $\sigma_{\min}\left(\boldsymbol{G}_{\text{R}}\right)$ are lower bounded by $3/4$.

Note, however, that the class of incoherent signals are far beyond the ones discussed above.

III-B Theoretical Guarantees

With the above incoherence measure, the main theoretical guarantees are provided in the following three theorems each accounting for a distinct data model: 1) noiseless measurements, 2) measurements contaminated by bounded noise, and 3) measurements corrupted by a constant proportion of arbitrary outliers.

Exact recovery is possible from a minimal number of noise-free samples, as asserted in the following theorem.

Let $\boldsymbol{X}$ be a data matrix of form (3), and $\Omega$ the random location set of size $m$. Suppose that the incoherence property (23) holds and that all measurements are noiseless. Then there exists a universal constant $c_{1}>0$ such that $\boldsymbol{X}$ is the unique solution to EMaC with probability exceeding $1-\left(n_{1}n_{2}\right)^{-2}$, provided that

Theorem 1 asserts that under some mild deterministic incoherence condition such that $\mu_{1}$ scales as a small constant, EMaC admits prefect recovery as soon as the number of measurements exceeds $\mathcal{O}(r\log^{4}\left(n_{1}n_{2}\right))$. Since there are $\Theta(r)$ degrees of freedom in total, the lower bound should be no smaller than $\Theta(r)$. This demonstrates the orderwise optimality of EMaC except for a logarithmic gap. We note, however, that the polylog factor might be further refined via finer tuning of concentration of measure inequalities.

It is worth emphasizing that while we assume random observation models, the data model is assumed deterministic. This differs significantly from , which relies on randomness in both the observation model and the data model. In particular, our theoretical performance guarantees rely solely on the frequency locations irrespective of the associated amplitudes. In contrast, the results in require the phases of all frequency spikes to be i.i.d. drawn in a uniform manner in addition to a separation condition.

III-B2 Stable Recovery in the Presence of Bounded Noise

Our method enables stable recovery even when the time domain samples are noisy copies of the true data. Here, we say the recovery is stable if the solution of Noisy-EMaC is close to the ground truth in proportion to the noise level. To this end, we provide the following theorem, which is a counterpart of Theorem 1 in the noisy setting, whose proof is inspired by .

with probability exceeding $1-(n_{1}n_{2})^{-2}$.

Theorem 2 reveals that the recovered enhanced matrix (which contains $\Theta(n_{1}^{2}n_{2}^{2})$ entries) is close to the true enhanced matrix at high SNR. In particular, the average entry inaccuracy of the enhanced matrix is bounded above by $\mathcal{O}(n_{1}^{3}n_{2}^{3}\delta)$, amplified by the subsampling factor. In practice, one is interested in an estimate of $\boldsymbol{X}$, which can be obtained naively by randomly selecting an entry in $\Omega_{\text{e}}(k,l)$ as $\hat{X}_{k,l}$, then we have

This yields that the per-entry noise of $\hat{\boldsymbol{X}}$ is about $\mathcal{O}(n_{1}^{2.5}n_{2}^{2.5}\delta)$, which is further amplified due to enhancement by a factor of $n_{1}n_{2}$. However, this factor arises from an analysis artifact due to our simple strategy to deduce $\hat{\boldsymbol{X}}$ from $\hat{\boldsymbol{X}_{\text{e}}}$, and may be elevated. We note that in numerical experiments, Noisy-EMaC usually generates much better estimates, usually by a polynomial factor. The practical applicability will be illustrated in Section V.

It is worth mentioning that to the best of our knowledge, our result is the first stability result with partially observed data for spectral compressed sensing off the grid. While the atomic norm approach is near-minimax with full data , it is not clear how it performs with partially observed data.

III-B3 Robust Recovery in the Presence of Sparse Outliers

Interestingly, Robust-EMaC can provably tolerate a constant portion of arbitrary outliers. The theoretical performance is formally summarized in the following theorem.

Let $\boldsymbol{X}$ be a data matrix with matrix form (3), and $\Omega$ a random location set of size $m$. Set $\lambda=\frac{1}{\sqrt{m\log\left(n_{1}n_{2}\right)}}$, and assume $\tau\leq 0.1$ is some small positive constant. Then there exist a numerical constant $c_{1}>0$ depending only on $\tau$ such that if (23) holds and

then Robust-EMaC is exact, i.e. the minimizer $(\hat{\boldsymbol{M}},\hat{\boldsymbol{S}})$ satisfies $\hat{\boldsymbol{M}}=\boldsymbol{X}$, with probability exceeding $1-(n_{1}n_{2})^{-2}$.

Note that $\tau\leq 0.1$ is not a critical threshold. In fact, one can prove the same theorem for a larger $\tau$ (e.g. $\tau\leq 0.25$) with a larger absolute constant $c_{1}$. However, to allow even larger $\tau$ (e.g. in the regime where $\tau\geq 50\%$), we need the sparse components exhibit random sign patterns.

Theorem 3 specifies a candidate choice of the regularization parameter $\lambda$ that allows recovery from a few samples, which only depends on the size of $\Omega$ but is otherwise parameter-free. In practice, however, $\lambda$ may better be selected via cross validation. Furthermore, Theorem 3 demonstrates the possibility of robust recovery under a constant proportion of sparse corruptions. Under the same mild incoherence condition as for Theorem 1, robust recovery is possible from $\mathcal{O}\left(r^{2}\log^{3}\left(n_{1}n_{2}\right)\right)$ samples, even when a constant proportion of the samples are arbitrarily corrupted. As far as we know, this provides the first theoretical guarantees for separating sparse measurement corruptions in the off-grid compressed sensing setting.

III-C Extension to Higher-Dimensional and Damping Frequency Models

By letting $n_{2}=1$ the above 2-D frequency model reverts to the line spectrum model. The EMaC algorithm and the main results immediately extend to higher dimensional frequency models without difficulty. In fact, for $K$-dimensional frequency models, one can arrange the original data into a $K$-fold Hankel matrix of rank at most $r$. For instance, consider a 3-D model such that

An enhanced form can be defined as a 3-fold Hankel matrix such that

where $\boldsymbol{X}_{i,\text{e}}$ denotes the 2-D enhanced form of the matrix consisting of all entries $X_{l_{1},l_{2},l_{3}}$ obeying $l_{3}=i$. One can verify that $\boldsymbol{X}_{\text{e}}$ is of rank at most $r$, and can thereby apply EMaC on the 3-D enhanced form. To summarize, for $K$-dimensional frequency models, EMaC (resp. Noisy-EMaC, Robust-EMaC) searches over all $K$-fold Hankel matrices that are consistent with the measurements. The theoretical performance guarantees can be similarly extended by defining the respective Dirichlet kernel in 3-D and the coherence measure. In fact, all our analyses can be extended to handle damping modes, when the frequencies are not of time-invariant amplitudes. We omit the details for conciseness.

IV Structured Matrix Completion

One problem closely related to our method is completion of multi-fold Hankel matrices from a small number of entries. While each spectrally sparse signal can be mapped to a low-rank multi-fold Hankel matrix, it is not clear whether all multi-fold Hankel matrices of rank $r$ can be written as the enhanced form of a signal with spectral sparsity $r$. Therefore, one can think of recovery of multi-fold Hankel matrices as a more general problem than the spectral compressed sensing problem. Indeed, Hankel matrix completion has found numerous applications in system identification , natural language processing , computer vision , magnetic resonance imaging , etc.

There has been several work concerning algorithms and numerical experiments for Hankel matrix completions . However, to the best of our knowledge, there has been little theoretical guarantee that addresses directly Hankel matrix completion. Our analysis framework can be straightforwardly adapted to the general $K$-fold Hankel matrix completions. Below we present the performance guarantee for the two-fold Hankel matrix completion without loss of generality. Notice that we need to modify the definition of $\mu_{1}$ as stated in the following theorem.

Condition (27) requires that the left and right singular vectors are sufficiently uncorrelated with the observation basis. In fact, condition (27) is a weaker assumption than (23).

It is worth mentioning that a low-rank Hankel matrix can often be converted to its low-rank Toeplitz counterpart, by reversely ordering all rows of the Hankel matrix. Both Hankel and Toeplitz matrices are effective forms that capture the underlying harmonic structures. Our results and analysis framework extend to low-rank Toeplitz matrix completion problem without difficulty.

V Numerical Experiments

In this section, we present numerical examples to evaluate the performance of EMaC and its variants under different scenarios. We further examine the application of EMaC in image super resolution. Finally, we propose an extension of singular value thresholding (SVT) developed by Cai et. al. that exploits the multi-fold Hankel structure to handle larger scale data sets.

To evaluate the practical ability of the EMaC algorithm, we conducted a series of numerical experiments to examine the phase transition for exact recovery. Let $n_{1}=n_{2}$, and we take $k_{1}=k_{2}=\lceil(n_{1}+1)/2\rceil$ which corresponds to the smallest $c_{\text{s}}$. For each $(r,m)$ pair, 100 Monte Carlo trials were conducted. We generated a spectrally sparse data matrix $\boldsymbol{X}$ by randomly generating $r$ frequency spikes in $[0,1)\times[0,1)$, and sampled a subset $\Omega$ of size $m$ entries uniformly at random. The EMaC algorithm was conducted using the convex programming modeling software CVX with the interior-point solver SDPT3 . Each trial is declared successful if the normalized mean squared error (NMSE) satisfies $\|\hat{\boldsymbol{X}}-\boldsymbol{X}\|_{\text{F}}/\|\boldsymbol{X}\|_{\text{F}}\leq 10^{-3}$, where $\hat{\boldsymbol{X}}$ denotes the estimate returned by EMaC. The empirical success rate is calculated by averaging over 100 Monte Carlo trials.

V-B Stable Recovery from Noisy Data

Fig. 3 further examines the stability of the proposed algorithm by performing Noisy-EMaC with respect to different parameter $\delta$ on a noise-free dataset of $r=4$ complex sinusoids with $n_{1}=n_{2}=11$. The number of random samples is $m=50$. The reconstructed NMSE grows approximately linear with respect to $\delta$, validating the stability of the proposed algorithm.

V-C Comparison with Existing Approaches for Line Spectrum Estimation

Suppose that we randomly observe $64$ entries of an $n$-dimensional vector ($n=127$) composed of $r=4$ modes. For such 1-D signals, we compare EMaC with the atomic norm approach as well as basis pursuit assuming a grid of size $2^{12}$. For the atomic norm and the EMaC algorithm, the modes are recovered via linear prediction using the recovered data . Fig. 4 demonstrates the recovery of mode locations for three cases, namely when (a) all the modes are on the DFT grid along the unit circle; (b) all the modes are on the unit circle except two closely located modes that are off the presumed grid; (c) all the modes are on the unit circle except that one of the two closely located modes is a damping mode with amplitude $0.99$. In all cases, the EMaC algorithm successfully recovers the underlying modes, while the atomic norm approach fails to recover damping modes, and basis pursuit fails with both off-the-grid modes and damping modes.

We further compare the phase transition of the EMaC algorithm and the atomic norm approach in for line spectrum estimation. We assume a 1-D signal of length $n=n_{1}=127$ and the pencil parameter $k_{1}$ of EMaC is chosen to be $64$. The phase transition experiments are conducted in the same manner as Fig. 2. In the first case, the spikes are generated randomly as Fig. 2 on a unit circle; in the second case, the spikes are generated until a separation condition is satisfied $\Delta:=\min_{i_{1}\neq i_{2}}|f_{i_{1}}-f_{i_{2}}|\geq 1.5/n$. Fig. 5 (a) and (b) illustrate the phase transition of EMaC and the atomic norm approach when the frequencies are randomly generated without imposing the separation condition. The performance of the atomic norm approach degenerates severely when the separation condition is not met; on the other hand, the EMaC gives a sharp phase transition similar to the 2D case. When the separation condition is imposed, the phase transition of the atomic norm approach greatly improves as shown in Fig. 5 (c), while the phase transition of EMaC still gives similar performance as in Fig. 5 (a) (We omit the actual phase transition in this case.) However, it is worth mentioning that when the sparsity level is relatively high, the required separation condition is in general difficult to be satisfied in practice. In comparison, EMaC is less sensitive to the separation requirement.

V-D Robust Line Spectrum Estimation

Consider the problem of line spectrum estimation, where the time domain measurements are contaminated by a constant portion of outliers. We conducted a series of Monte Carlo trials to illustrate the phase transition for perfect recovery of the ground truth. The true data $\boldsymbol{X}$ is assumed to be a $125$-dimensional vector, where the locations of the underlying frequencies are randomly generated. The simulations were carried out again using CVX with SDPT3.

Fig. 6(a) illustrates the phase transition for robust line spectrum estimation when $10\%$ of the entries are corrupted, which showcases the tradeoff between the number $m$ of measurements and the recoverable spectral sparsity level $r$. One can see from the plot that $m$ is approximately linear in $r$ on the phase transition curve even when 10% of the measurements are corrupted, which validates our finding in Theorem 3. Fig. 6(b) illustrates the success rate of exact recovery when we obtain samples for all entry locations. This plot illustrates the tradeoff between the spectral sparsity level and the number of outliers when all entries of the corrupted $\boldsymbol{X}^{o}$ are observed. It can be seen that there is a large region where exact recovery can be guaranteed, demonstrating the power of our algorithms in the presence of sparse outliers.

V-E Synthetic Super Resolution

V-F Singular Value Thresholding for EMaC

The above Monte Carlo experiments were conducted using the advanced SDP solver SDPT3. This solver and many other popular ones (e.g. SeDuMi) are based on interior point methods, which are typically inapplicable to large-scale data. In fact, SDPT3 fails to handle an $n\times n$ data matrix when $n$ exceeds 19, which corresponds to a $100\times 100$ enhanced matrix.

One alternative for large-scale data is the first-order algorithms tailored to matrix completion problems, e.g. the singular value thresholding (SVT) algorithm . We propose a modified SVT algorithm in Algorithm 1 to exploit the Hankel structure.

In particular, two operators are defined as follows:

$\mathcal{D}_{\tau_{t}}(\cdot)$ in Algorithm 1 denotes the singular value shrinkage operator. Specifically, if the SVD of $\boldsymbol{X}$ is given by $\boldsymbol{X}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{*}$ with $\boldsymbol{\Sigma}=\text{diag}\left(\left\{\sigma_{i}\right\}\right)$, then

where $\tau_{t}>0$ is the soft-thresholding level.

In the $K$-dimensional frequency model, $\mathcal{H}_{\boldsymbol{X}^{\text{o}}}(\boldsymbol{Q}_{t})$ denotes the projection of $\boldsymbol{Q}_{t}$ onto the subspace of enhanced matrices (i.e. $K$-fold Hankel matrices) that are consistent with the observed entries.

Consequently, at each iteration, a pair $\left(\boldsymbol{Q}_{t},\boldsymbol{M}_{t}\right)$ is produced by first performing singular value shrinkage and then projecting the outcome onto the space of $K$-fold Hankel matrices that are consistent with observed entries.

The key parameter that one needs to tune is the threshold $\tau_{t}$. Unfortunately, there is no universal consensus regarding how to tweak the threshold for SVT type of algorithms. One suggested choice is $\tau_{t}=0.1\sigma_{\max}\left(\boldsymbol{M}_{t}\right)/\left\lceil\frac{t}{10}\right\rceil$, which works well based on our empirical experiments.

Fig. 8 illustrates the performance of Algorithm 1. We generated a true $101\times 101$ data matrix $\boldsymbol{X}$ through a superposition of $30$ random complex sinusoids, and revealed 5.8% of the total entries (i.e. $m=600$) uniformly at random. The noise was i.i.d. Gaussian giving a signal-to-noise amplitude ratio of $10$. The reconstructed vectorized signal is superimposed on the ground truth in Fig. 8. The normalized reconstruction error was $\|\hat{\boldsymbol{X}}-\boldsymbol{X}\|_{\text{F}}/\left\|\boldsymbol{X}\right\|_{\text{F}}=0.1098$, validating the stability of our algorithm in the presence of noise.

VI Proof of Theorems 1 and 4

EMaC has similar spirit as the well-known matrix completion algorithms , except that we impose Hankel and multi-fold Hankel structures on the matrices. While has presented a general sufficient condition for exact recovery (see [31, Theorem 3]), the basis in our case does not exhibit desired coherence properties as required in , and hence these results cannot deliver informative estimates when applied to our problem. Nevertheless, the beautiful golfing scheme introduced in lays the foundation of our analysis in the sequel. We also note that the analyses adopted in rely on a desired joint incoherence property on $\boldsymbol{U}\boldsymbol{V}^{*}$, which has been shown to be unnecessary .

For concreteness, the analyses in this paper focus on recovering harmonically sparse signals as stated in Theorem 1, since proving Theorem 1 is slightly more involved than proving Theorem 4. We note, however, that our analysis already entails all reasoning required for establishing Theorem 4.

Denote by $\mathcal{A}_{\left(k,l\right)}\left(\boldsymbol{M}\right)$ the projection of $\boldsymbol{M}$ onto the subspace spanned by $\boldsymbol{A}_{(k,l)}$, and define the projection operator onto the space spanned by all $\boldsymbol{A}_{(k,l)}$ and its orthogonal complement as

There are two common ways to describe the randomness of $\Omega$: one corresponds to sampling without replacement, and another concerns sampling with replacement (i.e. $\Omega$ contains $m$ indices $\left\{\boldsymbol{a}_{i}\in[n_{1}]\times[n_{2}]:1\leq i\leq m\right\}$ that are i.i.d. generated). As discussed in [31, Section II.A], while both situations result in the same order-wide bounds, the latter situation admits simpler analysis due to independence. Therefore, we will assume that $\Omega$ is a multi-set (possibly with repeated elements) and $a_{i}$’s are independently and uniformly distributed throughout the proofs of this paper, and define the associated operators as

We also define another projection operator $\mathcal{A}^{\prime}_{\Omega}$ similar to (29), but with the sum extending only over distinct samples. Its complement operator is defined as $\mathcal{A}^{\prime}_{\Omega^{\perp}}:=\mathcal{A}-\mathcal{A}^{\prime}_{\Omega}$. Note that $\mathcal{A}_{\Omega}\left(\boldsymbol{M}\right)=0$ is equivalent to $\mathcal{A}^{\prime}_{\Omega}(\boldsymbol{M})=0$. With these definitions, EMaC can be rewritten as the following general matrix completion problem:

To prove exact recovery of convex optimization, it suffices to produce an appropriate dual certificate, as stated in the following lemma.

Consider a multi-set $\Omega$ that contains $m$ random indices. Suppose that the sampling operator $\mathcal{A}_{\Omega}$ obeys

If there exists a matrix $\boldsymbol{W}$ satisfying

Condition (31) will be analyzed in Section VI-B, while a dual certificate $\boldsymbol{W}$ will be constructed in Section VI-C. The validity of $\boldsymbol{W}$ as a dual certificate will be established in Sections VI-C - VI-E. These are the focus of the remaining section.

Lemma 1 requires that $\mathcal{A}_{\Omega}$ be sufficiently incoherent with respect to the tangent space $T$. The following lemma quantifies the projection of each $\boldsymbol{A}_{(k,l)}$ onto the subspace $T$.

for all $\left(k,l\right)\in[n_{1}]\times[n_{2}]$. For any $\boldsymbol{a},\boldsymbol{b}\in[n_{1}]\times[n_{2}]$, one has

Recognizing that (35) is the same as (27), the following proof also establishes Theorem 4. Note that Lemma 2 immediately leads to

As long as (37) holds, the fluctuation of $\mathcal{P}_{T}\mathcal{A}_{\Omega}\mathcal{P}_{T}$ can be controlled reasonably well, as stated in the following lemma. This justifies Condition (31) as required by Lemma 1.

Suppose that (37) holds. Then for any small constant $0<\epsilon\leq\frac{1}{2}$, one has

VI-C Construction of Dual Certificates

Now we are in a position to construct the dual certificate, for which we will employ the golfing scheme introduced in . Suppose that we generate $j_{0}$ independent random location multi-sets $\Omega_{i}$ ($1\leq i\leq j_{0}$), each containing $\frac{m}{j_{0}}$ i.i.d. samples. This way the distribution of $\Omega$ is the same as $\Omega_{1}\cup\Omega_{2}\cup\cdots\cup\Omega_{j_{0}}$ . Note that $\Omega_{i}$’s correspond to sampling with replacement. Let

represent the undersampling factors of $\Omega$ and $\Omega_{i}$, respectively.

Consider a small constant $\epsilon<\frac{1}{e}$, and pick $j_{0}:=3\log_{\frac{1}{\epsilon}}n_{1}n_{2}$. The construction of the dual matrix $\boldsymbol{W}$ then proceeds as follows:

We will establish that $\boldsymbol{W}$ is a valid dual certificate by showing that $\boldsymbol{W}$ satisfies the conditions stated in Lemma 1, which we now proceed step by step.

lie within the subspace of matrices supported on $\Omega$ or the subspace $\mathcal{A}^{\perp}$. This validates that $\mathcal{A}^{\prime}_{\Omega^{\perp}}\left(\boldsymbol{W}\right)=0$, as required in (32).

Secondly, the recursive construction procedure of $\boldsymbol{F}_{i}$ allows us to write

This allows us to bound $\left\|\mathcal{P}_{T}\left(\boldsymbol{F}_{i}\right)\right\|_{\text{F}}$ as

Finally, it remains to be shown that $\left\|\mathcal{P}_{T^{\perp}}\left(\boldsymbol{W}\right)\right\|\leq\frac{1}{2}$, which we will establish in the next two subsections. In particular, we first introduce two key metrics and characterize their relationships in Section VI-D. These metrics are crucial in bounding $\left\|\mathcal{P}_{T^{\perp}}\left(\boldsymbol{W}\right)\right\|$, which will be the focus of Section VI-E.

VI-D Two Metrics and Key Lemmas

In this subsection, we introduce the following two norms

Based on these two metrics, we can derive several technical lemmas which, taken collectively, allow us to control $\left\|\mathcal{P}_{T^{\perp}}\left(\boldsymbol{W}\right)\right\|$. Specifically, these lemmas characterize the mutual dependence of three norms $\left\|\cdot\right\|$, $\left\|\cdot\right\|_{\mathcal{A},2}$ and $\left\|\cdot\right\|_{\mathcal{A},\infty}$.

For any given matrix $\boldsymbol{M}$, there exists some numerical constant $c_{2}>0$ such that

with probability at least $1-\left(n_{1}n_{2}\right)^{-10}$.

Assume that there exists a quantity $\mu_{5}$ such that

For any given matrix $\boldsymbol{M}$, with probability exceeding $1-\left(n_{1}n_{2}\right)^{-10}$,

For any given matrix $\boldsymbol{M}\in T$, there is some absolute constant $c_{4}>0$ such that

with probability exceeding $1-\left(n_{1}n_{2}\right)^{-10}$.

Lemma 5 combined with Lemma 6 gives rise to the following inequality. Consider any given matrix $\boldsymbol{M}\in T$. Applying the bounds (46) and (47), one can derive

with probability exceeding $1-\left(n_{1}n_{2}\right)^{-10}$, where $c_{5}=\max\left\{c_{3},c_{4}\right\}$. This holds under the hypothesis (45).

Now we are ready to show how we may combine the above lemmas to develop an upper bound on $\left\|\mathcal{P}_{T^{\perp}}\left(\boldsymbol{W}\right)\right\|$. By construction, one has

Each summand can be bounded above as follows

Since $\boldsymbol{F}_{0}=\boldsymbol{U}\boldsymbol{V}^{*}$, it remains to control $\left\|\boldsymbol{U}\boldsymbol{V}^{*}\right\|_{\mathcal{A},\infty}$ and $\left\|\boldsymbol{U}\boldsymbol{V}^{*}\right\|_{\mathcal{A},2}$. We have the following lemma.

With the incoherence measure $\mu_{1}$, one can bound

and for any $(\alpha,\beta)\in\left[n_{1}\right]\times\left[n_{2}\right]$,

In particular, the bound (55) translates into

Substituting (53) and (54) into (52) gives

as required. So far, we have successfully verified that with high probability, $\boldsymbol{W}$ is a valid dual certificate, and hence by Lemma 1 the solution to EMaC is exact and unique.

VII Proof of Theorem 3

The algorithm Robust-EMaC is inspired by the well-known robust principal component analysis that seeks a decomposition of low-rank plus sparse matrices, except that we impose multi-fold Hankel structures on both the low-rank and sparse matrices. Following similar spirit as to the proof of Theorem 1, the proof here is based on duality analysis, and relies on the golfing scheme to construct a valid dual certificate.

In this section, we prove the results for a slightly different sampling model as follows.

The location multi-set $\Omega^{\text{clean}}$ of observed uncorrupted entries is generated by sampling $\left(1-\tau\right)\rho n_{1}n_{2}$ i.i.d. entries uniformly at random.

The location multi-set $\Omega$ of observed entries is generated by sampling $\rho n_{1}n_{2}$ i.i.d. entries uniformly at random, with the first $\left(1-\tau\right)\rho n_{1}n_{2}$ entries coming from $\Omega^{\text{clean}}$.

The location set $\Omega^{\text{dirty}}$ of observed corrupted entries is given by $\Omega^{\prime}\backslash\Omega^{\text{clean}^{\prime}}$, where $\Omega^{\prime}$ and $\Omega^{\text{clean}^{\prime}}$ denote the sets of distinct entry locations in $\Omega$ and $\Omega^{\text{clean}}$, respectively.

As mentioned in the proof of Theorem 1, this slightly different sampling model, while resulting in the same order-wise bounds, significantly simplifies the analysis due to the independence assumptions.

We will prove Theorem 3 under an additional random sign condition, that is, the signs of all non-zero entries of $\boldsymbol{S}$ are independent zero-mean random variables. Specifically, we will prove the following theorem.

Suppose that $\boldsymbol{X}$ obeys the incoherence condition with parameter $\mu_{1}$, and let $\lambda=\frac{1}{\sqrt{m\log\left(n_{1}n_{2}\right)}}$. Assume that $\tau\leq 0.2$ is some small positive constant, and that the signs of nonzero entries of $\boldsymbol{S}$ are independently generated with zero mean. If

then Robust-EMaC succeeds in recovering $\boldsymbol{X}$ with probability exceeding $1-(n_{1}n_{2})^{-2}$.

In fact, a simple derandomization argument introduced in [33, Section 2.2] immediately suggests that the performance of Robust-EMaC under the fixed-sign pattern is no worse than that under the random-sign pattern with sparsity parameter $2\tau$, i.e. the condition on the signs pattern of $\boldsymbol{S}$ is unnecessary and Theorem 3 follows after we establish Theorem 5. As a result, the section will focus on Theorem 3 with random sign patterns, which are much easier to analyze.

We adopt similar notations as in Section VI-A. That said, if we generate $\rho n_{1}n_{2}$ i.i.d. entry locations $\boldsymbol{a}_{i}$’s uniformly at random, and let the multi-sets $\Omega$ and $\Omega^{\text{clean}}$ contain respectively $\{\boldsymbol{a}_{i}|1\leq i\leq\rho n_{1}n_{2}\}$ and $\{\boldsymbol{a}_{i}|1\leq i\leq\rho(1-\tau)n_{1}n_{2}\}$), then

corresponding to sampling with replacement. Besides, $\mathcal{A}^{\prime}_{\Omega}$ (resp. $\mathcal{A}^{\prime}_{\Omega^{\text{clean}}}$) is defined similar to $\mathcal{A}_{\Omega}$ (resp. $\mathcal{A}_{\Omega^{\text{clean}}}$), but with the sum extending only over distinct samples.

We will establish that exact recovery can be guaranteed, if we can produce a valid dual certificate as follows.

for any matrix $\boldsymbol{M}$. If there exist a regularization parameter $\lambda$ $\left(0<\lambda<1\right)$ and a matrix $\boldsymbol{W}$ obeying

then Robust-EMaC is exact, i.e. the minimizer $\left(\hat{\boldsymbol{M}},\hat{\boldsymbol{S}}\right)$ satisfies $\hat{\boldsymbol{M}}=\boldsymbol{X}$.

We note that a reasonably tight bound on $\left\|\mathcal{P}_{T}\mathcal{A}\mathcal{P}_{T}-\frac{1}{\rho\left(1-\tau\right)}\mathcal{P}_{T}\mathcal{A}_{\Omega^{\text{clean}}}\mathcal{P}_{T}\right\|$ has been developed by Lemma 3. Specifically, there exists some constant $c_{1}>0$ such that if $\rho\left(1-\tau\right)n_{1}n_{2}>c_{1}\mu_{1}c_{\text{s}}r\log\left(n_{1}n_{2}\right)$, then one has

with probability exceeding $1-\left(n_{1}n_{2}\right)^{-4}$. Besides, Chernoff bound indicates that with probability exceeding $1-\left(n_{1}n_{2}\right)^{-3}$, none of the entries is sampled more than $10\log\left(n_{1}n_{2}\right)$ times. Equivalently,

Our objective in the remainder of this section is to produce a dual matrix $\boldsymbol{W}$ satisfying Condition (59).

VII-B Construction of Dual Certificate

Suppose that we generate $j_{0}$ independent random location multi-sets $\Omega_{j}^{\text{clean}}$, where $\Omega_{j}^{\text{clean}}$ contains $qn_{1}n_{2}$ i.i.d. samples uniformly at random. Here, we set $q:=\frac{\left(1-\tau\right)\rho}{j_{0}}$ and $\epsilon<\frac{1}{e}$. This way the distribution of the multi-set $\Omega$ is the same as $\Omega_{1}^{\text{clean}}\cup\Omega_{2}^{\text{clean}}\cup\cdots\cup\Omega_{j_{0}}^{\text{clean}}$.

We now propose constructing a dual certificate $\boldsymbol{W}$ as follows:

Take $\lambda=\frac{1}{\sqrt{m\log\left(n_{1}n_{2}\right)}}$. Note that the construction of $\boldsymbol{W}$ proceeds with a similar procedure as in Section VI-C, except that $\boldsymbol{F}_{0}$ and $\Omega_{i}$ are replaced by $\mathcal{P}_{T}\left(\boldsymbol{U}\boldsymbol{V}^{*}-\lambda\mbox{sgn}\left(\boldsymbol{S}_{\text{e}}\right)\right)$ and $\Omega_{i}^{\text{clean}}$, respectively.

We will justify that $\boldsymbol{W}$ is a valid dual certificate, by examining the conditions in (59) step by step.

with probability exceeding $1-(n_{1}n_{2})^{-3}$. Apply the same argument as for (40) to derive

(2) The second condition relies on an upper bound on $\left\|\mathcal{P}_{T^{\perp}}\left(\boldsymbol{W}+\lambda\mbox{sgn}\left(\boldsymbol{S}_{\text{e}}\right)\right)\right\|$. To this end, we proceed by controlling $\left\|\mathcal{P}_{T^{\perp}}\left(\boldsymbol{W}\right)\right\|$ and $\left\|\mathcal{P}_{T^{\perp}}\left(\lambda\mbox{sgn}\left(\boldsymbol{S}_{\text{e}}\right)\right)\right\|$ separately. Applying the same argument as for (52) suggests

where the second inequality follows since $\left\|\boldsymbol{M}\right\|_{\mathcal{A},2}\leq\sqrt{n_{1}n_{2}}\left\|\boldsymbol{M}\right\|_{\mathcal{A},\infty}$, and the last inequality arises from the fact that

Suppose that $s$ is a positive constant. then one has

for some constant $c_{9}>0$ with probability at least $1-(n_{1}n_{2})^{-4}$.

with probability exceeding $1-(n_{1}n_{2})^{-4}$.

It remains to control the term $\left\|\mathcal{P}_{T^{\perp}}\left(\lambda\mbox{sgn}\left(\boldsymbol{S}_{\text{e}}\right)\right)\right\|$, which is supplied in the following lemma.

Suppose that $\tau$ is a small positive constant, then one has

with probability at least $1-(n_{1}n_{2})^{-5}$.

(3) By construction, one has $\mathcal{A}^{\prime}_{\left(\Omega^{\text{clean}}\right)^{\perp}}\left(\boldsymbol{W}\right)=0$.

(4) The last step is to bound $\left\|\mathcal{A}^{\prime}_{\Omega^{\text{clean}}}\left(\boldsymbol{W}\right)\right\|_{\infty}$, which is apparently bounded above by $\left\|\mathcal{A}_{\Omega^{\text{clean}}}\left(\boldsymbol{W}\right)\right\|_{\infty}$. The construction procedure together with Lemma 6 allows us to bound

where the last inequality follows from (63). As a result, one can deduce

where the last inequality is obtained by setting $m>c_{12}\mu_{1}^{2}c_{\text{s}}^{2}r^{2}\log^{3}\left(n_{1}n_{2}\right)$ for some constant $c_{12}>0$.

To sum up, we have verified that $\boldsymbol{W}$ satisfies the four conditions required in (59), and is hence a valid dual certificate. This concludes the proof.

VIII Concluding Remarks

We present an efficient algorithm to estimate a spectrally sparse signal from its partial time-domain samples that does not require prior knowledge on the model order, which poses spectral compressed sensing as a low-rank Hankel structured matrix completion problem. Under mild incoherence conditions, our algorithm enables recovery of the multi-dimensional unknown frequencies with infinite precision, which remedies the basis mismatch issue that arises in conventional CS paradigms. We have shown both theoretically and numerically that our algorithm is stable against bounded noise and a constant proportion of arbitrary corruptions, and can be extended numerically to tasks such as super resolution. To the best of our knowledge, our result on Hankel matrix completion is also the first theoretical guarantee that is close to the information-theoretical limit (up to some logarithmic factor).

Our results are based on uniform random observation models. In particular, this paper considers directly taking a random subset of the time domain samples, it is also possible to take a random set of linear mixtures of the time domain samples, as in the renowned CS setting . This again can be translated into taking linear measurements of the low-rank $K$-fold Hankel matrix, given as $\boldsymbol{y}=\mathcal{B}(\boldsymbol{X}_{\text{e}})$. Unfortunately, due to the Hankel structures, it is not clear whether $\mathcal{B}$ exhibits approximate isometry property. Nonetheless, the technique developed in this paper can be extended without difficulty to analyze linear measurements, in a similar flavor of a golfing scheme developed for CS in .

It remains to be seen whether it is possible to obtain performance guarantees of the proposed EMaC algorithm similar to that in for super resolution. It is also of great interest to develop efficient numerical methods to solve the EMaC algorithm in order to accommodate large datasets.

IX Acknowledgement

This work was supported in part by the startup grant of The Ohio State University to Y. Chi. The authors thank Mr. Yuanxin Li for preparing Fig. 4 and Fig. 5.

Appendix A Bernstein Inequality

Our analysis relies heavily on the Bernstein inequality. To simplify presentation, we state below a user-friendly version of Bernstein inequality, which is an immediate consequence of [58, Theorem 1.6].

Then there exists a universal constant $c_{0}>0$ such that for any integer $a\geq 2$,

with probability at least $1-(d_{1}+d_{2})^{-a}$.

Appendix B Proof of Lemma 1

Consider any valid perturbation $\boldsymbol{H}$ obeying $\mathcal{P}_{\Omega}\left(\boldsymbol{X}+\boldsymbol{H}\right)=\mathcal{P}_{\Omega}\left(\boldsymbol{X}\right)$, and denote by $\boldsymbol{H}_{\text{e}}$ the enhanced form of $\boldsymbol{H}$. We note that the constraint requires $\mathcal{A}^{\prime}_{\Omega}\left(\boldsymbol{H}_{\text{e}}\right)=0$ (or $\mathcal{A}{}_{\Omega}\left(\boldsymbol{H}_{\text{e}}\right)=0$) and $\mathcal{A}^{\perp}\left(\boldsymbol{H}_{\text{e}}\right)=0$. In addition, set $\boldsymbol{Z}_{0}=\mathcal{P}_{T^{\perp}}\left(\boldsymbol{B}\right)$ for any $\boldsymbol{B}$ that satisfies $\left\langle\boldsymbol{B},\mathcal{P}_{T^{\perp}}\left(\boldsymbol{H}_{\text{e}}\right)\right\rangle=\left\|\mathcal{P}_{T^{\perp}}\left(\boldsymbol{H}_{\text{e}}\right)\right\|_{*}$ and $\left\|\boldsymbol{B}\right\|\leq 1$. Therefore, $\boldsymbol{Z}_{0}\in T^{\perp}$ and $\left\|\boldsymbol{Z}_{0}\right\|\leq 1$, and hence $\boldsymbol{U}\boldsymbol{V}^{*}+\boldsymbol{Z}_{0}$ is a sub-gradient of the nuclear norm at $\boldsymbol{X}_{\text{e}}$. We will establish this lemma by considering two scenarios separately.

(1) Consider first the case in which $\boldsymbol{H}_{\text{e}}$ satisfies

Since $\boldsymbol{U}\boldsymbol{V}^{*}+\boldsymbol{Z}_{0}$ is a sub-gradient of the nuclear norm at $\boldsymbol{X}_{\text{e}}$, it follows that

where (68) holds from (32), and (69) follows from the property of $\boldsymbol{Z}_{0}$ and the fact that $\left(\mathcal{A}^{\prime}_{\Omega}+\mathcal{A}^{\perp}\right)\left(\boldsymbol{H}_{\text{e}}\right)=0$. The last term of (69) can be bounded as

where the last inequality follows from the assumptions (33) and (34). Plugging this into (69) yields

where (71) follows from the inequality $\left\|\boldsymbol{M}\right\|_{*}\geq\left\|\boldsymbol{M}\right\|_{\text{F}}$ and (67). Therefore, $\boldsymbol{X}_{\text{e}}$ is the minimizer of EMaC.

We still need to prove the uniqueness of the minimizer. The inequality (71) implies that $\left\|\boldsymbol{X}_{\text{e}}+\boldsymbol{H}_{\text{e}}\right\|_{*}=\left\|\boldsymbol{X}_{\text{e}}\right\|_{*}$ holds only when $\left\|\mathcal{P}_{T^{\perp}}\left(\boldsymbol{H}_{\text{e}}\right)\right\|_{\text{F}}=0$. If $\left\|\mathcal{P}_{T^{\perp}}\left(\boldsymbol{H}_{\text{e}}\right)\right\|_{\text{F}}=0$, then $\left\|\mathcal{P}_{T}\left(\boldsymbol{H}_{\text{e}}\right)\right\|_{\text{F}}\leq\frac{n_{1}^{2}n_{2}^{2}}{2}\left\|\mathcal{P}_{T^{\perp}}\left(\boldsymbol{H}_{\text{e}}\right)\right\|_{\text{F}}=0$, and hence $\mathcal{P}_{T^{\perp}}\left(\boldsymbol{H}_{\text{e}}\right)=\mathcal{P}_{T}\left(\boldsymbol{H}_{\text{e}}\right)=0$, which only occurs when $\boldsymbol{H}_{\text{e}}=0$. Hence, $\boldsymbol{X}_{\text{e}}$ is the unique minimizer in this situation.

(2) On the other hand, consider the complement scenario where the following holds

We would first like to bound $\left\|\left(\frac{n_{1}n_{2}}{m}\mathcal{A}_{\Omega}+\mathcal{A}^{\perp}\right)\mathcal{P}_{T}\left(\boldsymbol{H}_{\text{e}}\right)\right\|_{\text{F}}$ and $\left\|\left(\frac{n_{1}n_{2}}{m}\mathcal{A}_{\Omega}+\mathcal{A}^{\perp}\right)\mathcal{P}_{T^{\perp}}\left(\boldsymbol{H}_{\text{e}}\right)\right\|_{\text{F}}$. The former term can be lower bounded by

On the other hand, since the operator norm of any projection operator is bounded above by $1$, one can verify that

where $a_{i}$ ($1\leq i\leq m$) are $m$ uniform random indices that form $\Omega$. This implies the following bound:

where the last inequality arises from our assumption. Combining this with the above two bounds yields

which immediately indicates $\mathcal{P}_{T^{\perp}}\left(\boldsymbol{H}_{\text{e}}\right)=0$ and $\mathcal{P}_{T}\left(\boldsymbol{H}_{\text{e}}\right)=0$. Hence, (72) can only hold when $\boldsymbol{H}_{\text{e}}=0$.

Appendix C Proof of Lemma 2

Since $\boldsymbol{U}$ (resp. $\boldsymbol{V}$) and $\boldsymbol{E}_{\text{L}}$ (resp. $\boldsymbol{E}_{\text{R}}$) determine the same column (resp. row) space, we can write

Note that $\sqrt{\omega_{k,l}}\boldsymbol{E}_{\text{L}}^{*}\boldsymbol{A}_{(k,l)}$ consists of $\omega_{k,l}$ columns of $\boldsymbol{E}_{\text{L}}^{*}$ (and hence it contains $r\omega_{k,l}$ nonzero entries in total). Owing to the fact that each entry of $\boldsymbol{E}_{\text{L}}^{*}$ has magnitude $\frac{1}{\sqrt{k_{1}k_{2}}}$, one can derive

A similar argument yields $\left\|\boldsymbol{A}_{(k,l)}\boldsymbol{E}_{\text{R}}^{*}\right\|_{\text{F}}^{2}\leq\frac{c_{\text{s}}r}{n_{1}n_{2}}$. Combining $\sigma_{\min}\left(\boldsymbol{E}_{\text{L}}^{*}\boldsymbol{E}_{\text{L}}\right)\geq\frac{1}{\mu_{1}}$ and $\sigma_{\min}\left(\boldsymbol{E}_{\text{R}}\boldsymbol{E}_{\text{R}}^{*}\right)\geq\frac{1}{\mu_{1}}$, (35) follows by plugging these facts into the above equations.

To show (36), since $\left|\left\langle\boldsymbol{A}_{\boldsymbol{b}},\mathcal{P}_{T}\left(\boldsymbol{A}_{\boldsymbol{a}}\right)\right\rangle\right|=\left|\left\langle\mathcal{P}_{T}\left(\boldsymbol{A}_{\boldsymbol{b}}\right),\boldsymbol{A}_{\boldsymbol{a}}\right\rangle\right|$, we only need to examine the situation where $\omega_{\boldsymbol{b}}<\omega_{\boldsymbol{a}}$. Observe that

Owing to the multi-fold Hankel structure of $\boldsymbol{A}_{\boldsymbol{a}}$, the matrix $\boldsymbol{U}\boldsymbol{U}^{*}\sqrt{\omega_{\boldsymbol{a}}}\boldsymbol{A}_{\boldsymbol{a}}$ consists of $\omega_{\boldsymbol{a}}$ columns of $\boldsymbol{U}\boldsymbol{U}^{*}$. Since there are only $\omega_{\boldsymbol{b}}$ nonzero entries in $\boldsymbol{A}_{\boldsymbol{b}}$ each of magnitude $\frac{1}{\sqrt{\omega_{\boldsymbol{b}}}}$, we can derive

Each entry of $\boldsymbol{U}\boldsymbol{U}^{*}$ is bounded in magnitude by

We still need to bound the magnitude of $\left\langle\boldsymbol{U}\boldsymbol{U}^{*}\boldsymbol{A}_{\boldsymbol{a}}\boldsymbol{V}\boldsymbol{V}^{*},\boldsymbol{A}_{\boldsymbol{b}}\right\rangle$. One can observe that for the $k$th row of $\boldsymbol{U}\boldsymbol{U}^{*}$:

Similarly, for the $l$th column of $\boldsymbol{V}\boldsymbol{V}^{*}$, one has $\left\|\boldsymbol{V}\boldsymbol{V}^{*}\boldsymbol{e}_{l}\right\|_{\text{F}}\leq\sqrt{\frac{\mu_{1}c_{\text{s}}r}{n_{1}n_{2}}}$. The magnitude of the entries of $\boldsymbol{U}\boldsymbol{U}^{*}\boldsymbol{A}_{\boldsymbol{a}}\boldsymbol{V}\boldsymbol{V}^{*}$ can now be bounded by

where we used $\left\|\boldsymbol{A}_{\boldsymbol{a}}\right\|=1/\sqrt{\omega_{\boldsymbol{a}}}$. Since $\boldsymbol{A}_{\boldsymbol{b}}$ has only $\omega_{\boldsymbol{b}}$ nonzero entries each has magnitude $\frac{1}{\sqrt{\omega_{\boldsymbol{b}}}}$, one can verify that

The above bounds (75), (76) and (77) taken together lead to (36).

Appendix D Proof of Lemma 3

for any $(k,l)\in[n_{1}]\times[n_{2}]$. For any matrix $\boldsymbol{M}$, we can compute

where the last inequality follows from (37). This further gives

where (80) uses (79). Applying Lemma 11 yields that there exists some constant $0<\epsilon\leq\frac{1}{2}$ such that

with probability exceeding $1-\left(n_{1}n_{2}\right)^{-4}$, provided that $m>c_{1}\mu_{1}c_{\text{s}}r\log\left(n_{1}n_{2}\right)$ for some universal constant $c_{1}>0$.

Appendix E Proof of Lemma 4

Suppose that $\mathcal{A}_{\Omega}=\sum_{i=1}^{m}\mathcal{A}_{\boldsymbol{a}_{i}}$, where $\boldsymbol{a}_{i}$, $1\leq i\leq m$, are $m$ independent indices drawn uniformly at random from $[n_{1}]\times[n_{2}]$. Define

where the first inequality follows since $\frac{1}{m}\sum_{k,l}\mathcal{A}_{(k,l)}\left(\boldsymbol{M}\right)=\frac{1}{m}\mathcal{A}\left(\boldsymbol{M}\right)$, and the last inequality arises from the fact that all non-zero entries of $\boldsymbol{A}_{(k,l)}\cdot\boldsymbol{A}_{(k,l)}^{\top}$ lie on its diagonal and are bounded in magnitude by $\frac{1}{\omega_{k,l}}$. This immediately suggests

On the other hand, the operator norm of each $\boldsymbol{S}_{(k,l)}$ can be bounded as follows

where (83) holds since $\left\|\boldsymbol{A}_{(k,l)}\right\|=\frac{1}{\sqrt{\omega_{k,l}}}$ and the last equality follows by applying the definition of $\left\|\cdot\right\|_{\mathcal{A},\infty}$.

Finally, we combine the above two bounds together with Bernstein inequality (Lemma 11) to obtain

with high probability, where $c_{2}>0$ is some absolute constant.

Appendix F Proof of Lemma 5

Write $\mathcal{A}_{\Omega}=\sum_{i=1}^{m}\mathcal{A}_{\boldsymbol{a}_{i}}$, where $\boldsymbol{a}_{i}$ ($1\leq i\leq m$) are $m$ independent indices uniformly drawn from $[n_{1}]\times[n_{2}]$. By the definition of $\left\|\boldsymbol{M}\right\|_{\mathcal{A},2}$, we need to examine the components

Define a set of variables $z_{(\alpha,\beta)}$’s to be

The definition of $\left\|\boldsymbol{M}\right\|_{\mathcal{A},2}$ allows us to express

where $\boldsymbol{z}_{(\alpha,\beta)}$’s are defined to be $n_{1}n_{2}$-dimensional vectors

where (86) follows from the definition of $\mu_{5}$ in (45). Now it follows that

where (87) follows from (42). On the other hand,

with high probability for some numerical constant $c_{3}>0$, which completes the proof.

Appendix G Proof of Lemma 6

From Appendix F, it is straightforward that

where $z_{\boldsymbol{a}_{i}}^{\left(k,l\right)}$’s are defined as (84). Using similar techniques as (86), we can obtain

where we have made use of the fact (36). As a result, one has

The Bernstein inequality in Lemma 11 taken collectively with the union bound yields that

with high probability for some constant $c_{4}>0$, completing the proof.

Appendix H Proof of Lemma 7

To bound $\left\|\boldsymbol{U}\boldsymbol{V}^{*}\right\|_{\mathcal{A},\infty}$, observe that there exists a unitary matrix $\boldsymbol{B}$ such that

For any $\left(k,l\right)\in[n_{1}]\times[n_{2}]$, we can then bound

Since $\boldsymbol{A}_{(k,l)}$ has only $\omega_{k,l}$ nonzero entries each of magnitude $\frac{1}{\sqrt{\omega_{k,l}}}$, this leads to

The rest is to bound $\left\|\boldsymbol{U}\boldsymbol{V}^{*}\right\|_{\mathcal{A},2}$ and $\left\|\mathcal{P}_{T}\left(\sqrt{\omega_{k,l}}\boldsymbol{A}_{\left(k,l\right)}\right)\right\|_{\mathcal{A},2}$. Observe that the $i$th row of $\boldsymbol{U}\boldsymbol{V}^{*}$ obeys

Moreover, the matrix $\mathcal{P}_{T}\left(\sqrt{\omega_{\alpha,\beta}}\boldsymbol{A}_{(\alpha,\beta)}\right)$ enjoys similar properties as well, which we briefly reason as follows. First, the matrix $\boldsymbol{U}\boldsymbol{U}^{*}\left(\sqrt{\omega_{\alpha,\beta}}\boldsymbol{A}_{(\alpha,\beta)}\right)$ obeys

since the operator norm of $\boldsymbol{U}$ and $\sqrt{\omega_{\alpha,\beta}}\boldsymbol{A}_{(\alpha,\beta)}$ are both bounded by 1. The same bound for $\sqrt{\omega_{\alpha,\beta}}\boldsymbol{A}_{(\alpha,\beta)}\boldsymbol{V}\boldsymbol{V}^{*}$ can be demonstrated via the same argument as for $\boldsymbol{U}\boldsymbol{U}^{*}\left(\sqrt{\omega_{\alpha,\beta}}\boldsymbol{A}_{(\alpha,\beta)}\right)$. Additionally, for $\boldsymbol{U}\boldsymbol{U}^{*}\left(\sqrt{\omega_{\alpha,\beta}}\boldsymbol{A}_{(\alpha,\beta)}\right)\boldsymbol{V}\boldsymbol{V}^{*}$ one has

Now our task boils down to bounding $\left\|\boldsymbol{M}\right\|_{\mathcal{A},2}$ for some matrix $\boldsymbol{M}$ satisfying some energy constraints per row, which subsumes $\left\|\boldsymbol{U}\boldsymbol{V}^{*}\right\|_{\mathcal{A},2}$ and $\left\|\mathcal{P}_{T}\left(\sqrt{\omega_{k,l}}\boldsymbol{A}_{\left(k,l\right)}\right)\right\|_{\mathcal{A},2}$ as special cases. We can then conclude the proof by applying the following lemma.

Denote by the set $\mathcal{M}$ of feasible matrices satisfying

Then there exists some universal constant $c_{3}>0$ such that

For ease of presentation, we split any matrix $\boldsymbol{M}$ into 4 parts, which are defined as follows

$\boldsymbol{M}^{\left(1\right)}$: the matrix containing all upper triangular components of all upper triangular blocks of $\boldsymbol{M}$;

$\boldsymbol{M}^{\left(2\right)}$: the matrix containing all lower triangular components of all upper triangular blocks of $\boldsymbol{M}$;

$\boldsymbol{M}^{\left(3\right)}$: the matrix containing all upper triangular components of all lower triangular blocks of $\boldsymbol{M}$;

$\boldsymbol{M}^{\left(4\right)}$: the matrix containing all lower triangular components of all lower triangular blocks of $\boldsymbol{M}$.

Here, we use the term “upper triangular” and “lower triangular” in short for “left upper triangular” and “right lower triangular”, which are more natural for Hankel matrices. Instead of maximizing $\left\|\boldsymbol{M}\right\|_{\mathcal{A},2}$ directly, we will handle $\max_{\boldsymbol{M}\in\mathcal{M}}\|\boldsymbol{M}^{(l)}\|_{\mathcal{A},2}^{2}$ for each $1\leq l\leq 4$ separately, owing to the fact that

In the sequel, we only demonstrate how to control $\|\boldsymbol{M}^{(1)}\|_{\mathcal{A},2}$. Similar bounds can be derived for $\|\boldsymbol{M}^{(l)}\|_{\mathcal{A},2}$ ($2\leq l\leq 4$) via very similar argument.