A Fast and Accurate Sparse Continuous Signal Reconstruction by Homotopy DCD with Non-Convex Regularization

Abstract

In recent years, various applications regarding sparse continuous signal recovery such as source localization, radar imaging, communication channel estimation, etc., have been addressed from the perspective of compressive sensing (CS) theory. However, there are two major defects that need to be tackled when considering any practical utilization. The first issue is off-grid problem caused by the basis mismatch between arbitrary located unknowns and the pre-specified dictionary, which would make conventional CS reconstruction methods degrade considerably. The second important issue is the urgent demand for low-complexity algorithms, especially when faced with the requirement of real-time implementation. In this paper, to deal with these two problems, we have presented three fast and accurate sparse reconstruction algorithms, termed as HR-DCD, Hlog-DCD and Hl_p-DCD, which are based on homotopy, dichotomous coordinate descent (DCD) iterations and non-convex regularizations, by combining with the grid refinement technique. Experimental results are provided to demonstrate the effectiveness of the proposed algorithms and related analysis.

1. Introduction

Compressive sensing (CS) and sparse signal representation have received widespread attention and increasing interest over the past few years [1,2], which are motivated by the sparse nature of the real world data and the advantages of the CS theory. The applications of CS in numerous areas have been widely investigated in the literature, such as magnetic resonance imaging (MRI) [3], synthetic aperture radar (SAR) imaging [4], inverse synthetic aperture radar (ISAR) imaging [5], passive radar imaging [6], direction-of-arrival (DOA) estimation [7], communication channel estimation [8], seismic signal processing [9], spectral estimation [10], image processing [11], speech enhancement [12], etc. Generally speaking, those disciplines explore the following linear signal model [13]:

$$\mathbf{y}=\mathbf{H}\mathbf{w}+\mathbf{n}$$

(1)

where y ∈ ℂ^M×1 is a M × 1complex vector, H ∈ ℂ^M×N is a M × N complex measurement matrix, w ∈ ℂ^N×1 represents the N × 1 unknown complex signals of interest, and n ∈ ℂ^M×1 denotes a M × 1complex noise vector. Suppose w is sparse or compressible in a known dictionary Ψ ∈ ℂ^N×N, i.e., w = Ψx, where x ∈ ℂ^N×1 is a K-sparse vector, namely it can be approximated by its K largest coefficients or its coefficients satisfy a power decay law with K strongest coefficients. Therefore, linear measurements are obtained in CS as:

$$\mathbf{y}=\mathbf{A}\mathbf{x}+\mathbf{n}$$

(2)

where A = HΨ ∈ ℂ^M×N represents the sensing matrix. Although the above equation is usually ill-posed, the CS theory has shown that if A satisfies some certain conditions, we can construct x and w stably from highly undersampled measurements y [14].

1.1. Off-Grid Problem in CS-Based Methods

In the CS processing procedure, the first necessary step is to design a dictionary through a discretization operation with the assumption that the elements of unknown x lie exactly on those pre-defined grids corresponding to A. Obviously, this is practically impossible since the continuous nature of the unknowns of x, such as the unknown directions, may not fall into the predefined angular grids in DOA, and the scatterers of the target may not locate exactly on the pre-discretized scene grids in radar imaging. Hence, once faced with a continuous signal, the off-grid problem in using the conventional sparse recovery techniques is inevitable, no matter how densely we grid x. Previous researches [15–17] have demonstrated that the traditional CS-based methods would be severely affected when the off-grid problem emerges, and [18] also claimed that the off-grid problem is one of the major constraints in popularizing CS-based methods in practice. It is worth noting that there may be other factors leading to basis mismatch, for example, in the radar imaging field, unsatisfactory system parameters (e.g., position errors of transceivers [19] and phase synchronization mismatch [20]) are also likely to degrade the performance of conventional CS-based methods from our previous researches, however this paper only focuses on the off-grid problem, and assumes the system errors are small enough.

The solutions for off-grid problem have been broadly studied in previous literatures [21–30]. So far, main research topics include three types as follows:

(a)

Direct method according to theories of Xampling and the finite rate of innovation (FRI). This scheme is first introduced by Michaeli et al. as a general framework within various solutions for analog signals [21]. The significant advantage of this method is its direct solution, without any pre-discretization of the continuous signal at first. Thus it is less sensitive to off-grid problem compared with traditional CS-based methods. However, this method is mainly based on spectral estimation algorithms (e.g., ESPIRT [22], MUSIC [6], Matrix Pencil [23]), which needs special signal expression in linear measurement model and its performance may suffer from low SNR and few snapshots [24].

(b)

Off-grid CS method. This kind of method is under a first order Taylor expansion model, and thus is quite sensitively depended on the model's accuracy, and cannot give a thorough solution when higher order approximations are significant [23]. Furthermore, it is likely to involve highly-computational burden [25,26] for alternatively finding the sparse solution and estimating the off-grid error, especially in large scale problem.

(c)

Grid refinement approach. The idea of the grid refinement technique was firstly introduced by Malioutov et al. [27] to mitigate the effect of limiting estimates to a grid of spatial locations in source localization problem. Then it is generalized by Liu et al. [28] for locating two-dimensional multiple underwater acoustic sources. Recently, Guldogan et al. [29] have proposed a novel grid refinement algorithm to alleviate off-grid problem by using a particle swarm optimization (PSO) and orthogonal matching pursuit (OMP), which makes use of the PSO to perturb the location of each grid point.

Since the “grid refinement approach” is iteratively refined to match with the desired resolution of the off-grid components by “coarse and fine grid partition” operations, and does not have the drawbacks of (a) and (b), this paper follows the idea of “grid refinement approach” which makes a coarse grid first instead of having an universally fine grid to reduce the complexity, then achieves fine grids around the peaks using more refinement levels. In this way, after a few iterations of refining process, it becomes fine enough that off-grid problem effect is negligible.

1.2. Fast and Efficient Algorithms for Real-Time Implementation

Furthermore, there exists another common challenge by utilizing the CS-based methods, i.e., we need fast and efficient algorithms for real-time system implementation, particularly for digital electronic circuits (e.g., ARM, FPGA, DSP) [30]. From previous researches, sparse recovery techniques can be roughly divided into two families, that are greedy methods (e.g., MP [31], OMP [32], GP [33]) and optimization based methods (e.g., l₁ optimization [34], smoothed l₀ optimization [35], non-convex optimization [36]). Generally speaking, greedy methods have the advantages of lower complexity, faster speed, less storage requirement, and flexible implementation compared with optimization based methods, and are considered the most suitable candidates for hardware implementation. However, their performances are inferior to those of the optimization based methods, such as l₁ norm minimization basis pursuit de-noising (BPDN) [34].

As the coordinate descent (CD) search has an inherent property of being low complexity when signals are sparse, Zakharov et al. have successfully exploited dichotomous CD (DCD) iterations for solving LS [37], RLS [38] and MVDR beamforming problems [39] using FPGAs for real-time implementation. In addition, to solve the reweighted l₁ minimization problem, they have developed a greedy algorithm called Hl1-DCD [40], i.e., homotopy DCD method with reweighted l₁ penalty, which has shown a high recovery performance with relatively low complexity. It turns out that the overall complexity of the algorithm is comparable to that of MP. Moreover, homotopy iterations result in high accuracy of the solution, even higher than that of the YALLl algorithm. Meanwhile, Liu and Zakharov et al. [28] have utilized the low complexity homotopy approach combined with CD search for locating underwater acoustic sources by solving the multi-frequency BPDN problem. The proposed method is evaluated by applying to simulated and real experimental data. As far as we know, previous studies [28,37–40] are mainly focusing on homotopy DCD with convex regularizations, while there are seldom researches linking to homotopy DCD with non-convex regularizations.

1.3. Our Contribution

From the above discussion, we have addressed two major issues (i.e., off-grid problem and efficient algorithms for real-time implementation) of applying CS to sparse continuous signal reconstruction. In this paper, motivated by the idea of Hl1-DCD [40] and grid refinement approach [27], we present three fast and accurate sparse reconstruction algorithms (i.e., HR-DCD, Hlog-DCD and Hl_p-DCD) which are based on homotopy, dichotomous coordinate descent (DCD) iterations and non-convex regularization, combining with the grid refinement technique to deal with the aforementioned issues. The main contributions of this paper are as follows: (1) we formulate the sparse recovery problem by homotopy DCD method with three typical classes of non-convex penalties, which are proved to recover sparsity in a more efficient way than homotopy DCD method with convex penalties as shown in Zakharov's previous researches [28,37–40]; (2) further, the grid refinement technique [27] is utilized to combine with our algorithms to alleviate the effect of off-grid problem and reduce the complexity simultaneously; (3) experimental results of three representative applications (DOA, passive radar imaging, ISAR imaging) are carried out to verify the effectiveness of the proposed methods.

The outline of this paper is as follows: in Section 2, two sparse recovery problems are discussed. Section 3 describes the proposed methods and related analysis. In Section 4, extensive experimental results are presented to verify our methods. Finally, Section 5 draws the conclusions.

Throughout this paper we shall make use of the following notation: bold-case letters are reserved for vectors and matrices, e.g., x is a vector, A is a matrix; Elements of A and x are represented as A_n,p and x_n, respectively. I and I^c denote the support and its complement, respectively. Further, we represent by R^(q) the q-th column of R; A_I a matrix obtained from A keeping only columns corresponding to I; x_I the subset of x that contains entries from x corresponding to I; ‖x‖_p denotes the l_p norm of a vector, ‖A‖₂ is the spectral norm of the matrix A; 〈,〉 denotes the inner product; (·)^H, (·)*, Re{·}denote the conjugate transpose, conjugate and real part operations, respectively.

2. Problem Formulation and Motivations

This section briefly introduces the general sparse recovery framework, upon which we develop our algorithms. As stated before, various applications can be represented by the linear model as Equation (2) shows.

When considering the practical applications, sparse recovery techniques require low-complexity algorithms just as some kinds of greedy methods which are suitable for real-time implementation. However, the performances of greedy methods are compared unfavorably with optimization based methods (e.g., using l₁ norm, reweighted l₁ norm, l_p norm, smoothed l₀ regularization). Recently, Zakharov et al. have proposed the homotopy DCD method with convex regularizations, including l₁ and reweighted l₁ norms, which have shown a high recovery performance and relatively low complexity. Moreover, most of the operations in the algorithms are additions, therefore, they are very suited for the hardware implementation of real-time operating systems. Motivated by the idea of homotopy DCD method and the fact that non-convex regularization usually yields a sparser solution than any convex penalty for a given residual energy (see Figure 1 for example), this paper proposes a fast and accurate sparse signal reconstruction by homotopy DCD technique with non-convex regularizations. Three typical non-convex penalty functions [12,41] are considered, i.e., the first order rational function penalty, the logarithmic penalty and the l_p penalty, respectively. With the corresponding penalties, we have achieved the derivation of the HR-DCD, Hlog-DCD and Hl_p-DCD algorithms in the next section.

Since x is usually distributed continuously in the corresponding space in many applications, the off-grid problem is likely to exist. For example, in the radar imaging field, the reflectivity centers of the scatterers are generally not located at exact on-grid spatial positions illustrated in Figure 2, which means the measurement matrix should cover the basis vectors corresponding to off-grid scatterers. As conventional CS-based methods do not consider the off-grid effect, therefore the related sparse reconstruction techniques suffer unacceptable degradation in image quality. This paper has utilized the grid refinement approach [27] shown in Figure 3 combined with the HR-DCD, Hlog-DCD and Hl_p-DCD algorithms to alleviate the impact of off-grid problem (see details in Section 3.2). Experimental results in Section 4 are provided to demonstrate the performance improvement of the proposed algorithms.

3. Homotopy DCD with Non-Convex Regularization Algorithms

We consider minimization of the following cost function to solve the problem Equation (2):

$$J\left[\mathbf{x}\right]=\frac{1}{2}{\Vert \mathbf{y}-\mathbf{A}\mathbf{x}\Vert }_2^2+\tau \rho \left(\mathbf{x}\right)$$

(1)

where the sparsity measurement function ρ(x) is generally separable, or ρ(x)= Σ_i φ (x_i) where φ(·) denotes the parameterized form and satisfies such conditions that ρ(x) is sparsity-preserving, and τ ≥ 0 is a regularization parameter. As mentioned in previous researches, we can use DCD iterations for minimizing the cost function, which is represented in a fixed point. The DCD algorithm is appealing for practical designs as it operates at the bit level, resulting in stable hardware implementations, which is shown in

Table 1. DCD Algorithm.

Step	Equation
	Initialization: x = 0, r = y, R = AHA, m = 0, δ = H
1	for m = 1 to Mb, repeat:
2	δ = δ/2, α = [δ,− δ, jδ, − jδ]
3	for i = 1 to N, repeat:
4	for k = 1 to 4, repeat:
5	ΔJ = δ2Ri,i/2 − Re{ ${\alpha }_k^{*}$ (A(i))H r}+τ(φ(xi+αk) − φ(xi))
6	if ΔJ < 0, do:
7	xi ← xi + αk, r =r − αkA(i)
8	end if
9	end for
10	end for
11	end for

. Different from Zakharov's studies [28,37–40], this paper utilizes the non-convex regularizations instead of convex regularizations embedding with the DCD algorithm, and chooses three typical penalty functions as follows:

(a)

The first order rational function penalty:

$$\varphi \left(x_i\right)=\frac{\left|x_i\right|}{1+\frac{a}{2}\left|x_i\right|}$$

(2)

(b)

The logarithmic penalty:

$$\varphi \left(x_i\right)=\frac{1}{a}log\left(1+a\left|x_i\right|\right)$$

(3)

(c)

The l_p penalty:

$$\varphi \left(x_i\right)={\left|x_i\right|}^p$$

(4)

In order to reduce the complexity of the DCD algorithm, we develop a greedy algorithm that is based on homotopy method with respect to a set of the parameter τ. Besides, the updates of our DCD algorithm are only done within the support instead of all elements. As the support is usually much smaller than N, therefore, the complexity is further reduced. We consider the first order rational function penalty at first, and develop HR-DCD algorithm. Then Hlog-DCD algorithm and Hl_p-DCD algorithm can be similarly obtained.

The following two propositions define the rules for adding/removing elements into/from the support of HR-DCD algorithm.

Let I be the support at some homotopy iteration, and I^c be its complementary set. Denote r = y − Ax, R = A^H A and b = A^H r.

Proposition 1: Add the t-th (t ∈ I^c) element into the support I according to the rule:

$$\begin{array}{cc}t=arg\underset{k\in I^c}{max}\frac{a\left|b_k\right|-R_{k,k}}{a{R}_{k,k}^2}\left\{{\left(R_{k,k}+a\left|b_k\right|\right)}^2-4a{R}_{k,k}\tau \right\}& s.t.{\left(R_{t,t}+a\left|b_t\right|\right)}^2>4a{R}_{t,t}\end{array}\tau$$

(5)

Proposition 2: Remove the t-th (t ∈ I) element from the support I according to the rule

$$\begin{array}{cc}t=arg\underset{k\in I}{min}\frac{1}{2}{\left|x_k\right|}^2{R}_{k,k}+\Re \left\{x_k^{*}{b}_k\right\}-\tau \frac{\left|x_k\right|}{1+\frac{a}{2}\left|x_k\right|}& s.t.\frac{1}{2}{\left|x_t\right|}^2{R}_{t,t}+\Re \left\{x_t^{*}{b}_t\right\}-\tau \frac{\left|x_t\right|}{1+\frac{a}{2}\left|x_t\right|}\end{array}<0$$

(6)

We now prove Proposition 1.

Proof: Let the t-th element x_t = 0 be activated as x̂ = α = |α|e ^jarg{α}, and the updated solution vector is denoted as x̂ The update of the cost function Equation (3) is then given by:

$$\mathrm{\Delta }J=J\left[\widehat{\mathbf{x}}\right]-J\left[\mathbf{x}\right]=\frac{1}{2}{\left|\alpha \right|}^2{R}_{t,t}-\Re \left\{{\alpha }^{*}{b}_t\right\}+\tau \frac{\left|\alpha \right|}{1+\frac{a}{2}\left|\alpha \right|}$$

(7)

The cost function is reduced if ΔJ < 0. For a fixed |α|, ΔJ achieves a minimum value if arg{α} = arg{b_t}, and in this case:

$$\mathrm{\Delta }J=\frac{1}{2}{\left|\alpha \right|}^2{R}_{t,t}-\left|\alpha \right|\left|b_t\right|+\tau \frac{\left|\alpha \right|}{1+\frac{a}{2}\left|\alpha \right|}$$

(8)

Let:

$$\begin{array}{ll}\mathrm{\Delta }J\prime \hfill & =1\left(1+\frac{a}{2}\left|\alpha \right|\right)\mathrm{\Delta }J\hfill \\ & =\left|\alpha \right|\left\{a{R}_{t,t}{\left|\alpha \right|}^2\left(R_{t,t}-a\left|b_t\right|\right)\left|\alpha \right|+4\left(\tau -\left|b_t\right|\right)\right\}\hfill \\ & =a{R}_{t,t}\left|\alpha \right|\left\{{\left(\left|\alpha \right|+\frac{R_{t,t}-a\left|b_t\right|}{a{R}_{t,t}}\right)}^2-{\left(\frac{R_{t,t}-a\left|b_t\right|}{a{R}_{t,t}}\right)}^2+\frac{1\left(\tau -\left|b_t\right|\right)}{a{R}_{t,t}}\right\}\hfill \end{array}$$

(9)

For |α| > 0, if $\left|\alpha \right|=-\frac{R_{t,t}-a\left|b_t\right|}{a{R}_{t,t}}$, we have:

$$\mathrm{\Delta }J\prime =\frac{1}{a{R}_{t,t}}\left|\alpha \right|\left\{4a{R}_{t,t}\tau -{\left(R_{t,t}+a\left|b_t\right|\right)}^2\right\}$$

(10)

If ΔJ′ < 0, then 4aR_t,t τ − (R_t,t + a|b_t|)² < 0. Thus, there exists a value of α that results in ΔJ < 0, i.e., in reducing the cost function. It is seen from the last expression that if we want to add a new element to the solution vector, the index t of the element should correspond to the maximum of $\frac{a\left|b_t\right|-R_{t,t}}{a{R}_{t,t}^2}\left\{{\left(R_{t,t}+a\left|b_t\right|\right)}^2-4a{R}_{t,t}\tau \right\}$ over t ∈ I^c. In this case, we will obtain the largest decrement of the cost function.

Thus, the value of α that for a fixed t such that 4aR_t,t τ − (R_t,t + a|b_t|)² < 0 results in the decrement of the cost function is given by:

$$\alpha =\frac{a\left|b_t\right|-R_{t,t}}{a{R}_{t,t}}{e}^{jarg\left\{b_t\right\}}$$

(11)

According to the above statements, the (n + 1)-th nonzero element to be activated in x̂ should satisfy:

$$\begin{array}{cc}t=arg\underset{k\in I^c}{max}\frac{a\left|b_k\right|-R_{k,k}}{a{R}_{k,k}^2}\left\{{\left(R_{k,k}+a\left|b_k\right|\right)}^2-4a{R}_{k,k}\tau \right\}& s.t.{\left(R_{t,t}+a\left|b_t\right|\right)}^2>4a{R}_{t,t}\tau \end{array}$$

(12)

For |α| > 0, we have:

$$\begin{array}{ll}\mathrm{\Delta }J\hfill & \ge \frac{1}{2}{\left|\alpha \right|}^2{R}_{t,t}-\left|\alpha \right|\left|b_t\right|+\tau \frac{\left|\alpha \right|}{1+\frac{a}{2}\left|\alpha \right|}\hfill \\ & \ge \frac{\left|\alpha \right|}{4\left(1+\frac{a}{2}\left|\alpha \right|\right)}\left\{\left(\left|\alpha \right|R_{t,t}-2\left|b_t\right|\right)\left(2+a\left|\alpha \right|\right)+4\tau \right\}\hfill \end{array}$$

(13)

If (R_t,t + a|b_t|)² < 4aR_t,t τ, then in this case:

$$\begin{array}{ll}\mathrm{\Delta }J\hfill & >\frac{\left|\alpha \right|}{4\left(1+\frac{a}{2}\left|\alpha \right|\right)}\left\{\left(\left|\alpha \right|R_{t,t}-2\left|b_t\right|\right)\left(2+a\left|\alpha \right|\right)+\frac{1}{a{R}_{t,t}}{\left(R_{t,t}+a\left|b_t\right|\right)}^2\right\}\hfill \\ & >\frac{\left|\alpha \right|}{4a{R}_{t,t}\left(1+\frac{a}{2}\left|\alpha \right|\right)}{\left(a{R}_{t,t}\left|\alpha \right|+R_{t,t}-a\left|b_t\right|\right)}^2\ge 0\hfill \end{array}$$

(14)

Therefore, if |α|> 0 then for (R_t,t + a|b_t|)² < 4aR_t,t τ, we have ΔJ > 0, i.e., the cost function increases.

We now prove the Proposition 2:

Proof: Let the t-th element x_t ≠ 0 be activated as x̂_t = 0 and the updated solution vector is denoted as x̂. The update of the cost function is then given by:

$$\mathrm{\Delta }J=J\left[\widehat{x}\right]-J\left[x\right]=\frac{1}{2}{\left|x_t\right|}^2{R}_{t,t}+\Re \left\{x_t^{*}{b}_t\right\}-\tau \frac{\left|x_t\right|}{1+\frac{a}{2}\left|x_t\right|}$$

(15)

Thus, if ${\left|x_t\right|}^2{R}_{t,t}/2+Re\left\{x_t^{*}{b}_t\right\}-\tau \left|x_t\right|/(1+\frac{a}{2}\left|x_t\right|)<0$, then ΔJ < 0, there exists a nonzero value of the t-th element that decreases the cost function that should be removed from the support.

Combining the DCD algorithm in Table 1 and the above two propositions, we arrive at the HR-DCD algorithm presented in

Table 2. HR-DCD Algorithm.

Step	Equation
	Initialization: x = 0, I = Ø, r = y, b = AH r, R = AH A
1	Choose the first (t-th) element into the support according to:
	t = arg maxk |bk|2; I ={t}, τ = |bt|
2	Repeat until a termination condition is met:
3	Solve ${min}_{x_I}{\Vert y-A_I{x}_I\Vert }_2^2/2+\tau {\sum }_{n\in I}\left|x_n\right|/\left(1+a\left|x_n\right|/2\right)$ on the support I and update r using DCD iterations
4	Update the regularization parameter: τ ← γτ, 0 < γ < 1
5	Remove the t-th element from I according to the rule:
	$\begin{array}{cc}t=arg{min}_{k\in I}{\left|x_k\right|}^2{R}_{k,k}/2+Re\left\{x_k^{*}{b}_k\right\}-\tau \left|x_k\right|/\left(1+a\left|x_k\right|/2\right)& s.t.{\left|x_t\right|}^2{R}_{t,t}/2+Re\end{array}\left\{x_t^{*}{b}_t\right\}-\tau \left|x_t\right|/\left(1+a\left|x_t\right|/2\right)<0$
	If the t-th element is removed, update r: r ← r + xtA(t), I ← I/t
6	Add the t-th element into I according to the rule:
	$\begin{array}{cc}t=arg{max}_{k\in I^c}\left(a\left|b_k\right|-R_{k,k}\right)/a{R}_{k,k}^2\cdot \left\{{\left(R_{k,k}+a\left|b_k\right|\right)}^2-4a{R}_{k,k}\tau \right\}& s.t.{\left(R_{t,t}+a\left|b_t\right|\right)}^2\end{array}>4a{R}_{t,t}\tau$
	If the t-th element is added, update: I ← I ∪ t

.

The HR-DCD algorithm starts with the zero support I = Ø, and τ = |b_t|. For each τ, it updates the support and minimizes the cost function according to the corresponding propositions. Besides, the algorithm will terminate if the iteration times reach the preset parameter L or τ ≤ τ_min (τ_min is a predefined parameter).

Remark 1: When the penalty is logarithmic function, we give two similar propositions as follows:

Proposition 3: Add the t-th (t ∈ I^c) element into the support according to the rule:

$$\begin{array}{cc}t=arg\underset{k\in I^c}{max}\left|b_k\right|& s.t.\left|b_t\right|>\tau \end{array}$$

(16)

Proposition 4: Remove the t-th (t ∈ I) element from the support according to the rule:

$$\begin{array}{cc}t=arg\underset{k\in I}{min}\frac{1}{2}{\left|x_k\right|}^2{R}_{k,k}\Re \left\{x_k^{*}{b}_k\right\}-\tau \frac{1}{a}log\left(1+a\left|x_k\right|\right)& s.t.\frac{1}{2}{\left|x_t\right|}^2{R}_{t,t}\Re \left\{x_t^{*}{b}_t\right\}-\tau \frac{1}{a}log\left(1+a\left|x_t\right|\right)\end{array}<0$$

(17)

Remark 2: When the penalty is l_p norm, we can also give two similar propositions as follows:

Proposition 5: Add the t-th (t ∈ I^c) element into the support according to the rule:

$$\begin{array}{cc}t=arg\underset{k\in I^c}{max}\left|b_k\right|& s.t.{\left|b_t\right|}^{2-p}>2\tau R_{t,t}^{1-p}\end{array}$$

(18)

Proposition 6: Remove the t-th (t ∈ I) element from the support according to the rule:

$$\begin{array}{cc}t=arg\underset{k\in I}{min}\frac{1}{2}{\left|x_k\right|}^2{R}_{k,k}+\Re \left\{x_k^{*}{b}_k\right\}-\tau {\left|x_k\right|}^p& s.t.\end{array}\frac{1}{2}{\left|x_t\right|}^2{R}_{t,t}+\Re \left\{x_t^{*}{b}_t\right\}-\tau {\left|x_t\right|}^p<0$$

(19)

The proofs for Equations (18)–(21) will be shown in the Appendix. By combining the corresponding propositions with related DCD algorithms, it is easy to obtain Hlog-DCD algorithm and Hl_p-DCD algorithm, respectively, and we omit the tabulated expressions here for simplicity.

3.1. Complexity Analysis

The complexity is given by P ≅ 8MN + 4NL + 2C_uN +C_i +2N_debL_g The term 8MN is for computing the initial b. The term 4NL is for selection of elements in the support. The term 2C_uN is for updating b in the total number C_u of successful DCD iterations, and each update requires only 2N real valued additions, as multiplications by α_t are bit-shifts. The term C_i takes into account C_i (total number of DCD iterations) tests to decide if the DCD iteration is successful. The debiasing (2N_debL_g) is now done by using extra N_deb DCD iterations on the finally identified support with size L_g.

3.2. Off-Grid Problem Solution

As stated before, we explore the idea of grid refinement technique [27] to reduce the number of the arbitrarily located potential positions of x. Two-dimensional passive radar imaging is selected as an example to illustrate the procedures of solution, which is realized as follows:

(1)

Under a coarse resolution Δx_q, Δy_q, where Δx_q, Δy_q represent the range step and azimuth step respectively, calculate the coarse localizations (x_n, y_n), = 1, ⋯, S results by finding S largest amplitude peaks using the proposed algorithms, i.e., HR-DCD, Hlog-DCD and Hl_p-DCD.

(2)

Build a denser grid around the estimated locations (x_n, y_n) with a finer resolution (Δx_q, Δy_q) as Figure 2 mentioned in Section 2.

(3)

Set q = q + 1 and $\mathrm{\Delta }{x}_q=\frac{1}{2}\mathrm{\Delta }{x}_q$, $\mathrm{\Delta }{y}_q=\frac{1}{2}\mathrm{\Delta }{y}_q$, repeat for Q times from Step (1) until the desired resolution is achieved.

In order to reduce the multiplications, we use dichotomy for grid refinement, which can use bit-shifts instead of multiplications to be implement in hardware system. According to our experimental results in Section 4, no more than 10-step grid refinement process is needed.

3.3. Extension to Multiple Measurement Vectors (MMV) Case

Although we have derived HR-DCD, Hlog-DCD and Hl_p-DCD algorithms from a single measurement vector (SMV), the results can be generalized to MMV case. In some kinds of applications, such as wideband source localization [42], multiple frequency bins are explored. As we all know that single frequency results can be easily disturbed by noise and interference, and in this issue, we can use a frequency diversity to achieve better performance. A traditional method of modifying the above approach to multi-frequency signals is averaging the results of all the frequencies. But due to the presence of noise, environmental mismatch and interference, different frequencies may give different results, thus this combining method may not work well.

By the similar approach used in [28], our proposed methods require that for all frequencies to have the same support, which utilize the joint sparsity pattern, and choose an element to add to the support according to the corresponding propositions to achieve greater frequency diversity and avoid the possible wrong results. And the final result can be given by averaging all the results of all the frequencies through this way.

4. Experimental Results

In this section, we will present several simulation results which illustrate the effectiveness of the proposed methods. All experiments are performed by using MATLAB R2013a on a PC equipped with an Inter Core i7 3770k CPU, 3.5GHz and 32 GB memory. The state-of-the-art sparse recovery methods such as MP, BP, FOCUSS, SBL, Hl1-DCD, etc. are selected for comparison. As a benchmark we will utilize an oracle sparse recovery (OSR) to represent the best inversion performance, which knows the true support of x. In addition, we set parameter a = 10 in Equation (4) for HR-DCD, in Equation (5) for Hlog-DCD, and p = 0.8 in Equation (6) for Hl_p-DCD, and M_b = 10, τ_min = 0.02, γ = 0.9 by experience from extensive simulation results. Other predefined parameters such as H, L, S, and the grid refinement times Q, etc., are decided according to the actual conditions, which will be given in details in the next text.

4.1. DOA Estimation

As described by Malioutov et al., for one snapshot case, the DOA problem of Equation (5) in [27] is equivalent to Equation (2) in this paper. With the assumption that the sources can be viewed as point sources and their number is small, the underling spatial spectrum is sparse and it can be solved via sparse recovery methods mentioned above. The first simulation considers the scenario that there are five uncorrelated signals impinging from [−47.7°, −28.5°, −9.2°, 11.6°, 30.1°]. The grid is set to be within the range of −90° to 90° with 1° spacing. A 30-element uniform linear array (ULA) spaced in half-wavelength units is used, and the signal-to-noise ratio (SNR) is set to be 20 dB. Here we only consider one snapshot case, and set H = 1, L = S = 8, Q = 10. Figure 4 depicts the solutions solved by different methods, and the vertical dashed lines mark the true directions. Obviously, our methods find the positions and the amplitudes exactly and achieve much sparser solutions than MP and BP. Moreover, they have a better amplitude estimation than Hl1-DCD [40].

We give a time cost for the aforementioned methods in

Table 3. Time consumption of different methods.

	MP	BP	Hl1-DCD	HR-DCD	Hlog-DCD	Hlp-DCD
Time(s)	0.013	0.721	0.130	0.122	0.129	0.126

. From the results, we can see that our methods have the same time cost as the Hl₁-DCD. Although they are a little longer than MP, but they are much less costly than BP. However, the significant advantages of our methods are that they are very easy to implement in a hardware platform because only bit-shifts are needed instead of multiplications similar as [37–40], besides, they can achieve the similar reconstruction performance as BP without off-grid problem.

Next we compare the mean squared error (MSE) of position and amplitude estimation performances by applying the proposed methods and several popular methods (MP, BP, FOCUSS, SBL, Hl1-DCD, OSR), which are defined as follows:

$$\mathit{\text{Position}}\_\mathit{\text{MSE}}=\frac{{\Vert \mathit{\text{Position}}\left(\widehat{\mathbf{x}}\right)-\mathit{\text{Position}}\left(\mathbf{x}\right)\Vert }_2^2}{{\Vert \mathit{\text{Position}}\left(\mathbf{x}\right)\Vert }_2^2}$$

(20)

$$\mathit{\text{Amplitude}}\_\mathit{\text{MSE}}=\frac{{\Vert \widehat{\mathbf{x}}-\mathbf{x}\Vert }_2^2}{{\Vert \mathbf{x}\Vert }_2^2}$$

(21)

where x̂ is an estimate of x.

Figure 5a compares the MSEs of DOA estimation results of different methods under varying SNR, which are averaged by 100 Monte Carlo trials. From the curves of MSE of position versus SNR shown in Figure 5a, it can be seen that Hl₁-DCD and the proposed methods have similar results in position estimation, and outperform their counterparts. However, our methods have better amplitude estimation accuracy than Hl₁-DCD from Figure 5b, and are even very close to the OSR performance.

The MSEs of DOA estimation versus number of signals shown in Figure 6 are obtained over 100 Monte Carlo trials with SNR equals to 20 dB. We can clearly see that our methods perform best compared to other methods with the same parameters, especially on the amplitude estimation results.

4.2. Passive Radar Imaging

In this part, we aim to test the performance of the proposed methods in a passive radar imaging system. Based on the final echo Equation (17) stated in [6] which is similar to Equation (2) in this paper, sparse recovery techniques can be utilized to solve the problem of passive radar imaging. In the simulation, we choose seven digital video broadcasting (DVB) satellites as opportunity transmitters and 10 receivers, and the initialization grids are 21 × 41 represented for the imaging scene 20 m × 20 m (here we use the axis function in matlab to make the scene limited to be 12 m × 10 m so as to achieve a better visual effect). There are nine off-grid scatterers with different reflection coefficients in the scene of interest as the red circles illustrated in Figure 7. The number of frequency samplings corresponding to each transmitter is 10, and the SNR equals to 20 dB. Other simulation parameters are the same as displayed in [6]. Here we set H = 2, L = S = 10 and grid refinement times Q = 10.

Figure 7a shows the poor imaging result by BP due to the reason that in practice the scatterers are not located exactly on the pre-discretized grids, and the echoes are mismatched with the measurement matrix, therefore the performances of CS-based reconstructions are generally unsatisfactory. Figure 7b is the imaging result by Hl₁-DCD, and we can see that it can seek the exact locations of the target, however the amplitude estimation of the reflection coefficients is not precise.

Figure 7c–e represents the imaging results by the proposed methods. As expected, the proposed methods do not have the off-grid problem from the perspective of both the position and amplitude estimation results, which show the potential of our methods to be applied in practical passive radar system.

Figure 8 shows the MSEs of position and amplitude estimation results versus SNR by applying different methods, which are averaged over 100 Monte Carlo trials. In the simulation, we have assumed that there are nine off-grid scatterers with unit reflection coefficients distributed randomly in imaging region. It is seen that our methods outperform the conventional methods (MP, BP, FOCUSS, SBL) in further improvement of MSEs of position and amplitude estimates when off-grid target emerges. In the meantime, they have better amplitude estimation result than Hl₁-DCD with the same parameters shown in Figure 8b. Moreover, the performances of our methods approach the OSR very closely under varying SNR, both in position and amplitude estimates.

4.3. ISAR Imaging

We utilize the CS ISAR imaging model as discussed in [43], and assume that the translational motion compensation has been already achieved and thus only the rotational motion is considered. The relation between the received signal and the complex reflection coefficients can be written as Equation (10) in [43], which is equivalent to Equation (2) in this paper.

We use the quasi real data of an airplane (B-727) provided by the U.S. Naval Research Laboratory to test the feasibility and performance of the proposed methods, which is available on the website http://airborne.nrl.navy.mil/∼vchen/tftsa.html. The stepped frequency radar operates at 9 GHz with the equivalent pulse repetition interval (PRI) of 3.2 ms and has a bandwidth of 150 MHz. For each pulse train, 64 complex range samples were saved, and the file contains 200 successive pulse trains within the long coherent processing interval (CPI). An additive noise is added to the original B-727s data, and the SNR is set to 10 dB. Herein we choose 40 range cell numbers and 20 cross range cell numbers for short CPI test, and the reconstruction of target spatial domain is discretized with 64 range bins and 32 cross range bins, besides, we set H = 2, L = S = 60 and Q = 5. Since in radar imaging field, the position estimates are very important when off-grid problem exists.

Figure 9 are the position recovery results by MP, SBL, FOCUSS, BP, Matrix-Pencil, ESPRIT, Hl₁-DCD and the proposed methods (HR-DCD, Hlog-DCD and Hl_p-DCD) shown as the red circles, and we use the conventional FFT-based ISAR image as the background for fair comparison. As stated above, the imaging performances of the classical methods (MP, BP, FOCUSS, SBL) are not good because of the off-grid errors. Matrix Pencil and ESPRIT can direct solve the positions of the plane, but their performances suffer from low SNR and few snapshots, therefore the imaging results are not good enough. In contrast, our methods can deal with the off-grid target, and achieve satisfying imaging performances which have better visual effect for the airplane shape. Moreover, it is obvious that the imaging results of our methods are better than Hl₁-DCD with the same parameters, especially on the details of airplane wings and tails as illustrated in zoomed-in regions of Figure 9h–j, which demonstrate the advantage of homotopy DCD method with non-convex regularization compared with convex regularization under the same measurements.

In order to quantitatively evaluate the amplitude estimation performances of the obtained ISAR images via different methods, we use the correlation value [44] to evaluate the similarity between the recovered images and the reference image, and the image entropy [44] to measure the focusing quality of the recovered images. They are defined as follows:

$$\mathit{\text{Cor}}\left({\mathbf{x}}_r,{\mathbf{x}}_o\right)=\frac{\left|\langle {\mathbf{x}}_r,{\mathbf{x}}_o\rangle \right|}{{\Vert {\mathbf{x}}_r\Vert }_2{\Vert {\mathbf{x}}_o\Vert }_2}$$

(22)

$$\mathit{\text{Ent}}\left({\mathbf{x}}_o\right)=-\sum _i{p}_{i}log{p}_i$$

(23)

where x_r and x_o denote the reference image vector using full measurements and the recovered image vector by different methods, and p_i is the histogram of the recovered gray level (0∼255) image.

The correlation and entropy values of the recovered ISAR images under different methods are summarized in

Table 4. Correlation and entropy values by different methods.

	MP	SBL	FOCUSS	BP	Matrix Pencil	ESPRIT	Hl1-DCD	HR-DCD	Hlog-DCD	Hlp-DCD
Cor	0.817	0.890	0.887	0.913	0.816	0.798	0.925	0.952	0.957	0.965
Ent	0.826	0.811	0.836	0.859	0.810	0.802	0.713	0.689	0.683	0.613

. It can be seen that the proposed methods have a higher correction with the reference image and a lower entropy than their counterparts, and thereby exhibit a better capability of target-information extraction.

5. Conclusions

Two major problems of applying CS to sparse continuous signal reconstruction are discussed in this paper, and we have presented computationally efficient and accurate methods to overcome the difficulties. For solving the off-grid problem, we utilized the grid refinement technique. In the meantime, in order to obtain high recovery performance and keep low calculation complexity, we propose a fast and accurate homotopy DCD reconstruction combined with three typical non-convex regularizations, which promotes sparsity more strongly than any convex penalty function can. Extensive experiments have been conducted to validate and compare the performances of the proposed methods with several popular solvers. Our future work will try to synthesize with parallel sparse optimization technique using multi-core CPUs/GPUs [45], which may provide a viable solution to real-time potential applications.