Spectrally Sparse Signal Recovery via Hankel Matrix Completion with Prior Information

Spectrally sparse signal recovery refers to recovering a spectrally sparse signal from a small number of time domain samples, which is fundamental in various applications, such as medical imaging[1], radar imaging[2], analog-to-digital conversion[3] and channel estimation[4]. Let $\bm{x}=[x_{0},\ldots,x_{n-1}]^{T}\in\mathbb{C}^{n}$ denote the one-dimensional spectrally sparse signal to be estimated. Each entry of the desired signal $\bm{x}$ is a weighted superposition of $r$ complex sinusoids

where $k=0,\ldots,n-1$, $\{f_{1},\ldots,f_{r}\}$ and $\{w_{1},\ldots,w_{r}\}$ denote the normalized frequencies and amplitudes for the $r$ sinusoids, respectively, and $f_{l}\in[0,1)$ for $l=1,\ldots,r$.

In many practical applications, we only have access to a small subset of signal samples. For example, in the field of computed tomography (CT), only part of the desired signals can be observed to protect the patients from high-dose radiation[5]; in wideband signal sampling, it’s challenging to build analog-to-digital converter according to Shannon sampling theorem, and hence only a small number of samples of the wideband signals can be acquired for reconstruction[3]. Therefore, we have to figure out a way to recover the original signal $\bm{x}$ from its random undersampled observations

where $\Omega\in\{0,\ldots,n-1\}$ denote the index set of the entries we observe, $\bm{e}_{k}$ denotes the $k$-th canonical basis of $\mathbb{R}^{n}$, and $\mathcal{P}_{\Omega}(\cdot)$ denotes the projection operator on the sampling index set $\Omega$, i.e., $\mathcal{P}_{\Omega}(\bm{z})=\sum_{k\in\Omega}\left\langle\bm{z},\bm{e}_{k}\right\rangle\bm{e}_{k}$ for $\bm{z}\in\mathbb{C}^{n}$.

In order to reconstruct $\bm{x}$, structured low-rank completion methods have been proposed by using the low-rank Hankel structure of the spectrally sparse signals in the lifting domain

where $\text{Rank}(\cdot)$ returns the rank of matrix, and $\mathcal{H}:\mathbb{C}^{n}\to\mathbb{C}^{n_{1}\times n_{2}}$ is a linear lifting operator to generate the Hankel low-rank structure. In particular, for a vector $\bm{x}=[x_{0},x_{1},\ldots,x_{n-1}]^{T}\in\mathbb{C}^{n}$, the Hankel matrix $\mathcal{H}(\bm{x})$ is defined as

where $d$ denotes the matrix pencil parameter, $n_{1}=n-d+1$ and $n_{2}=d$.

By using Vandermonde decomposition, the Hankel matrix $\mathcal{H}(\bm{x})$ can be decomposed as

where $\bm{y}_{l}=[1,e^{i2\pi f_{l}},\ldots,e^{i2\pi(n_{1}-1)f_{l}}]^{T}$ and $\bm{z}_{l}=[1,e^{i2\pi f_{l}},\ldots,e^{i2\pi(n_{2}-1)f_{l}}]^{T}$, $l=1,\ldots,r$. When the frequencies are all distinct and $r\ll\min\{n_{1},n_{2}\}$, $\mathcal{H}(\bm{x})$ is a low-rank matrix with $\text{Rank}(\mathcal{H}(\bm{x}))\leq r$.

Since Eq. (1) is a non-convex problem and solving it is NP-hard, an alternative approach based on convex relaxation is proposed to complete the low rank matrix, that is, Hankel matrix completion program [6, 7]

where $\left\lVert\cdot\right\rVert_{*}$ denotes the nuclear norm. Theoretical analysis was given to show that $O(r\log^{4}n)$ samples are enough to recover the desired signal with high probability [6].

Apart from the sparsity constraint in the spectral domain, a reference signal $\bm{\phi}$ that is similar to the original signal $\bm{x}$ sometimes is available to us. There are two main sources of this kind of prior information. The first source comes from natural (non-self-constructed) signals. In high resolution MRI [8, 9, 10], adjacent slices show good similarity with each other; in multiple-contrast MRI [11, 12, 13], different contrasts in the same scan are similar in structure; in dynamic CT [14], the scans for the same slice in different time have similar characteristics. The second source comes from self-constructed signals. One way is to use classical method to construct a similar signal. For example, filtered backprojection reconstruction algorithm from the dynamic scans was used to construct the prior information in dynamic CT [15]; smooth method was used to generate prior information in sparse-view CT [16]; the standard spectrum of dot object was used as the prior information in radar imaging [17]. The other way is to use machine learning to generate a similar signal. In [18], the authors generated the reference image by using a CNN network; similarly, other algorithms from deep learning can be used to create reference signals, see e.g. in [19, 20].

In this paper, we propose a convex approach to integrate prior information into the reconstruction of spectrally sparse signals by maximizing the correlation between signal $\bm{z}$ and prior information $\bm{\phi}$ in the lifting domain

where $\lambda>0$ is a tradeoff parameter, $\mathcal{G}(\cdot)=\mathcal{\mathcal{F}(\mathcal{H}(\cdot))}$ is a composition operator, $\left\langle\bm{X},\bm{Y}\right\rangle=(\mbox{vec}(\bm{Y}))^{H}\mbox{vec}(\bm{X})$ is the inner product and $\text{Re}(\cdot)$ returns the real part of a complex number. Here, $\mathcal{F}:\mathbb{C}^{n_{1}\times n_{2}}\to\mathbb{C}^{n_{1}\times n_{2}}$ is a suitable operator to be designed in the sequel. Theoretical guarantees are provided to show that our method has better performance than vanilla Hankel matrix completion when the prior information is reliable. In addition, we propose an Alternating Direction Method of Multipliers (ADMM)-based optimization algorithm for efficient reconstruction of the desired signals.

In this section, we introduce some important notation and definitions, which will be used in the sequel.

Let $\{\bm{A}_{k}\}_{k=0}^{n-1}\in\mathbb{C}^{n_{1}\times n_{2}}$ be an orthonormal basis of Hankel matrices [7, 34], which is defined as

where

and $|S|$ returns the cardinality of the set $S$. Then $\mathcal{H}(\bm{x})$ can be expressed as

where $\mathcal{A}_{k}(\bm{X})\triangleq\left\langle\bm{X},\bm{A}_{k}\right\rangle\bm{A}_{k}$ for $\bm{X}\in\mathbb{C}^{n_{1}\times n_{2}}$.

Let $\mathcal{H}(\bm{x})=\bm{U}\bm{\Sigma}\bm{V}^{H}$ denote the compact singular value decomposition (SVD) of $\mathcal{H}(\bm{x})$ with $\bm{U}\in\mathbb{C}^{n_{1}\times r}$, $\bm{\Sigma}\in\mathbb{R}^{r\times r}$ and $\bm{V}\in\mathbb{C}^{r\times n_{2}}$. Let the subspace $\mathcal{\mathcal{T}}$ denote the support of $\mathcal{H}(\bm{x})$ and $\mathcal{\mathcal{T}}^{\perp}$ be its orthogonal complement. Let $\operatorname{sgn}(\widetilde{\bm{X}})=\widetilde{\bm{U}}\widetilde{\bm{V}}^{H}$ denote the sign matrix of $\widetilde{\bm{X}}$, where $\widetilde{\bm{X}}=\widetilde{\bm{U}}\widetilde{\bm{\Sigma}}\widetilde{\bm{V}}^{H}$ denotes the compact SVD of $\widetilde{\bm{X}}$.

In order to analyses the matrix completion problem theoretically, we need to introduce the standard incoherence condition as follows

We also need to introduce the following norms which, respectively, measure the largest spectral norm among matrices $\{\mathcal{A}_{k}(\bm{X})\}_{k=1}^{n}$ and the $\ell_{2}$ norm of $\{\left\lVert\mathcal{A}_{k}(\bm{X})\right\rVert\}_{k=1}^{n}$ [6, 7]

and

In this section, we start by giving the theoretical guarantees for the proposed method. Then we extend the analysis to noisy circumstance. Our main result shows that when the prior information is reliable, the proposed approach (3) can outperform previous approach (2) by a logarithmic factor.

By straightly following [35, Theorem 7], an extension to the noisy version with bounded noise can be shown as follows.

In this section, we extend the analysis from one-dimensional signal to multi-dimensional signal. Consider a $d$-way tensor $\mathcal{X}^{d}\in\mathbb{C}^{N_{1}\times\ldots\times N_{d}}$, each of whose entries can be denoted as

for $(k_{1},\ldots,k_{d})\in\{1,\ldots,N_{1}\}\times\ldots\times\{1,\ldots,N_{d}\}$. Denote $\bm{f}_{l}=(f_{1,l},\ldots,f_{d,l})\in[0,1)^{d}$ as the frequency vector for $l=1,\ldots,r$.

Similar to the one-dimensional case, we use multi-level Hankel operator to lift $\mathcal{X}^{d}$ to a low-rank matrix. See Eq. (12), where $\mathcal{X}^{d-1}(i)=\mathcal{X}^{d}(:,\ldots,:,i),\,i=0,\ldots,\max\{n_{d}-1,N_{d}-n_{d}\}$. When $d=1$, the above operator degrades to normal Hankel operator. According to [6], the rank of the lifted matrix satisfies $\mbox{Rank}(\mathcal{H}^{d}\left(\mathcal{X}^{d}\right))\leq r$ by using high-dimensional Vandermonde decomposition. Let $\Phi^{d}$ denote the prior information of $\mathcal{X}^{d}$. We can complete the $d$-way tensor by using the following nuclear norm minimization

where $\Omega=\{(k_{1},\ldots,k_{d})\in\{0,\ldots,N_{1}-1\}\times\ldots\times\{0,\ldots,N_{d}-1\}\}$ denotes the index set of known entries of $\mathcal{X}^{d}$ and $\left\langle\mathcal{X}^{d},\mathcal{Y}^{d}\right\rangle=(\mbox{vec}(\mathcal{Y}^{d}))^{H}\mbox{vec}(\mathcal{X}^{d})$ denotes the inner product of the $d$-way tensors.

Let $\bm{E}(k_{1},\ldots,k_{d})$ be the canonical basis in the domain $\mathbb{C}^{N_{1}\times\ldots\times N_{d}}$, and define the following orthonormal basis,

where

and

So each $d$-way tensor $\mathcal{X}^{d}$ can be rewritten as

It’s straightforward to extend the theoretical guarantee from one-dimensional case to multi-dimensional case.

Due to the high computational complexity of low rank methods, we decide to use the non-convex method to solve the problem. First we use matrix factorization to decompose $\mathcal{H}(\bm{z})$ to two low complexity matrices, i.e. $\mathcal{H}(\bm{z})=\bm{U}\bm{V}^{H}$ with $\bm{U}\in\mathbb{C}^{n_{1}\times r}$ and $\bm{V}\in\mathbb{C}^{n_{2}\times r}$. Then we use Alternating Direction Method of Multipliers (ADMM) [36] to solve the problem.

First of all, we denote the nuclear norm as follows [37, Lemma 8]

Incorporating (14) into the problem (3) yields

Then, we start ADMM by the following argumented Lagrange function

where $\mu>0$ is an absolute constant, and $\Pi(\bm{z})$ is an indicator function

Next, we decompose (15) into three subproblems to get $\bm{z}^{(n+1)}$, $\bm{U}^{(n+1)}$ and $\bm{V}^{(n+1)}$

and the Lagrangian update is

By simple calculations, we can obtain

where $\mathcal{P}_{\Omega^{c}}$ denotes the projection operator on $\Omega^{c}$ and $\mathcal{H}^{\dagger}(\cdot)$ denotes the Penrose-Moore pseudo-inverse mapping corresponding to $\mathcal{H}(\cdot)$.

Then, by taking the derivative of the other two problems and setting them to zero, we can otain

and

The last question is how to initialize $\bm{U}$ and $\bm{V}$. In order to converge quickly, the authors in [38] uses an algorithm named LMaFit [39], which is

However, LMaFit only uses the undersampled measurements and cannot guarantee that $\bm{Z}$ has the lifting matrix structure. Instead, we can take advantage of the reference image to initialize $\bm{U}$ and $\bm{V}$ by using truncated SVD when the reference image is reliable. Here, we use the value of $\left\lVert\mathcal{P}_{\Omega}(\bm{x}-\bm{\phi})\right\rVert$ as a criterion to choose the suitable initialization strategy.

Then we can give the corresponding algorithm in Algorithm 1. Here, $r\text{-SVD}(\cdot)$ returns the results of truncated SVD. And $\text{LMaFit}(\mathcal{H},\mathcal{P}_{\Omega}({\bm{x}}))$ denotes the algorithm in (16).

In this section, we carry on numerical simulations to show the improvement of the proposed method (3) compared to standard Hankel matrix completion (2). Besides, we compare the performance under two different solvers: CVX solver and ADMM-solver. Here, we use CVX package[40, 41] to get the convex results and use Algorithm 1 to get the ADMM results.

In this paper, we have integrated prior information to improve the performance of spectrally sparse signal recovery via structured matrix completion problem and have provided the related performance guarantees. Furthermore, we have designed corresponding ADMM algorithm to reduce the computational complexity. Both the theoretical and experimental results show that the proposed scheme outperforms standard Hankel matrix completion.