Large Scale 2D Spectral Compressed Sensing in Continuous Domain

TL;DR

This paper introduces a semidefinite programming approach for large-scale 2D spectral compressed sensing in the continuous domain, enabling recovery of high-resolution signals from limited samples.

Contribution

It presents a novel semidefinite program that scales to 500x500 signals, significantly surpassing traditional methods in size handling.

Findings

Successfully recovers 2D spectrally sparse signals from partial observations.

Handles large-scale signals of size up to 500x500.

Outperforms traditional methods limited to 20x20 signals.

Abstract

We consider the problem of spectral compressed sensing in continuous domain, which aims to recover a 2-dimensional spectrally sparse signal from partially observed time samples. The signal is assumed to be a superposition of s complex sinusoids. We propose a semidefinite program for the 2D signal recovery problem. Our model is able to handle large scale 2D signals of size 500*500, whereas traditional approaches only handle signals of size around 20*20.

Tables

Table 1. Table 1 : Comparison of running times (in seconds) between existing approach [ 12 ] and the proposed SDP formulation ( 2.4 ) for different sizes n 𝑛 n and observations m 𝑚 m

$n$	$m$	[12] via CVX	(2.4) via CVX	(2.4) via ADMM
15	80	38.52	0.55	0.050
16	90	60.57	0.26	0.059
17	100	103.47	0.29	0.060
18	110	161.27	0.49	0.063
19	120	239.09	0.39	0.064
20	130	363.61	0.43	0.074
21	140	534.03	0.68	0.072
22	150	850.19	0.54	0.076
23	160	crashed	0.60	0.087

Equations

$$\begin{array}{rl}{X}^{★}(j,k)& =\sum _{p=1}^s{c}_p{e}^{i2\pi (f_{p1}j+f_{p2}k)},\\ (j,k)& \in T,\end{array}$$

$${\mathbf{x}}^{★}(j)=\sum _{p=1}^s{c}_p{e}^{i2\pi f_{p}j},j=0,\dots ,n-1,$$

$$X^{★}$$

$$=\sum _{p=1}^s{c}_p{{\mathbf{v}}_{\mathbf{p}1}{\mathbf{v}}_{\mathbf{p}2}}^H$$

$$C$$

$$V_1$$

$$\displaystyle=\left[\begin{array}[]{ccc}e^{i2\pi f_{11}0}&\cdots&e^{i2\pi f_{s1}0}\\ \vdots&\ddots&\vdots\\ e^{i2\pi f_{11}(n-1)}&\cdots&e^{i2\pi f_{s1}(n-1)}\end{array}\right],$$

$$V_2$$

$$\displaystyle=\left[\begin{array}[]{ccc}e^{-i2\pi f_{12}0}&\cdots&e^{-i2\pi f_{s2}0}\\ \vdots&\ddots&\vdots\\ e^{-i2\pi f_{12}(n-1)}&\cdots&e^{-i2\pi f_{s2}(n-1)}\end{array}\right].$$

$$\underset{X,V_1,C,V_2}{\min }$$

$$X=V_{1}C{V}_2^H$$

$$\text{where}{V}_1,V_2\text{are full rank Vandermonde matrices},$$

$$\displaystyle C~{}\text{diagonal as in (\ref{form_matrix})},$$

$$X(j,k)=X^{★}(j,k),(j,k)\in T,$$

$$\underset{T_1,T_2,X}{\min }$$

$$\displaystyle\left[\begin{array}[]{cc}T_{1}&X\\ X^{H}&T_{2}\end{array}\right]\succeq 0,$$

$$T_1,T_2\text{are Toeplitz},$$

$$X(j,k)=X^{★}(j,k),(j,k)\in T,$$

$$\underset{W_1,W_2,X}{\min }$$

$$\displaystyle\left[\begin{array}[]{cc}W_{1}&X\\ X^{H}&W_{2}\end{array}\right]\succeq 0,$$

$$X(j,k)=X^{★}(j,k),(j,k)\in T,$$

$$\underset{X}{\min }||X|{|}_{\ast }$$

$$X(j,k)=X^{★}(j,k),(j,k)\in T.$$

$$T_1=U_1{D}_1{U}_1^H,T_2=U_2{D}_2{U}_2^H$$

$$X=U_{1}A{U}_2^H$$

$$\underset{M,N}{\min }$$

$$M⪰0,$$

$$\displaystyle N=\left[\begin{array}[]{cc}T_{1}&X\\ X^{H}&T_{2}\end{array}\right],$$

$$T_1,T_2\text{are Toeplitz}$$

$$X(j,k)=X^{★}(j,k),(j,k)\in T,$$

$$M-N=0.$$

$$M^{k+1}$$

$$\text{s.t.}M⪰0,$$

$$N^{k+1}$$

$$\displaystyle\text{s.t.}~{}~{}N=\left[\begin{array}[]{cc}T_{1}&X\\ X^{H}&T_{2}\end{array}\right],$$

$$T_1,T_2\text{are Toeplitz},$$

$$X(j,k)=X^{★}(j,k),(j,k)\in T,$$

$$U^{k+1}$$

$$M^{k+1}=V{D}_{+}{V}^H,$$

$$T_1^{k+1}(i+l,i)=\sum _{i=\max \{1,1-l\}}^{\min \{n-l,n\}}[M_1^{k+1}+U_1^k](i+l,i)/(n-l),$$

$$T_2^{k+1}(i+l,i)=\sum _{i=\max \{1,1-l\}}^{\min \{n-l,n\}}[M_2^{k+1}+U_2^k](i+l,i)/(n-l),$$

$$\text{for each }l=-(n-1),\dots ,0,\dots ,n-1$$

$$\text{and }\forall i=\max \{1,1-l\},\dots ,\min \{n-l,n\}.$$

$$\displaystyle X^{k+1}(j,k)=\left\{\begin{array}[]{ll}X^{\bigstar}(j,k),&(j,k)\in T,\\ \\ \left[M_{3}^{k+1}+U_{3}^{k}\right](j,k),&\text{otherwise.}\end{array}\right.$$

$$\Vert r^k{\Vert }_F$$

Full text

**Index Terms— ** Compressed sensing, sparse recovery, Toeplitz matrices, matrix completion

1 Introduction

Spectral compressed sensing aims to recover a spectrally sparse signal from sub-Nyquist sampling of time samples. Consider a frequency-sparse signal $X^{\bigstar}$ which consists of a superposition of $s$ complex sinusoids with 2-dimensional frequencies $\bm{f}_{p}=(f_{p1},f_{p2})$, where $f_{pi}\in[0,1),i=1,2$. Suppose we observe $m$ regularly-spaced time samples at times $(j,k)$ which form an index set $T\subset\{0,\ldots,n-1\}\times\{0,\ldots,n-1\}$. We would like to recover the signal of the following form:111Here we consider square matrices for simplicity, non-square cases can be readily generalized.

[TABLE]

where $c_{p}$ is the coefficient of each complex sinusoid.

In the foundational works of Candès et al. [1] and Donoho [2], successful signal recoveries from very few samples were guaranteed with high probabilities under the assumption that the unknown frequencies lie on a grid. However, in real-world applications, frequencies can possibly lie anywhere on the continuous domain [3, 4]. If they are still assumed to lie on a certain grid, the recovery might suffer from the “basis mismatch” issue and the recovery might be far from the original true signal [5]. Other approaches for frequency-sparse signal recovery have also been proposed, for example, finite rate of innovation (FROI) sampling [6, 7]. Recently, people have been working on the recovery based on continuous frequency domain. Tang et al. [8] showed that $\mathcal{O}(s\log s\log n)$ random samples guarantee exact frequency recovery of a 1D spectrally sparse signal

[TABLE]

with high probability over continuous frequency domain, given the condition that the frequencies are well separated. An atomic norm minimization approach was proposed therein, which was reformulated as an equivalent solvable semidefinite program (SDP). Despite the 1D case, there are more applications involving the recovery of $d$-dimensional signals ($d\geq 2$). For example, super-resolution imaging [9] and nuclear magnetic resonance (NMR) spectroscopy [10]. Chi and Chen [11, 12] extended Tang’s approach to 2D frequency models. However, the equivalence between the atomic norm minimization and the proposed SDP was not guaranteed. Xu et al. [13] proposed equivalent SDP formulations of atomic norm minimization to recover signals with 2D off-the-grid frequencies using theories from trigonometric polynomials. Moreover, the SDP formulations were generalized to higher dimensions as well.

In Chi’s and Xu’s approaches, the $n\times n$ signal matrix $X$ appeared as an $n^{2}\times 1$ vector in the positive semidefinite constraints. The resulting SDP was on a space of dimension at least $\mathcal{O}(n^{2})\times\mathcal{O}(n^{2})$. For example, to recover a $100\times 100$ signal $X$, we would have to solve an SDP of size at least $10001\times 10001$, which is computationally prohibitive.

In this paper, we propose a model to recover 2D spectrally sparse signals, which reduces the size of the SDP to $\mathcal{O}(n)\times\mathcal{O}(n)$, while at the same time requires a small number of observed samples to guarantee exact recovery with high probability, which means the model is able to handle large scale 2D signals compared to existing approaches.

The remainder of this paper is organized as follows. In Section 2, we introduce our proposed SDP formulation for large scale 2D signal recovery. In Section 3, we give numerical experiments to validate our model.

2 Semidefinite Formulation for 2D Signal Recovery

Notice that the 2D signal (1.1) has the matrix form

[TABLE]

where $C$ is a diagonal matrix and $V_{1}$ and $V_{2}$ are Vandermonde matrices of the following form222In (2.2), the negative exponents in entries of $V_{2}$ arise from the conjugate transpose of $\bm{v_{p2}}$ in (2.1) in order to be consistent with the original form of the signal (1.1) .:

[TABLE]

We assume that $s\ll n$. The $i$th $(i=1,2)$ frequency components $f_{pi}$ are distinct across $p$ $(p=1,\ldots s)$, so that $\text{rank}~{}(X^{\bigstar})=s$ and can be considered low-rank. The recovery problem can thus be formulated as the following rank minimization problem:

[TABLE]

where $T$ is the set of indices of the observed entries and $|T|=m$.

We propose the following semidefinite program for Problem (2)

[TABLE]

The advantage of this SDP over the ones proposed in [12, 13] is that the size of the semidefinite constraint is $2n\times 2n$ instead of $(n^{2}+1)\times(n^{2}+1)$, which significantly reduces the memory needed for computation.

2.1 Why it works

Our SDP model (2.4) is inspired from the SDP model in [14]

[TABLE]

which is equivalent to the problem of minimizing the nuclear norm of a matrix $X$:

[TABLE]

It has been proved to be an effective approach in the low-rank matrix completion problem [15].

Our model (2.4) for 2D signal recovery can be viewed as an extension of the low-rank matrix completion problem in the following sense:

First, model (2.4) will yield a low-rank solution. Notice that it minimizes the trace of a positive semidefinite matrix \left[\begin{array}[]{cc}T_{1}&X\\ X^{H}&T_{2}\end{array}\right], which is equivalent to minimizing the sum of all the (nonnegative) eigenvalues. It can be viewed as a relaxation of minimizing the number of nonzero eigenvalues, therefore minimizing rank\left(\left[\begin{array}[]{cc}T_{1}&X\\ X^{H}&T_{2}\end{array}\right]\right) and hence rank$(X)$.

Second, we would like to incorporate the additional structure of the true signal matrix $X^{\bigstar}=V_{1}CV_{2}^{H}$ into the recovered solution $X_{\text{rec}}$, therefore imposing the Toeplitz constraints on matrices $T_{1}$ and $T_{2}$ so that the problem becomes our model (2.4). This idea works due to the Caratheodory-Toeplitz Lemma [8], which states that for Toeplitz positive semidefinite matrices $T_{1}$ and $T_{2}$, we have the following decompositions:

[TABLE]

where $D_{1},D_{2}$ are $r\times r$ diagonal matrices and $U_{1},U_{2}$ are full rank $n\times r$ Vandermonde matrices. Following the similar argument in [8], the Vandermonde decomposition (2.7) together with the positive semidefinite constraint in (2.4) imply that $X$ must have the following form

[TABLE]

Notice that $U_{l}$ in (2.8) have exactly the same form as $V_{l}$ in (2.2) $(l=1,2)$ (they are all Vandermonde matrices). The only difference between the true signal $X^{\bigstar}=V_{1}CV_{2}^{H}$ and the recovered signal $X_{\text{rec}}=U_{1}AU_{2}^{H}$ is that the matrix $C$ is required to be diagonal in the true signal, whereas in the recovered signal, there is no such restriction on the matrix $A$, which means that the recovered signal is a structural approximation of the true signal. Therefore, by solving Problem (2.4), we expect to get a low-rank solution $X_{\text{rec}}$ that resembles the form of the true solution $X^{\bigstar}$.

In a nutshell, (2.4) can be viewed as a low-rank matrix completion problem where the matrix has the additional structure of Vandermonde factors. We will demonstrate through numerical experiments in Section 3 that the solution $X_{\text{rec}}$ is exactly $X^{\bigstar}$ with high probability given only a small number of randomly observed time samples.

3 Numerical Experiments

We solve (2.4) by using the alternating direction method of multipliers (ADMM). To synthesize the $n\times n$ true signal matrix $X^{\bigstar}$, the $s$ frequency pairs are randomly drawn from $[0,1]^{2}$ with a minimum separation condition [8] $\Delta_{\text{min}}\triangleq\min_{p\neq q}\{|f_{p1}-f_{q1}|,|f_{p2}-f_{q2}|\}\geq 1/n$, where $p,q=1,\ldots,s$. The amplitudes $|c_{p}|$ are drawn randomly from the distribution $0.5+w^{2}$ with $w$ a zero mean unit variance Gaussian random variable. The phases $e^{i\phi_{p}}$ are drawn uniformly at random in $[0,2\pi)$. A total of $m$ observed entries are randomly chosen from $X^{\bigstar}$. Different values of parameters $m$ and $s$ are used in the implementations of ADMM below to investigate the dependence of $m$ on $s$. The recovery is considered successful if the recovered signal $X_{\text{rec}}$ and the true synthetic signal $X^{\bigstar}$ satisfy $\frac{\|X_{\text{rec}}-X^{\bigstar}\|_{F}}{\|X^{\bigstar}\|_{F}}\leq 10^{-3}$.

3.1 ADMM Details

To apply ADMM, we first reformulate (2.4) as

[TABLE]

This new formulation splits the semidefinite constraint and the Toeplitz constraint so that they contain $M$ and $N$ separately.

ADMM consists of the following three updates:

[TABLE]

where $\rho$ is a user-specified parameter in ADMM [16] . We pick $\rho=0.1$ in the experiments.

M-update can be simplified to

[TABLE]

where $D_{+}$ is $D$ with all negative diagonal entries replaced with 0. And $VDV^{H}$ is the eigendecomposition of $N^{k}-U^{k}-\frac{1}{\rho}I$.

N-update contains 4 blocks that can be simplified to

[TABLE]

ADMM terminates when both the primal residual $r^{k}=M^{k}-N^{k}$ and dual residual $s^{k}=\rho(N^{k-1}-N^{k})$ satisfy the following conditions, as detailed in [16].

[TABLE]

where $\epsilon^{\text{abs}}$ and $\epsilon^{\text{rel}}$ are user-specified parameters and here we set both of them to be $10^{-5}$. Once the algorithm terminates, the upper right block of $N$ is the output $X_{\text{rec}}$. The entire ADMM algorithm can be summarized as:

3.2 Phase Transition

We investigate the dependence of $m$ on $s$. Here we fix $n=50$ and implement ADMM with different $m$ and $s$ values. For each $(m,s)$ pair, the simulation is repeated 20 times. The gray level indicates the number of successful recoveries among all 20 repetitions. In Figure 1, we compare between the phase transition plots of the proposed model (2.4) and the nuclear norm minimization in [14] (i.e., after removing the Toeplitz constraint in (2.4)).

We can see that the Toeplitz constraint in our model improves the phase transition because it captures the structure of the signal (i.e., the Vandermonde factors). Also, Figure 1 (left) indicates that the dependence of $m$ on $s$ is linear.

3.3 Comparison with Existing Approach

Here we compare our approach with Chi’s approach (i.e., solving SDP (16) in [12]). We implemented all the experiments on a PC with a 8-core Intel i7-4790 3.60 GHz CPU and 16 GB memory. First, we synthesized 2D signal matrices with different sizes $n$ from $s=5$ frequency components and randomly drew $m$ time samples. Then we implemented Chi’s approach via CVX, and our proposed SDP (2.4) via both CVX and ADMM. The running times are listed in Table 1.

We can see that as the size of the SDP slowly increases, the time required to solve the SDP in [12] grows quickly, whereas the time required for our proposed SDP remains small. Actually, as $n$ goes beyond 23, the SDP in [12] crashes since the size of the SDP is too large to compute. On the contrary, our proposed SDP is still capable of handling large scale problems in a reasonable amount of time, which will be demonstrated in Section 3.4.

3.4 A Large Scale Example

We implemented ADMM on SDP (2.4) to recover a $500\times 500$ signal generated from $s=10$ random frequencies on $[0,1)\times[0,1)$ with $m=5000$ ($2\%$) randomly observed samples. The ADMM implementation details are demonstrated in Figure  2. The relative error $\frac{\|X^{\bigstar}-X_{\text{rec}}\|_{F}}{\|X^{\bigstar}\|_{F}}=1.4633\times 10^{-4}$, indicating successful recovery. The Fast ADMM by Goldstein et al.[17] was applied to accelerate convergence. The running time is around 24 min on the same PC.

The size of SDP (2.4) for this problem is $1000\times 1000$. If the approaches in [12, 13] were adopted, it would become $250,001\times 250,001$, which would be impossible to compute via existing solvers. Therefore, our approach achieved better performance in terms of memory savings.

4 Conclusion

In this paper, we proposed a model for the recovery of large scale 2D spectrally sparse signal in continuous domain. This model is able to handle 2D signal matrices with much larger size ($500\times 500$) compared with existing approaches (around $20\times 20$). We demonstrated through numerical experiment that the proposed SDP formulation indeed outperforms existing approaches in terms of running time and memory.