Spectral Compressed Sensing via Projected Gradient Descent

Introduction

In this paper, we are interested in the problem of reconstructing a spectrally sparse signal with or without damping from its nonuniform time-domain samples. Let $x(t)$ be a one-dimensional signal. We say that $x(t)$ is spectrally sparse if it is superposition of a few complex sinusoids, namely

where $\imath=\sqrt{-1}$, $r$ is the model order, $f_{k}$ is the frequency of each sinusoid, $d_{k}$ is the weight of each sinusoid, and $\tau_{k}\geq 0$ is a damping factor. Let $n>0$ be a natural number. Without loss of generality, we assume $f_{k}\in[0,1)$ and consider the samples of $x(t)$ at all the integer values from to $n-1$, denoted $\bm{x}$. That is,

Spectrally sparse signals of the form (1) and the corresponding sampling model in (2) arise in many areas of science and engineering including magnetic resonance imaging , fluorescence microscopy , radar imaging , nuclear magnetic resonance spectroscopy , and analog-to-digital conversion . However, in those real-world applications, full sampling at all the points on a uniform grid is either time-consuming or technically prohibited. In addition, the signal may become too weak to be detected after a certain period of time when $\tau_{k}>0$. Therefore, for the purpose of more efficient data acquisition, nonuniform sampling is typically used in practice. When restricted to the sampling model in (2), this means that only partial entries of $\bm{x}$ are known and we need to estimate the missing ones. Let $\Omega$ be subset of $\{0,\cdots,n-1\}$ corresponding to the observed entries, and let $\mathcal{P}_{\Omega}$ be the associated sampling operator which acquires only the entries indexed by $\Omega$. Then the task can be formally expressed as:

2 Prior Art and Main Contributions

It is clear that (3) is a task that cannot be achieved if $\bm{x}$ does not have any intrinsic simple structures. Fortunately, the signal of interest in this paper is spectrally sparse. Moreover, the number of degrees of freedom in $\bm{x}$ is completely determined by the number of Fourier modes in $x(t)$, which is proportional to $r$ and independent of $n$. This key observation suggests the possibility of reconstructing $\bm{x}$ from its partial revealed entries, which can be further achieved by exploiting the simplicity of $\bm{x}$ in different ways.

Note that we are mainly interested in the scenario where $\bm{x}$ only has a few Fourier components (i.e., $r$ is small). Thus, one can utilize the sparsity of $\bm{x}$ in the frequency domain to design reconstruction algorithms. In particular, if there is no damping in $\bm{x}$, spectral compressed sensing can be recast as a conventional compressed sensing problem after discretization of the Fourier domain; so many existing algorithms for compressed sensing are available, such as Basis Pursuit , IHT , CoSaMP and SP . However, the performance of the compressed sensing approach for spectrally sparse signal recovery suffers from the mismatch error between the true frequencies and the discrete frequencies . A grid-free approach was developed in which exploited the frequency sparsity of $\bm{x}$ in a continuous manner via the atomic norm minimization (ANM). It was shown in that ANM could achieve exact recovery from $O(r\log(r)\log(n))$ random time-domain samples under some mild conditions.

By the Vandermonde decomposition, one may easily see that the Hankel matrix computed from a spectrally sparse signal is low rank when $r$ is small relative to $n$. Consequently, spectral compressed sensing can be reformulated as a low rank Hankel matrix completion problemSee Section 2.1 for details.. Inspired by low rank matrix completion , another grid-fee method known as enhanced matrix completion (EMaC) was developed in by reformulating the non-convex low rank Hankel matrix completion problem into a convex Hankel matrix nuclear norm minimization problem. EMaC was shown to be able to reconstruct a spectrally sparse signal with high probability provided the number of observed entries is $O(r\log^{4}(n))$. The same approach was studied in under the Gaussian random sampling model, and various first-order methods were discussed in for the regularized Hankel matrix nuclear norm minimization problem. Alternative to EMaC, there have been several non-convex algorithms which were designed to solve the low rank Hankel matrix completion directly. Examples include PWGD , IHT and FIHT . Compared to the convex approaches such as ANM and EMaC, those non-convex algorithms are typically much more efficient, especially for higher dimensional problems. Moreover, inspired by the guarantee analysis of Riemannian optimization for low rank matrix reconstruction , it was shown in that FIHT with a proper initial guess was able to reconstruct a spectrally sparse signal with high probability from $O(r^{2}\log^{2}(n))$ random observations. For multi-dimensional spectrally sparse signal recovery problems, we can also exploit the low rank tensor structure of the signal when developing recovery algorithms, see for example and references therein.

The main contributions of this work are two-fold. Firstly, we present a new non-convex algorithm for spectral compressed sensing via low rank Hankel matrix completion, which we refer to as Projected Gradient Descent (PGD). Extensive empirical performance comparisons show that PGD is competitive with other state-of-the-art spectral compressed sensing algorithms both in terms of the problem size that can be solved and in terms of overall computation time. Secondly, exact recovery guarantee has been established for PGD, showing that PGD can successfully recover a spectrally sparse signal from $O(r^{2}\log(n))$ random observed entries.

Although we focus on spectrally sparse signal recovery in this paper, the proposed PGD algorithm can be easily extended to the general low rank Hankel matrix completion problem. Moreover, the recovery guarantee analysis equally applies provided the underlying target matrix is incoherentSee Definition 2.1.. Low-rank Toeplitz matrices can also be provably recovered from partial revealed entries by a slightly modified version of PGD.

3 Outline and Notation

The remainder of this paper is organized as follows. We present the details of PGD, along with its recovery guarantee in Section 2. In Section 3 we evaluate the empirical performance of PGD with a set of numerical experiments. The proof of the exact recovery guarantee is presented in Section 4. We conclude the paper with some potential future directions in Section 5.

Algorithm and Main Result

Thus, one has $[\mathcal{H}\bm{z}]^{(i,j)}=z_{i+j}$ for $i\in[n_{1}]\mbox{ and }j\in[n_{2}].$ In particular, the $(i,j)$-th entry of the Hankel matrix formed from the spectrally sparse signal $\bm{x}$ is given by

If we let $w_{k}=e^{(2\pi\imath f_{k}-\tau_{k})}$ for $k=1,\cdots,r$, it follows immediately that $\mathcal{H}\bm{x}$ admits the following Vandermonde decomposition:

and $\bm{D}=\operatorname{diag}(d_{1},\cdots,d_{r})$. Moreover, one has $\operatorname{rank}(\mathcal{H}\bm{x})=r$ provided the frequencies $\{f_{k}\}_{k=1}^{r}$ are different with each other and the diagonal entries of $\bm{D}$ are all nonzeros.

where $w_{a}$ in the second line is the number of entries in the $a$-th skew-diagonal of an $n_{1}\times n_{2}$ matrix, and $\mathcal{D}$ in the last line is a linear map which scales the $a$-th entry of a vector by a factor of $\sqrt{w_{a}}$ for all $a=0,\cdots,n-1$. We have seen that $\mathcal{H}\bm{x}$ is a rank $r$ matrix. Thus, to reconstruct $\bm{x}$, we may seek a signal $\bm{z}$ such that $\operatorname{rank}(\mathcal{H}\bm{z})=r$ and $\mathcal{H}\bm{z}$ fits the revealed skew-diagonals of $\mathcal{H}\bm{x}$ as well as possible by solving a rank constraint weighted least square problem:

which will be our primary focus in this paper. A more direct interpretation of (5) is as follows. Since $\bm{y}=\mathcal{D}\bm{x}$, $\mathcal{P}_{\Omega}(\bm{y})=\mathcal{P}_{\Omega}(\mathcal{D}\bm{x})=\mathcal{D}\mathcal{P}_{\Omega}(\bm{x})$, $\operatorname{rank}(\mathcal{G}\bm{y})=\operatorname{rank}(\mathcal{H}\bm{x})=r$, and $\mathcal{D}$ is invertible, one can instead attempt to reconstruct $\bm{y}$ from $\mathcal{P}_{\Omega}(\bm{y})$ by seeking a signal that corresponds to a low rank Hankel matrix and fits the observations as well as possible.

2 Algorithm: Projected Gradient Descent

Thus, by further noting that $\bm{z}=\mathcal{G}^{*}(\mathcal{G}\bm{z})=\mathcal{G}^{*}(\bm{Z}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}\bm{Z}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}^{*})$, we can rewrite (5) using $\bm{Z}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}$ and $\bm{Z}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}$ as

which is an equality constraint minimization problem. Alternatively, (6) can be interpreted as follows: we estimate the rank $r$ matrix $\mathcal{G}\bm{y}$ by a Hankel matrix of the form $\bm{Z}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}\bm{Z}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}^{*}$ that minimizes the mismatch in the measurement domain. Once $\mathcal{G}\bm{y}$ is reconstructed, one can recover $\bm{y}$ via $\bm{y}=\mathcal{G}^{*}(\mathcal{G}\bm{y})$.

Putting the constraint and the objective function in (6) together allows us to consider an optimization problem without the equality constraint by minimizing

denotes the concatenation of $\bm{Z}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}$ and $\bm{Z}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}$, and the weight $p=m/n$ is the sampling ratio. Let $\mathcal{G}\bm{y}=\bm{U}\bm{\Sigma}\bm{V}^{*}$ be the reduced singular value decomposition (SVD) of $\mathcal{G}\bm{y}$. Define

where $\bm{M}_{\scaleto{\bm{U}\mathstrut}{6pt}}=\bm{U}\bm{\Sigma}^{1/2}$ and $\bm{M}_{\scaleto{\bm{V}\mathstrut}{6pt}}=\bm{V}\bm{\Sigma}^{1/2}$. It is easily shown that $f(\bm{Z})=0$ and thus achieves its minimum for the set of matrices

Note that (8) is also a set of solutions for the equality constrained problem (6). Among this set of solutions, there are ones which are highly unbalanced, i.e., these having $\left\|\bm{Z}_{\scaleto{\bm{U}\mathstrut}{6pt}}\right\|_{F}\rightarrow 0$ and $\left\|\bm{Z}_{\scaleto{\bm{V}\mathstrut}{6pt}}\right\|_{F}\rightarrow\infty$, or vice versa. For example, let $\bm{Z}_{\scaleto{\bm{U}\mathstrut}{6pt}}=\alpha\bm{M}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}$ and $\bm{Z}_{\scaleto{\bm{V}\mathstrut}{6pt}}=\alpha^{-1}\bm{M}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}$ for $\alpha$ being a real number that approaches either zero or infinity. Those solutions are unfavorable for the purpose of both computation and analysis. In order to reduce the solution space and avoid the occurrence of the pathological solutions, we add the regularizer function

to $f(\bm{Z})$ and instead consider the minimization problem with respect to

where $\lambda>0$ is to be determined. Here, $g(\bm{Z})$ in some sense penalizes the mismatch between the sizes of $\bm{Z}_{\scaleto{\bm{U}\mathstrut}{6pt}}$ and $\bm{Z}_{\scaleto{\bm{V}\mathstrut}{6pt}}$, and it was also used in rectangular low rank matrix recovery, see .

Now, the set of solutions that minimizes $F(\bm{Z})$ or at which $F(\bm{Z})=0$ is given by

Let $\bm{M}^{*}\bm{Z}=\bm{Q}_{1}\bm{\Lambda}\bm{Q}_{2}^{*}$ be the SVD of $\bm{M}^{*}\bm{Z}$. By the Von Neumann’s trace inequality , the above minimum is achieved at the unitary matrix $\bm{Q}_{\scaleto{\bm{Z}\mathstrut}{6pt}}$ given by

2.2 Which Feasible Set?

As we have already seen, the goal in spectrally sparse signal recovery is in fact to reconstruct a low rank Hankel matrix matrix $\mathcal{G}\bm{y}$ from its partial revealed skew-diagonals. In general, it is impossible to reconstruct a low rank matrix from entry-wise sampling unless its singular vectors are weakly correlated with the sampling basis. Here, we are interested in $\mu_{0}$-incoherent matrix which was first introduced in for low rank matrix completion.

With $\mathcal{G}\bm{y}=\bm{U}\bm{\Sigma}\bm{V}^{*}$ being the SVD of $\mathcal{G}\bm{y}$, we say $\mathcal{G}\bm{y}$ is $\mu_{0}$-incoherent if there exists an absolute numerical constant $\mu_{0}>0$ such that

A sufficient condition for $\mathcal{G}\bm{y}$ to be $\mu_{0}$-incoherent can be derived based on the Vandermonde decomposition of $\mathcal{G}\bm{y}$. Assume that

which implies $\mathcal{G}\bm{y}$ is $\mu_{0}$-incoherent. Moreover, [31, Thm. 2] says that (12) holds for undamping signals provided the minimum wrap-around distance between each pair of the frequencies of the spectrally sparse signal is greater than about $2/n$.

Let $\mu$ and $\sigma$ be two numerical constants such that $\mu\geq\mu_{0}$ and $\sigma\geq\sigma_{1}(\mathcal{G}\bm{y})$. When $\mathcal{G}\bm{y}$ is $\mu_{0}$-incoherent, the matrix $\bm{M}$ constructed in (7) satisfies $\left\|\bm{M}\right\|_{2,\infty}\leq\sqrt{\mu c_{s}r\sigma/n}$. Moreover, letting $\mathcal{C}$ be a convex set defined as

it is evident that $\mathcal{S}\subset\mathcal{C}$. Therefore, we can restrict our search on the feasible set $\mathcal{C}$ when computing the minimum or zero value of $F(\bm{Z})$.

2.3 Algorithm

The discussion above tells us that we can reconstruct the low rank factors $\bm{M}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}$ and $\bm{M}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}$ of the ground truth matrix $\mathcal{G}\bm{y}$ by minimizing the function $F(\bm{Z})$ on the feasible set $\mathcal{C}$, namely

where $F(\bm{Z})$ is defined in (9) and $\mathcal{C}$ is defined in (13). We present a simple projected gradient descent algorithm for this problem, see Algorithm 1.

In each iteration of the algorithm, the current estimate $\bm{Z}^{k}$ is updated along the negative gradient descent direction $-\nabla F(\bm{Z}^{k})$, using a stepsize $\eta$, followed by projection onto the convex set $\mathcal{C}$. Since we are working with complex matrices, the gradient $F(\bm{Z})$ of a matrix $\bm{Z}$ is calculated under the Wirtinger calculus, given by

PGD can be implemented very efficiently and the main computational cost per iteration is $O(r^{2}n+rn\log(n))$ flops, which lies in the computation of $\nabla F(\bm{Z})$ in each iteration. Taking the computation of $\nabla F_{{\scaleto{\bm{U}\mathstrut}{6pt}}}(\bm{Z})$ as an example, we note that

Clearly, the second term can be computed using $O(r^{2}n)$ flops. Let $\bm{w}=p^{-1}\mathcal{P}_{\Omega}(\mathcal{G}^{*}(\bm{Z}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}\bm{Z}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}^{*})-\bm{y})-\mathcal{G}^{*}(\bm{Z}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}\bm{Z}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}^{*})$. Since we can compute $\mathcal{G}^{*}(\bm{Z}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}\bm{Z}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}^{*})$ by $r$ fast convolutions, $\bm{w}$ can be obtained using $O(rn\log(n))$ flops. Moreover, $(\mathcal{G}\bm{w})\bm{Z}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}$ can be computed via $r$ fast Hankel matrix-vector multiplications that also cost $O(rn\log(n))$ flops.

Before proceeding, it is worth noting that non-convex (projected) gradient decent methods have received intensive investigations for other low rank matrix recovery problems, such as unstructured low rank matrix recovery and matrix completion , phase retrieval , robust principle component analysis , and blind deconvolution . In those papers, lower bounds on the sampling complexity have been established under different random measurement models, showing that the number of measurements needed for the successful recovery of the target matrices is essentially determined by the number of degrees of freedom in the matrices. In particular, a projected gradient descent algorithm was studied in for unstructured rectangular low rank matrix completion. The convergence analysis of PGD in this paper is directly inspired by , though the technical details are substantially different.

3 Main Result

In the guarantee analysis of PGD, we assume $\mu$ and $\sigma$ in (13) are two tuning parameters obeying $\mu\geq\mu_{0}$ and $\sigma\geq\sigma_{1}(\mathcal{G}\bm{y})$ so that $\bm{M}\in\mathcal{C}$. For conciseness, we take $\sigma=\sigma_{1}(\bm{L}_{0})/(1-\varepsilon_{0})$ for some $0<\varepsilon_{0}<1$ and will later show that $\sigma\geq\sigma_{1}(\mathcal{G}\bm{y})$ with high probability.

Assume $\mathcal{G}\bm{y}$ is $\mu_{0}$-incoherent. Let $\varepsilon_{0}$ be a absolute constant obeying $0<\varepsilon_{0}\leq 1/11$. Let $\mu\geq\mu_{0}$ and $\sigma=\sigma_{1}(\bm{L}_{0})/(1-\varepsilon_{0})$. If we take $\lambda=1/4$ in (9), then with probability at least $1-c_{1}\cdot n^{-2}$, the sequence $\left\{\bm{Z}^{k}\right\}_{k\geq 1}$ returned by Algorithm 1 obeys

provided $m\geq c_{2}\hskip 1.42271pt\varepsilon_{0}^{-2}\mu^{2}c_{s}^{2}\kappa^{2}r^{2}\log(n)$, where $\kappa=\sigma_{1}(\mathcal{G}\bm{y})/\sigma_{r}(\mathcal{G}\bm{y})$.

1). After an approximation of $\mathcal{G}\bm{y}$, given by $\bm{Z}^{k}_{\scaleto{\bm{U}\mathstrut}{6pt}}(\bm{Z}^{k}_{\scaleto{\bm{V}\mathstrut}{6pt}})^{*}$, is obtained from PGD, we can estimate $\bm{y}$ by $\bm{y}^{k}=\mathcal{G}^{*}(\bm{Z}^{k}_{\scaleto{\bm{U}\mathstrut}{6pt}}(\bm{Z}^{k}_{\scaleto{\bm{V}\mathstrut}{6pt}})^{*})$, and in turn estimate $\bm{x}$ by $\mathcal{D}^{-1}\bm{y}^{k}$. Recall from (11) that $\bm{Q}_{{\scaleto{\bm{Z}\mathstrut}{6pt}}^{k}}$ is a unitary matrix which obeys ${\mbox{dist}}(\bm{Z}^{k},\bm{M})=\left\|\bm{Z}^{k}-\bm{M}\bm{Q}_{{\scaleto{\bm{Z}\mathstrut}{6pt}}^{k}}\right\|_{F}$. A simple calculation yields

2). After each iteration, Theorem 2.1 implies that the distance between the estimate given by PGD and $\bm{M}$ is reduced by at least of a factor of $1-O(1/(\mu c_{s}r\kappa)^{2})$. Thus, after $k\approx O((\mu c_{s}r\kappa)^{2}\log(1/\epsilon))$ iterations, one has ${\mbox{dist}^{2}}(\bm{Z}^{k},\bm{M})\leq\epsilon\cdot{\mbox{dist}^{2}}(\bm{Z}^{0},\bm{M})$.

3). It was shown in that FIHT can achieve exact recovery when the number of revealed entries is of order $O(\kappa^{6}r^{2}\log^{2}(n))$. In contrast, the sampling complexity of PGD is only a quadratic function of $\kappa$ and a linear function of $\log(n)$. Moreover, the exact recovery guarantee of FIHT relies on a more complicated initialization scheme which requires a partition of the observed entries into $O(\log(n))$ groups, while the initial guess constructed for the exact recovery guarantee of PGD can be computed much more easily.

4 Extension to Higher Dimension

So far we have restricted our attention to one-dimensional spectrally sparse signal reconstruction problem. Our algorithm and results can be extended to higher dimensions based on the Hankel structures of multi-dimensional spectrally sparse signals. Without loss of generality, we discuss the two-dimensional setting but emphasize that the situation in general $d$-dimensions is similar.

The two-fold Hankel matrix of $\bm{X}$ is given by

where each block is an $n_{1}\times(N_{1}-n_{1}+1)$ Hankel matrix corresponding to a column of $\bm{X}$,

Clearly, $\mathcal{H}\bm{X}$ is an $(n_{1}n_{2})\times(N_{1}-n_{1}+1)(N_{2}-n_{2}+1)$ matrix. Letting $i=i_{1}+i_{2}\cdot n_{1}$ and $j=j_{1}+j_{2}\cdot(N_{1}-n_{1}+1)$, the $(i,j)$-th entry of $\mathcal{H}\bm{X}$ is given by

For $k=1,\cdots,r$, we define the four vectors $\bm{w}_{k}^{[n_{1}]}$, $\bm{w}_{k}^{[N_{1}-n_{1}+1]}$, $\bm{z}_{k}^{[n_{2}]}$, and $\bm{z}_{k}^{[N_{2}-n_{2}+1]}$ as

Let $\bm{E}_{L}$ be an $(n_{1}n_{2})\times r$ matrix with the $k$-th column being given by $\bm{z}_{k}^{[n_{2}]}\otimes\bm{w}_{k}^{[n_{1}]}$, and let $\bm{E}_{R}$ be an $(N_{1}-n_{1}+1)(N_{2}-n_{2}+1)\times r$ matrix with the $k$-th column being given by $\bm{z}_{k}^{[N_{2}-n_{2}+1]}\otimes\bm{w}_{k}^{[N_{1}-n_{1}+1]}$. Then it follows from (17) that $\mathcal{H}\bm{X}$ admits the Vandermonde decomposition

where $\bm{D}=\operatorname{diag}(d_{1},\cdots,d_{r})$. Thus, it is self-evident that $\mathcal{H}\bm{X}$ is a rank $r$ matrix.

subject to a feasible set $\mathcal{C}$, where $\mathcal{G}^{*}$ is the adjoint of $\mathcal{G}$ which obeys $\mathcal{G}^{*}\mathcal{G}=\mathcal{I}$,

is an $(n_{1}n_{2}+(N_{1}-n_{1}+1)(N_{2}-n_{2}+1))\times r$ matrix, and $\mathcal{C}$ is a convex set similar to the one defined in (13) but the size of $\bm{Z}$ is different.

Therefore, a projected gradient descent algorithm can also be developed for the two-dimensional spectrally sparse signal reconstruction problem. Let $\mathcal{G}\bm{Y}=\bm{U}\bm{\Sigma}\bm{V}^{T}$ be the SVD of $\mathcal{G}\bm{Y}$. We say $\mathcal{G}\bm{Y}$ is $\mu_{0}$-incoherent if there exists a numerical constant $\mu_{0}>0$ such that

where $c_{s}=\max\{N_{1}N_{2}/(n_{1}n_{2}),N_{1}N_{2}/((N_{1}-n_{1}+1)(N_{2}-n_{2}+1))\}$. Based on [30, Theorem 1], one can show that $\mathcal{G}\bm{Y}$ ($=\mathcal{H}\bm{X}$) is $\mu_{0}$-incoherent if there is no damping in $\bm{X}$ and the minimum wrap-around distance between the underlying frequencies $\{f_{ik}\}_{k=1}^{r}$ is greater than about ${2}/{N_{i}}$ for $i=1,2$. Let

where $\bm{M}_{\scaleto{\bm{U}\mathstrut}{6pt}}=\bm{U}\bm{\Sigma}^{1/2}$ and $\bm{M}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}=\bm{V}\bm{\Sigma}^{1/2}$. If we assume $\mathcal{G}\bm{Y}$ is $\mu_{0}$-incoherent and $\mu$ and $\sigma$ in $\mathcal{C}$ are properly tuned such that $\bm{M}\in\mathcal{C}$, then the exact guarantee analysis of PGD for the one-dimensional case can be extended immediately to the two-dimensional case. It can be established that $O(\mu^{2}c_{s}^{2}\kappa^{2}r^{2}\log(N_{1}N_{2}))$ number of measurements are sufficient for PGD to achieve the successful recovery of a two-dimensional spectrally sparse signal.

Numerical Experiments

In this section, we conduct numerical experiments to evaluate the performance of PGDIn our random simulations, we didn’t find much difference between the performance of PGD and the performance of the gradient descent algorithm applied to $f(\bm{Z})$ directly. However, since the extra cost incurred by computing the gradient of $g(\bm{Z})$ and the projection $\mathcal{P}_{\mathcal{C}}(\bm{Z})$ is marginal, it is appealing to run PGD for its recovery guarantee.. The experiments are executed from MATLAB R2017a on a 64-bit Linux machine with multi-core Intel Xeon CPU E5-2667 v3 at 3.20GHz and 64GB of RAM. In Section 3.1, we investigate the largest number of Fourier components that can be successfully recovered by PGD. The tests are conducted on one-dimensional signals in large part due to the high computational cost of this type of simulations. Then we evaluate PGD against computational efficiency, robustness to additive noise, and sensitivity to mis-specification of model order on three-dimensional signals in Sections 3.2, 3.3, and 3.4, respectively. The initial guess of PGD is computed using the PROPACK package , and the parameters $\mu$ and $\sigma$ used in the projection are estimated from the initialization. Instead of using the constant stepsize suggested in the main result which appears to be conservative, we choose the stepsize via a backtracking line search in the implementation.

We evaluate the recovery ability of PGD in the framework of phase transition and compare it with ANM , EMaC and FIHT . ANM and EMaC are implemented using CVX with default parameters. The test spectrally sparse signals of length $n$ with $r$ frequency components are formed in the following way: each frequency $f_{k}$ is randomly generated from $[0,1)$, and the argument of each complex coefficient $d_{k}$ is uniformly sampled from $[0,2\pi)$ while the amplitude is selected to be $1+10^{0.5c_{k}}$ with $c_{k}$ being uniformly distributed on $$. We test two different settings for the frequencies: a) no separation condition is imposed on$\{f_{k}\}_{k=1}^{r}$, and b) the wrap-around distances between each pair of the randomly drawn frequencies are guaranteed to be greater than$1.5/n$. After a signal is formed,$m$of its entries are sampled uniformly at random. For a given triple$(n,r,m)$,$50$random tests are conducted. We consider an algorithm to have successfully reconstructed a test signal if the root mean squared error (RMSE) is less than$10^{-3}$,

We plot in Figure 1 the empirical recovery phase transition curves that identify the 80% success rate for each tested algorithm under the two different frequency settings. When the frequencies are separated by at least $1.5/n$, the right plot shows that ANM has the highest phase transition curve, and the phase transition curve of PGD closely tracks that of ANM. The performance of ANM degrades severely when there is no frequency separation requirement. In both of the frequency settings, the recovery phase transition curves of PGD are overall higher than that of EMaC. In the region of greatest interest where $p\leq 0.5$, the recovery phase transition curves of PGD are substantially higher than that of FIHT.

2 Computational Efficiency

PGD has the same leading-order computational complexity as FIHT, and both of them are able to handle large and high-dimensional signals. We compare the computational performance of these two algorithms on undamped and damped three-dimensional spectrally sparse signals of size $n=64\times 128\times 512$. Tests are conducted with $r\in\{20,30\}$ and $m\approx 130\log(n)$ in the undamped setting while $m\approx 0.03n$ in the damped setting, and we test signals which obey the frequency separation condition as well as signals which are fully random. As to the damping factors, for $1\leq k\leq r$, $1/\tau_{1k}$ is uniformly sampled from $[8~{}16]$, $1/\tau_{2k}$ is uniformly sampled from $[16~{}32]$, and $1/\tau_{3k}$ is uniformly sampled from $[64~{}128]$. For each triple of $(r,\mbox{undamped/damped, with/without separation})$, $10$ random problem instances are tested. FIHT is terminated when $\|\bm{x}^{k+1}-\bm{x}^{k}\|_{2}/\|\bm{x}^{k}\|_{2}\leq 10^{-3}$ or $\|\bm{x}^{k+1}-\bm{x}^{k}\|_{2}/\|\bm{x}^{k}\|_{2}\geq 2$ which usually implies divergence. PGD is terminated when $\|\bm{x}^{k+1}-\bm{x}^{k}\|_{2}/\|\bm{x}^{k}\|_{2}\leq 2\times 10^{-4}$. The average computational time (referred to as TIME) and average number of iterations (referred to as ITER) of FIHT and PGD over tests of successful recovery are summarized in Tables 1 and 2 for the undamped and damped signals, respectively. For the sake of completeness, we also include the ratio of successful recovery out of the 10 random tests (referred to as SR) for each algorithm in the tables.

First it is worth noting that PGD succeeded in all the $10$ random tests under each test setting when $r=30$, whereas FIHT only succeeded in a small fraction of the tests. Thus, Tables 1 and 2 show that PGD is able to more reliably recover signals that consist of a larger number of Fourier components, which coincides with our observations on one-dimensional signals in Section 3.1. The tables also show that FIHT requires fewer number of iterations and less computational time than PGD to achieve convergence for easier problem instances when $r=20$, while PGD is faster when $r=30$ and the test signals are undamped.

3 Robustness to Additive Noise

We demonstrate the performance of PGD under additive noise by conducting tests on 3D signals of the same size as in Section 3.2 but with measurements corrupted by the vector

where $\bm{x}$ is a reshaped three-dimensional spectrally sparse signal to be reconstructed, the entries of $\bm{w}$ are i.i.d. standard complex Gaussian random variables, and $\theta$ is referred to as the noise level.

Tests are conducted with $7$ different values of $\theta$ from $10^{-3}$ to 1, corresponding to $7$ equispaced signal-to-noise ratios (SNR) from 60 to 0 dB. For each value of $\theta$, 10 random instances are tested. PGD is terminated when $\|\bm{x}^{k+1}-\bm{x}^{k}\|_{2}/\|\bm{x}^{k}\|_{2}\leq 10^{-5}$. In our simulations, we fix $r=20$ and choose $m\in\{130\log(n),195\log(n)\}$ in the undamped setting while $m\in\{0.03n,0.045n\}$ in the damped setting. The frequencies of the test signals are randomly generated from $[0,1)$ without the separation requirement and the damping factors are generated in the same fashion as in Section 3.2. The average RMSE of the reconstructed signals (measured in negative dB) plotted against the input SNR values of the samples is presented in Figure 2. The plots display a desirable linear scaling between the relative reconstruction error and the noise level for both the undamped and damped signals. Moreover, the relative reconstruction error decreases linearly on a log-log scale as the number of measurements increases.

4 Sensitivity to Model Order

In practice, we may not know the exact model order of a spectrally sparse signal but only have an estimation of it. Thus, it is of great interest to examine the performance of PGD when the model order is under- or over- estimated. The experiments are conducted for three-dimensional signals of the same size as in Section 3.2. Here the true model order is $r=20$, and we observe $m=130\log(n)$ entries for undamped signals while $m=0.03n$ entries for damped signals. The frequencies are generated randomly and the damping factors are generated in the same way as in Section 3.2. Three noise levels are investigated: SNR$=\infty$ (noise-free), SNR$=20$ (light noise) and SNR$=0$ (heavy noise), and tests are conducted under the same additive noise model as in Section 3.3. For a fixed noise level, we test PGD starting from $r=5$ and then increase the value of $r$ by $5$ each time until the maximum value $40$ is reached. For each pair of $(\mbox{SNR},~{}r)$, $10$ random problem instances are tested, and PGD is terminated when $\|\bm{x}^{k+1}-\bm{x}^{k}\|_{2}/\|\bm{x}^{k}\|_{2}\leq 10^{-5}$. The median values of ITER and SNR when convergence is attained are reported in Tables 3 and 4 for undamped and damped signals, respectively. As expected, PGD achieves the best SNR when the input value of $r$ is equal to $20$ (the true model order). The SNR of the estimation is usually very low when $r$ is smaller than $20$ due to the systematic truncation error. On the other hand, even when $r$ is twice as large as the true model order, the SNR of the estimation is still desirable though it requires dramatically more number of iterations for PGD to converge.

Next, we suggest a rank increasing heuristic for PGD when the underlying model order is not known a priori. Starting from a sufficiently small $r$, we run PGD until convergence is reached (i.e., when $\|\bm{x}^{k+1}-\bm{x}^{k}\|_{2}/\|\bm{x}^{k}\|_{2}\leq 10^{-5}$). Then we compute and compare the relative residuals over the observed entries for the two successive testing values of $r$. If the relative residual is improved significantly, we increase the value of $r$; otherwise the algorithm is terminated. To validate the potential effectiveness of this heuristic, we test PGD for problem instances with SNR$=20$ for both undamped and damped signals, and with the values of $r$ increasing from 1 to 40. The computational results are presented in Figure 3, where we show the relative residual plotted against the values of $r$, as well as the change of the relative residual when $r$ is increased by one. The figure shows that when $r$ is greater than $20$, the improvement of the relative residuals becomes very marginal for both undamped and damped signals.

Proof of Theorem 2.1

The structure of the proof for Theorem 2.1 follows the typical two-step strategy in the convergence analysis of non-convex optimization algorithms: a basin of attraction is firstly established, in which the algorithm converges linearly to the true solution; and then it can be shown that the initial guess constructed in the algorithm lies inside the basin of attraction. We begin our presentation of the proof with a proposition about the initialization.

Suppose $\mathcal{G}\bm{y}$ is $\mu_{0}$-incoherent. If $m\geq c\hskip 1.42271pt\varepsilon_{0}^{-2}\mu c_{s}\kappa^{2}r^{2}\log(n)$, then one has $\bm{M}\in\mathcal{C}$ and

with probability at least $1-n^{-2}$, where in the second inequality we use the assumption $\mu_{0}\leq\mu$. Together with the assumption on $m$, it follows immediately that

Consequently, one has $\bm{M}\in\mathcal{C}$ since $\left\|\bm{M}\right\|_{2,\infty}\leq\sqrt{\sigma_{1}(\mathcal{G}\bm{y})}\max\{\left\|\bm{U}\right\|_{2,\infty},\left\|\bm{V}\right\|_{2,\infty}\}.$ Moreover, one can easily see that $\bm{M}\bm{Q}\in\mathcal{C}$ for all $r$ by $r$ unitary matrices $\bm{Q}$.

Let $\bm{A}$, $\bm{B}$, $\bm{C}$, and $\bm{D}$ be four $s\times r$ complex matrices with $s\geq r$. A simple calculation yields

where $\bm{a}_{i}$, $\bm{b}_{i}$, $\bm{c}_{i}$ and $\bm{d}_{i}$ are the $i$-th columns of $\bm{A}$, $\bm{B}$, $\bm{C}$ and $\bm{D}$ respectively. Then it follows that

which can be easily verified using (22). Substituting (23) into (21) gives

With Proposition 4.1 in place, the proof of Theorem 2.1 is complete if we can establish the local contraction property of Algorithm 1, as stated in the following proposition.

Assume $\bm{M}\in\mathcal{C}$. Let $\varepsilon_{0}$ be an absolute constant obeying $0<\varepsilon_{0}\leq\frac{1}{11}$. For any matrix $\bm{Z}\in\mathcal{C}$, define

There exists a numerical constant $\nu=\frac{1}{10}\sigma_{r}(\mathcal{G}\bm{y})$ such that with probability at least $1-c_{1}\cdot n^{-2}$,

holds for all $\bm{Z}$ obeying ${\mbox{dist}^{2}}(\bm{Z},\bm{M})\leq 3\varepsilon_{0}^{2}\sigma_{r}(\mathcal{G}\bm{y})$ provided

holds for all matrices $\bm{Z}$ within a small neighborhood of $\bm{M}$. Let $\bm{H}=\bm{Z}-\bm{M}\bm{Q}_{{\scaleto{\bm{Z}\mathstrut}{6pt}}}$. We follow a similar route as in and instead establish the regularity condition

for all matrices $\bm{Z}$ that are sufficiently close to $\bm{M}$. The notation of regularity condition was first introduced in to show the convergence of a non-convex gradient descent algorithm for phase retrieval and since then has been extended to many other problems, see and references therein. Once (25) is established, a little algebra yields

The proof of the regularity condition will occupy the remainder of this section. Even though the proof follows a well-established route, especially that in , the details of the proof are nevertheless quite involved and technical. Firstly, our objective function involves a transformation from the matrix domain to the vector domain, and an extra regularizer is also included to preserve the Hankel structure of the matrix. Secondly, we need to establish a key lemma which is closely related to the second largest eigenvalue of a special random graph, as presented in the next subsection.

The following lemma will play a key role in the proof of the regularity condition.

holds with probability at least $1-2n^{-2}$ provided $m\geq\frac{8}{3}\log(n)$.

Let $\bm{H}_{a},~{}a=0,\cdots,n-1$, be an $n_{1}\times n_{2}$ matrix with the $a$-th skew-diagonal entries being equal to one and all the other entries being equal to zero. Notice that $p^{-1}\sum_{k=1}^{m}\sum_{i+j=a_{k}}z_{i}w_{j}$ can be written as

Thus, the application of the Bernstein’s inequality (see for example [42, Theorem 1.6]) yields

Letting $t=\sqrt{\frac{24n^{2}\log(n)}{m}}$ gives

provided $m\geq\frac{8}{3}\log(n)$. Substituting this result into (26) concludes the proof. ∎

2 Proof of the Regularity Condition

The goal of this subsection is to show that the regularity condition (25) holds with high probability. Before proceeding to the formal proof, we first consider the expectation of ${\operatorname{{Re}}\left\langle\nabla F(\bm{Z}),\bm{H}\right\rangle}$ and see what lower bound can be anticipated. With a slight abuse of notation, we denote $\bm{M}\bm{Q}_{{\scaleto{\bm{Z}\mathstrut}{6pt}}}$ by $\bm{M}$ throughout this subsection for ease of presentation. Since there exists a close solution for $\bm{Q}_{{\scaleto{\bm{Z}\mathstrut}{6pt}}}$, as presented in (11), one can easily verify that

where in the second line we use $\bm{Z}=\bm{M}+\bm{H}$, and in the third line we use the inequality $a^{2}-3ab+2b^{2}\geq\frac{1}{2}a^{2}-\frac{5}{2}b^{2}$.

where the last equality follows from the fact $\bm{M}_{\scaleto{\bm{U}\mathstrut}{6pt}}\bm{M}_{\scaleto{\bm{V}\mathstrut}{6pt}}^{*}=\bm{U}\bm{\Sigma}\bm{V}^{*}$. Since $\left\|\bm{H}\right\|_{F}^{2}=\delta^{2}\left\|\bm{M}\right\|_{F}^{2}=2\delta^{2}\left\|\bm{\Sigma}\right\|_{*}$, the regularity condition (25) cannot be true for $f(\bm{Z})$ without the regularization function $g(\bm{Z})$. In this case, one can observe that the mismatch between $\bm{Z}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}^{*}\bm{Z}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}$ and $\bm{Z}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}^{*}\bm{Z}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}$ increases compared with the mismatch between $\bm{M}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}^{*}\bm{M}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}$ and $\bm{M}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}^{*}\bm{M}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}$ which is equal to zero. Because $g(\bm{Z})$ penalizes the mismatch between $\bm{Z}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}^{*}\bm{Z}_{{\scaleto{\bm{U}\mathstrut}{6pt}}}$ and $\bm{Z}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}^{*}\bm{Z}_{{\scaleto{\bm{V}\mathstrut}{6pt}}}$, one may intuitively expect that it can control the occurrence of this case so that $F(\bm{Z})=f(\bm{Z})+\lambda g(\bm{Z})$ could obey the regularity condition.

Let $\bm{D}=\begin{bmatrix}\bm{I}_{n_{1}}&\bm{0}\\ \bm{0}&-\bm{I}_{n_{2}}\end{bmatrix}$. We can bound $\operatorname{{Re}}\left\langle\nabla g(\bm{Z}),\bm{H}\right\rangle$ from below as

where the third equality follows from $\bm{M}^{*}\bm{D}\bm{M}=\bm{0}$, the fourth equality follows from

and the inequality follows from $\bm{H}^{*}\bm{M}=\bm{M}^{*}\bm{H}$, see (27).

If we take $\lambda=\frac{1}{4}$, then combining (28) and (29) together implies

That is, we have established a lower bound for the expectation of $\operatorname{{Re}}\left\langle\nabla F(\bm{Z}),\bm{H}\right\rangle$. As we will show later, $\operatorname{{Re}}\left\langle\nabla F(\bm{Z}),\bm{H}\right\rangle$ obeys a similar lower bound with high probability. Moreover, the right hand side of (25) can be bounded from above by a similar bound. Therefore, $F(\bm{Z})$ obeys the regularity condition for sufficiently small $\bm{H}$. Specifically, we are going to show the following two bounds,

hold with high probability provided $\left\|\bm{H}\right\|_{F}^{2}\leq 3\varepsilon_{0}^{2}\sigma_{r}(\mathcal{G}\bm{y})$ and $m\gtrsim\varepsilon_{0}^{-2}\mu^{2}c_{s}^{2}\kappa^{2}r^{2}\log(n)$ for $\varepsilon_{0}\leq\frac{1}{11}$. The above two inequalities are typically referred to as the local curvature property and the local smooth property of the function $F(\bm{Z})$ in the literature, see for example . Once they are established, one can easily see that $F(\bm{Z})$ obeys the regularity condition (25) with

Since $\operatorname{{Re}}\left\langle\nabla g(\bm{Z}),\bm{H}\right\rangle$ is deterministic and we have already obtained its lower bound in (29), it only remains to work out the lower bound for $\operatorname{{Re}}\left\langle\nabla f(\bm{Z}),\bm{H}\right\rangle$ and then combine it together with that for $\operatorname{{Re}}\left\langle\nabla g(\bm{Z}),\bm{H}\right\rangle$. Note that

where the second equality follows from the fact $(\mathcal{I}-\mathcal{G}\mathcal{G}^{*})(\bm{M}_{\scaleto{\bm{U}\mathstrut}{6pt}}\bm{M}_{\scaleto{\bm{V}\mathstrut}{6pt}}^{*})=\bm{0}$.

Lower bound for $I_{1}$. The first term $I_{1}$ can be bounded directly as follows:

where the first equality follows from that $\mathcal{G}\mathcal{G}^{*}$ is a projection operator, and the second inequality follows from $a^{2}-3ab+2b^{2}\geq\frac{11}{20}a^{2}-3b^{2}$.

Lower bound for $I_{2}$. Recall from Section 2.1 that $w_{a}$, $a=0,\cdots,n-1$, denotes the number of entries in the skew-diagonal of an $n_{1}\times n_{2}$ matrix. Let $\bm{G}_{a},~{}a=0,\cdots,n-1$, be an $n_{1}\times n_{2}$ matrix with the $a$-th skew-diagonal entries being equal to $1/\sqrt{w_{a}}$ and all the other entries being equal to zero. Then,

where the third equality and the last equality follow from (16), the first inequality follows from the Hölder inequality, and the second inequality follows from $a^{2}-3ab+2b^{2}\geq\frac{11}{20}a^{2}-3b^{2}$. Consequently,

We can bound $p^{-1}\sum_{k=1}^{m}\left|\left\langle\bm{G}_{a_{k}},\bm{H}_{\scaleto{\bm{U}\mathstrut}{6pt}}\bm{H}_{\scaleto{\bm{V}\mathstrut}{6pt}}^{*}\right\rangle\right|^{2}$ from above by Lemma 4.1 as follows:

where the fourth line follows from Lemma 4.1, the sixth line follows from

and the last line follows from (19) and the assumptions on $\left\|\bm{H}\right\|_{F}^{2}$ and $m$.

Lower bound for $\operatorname{{Re}}\left\langle\nabla f(\bm{Z}),\bm{H}\right\rangle$. Before finally showing the lower bound for $\operatorname{{Re}}\left\langle\nabla f(\bm{Z}),\bm{H}\right\rangle$, we need to define the tangent space of the rank $r$ matrix manifold at $\mathcal{G}\bm{y}$, denoted $T$. Given the SVD $\mathcal{G}\bm{y}=\bm{U}\bm{\Sigma}\bm{V}^{*}$, we define $T$ as

One can easily see that $\bm{H}_{\scaleto{\bm{U}\mathstrut}{6pt}}\bm{M}_{\scaleto{\bm{V}\mathstrut}{6pt}}^{*}+\bm{M}_{\scaleto{\bm{U}\mathstrut}{6pt}}\bm{H}_{\scaleto{\bm{V}\mathstrut}{6pt}}^{*}\in T$. Substituting the bound for $p^{-1}\sum_{k=1}^{m}\left|\left\langle\bm{G}_{a_{k}},\bm{H}_{\scaleto{\bm{U}\mathstrut}{6pt}}\bm{H}_{\scaleto{\bm{V}\mathstrut}{6pt}}^{*}\right\rangle\right|^{2}$ into (34) and then combining the lower bounds for $I_{1}$ and $I_{2}$ together yields

where the second inequality follows from the fact $\mathcal{G}\mathcal{G}^{*}$ is a projection operator, the third inequality holds with probability at least $1-n^{-2}$ (see Lemma A.3) under the assumption on $m$ and $\left\|\bm{H}\right\|_{F}^{2}$, and the last inequality follows from $\left\|\bm{H}_{\scaleto{\bm{U}\mathstrut}{6pt}}\bm{M}_{\scaleto{\bm{V}\mathstrut}{6pt}}^{*}\right\|_{F}\geq\sigma_{r}(\bm{M}_{{\scaleto{\bm{V}\mathstrut}{6pt}}})\left\|\bm{H}_{\scaleto{\bm{U}\mathstrut}{6pt}}\right\|_{F}$, $\left\|\bm{M}_{\scaleto{\bm{U}\mathstrut}{6pt}}\bm{H}_{\scaleto{\bm{V}\mathstrut}{6pt}}^{*}\right\|_{F}\geq\sigma_{r}(\bm{M}_{{\scaleto{\bm{U}\mathstrut}{6pt}}})\left\|\bm{H}_{\scaleto{\bm{V}\mathstrut}{6pt}}\right\|_{F}$, and the assumption $\varepsilon_{0}\leq\frac{1}{11}$.

Lower bound for $\operatorname{{Re}}\left\langle\nabla F(\bm{Z}),\bm{H}\right\rangle$. Let $\lambda=\frac{1}{4}$. Combining the lower bound in (35) for $\operatorname{{Re}}\left\langle\nabla f(\bm{Z}),\bm{H}\right\rangle$ and the lower bound in (29) for $\operatorname{{Re}}\left\langle\nabla g(\bm{Z}),\bm{H}\right\rangle$ together gives

and the assumption $\varepsilon_{0}\leq\frac{1}{11}$. This concludes the proof of (31).

2.2 Proof of (32)

it suffices to bound $\left\|\nabla f(\bm{Z})\right\|_{F}^{2}$ and $\left\|\nabla g(\bm{Z})\right\|_{F}^{2}$ separately.

Upper bound for $\left\|\nabla g(\bm{Z})\right\|_{F}^{2}$. We begin with the upper bound for $\left\|\nabla g(\bm{Z})\right\|_{F}^{2}$, which can be obtained in a straightforward way,

where the third equality follows from $\bm{Z}=\bm{M}+\bm{H}$ and $\bm{M}^{*}\bm{D}\bm{M}=\bm{0}$, the fourth inequality follows from $\left\|\bm{M}\right\|_{2}=\sqrt{2\sigma_{1}(\mathcal{G}\bm{y})}$ and

and the last line follows from the assumption $\varepsilon_{0}\leq\frac{1}{11}$.

In order to bound $\left\|\nabla f(\bm{Z})\right\|_{F}^{2}$, we consider $\left|\left\langle\nabla f(\bm{Z}),\bm{X}\right\rangle\right|^{2}$ for matrices $\bm{X}=\begin{bmatrix}\bm{X}_{\scaleto{\bm{U}\mathstrut}{6pt}}&\bm{X}_{\scaleto{\bm{V}\mathstrut}{6pt}}\end{bmatrix}^{T}$ with unit Frobenius norm (i.e., $\left\|\bm{X}_{\scaleto{\bm{U}\mathstrut}{6pt}}\right\|_{F}^{2}+\left\|\bm{X}_{\scaleto{\bm{V}\mathstrut}{6pt}}\right\|_{F}^{2}=1$). Note that

where in the last line we have utilized $\left\|\bm{X}_{\scaleto{\bm{U}\mathstrut}{6pt}}\right\|_{F}^{2}+\left\|\bm{X}_{\scaleto{\bm{V}\mathstrut}{6pt}}\right\|_{F}^{2}=1$. Because $\mathcal{I}-\mathcal{G}\mathcal{G}^{*}$ is a projection operator, $I_{3}$ can be bounded as follows:

We can bound $\left|\left\langle p^{-1}\mathcal{G}\mathcal{P}_{\Omega}\mathcal{G}^{*}(\bm{Z}_{\scaleto{\bm{U}\mathstrut}{6pt}}\bm{H}_{\scaleto{\bm{V}\mathstrut}{6pt}}^{*}),\bm{X}_{\scaleto{\bm{U}\mathstrut}{6pt}}\bm{Z}_{\scaleto{\bm{V}\mathstrut}{6pt}}^{*}\right\rangle\right|$ as follows:

where in the last line, we utilize $\left\|\bm{Z}\right\|_{2,\infty}^{2}\leq\mu c_{s}r\sigma/n$. Similar upper bounds can be established for the other three terms in (39). That is,

Combining these four upper bounds together yields

where in the last line we have used the fact $\left\|\bm{X}_{\scaleto{\bm{U}\mathstrut}{6pt}}\right\|_{F}^{2}+\left\|\bm{X}_{\scaleto{\bm{V}\mathstrut}{6pt}}\right\|_{F}^{2}=1$.

Upper bound for $\left\|\nabla f(\bm{Z})\right\|_{F}^{2}$. Substituting the upper bounds for $I_{3}$ and $I_{4}$ into (38) give the upper bound for $\left\|\nabla f(\bm{Z})\right\|_{F}^{2}$,

Upper bound for $\left\|\nabla F(\bm{Z})\right\|_{F}^{2}$. Noting $\lambda=1/4$, $\sigma\leq(1+\varepsilon_{0})\sigma_{1}(\mathcal{G}\bm{y})/(1-\varepsilon_{0})$, and $\varepsilon_{0}\leq 1/11$, after substituting the upper bounds for $\left\|\nabla f(\bm{Z})\right\|_{F}^{2}$ and $\left\|\nabla g(\bm{Z})\right\|_{F}^{2}$ into (36), we get

Discussion

We have proposed a novel algorithm for spectral compressed sensing by applying projected gradient descent updates to a non-convex functional. Exact recovery guarantee has been established, showing that $O(r^{2}\log(n))$ random observations are sufficient for the algorithm to achieve the successful recovery. Additionally, empirical evaluation shows that our algorithm is competitive with other state-of-the-art algorithms. In particular, our algorithm is superior to FIHT, a non-convex algorithm for spectral compressed sensing with provable recovery guarantees, in terms of phase transitions when the number of observations is small.

For future work, recovery stability of the proposed algorithm to additive noise will be investigated. The proofs presented in this paper should extend easily to bounded noise with a small magnitude. It remains to address whether or not our algorithm can achieve some statistically optimal rates under a stochastic noise model.

Recently, a line of research work has been devoted to the geometric analysis of non-convex optimization problems including dictionary learning , phase retrieval , low rank matrix sensing and matrix completion , tensor completion and robust PCA . It has been shown that the non-convex functionals for those problems have well-behaved landscape: all local minima are also globally optimal. Preliminary numerical results show that our projected gradient descent algorithm works equally well with random initialization, which suggests the geometric landscape of the objective function $F(\bm{Z})$ introduced in this paper may be similarly well-behaved.

Appendix A Supplementary Lemmas

Here we list three technical lemmas from the literature that have been used in the analysis of PGD.

Assume $\mathcal{G}{\bm{y}}$ is $\mu_{0}$-incoherent and let $\bm{L}_{0}=\mathcal{T}_{r}\mathcal{G}(p^{-1}\mathcal{P}_{\Omega}{(\bm{y})})$. Then,

holds with probability at least $1-n^{-2}$.

Assume $\mathcal{G}{\bm{y}}$ is $\mu_{0}$-incoherent, and let $T$ be the tangent space of the rank $r$ matrix manifold at $\mathcal{G}\bm{y}$. Then,

holds with probability at least $1-n^{-2}$.

References