Projected Wirtinger Gradient Descent for Digital Waves Reconstruction

Signals contain valuable information on the attributes and actions of the emitter, and reflect the unique features of the relevant phenomena. In terms of economy, the periodic fluctuations of signals, which differ in magnitude, phase, and frequency, have attracted much attention from the academia. In many cases, however, it is difficult to extract all the information from signals, due to the limitations in sampling instruments or measuring methods. Take the extraction of frequency components for instance. Most scholars could only obtain a few discrete components from the target signal. Therefore, the signal recovery from multiple measured data is of great importance to signal processing in various fields, such as GNP cyclical fluctuations in economics [1, 2], in acceleration of medical imaging, analog-to-digital conversion, inverse scattering and in seismic imaging.

This work aims to recover a signal that linearly combines one-dimensional complex sine components with R different but continuous frequencies. To recover the entire signal, it is necessary to fully analyze the few samples of the sine components in the time domain, and measure their frequencies. This recovery task is similar to that in the signal processing of the following fields: variations in gross national product [1, 2], speeding up the capture of medical images [3], transform from analog signals to digital signals [4], and imaging of earthquake consequences [5]. The traditional strategies for the task include Prony’s approach [6], Estimation of Signal Parameters via Rotational Invariance Techniques [7], and matrix pencil [8]. In all these strategies, the sampling speed meets the Nyquist sampling theorem (NST).

As a novel technique for signal recovery, compressed sensing can recover the signal, which is sparsely or approximately sparsely distributed in a finite discrete domain, with a sample size smaller than the requirement of the NST [9] [10]. In the real-world, the signal values are usually distributed in a continuous manner. In our problem, the signal frequencies fall within the interval of [0, 1). To make compressed sensing feasible, the continuous values of the target signal could be divided into a limited number of equidistant points. Nevertheless, basis mismatch will occur if the division is not sufficiently refined [11]. In this case, the division error will be so large as to induce a huge error in signal recovery. To make matters worse, the ultrafine division of the continuous signal values will incur an unbearably high cost [12].

In recent years, more and more scholars have proposed strategies to recover continuous values based on a few discrete and unbalanced data in the time domain. For example, Candès et al. [13] minimized the total variation to pinpoint the continuous frequencies from equidistant data. Xu et al. [14] recovered the frequencies of the target signal from unbalanced data by minimizing the atomic norm. Candès et al. [13], Xu et al. [14] transformed the recovery of signal frequencies into the completion of low-rank Toeplitz matrix. Tang et al. [15] treated the signal frequency recovery of unbalanced data as a low-rank Hamiltonian Monte-Carlo problem. The above strategies can recover signals with a high robustness. But the computing efficiency is rather low in the completion of the low-rank matrices: the number of unknown factors in the optimization process is the square of the signal dimension.

Many efforts have been paid to complete the relevant matrices in an efficient manner. For instance, Candès and Recht [16] adopted interior point strategy to derive a Hessian matrix, whose scale is O(N⁴), in the Newtonian operation. In addition, many first-order approaches [17] were developed for matrix completion in signal recovery. Nonetheless, these approaches require an unstructured dual matrix, whose scale is O(N²). Overall, none of these approaches could efficiently recover large-dimensional signals.

To improve the efficiency of large-dimensional signal recovery, this work presents a feasible-point algorithm to complete the low-rank Hamiltonian Monte-Carlo matrix. The algorithm does not implement convex optimization of the matrix, but adopts a memory of the size O(NR), and recovers large-dimensional signal efficiently, for the number of sine components is way larger than N. Inspired by previous research [18, 19], the convergence of the algorithm to global optimal solution was analyzed in details. After that, the fast iterative shrinkage-thresholding (FIST) algorithm was referred to accelerate the convergence of our algorithm [1]. Finally, simulations confirm that the initial algorithm and the speed up strategy could recover large-dimensional signal, which linearly combines one-dimensional complex sine components with R different but continuous frequencies.

The rest of this work is arranged as follows: Chapter 2 introduces the details on the signal recovery problem; Chapter 3 presents the initial algorithm, analyzes its convergence performance, and speeds up its convergence with a self-designed algorithm; Chapter 4 verifies the feasibility of our algorithm through simulations; Chapter 5 puts forward the conclusions.

3.1 Benchmark

To solve the non-convex optimization problem, the collection of smaller-than-R-order matrices with complex values and the collection of Hankel matrices with complex values that meet the observable data can be respectively expressed as:

$\mathcal{R}_{\mathbb{C}}^{R}=\left\{L \in \mathbb{C}^{N \times N} \mid \operatorname{rank}(L) \leq R\right\}$.   (8)

$\mathcal{H}=\left\{\mathcal{H} x \mid x \in \mathbb{C}^{2 N-1}, x_{\Theta}=\widetilde{x_{\Theta}}\right\}$     (9)

Then, the signal recovery task and non-convex optimization can be described as:

$\min _{L \in \mathcal{R}_{\mathbb{C}, H \in \mathcal{H}}^{i n}} \frac{1}{2}\|L-H\|_{F}^{2}$     (10)

where, $\mathcal{R}_{\mathbb{C}}^{R}$ is a collection and a smooth manifold; $\mathcal{H}$ is an affine space. Formula (10) was solved by a feasible point algorithm. The objective function with the actual values of complex parameters can be expressed as: $F(L, H):=\frac{1}{2}\|L-H\|_{F}^{2}$. Complex calculus theory shows that the objective function cannot be differentiated. But this function can be differentiated into a real part $\Re$ of parameters, and an imaginary part $\Im$ of parameters. The two parts were subject to gradient flow, respectively. It is assumed that:

$Z=\left[\begin{array}{l}

\mathrm{L} \\

\mathrm{H}

\end{array}\right]=\Re+i \Im$     (11)

Then, the objective function can be rewritten as $F(\Re, \Im)$. Next, the real part was updated by $\frac{\partial F}{\partial \Re}$, while the imaginary part by $\frac{\partial F}{\partial \Im}$. That is, the initial algorithm updates Z by $\frac{\partial F}{\partial \Re}+i \frac{\partial F}{\partial \Im}$. Inspired by Fischer [27], there is:

$\frac{\partial F}{\partial \Re}+i \frac{\partial F}{\partial \Im}=2 \frac{\partial F}{\partial \bar{Z}}$   (12)

Thus:

$2 \frac{\partial F}{\partial \bar{Z}}=\left[\begin{array}{l}

2 \frac{\partial F}{\partial \bar{L}} \\

2 \frac{\partial F}{\partial \bar{H}}

\end{array}\right]=\left[\begin{array}{l}

L-H \\

H-L

\end{array}\right]$     (13)

Let $\delta_{1} \text { and } \delta_{2}$ are positive step lengths; $\mathcal{P}_{\mathcal{R}_{\mathbb{C}}^{R}} \text { and } \mathcal{P}_{\mathcal{H}}$ be the projections to $\mathcal{R}_{\mathbb{C}}^{R} \text { and } \mathcal{H}$, respectively. Then, the initial algorithm can be established as:

$\left\{\begin{array}{c}

\mathrm{L}_{t+1} \in \mathcal{P}_{\mathcal{R}_{\mathbb{C}}^{R}}\left(\mathrm{~L}_{t}-\delta_{1}\left(\mathrm{~L}_{t}-\mathrm{H}_{t}\right)\right), \\

\mathrm{H}_{t+1} \in \mathcal{P}_{\mathcal{H}}\left(\mathrm{H}_{t}-\delta_{2}\left(\mathrm{H}_{t}-\mathrm{L}_{t+1}\right)\right).

\end{array}\right.$     (14)

Then, the aim of the initial algorithm is to search for $\mathcal{P}_{\mathcal{R}_{\mathbb{C}}^{R}} \text { and } \mathcal{P}_{\mathcal{H}}$. Suppose the columns of U_R and V_R are the first R singular vectors of X from the left and the right, respectively; $\Sigma_{R}$ is the diagonal matrix whose elements correspond to these vectors. Drawing on Golub and Van Loan [28], $\mathcal{P}_{\mathcal{R}_{\mathbb{C}}^{R}}(\mathrm{X})$, as the optimal rank approximation to X, can be described as $\mathcal{P}_{\mathcal{R}_{\mathbb{C}}^{R}}(\mathrm{X})=\mathrm{U}_{R} \Sigma_{R} \mathrm{~V}_{R}^{*}$. Then, the following lemma was proposed to describe $\mathcal{P}_{\mathcal{H}}$ in its closed form Lemma 1:

$\mathcal{P}_{\mathcal{H}}(\mathrm{X})=\mathcal{H} \mathrm{z}$,

where $z_{j}=\left\{\begin{array}{cc}

\tilde{x}_{j}, & \text { if } j \in \Theta \\

\operatorname{mean}\left\{X_{k l} \mid k+l=j\right\}, & \text { otherwise }

\end{array}\right..$               (15)

Proof. $\mathcal{P}_{\mathcal{H}}(\mathrm{X})$ is the results of the LS (least squares) problem:

$\mathcal{P}_{\mathcal{H}}(Z)=\operatorname{argmin}_{Z}\left\{\|Z-X\|_{F}^{2}: W \in \mathcal{H}\right\}$ 

$=\mathcal{H} \cdot \operatorname{argmin}_{z}\left\{\|\mathcal{H} z-X\|_{F}^{2}: z_{\Theta}=\tilde{x}_{\Theta}\right\}$

$=\mathcal{H} \cdot \operatorname{argmin}_{z}\left\{\sum_{j=0}^{2 N-2}\left[\sum_{k+l=j}\left(z_{j}-X_{k l}\right)^{2}: z_{\Theta}=\tilde{x}_{\Theta}\right]\right\}$        (16)

The results in the last row are clearly the z_j in formula (15). Q.E.D. During the iteration of our feasible-point algorithm, L_t and H_t separately fall into the sets of feasible values $\mathcal{R}_{\mathbb{C}}^{R} \text { and } \mathcal{H}$, which minimizes the computing load and memory occupation, provided that R is less than N. Whereas $\mathrm{L}_{t} \in \mathcal{R}_{\mathbb{C}}^{R}$, the relevant results only occupy a memory of O(NR). Moreover, the parametric representation of H_t only occupies a space of O(N). In addition, the initial algorithm only has to calculate the first R singular values and the relevant vectors during the projection computation. The computation can be completed without calling anything other than the matrix vector product (MVP) of $L_{t}-\delta_{1}\left(L_{t}-H_{t}\right)$. The MVP of L_t, which is an R-rank factorized value, can be computed in only O(NR) steps. Through fast Fourier transform [28], the MVP of H_t can be computed in only O(NlogN) steps. The second step of the initial algorithm only need to average L_t+1 along anti-diagonals [29-34].

3.2 Convergence analysis

Based on the data of Attouch et al. [18], the convergence of the initial algorithm was analyzed. Let us denote proper and lower semi-continuous (PLSC) functions as$\theta: \mathbb{R}^{n} \mapsto \mathbb{R} \cup\{+\infty\} \text { and } \omega: \mathbb{R}^{m} \mapsto \mathbb{R} \cup\{+\infty\}$, and C¹ function as $\phi: \mathbb{R}^{n} \times \mathbb{R}^{m} \mapsto \mathbb{R}$. In general, the non-convex optimization can be described as:

$\min _{x, y} \psi(x, y):=\phi(x, y)+\theta(x)+\omega(y)$     (17)

Drawing on Attouch et al. [18], formula (17) was solved through proximal alternating minimization:

$\left\{\begin{array}{c}

\mathrm{x}_{k+1} \in \operatorname{argmin}_{\mathrm{x} \in \mathbb{B}^{n}} \psi\left(\mathrm{x}, \mathrm{y}_{k}\right)+\frac{1}{2 \lambda_{k}}\left\|\mathrm{x}-\mathrm{x}_{k}\right\|_{2}^{2}, \\

\mathrm{y}_{k+1} \in \operatorname{argmin}_{\mathrm{y} \in \mathbb{R}^{m}} \psi\left(\mathrm{x}_{k+1}, \mathrm{y}\right)+\frac{1}{2 \mu_{k}}\left\|\mathrm{y}-\mathrm{y}_{k}\right\|_{2}^{2}.

\end{array}\right.$      (18)

Suppose function $\psi$ meets the Kurdyka-Lojasiewicz condition, and that $\nabla \phi$ is Lipschitz on bounded sets. Then, the convergence of formula (17) could be proved by the method specified by Attouch et al. [18]. Since it is difficult to verify, the Kurdyka-Lojasiewicz condition could be guaranteed by the semi-algebraic SA property. A PLSC function has SA property, as long as its graph forms an SA set. The condition for a subset $S \subset \mathbb{R}^{d}$ to be a real SA set was defined as follows: there is a limited number of real polynomial functions $g_{i j}, h_{i j}: \mathbb{R}^{d} \mapsto \mathbb{R}$ so that:

S=$\bigcup_{j=1}^{p} \bigcap_{i=1}^{q}$$\left\{\mathrm{u} \in \mathbb{R}^{d} \mid g_{i j}(\mathrm{u})=0, h_{i j}(\mathrm{u})<0\right\}$     (19)

Then, the sets $C \in \mathbb{R}^{n} \text { and } D \in \mathbb{R}^{m}$ were respectively indicated by functions $\theta=\delta_{C} \text { and } \omega=\delta_{D}$. Then, set C could be indicated by $\delta_{C}(\mathrm{x})=\left\{\begin{array}{cl}

0, & \text { if } \mathrm{x} \in C, \\

+\infty, & \text { if } \mathrm{x} \notin C .

\end{array}\right.$ Under the condition that $\phi(x, y)=\frac{1}{2}\|x-y\|_{2}^{2}$, the non-convex optimization can be converted into an alternating projection:

$\left\{\begin{array}{c}

\mathrm{x}_{k+1} \in \mathcal{P}_{C}\left(\mathrm{x}_{k}-\frac{1}{1+\delta_{k}}\left(\mathrm{x}_{k}-\mathrm{y}_{k}\right)\right), \\

\mathrm{y}_{k+1} \in \mathcal{P}_{D}\left(\mathrm{y}_{k}-\frac{1}{1+\mu_{k}}\left(\mathrm{y}_{k}-\mathrm{x}_{k+1}\right)\right).

\end{array}\right.$     (20)

According to Bolte et al. [7] and Attouch et al. [18], the convergence of formula (18) can be described by the following theorem.

Theorem 1. It is assumed that $C \subset \mathbb{R}^{n} \text { and } D \subset \mathbb{R}^{m}$ have SA property. Then, the alternating projection can produce (x_k, y_k), where 0<a<δ_k, and μ_k<b for all k:

Case 1: $\left\|\left(\mathrm{x}_{k}, \mathrm{y}_{k}\right)\right\|_{2} \rightarrow \infty \text { as } k \rightarrow \infty$, or (x_k, y_k) converges to a critical point of $\psi$;

Case 2: If (x₀, y₀) is feasible and sufficiently close to a global minimizer of $\psi$, (x₀, y₀) converges to that minimizer.

Further, the above theorem was introduced to the initial algorithm to test its convergence performance. Although the initial algorithm takes the same shape as alternating projection, the above theorem targets real numbers, while the initial algorithm handles the matrix space with complex values. Fortunately, any matrix with complex values could be converted into a matrix with real number by merging the real part with imaginary part. As stated before, the initial algorithm was constructed based on the gradient of the two parts. Because the objective function remains the same through the identification, it is only necessary to judge if $\mathcal{R}_{\mathbb{C}}^{R} \text { and } \mathcal{H}$ have SA property, when they are treated as collections of the two parts. The judgment was completed by the following lemmas:

Lemma 2. The following set has SA property:

$\mathcal{S}_{R}=\left\{[X, Y] \mid(X+l Y) \in \mathcal{R}_{C}^{R}\right\}$    (21)

Proof. Suppose

$\mathcal{P}_{r}=\left\{[X, Y] \mid X, Y \in \mathbb{R}^{N \times N}, \operatorname{rank}(X+\imath Y)=r\right\}$     (22)

and

$Q_{r}=\left\{[X, Y] \mid X, Y \in \mathbb{R}^{N \times N}, \operatorname{rank}\left(\left[\begin{array}{cc}

X & -Y \\

Y & X

\end{array}\right]\right)=2 r\right\}$    (23)

The first step is to demonstrate that $\mathcal{P}_{\mathrm{r}}=Q_{\mathrm{r}}$ by proving $\mathcal{P}_{\mathrm{r}} \subset Q_{\mathrm{r}} \text { and } Q_{\mathrm{r}} \subset \mathcal{P}_{\mathrm{r}} \mathrm{y}$. Suppose $[X, Y] \in \mathcal{P}_{r}$, and a singular value decomposition (SVD) of X+iY is $(X+i Y)=\left(U_{R e}+i U_{I m}\right) \Sigma\left(V_{R e}+i V_{I m}\right)^{*}$, where $\mathrm{U}_{\mathrm{Re}}, \mathrm{U}_{\mathrm{Im}}, \mathrm{V}_{\mathrm{Re}}, \text { and } \mathrm{V}_{\mathrm{Im}} \in \mathbb{R}^{\mathrm{N} \times \mathrm{r}} \text { and } \Sigma \in \mathbb{R}^{\mathrm{r} \times \mathrm{r}}$. Then, an SVD of $\left[\begin{array}{cc}

X & -Y \\

Y & X

\end{array}\right]$ can be expressed as:

$\left[\begin{array}{cc}

X & -Y \\

Y & X

\end{array}\right]=\left[\begin{array}{cc}

U_{R e} & -U_{I m} \\

U_{I m} & U_{R e}

\end{array}\right]\left[\begin{array}{cc}

\Sigma & 0 \\

0 & \Sigma

\end{array}\right]$$\left[\begin{array}{cc}

V_{R e} & -V_{I m} \\

V_{I m} & V_{R e}

\end{array}\right]^{*}$    (24)

Thus, $\operatorname{rank}\left(\left[\begin{array}{cc}

X & -Y \\

Y & X

\end{array}\right]\right)=2 r$ indicates that $[X, Y] \in \mathcal{Q}_{r}, \text { and that } \mathcal{P}_{\mathrm{r}} \subset \mathcal{Q}_{\mathrm{r}}$.

Next, it is assumed that $[X, Y] \in \mathcal{Q}_{r}$ Suppose $\left(\sigma,\left[\begin{array}{l}

u_{1} \\

u_{2}

\end{array}\right],\left[\begin{array}{l}

v_{1} \\

v_{2}

\end{array}\right]\right)$ is a singular triplet of $\left[\begin{array}{cc}

X & -Y \\

Y & X

\end{array}\right]$, then it is obvious that $\left(\sigma,\left[\begin{array}{c}

-u_{2} \\

u_{1}

\end{array}\right],\left[\begin{array}{c}

-v_{2} \\

v_{1}

\end{array}\right]\right)$. As a result, every singular value has an even multiplicity, and the SVD of $\left[\begin{array}{cc}

X & -Y \\

Y & X

\end{array}\right]$ must take the shape of formula (24). Thus, $(X+i Y)=\left(U_{R e}+i U_{I m}\right) \Sigma\left(V_{R e}+i V_{I m}\right)^{*}$ is an SVD of X+iY, indicating that $\operatorname{rank}(X+i Y)=r$ Hence, $[X, Y] \in \mathcal{P}_{r}, \text { and thus } \mathcal{Q}_{\mathrm{r}} \subset \mathcal{P}_{\mathrm{r}}$.

Because $\mathcal{Q}_{\mathrm{r}}$ is the overlap between the collection of rank-2r matrices with real numbers, and the linear subspace of matrices, whose form is $\left[\begin{array}{cc}

X & -Y \\

Y & X

\end{array}\right]$, it must have the SA property, according to the example in Bolte’s research [7]. Whereas $\mathcal{P}_{\mathrm{r}}=\mathcal{Q}_{\mathrm{r}}$, it can be seen that $\mathcal{P}_{\mathrm{r}}$ also has SA property.

Since $\mathcal{S}_{\mathrm{R}}=\bigcup_{\mathrm{r}=0}^{\mathrm{R}} \mathcal{P}_{\mathrm{r}}, \mathcal{S}_{\mathrm{R}}$ is a collection with SA property. Q.E.D.

Lemma 3. The following set has SA property:

$\mathcal{K}=\{[X, Y] \mid(X+i Y) \in \mathcal{H}\}$    (25)

Proof. For the linear operator $\mathcal{H}$, we have:

$\mathcal{H} x=\mathcal{H} \Re(x)+i \mathcal{H} \Im(x)$    (26)

If x meets $x_{\Theta}$$=\widetilde{x_{\Theta}}$, then $\mathfrak{R}\left(x_{\Theta}\right)=\mathfrak{R}\left(\widetilde{x_{\Theta}}\right)$ and $\widetilde{\mathfrak{J}}\left(x_{\Theta}\right)=\mathfrak{I}\left(\widetilde{x_{\Theta}}\right)$. Thus, we have:

$\mathcal{H}=\mathfrak{R}(\mathcal{H})+\mathrm{l} \mathfrak{S}(\mathcal{H})=\mathcal{K}_{1}+1 \mathcal{K}_{2}$    (27)

where, $\mathcal{K}_{1}=\left\{\mathcal{H} \mathrm{r} \mid \mathrm{r} \in \mathbb{R}^{2 \mathrm{~N}-2}, \mathrm{r}_{\Theta}=\Re\left(\widetilde{\mathrm{x}_{\Theta}}\right)\right\}$, $\mathcal{K}_{2}=\left\{\mathcal{H} i \mid \mathrm{i} \in \mathbb{R}^{2 \mathrm{~N}-2}, \mathrm{i}_{\Theta}=\mathfrak{I}\left(\widetilde{\mathrm{x}_{\Theta}}\right)\right\}$.

The above results indicate that $\mathcal{K}=\mathcal{K}_{1} \times \mathcal{K}_{2}$. Whereas $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$ are affine spaces, the product $\mathcal{K}$ of these spaces must be an affine space, which has SA property. Based on the above theorem and lemmas, the convergence of the initial algorithm can be described as the following theorem.

Theorem 2. Suppose (L_t, H_t) is produced by the initial algorithm under 0<δ₁, and δ₂<1.

Case 1: $\left\|\left(\mathrm{L}_{t}, \mathrm{H}_{t}\right)\right\|_{F} \rightarrow \infty$ as $t \rightarrow \infty$, or (L_t, H_t) converges.

Case 2: Suppose (L₀, H₀) is feasible and sufficiently close to a global minimizer of $\left\{\frac{1}{2}\|\mathrm{~L}-\mathrm{H}\|_{F}^{2} \mid \mathrm{L} \in \mathcal{R}_{\mathbb{C}}^{R}, \mathrm{H} \in \mathcal{H}\right\}$, then (L₀, H₀) converges to that minimizer.

The unboundedness in Case 1 can be solved easily by adding a constraint to $\mathcal{H}$.

3.3 Speeding up convergence

The convergence of the initial algorithm was sped up referring to the FIST algorithm [1]. The FIST algorithm can efficiently minimize the sum of two convex functions, in which one of them has a Lipschitz continuous gradient, by linearly combine the iterative results of the two functions. This linear combination strategy was adopted for the initial algorithm, despite that our problem is non-convex optimization. The speed-up strategy is described as follows:

For k₀=1, $\left(\mathrm{L}_{t}, \mathrm{H}_{k \backslash t}\right)$ can be produced by

$\left\{\begin{array}{c}\mathrm{L}_{t+1} \in \mathcal{P}_{\mathcal{R}_{\mathbb{C}}^{R}}\left(\mathrm{~L}_{t}-\delta_{1}\left(\mathrm{~L}_{t}-\widetilde{\mathrm{H}_{\mathrm{t}}}\right)\right), \\ \mathrm{H}_{t+1} \in \mathcal{P}_{\mathcal{H}}\left(\mathrm{H}_{t}-\delta_{2}\left(\widetilde{\mathrm{H}_{\mathrm{t}}}-\mathrm{L}_{t+1}\right)\right), \\ k_{t+1}=\frac{\sqrt{1+4 k_{t}^{2}}+1}{2}, \\ \widetilde{\mathrm{H}_{\mathrm{t}+1}}=\mathrm{H}_{t+1}+\frac{k_{t}-1}{k_{t+1}}\left(\mathrm{H}_{t+1}-\mathrm{H}_{t}\right).\end{array}\right.$    (28)

where, H is an affine subspace, the linear combination does not affect the feasibility of $\widetilde{\mathrm{H}_{\mathrm{t}+1}}$. Hence, the first and second steps of the speed-up strategy have the same computing load and memory occupation as those in the proposed algorithm. In addition, the computing load of the third and fourth steps are so small as to be negligible. Therefore, the speed-up strategy is as simple as the initial algorithm in terms of computation.

This paper presents a feasible-point algorithm for the recovery of signals with sparse spectrum, and continuous frequencies in the interval [0, 1], using only a few data in the time domain. Unlike the conventional recovery methods, the initial algorithm takes root in the non-convex optimization, and handles large-dimensional signals effectively. The convergence of the algorithm was verified through repeated simulations.

Moreover, the convergence of the initial algorithm was accelerated by referring to the FIST algorithm. Simulations show that the initial algorithm and the speed up strategy could recover signals on medium and large scales successfully from a limited amount of data in the time domain, and outperform contrastive methods like ANM and EMaC.

The future research will focus on the following aspects: identifying the theoretical limit on the sample size for signal recovery by our algorithms; applying our methods to the recovery of damped signals, e.g., the fluctuation GNP data; solving the problems in the programming of the Hankel matrix.