Sparse Recovery Algorithm for Compressed Sensing Using Smoothed l₀ Norm and Randomized Coordinate Descent

Abstract

Compressed sensing theory is widely used in the field of fault signal diagnosis and image processing. Sparse recovery is one of the core concepts of this theory. In this paper, we proposed a sparse recovery algorithm using a smoothed l₀ norm and a randomized coordinate descent (RCD), then applied it to sparse signal recovery and image denoising. We adopted a new strategy to express the (P⁰) problem approximately and put forward a sparse recovery algorithm using RCD. In the computer simulation experiments, we compared the performance of this algorithm to other typical methods. The results show that our algorithm possesses higher precision in sparse signal recovery. Moreover, it achieves higher signal to noise ratio (SNR) and faster convergence speed in image denoising.

1. Introduction

In recent years, compressed sensing (CS) has become an essential mathematical tool in the field of information theory [1,2]. It has been widely used in the field of fault diagnosis [3], image processing [4], data compression [5,6], and so on. CS theory suggests that a signal can be recovered from fewer samples than suggested by the Shannon–Nyquist sampling rate given that the original signal is sparse or approximately sparse in certain representation domains.

As shown in Figure 1, assuming that signal f (n × 1) will be recovered from incomplete sampling values, the sampled value of the signal f is as follows:

$$\mathit{y}=\mathit{\Phi }\mathit{f}$$

(1)

where Φ (m × n, m << n) is the measurement matrix, and y (m × 1) is an measurement vector. When the dimension of y is much lower than that of f, Equation (1) becomes an underdetermined equation and f is unable to be calculated. However, if the number of non-zero elements in f is K (K ≤ m), f can be calculated, in which case K is called sparsity.

Most natural signals are not sparse in practical applications. The sparse representation of f is obtained by a specific mathematical transformation Ψ, as shown in Equation (2) and Figure 2.

$$\mathit{f}=\mathit{\Psi }\mathit{s}$$

(2)

where Ψ is the sparse dictionary, and s is the sparse representation coefficient, such that signal f can then be recovered indirectly by constructing the sparse coefficients, which is shown in Equation (3) and Figure 3.

$$\mathit{y}=\mathit{\Phi }\mathit{\Psi }\mathit{s}$$

(3)

Here we set A = ΦΨ, then Equation (3) changes to y = As. If A satisfies restricted isometry property (RIP), s can be accurately recovered [1].

The CS theory consists of two parts. The first is the sparse representation. Because a natural signal is not sparse in a usual case, we need to first make it sparse by mathematical transformation. The second is sparse recovery, which means to recover the sparse signal. Among the two parts, the sparse recovery method directly determines the quality of signal recovery [7]. In this paper, we focus on proposing a sparse recovery algorithm for CS. In the following sections, signal s is set to be a sparse signal vector.

Generally speaking, the purpose of sparse recovery is to solve the l₀ norm minimization problem, as shown in Equation (4):

$$\{\begin{array}{c}\min {\Vert \mathit{s}\Vert }_0\\ \mathrm{s}.\mathrm{t}. \mathit{y}=\mathit{A}\mathit{s}\end{array}$$

(4)

Equation (4) aims at seeking the minimal ${\Vert \mathit{s}\Vert }_0$ solution of $\mathit{y}=\mathit{A}\mathit{s}$, which is also called the (P₀) problem. However, if m << n, the (P⁰) problem becomes the NP-hard, which means that we cannot solve Equation (4) directly in that condition [8,9,10].

One way to solve such (P⁰) problems are greedy algorithms, which include matching pursuit (MP) [11], orthogonal matching pursuit (OMP) [12], stage-wise orthogonal matching pursuit (StOMP) [13], compressive sampling matching pursuit (CoSaMP) [14], subspace pursuit (SP) [15], gradient pursuit (GP) [16], and so on. These algorithms are widely used in sparse recovery, but there are strict requirements for the sparsity. For instance, the OMP algorithm, only when $K<\sqrt{n}$ or A satisfies certain conditions, s can be recovered successfully [10].

In addition, we considered the following optimization problem:

$$\{\begin{array}{c}\min {\Vert \mathit{s}\Vert }_p^p\\ \mathrm{s}.\mathrm{t}. \mathit{y}=\mathit{A}\mathit{s}\end{array}$$

(5)

According to the literature [10], there are notably different properties of l_p norms for different values of p. As is shown in Figure 4, when p > 1, Equation (5) is a convex optimization problem and its solution is not sparse (the elements of the intersection point are all not 0). While when 0 ≤ p < 1, it becomes a non-convex optimization problem and its solution is sparse (only one element of the intersection point is 0). In particular, when p = 1, Equation (5) is called a basis pursuit (BP) algorithm and under certain conditions its solution is also sparse, but may not result in the sparsest solution compared to the case when 0 ≤ p < 1 [9].

In a term, if 0 ≤ p < 1, the solution of Equation (5) is more sparse. Unfortunately, this results in a non-convex optimization problem. In order to resolve this problem, researchers adopt the sum of a set of smooth concave functions to express the (P₀) problem approximately, which is called the smooth l₀ norm approximation algorithm.

For instance, literature [17] adopted an arctangent function as

$$f\left(s_i\right)=\frac{\pi }{2}\arctan \left(\frac{s_i^2}{2\sigma }\right)$$

(6)

where s_i is the elements in vector s. Then the (P₀) problem can be approximately expressed as

$$\{\begin{array}{c}\min {\displaystyle {\displaystyle \sum }_{i=1}^n}\frac{\pi }{2}\arctan \left(\frac{s_i^2}{2\sigma }\right)\\ \mathrm{s}.\mathrm{t}. \mathit{y}=\mathit{A}\mathit{s}\end{array}$$

(7)

Similarly, some research [18] adopted a hyperbolic tangent function as

$$f\left(s_i\right)=\frac{\exp \left(\frac{s_i^2}{2{\sigma }^2}\right)-\exp \left(-\frac{s_i^2}{2{\sigma }^2}\right)}{\exp \left(\frac{s_i^2}{2{\sigma }^2}\right)+\exp \left(-\frac{s_i^2}{2{\sigma }^2}\right)}$$

(8)

The recovery algorithm based on a smoothed l₀ norm and modified Newton method (L0MN) in the literature [19] adopted a function as

$$f\left(s_i\right)=\frac{\sigma }{s_i+\sigma }$$

(9)

The recovery algorithm based on a smoothed l₀ norm and gradient projection (L0GP) in the literature [20] adopted a function as

$$f\left(s_i\right)=\frac{{\left|s_i\right|}^2}{{\left|s_i\right|}^2+{\sigma }^2}$$

(10)

The smooth l₀ norm approximation algorithms can be easily solved by Lagrange multiplier technique. However, because f(s_i) is non-convex, these algorithms have the possibility of falling into a local optimum.

In this paper, we have improved Equation (9) to establish a new strategy to express the (P₀) problem approximately. However, the resulting optimization problem is still not a convex problem, hence finding a global optimal solution is not a trivial. For this, we selected a randomized coordinate descent (RCD) algorithm for solving the optimization, which does not guarantee one to find the global optimum, but it can avoid many local minima, and probably can land close to the global optimum.

2. Sparse Recovery Algorithm

2.1. Cost Function

Observing Equation (9), we find that it is a function with very simple structure. According to the literature [19], this equation can reduce the computational complexity while ensuring precision. However, as is shown in Figure 5, in very few cases, $s_i+\sigma =0$ will occur. This case can lead to the discontinuity of functions, which will affect the calculation results to a certain extent.

In order to avoid this defect, we give a continuous smooth fractional function as Equation (11),

$$f_{\sigma }\left(s_i\right)=-\frac{\sigma }{s_i^2+\sigma }$$

(11)

where $\sigma >0$. In this way, we guarantee the continuity of this function, which is shown in Figure 6.

As shown in Figure 7, when the value of σ decreases gradually, $f_{\sigma }\left(s_i\right)$ gradually approaches the form of l₀ norm, that is:

$$\underset{\sigma \to 0}{\lim }{f}_{\sigma }\left(s_i\right)=\{\begin{array}{c}-1, s_i=0\\ 0, s_i\ne 0\end{array}$$

(12)

It should be noted that when $\sigma \to 0$, although $f_{\sigma }\left(s_i\right)$ and l₀ norm look similar, they are essentially different. In order to use Equation (12) to express the (P₀) problem, we make the following attempts:

$$\{\begin{array}{c}\min {\displaystyle {\displaystyle \sum }_{i=1}^n}{f}_{\sigma }\left(s_i\right)\\ \mathrm{s}.\mathrm{t}. \mathit{y}=\mathit{A}\mathit{s}\end{array}$$

(13)

As is shown in Equation (11), when $\sigma \to 0$, the value $f_{\sigma }\left(s_i\right)$ can only be 0 or −1. If we want to minimize ${{\displaystyle \sum }}_{i=1}^n{f}_{\sigma }\left(s_i\right)$, we need to make as many terms as possible in this polynomial equal to −1. Then, according to Equation (11), this is equivalent to that which we need to make as many elements as possible in vector s equal to 0. Thus, this is consistent with the meaning of l₀ norm minimization.

In a word, when $\sigma \to 0$, the optimal solution of Equations (4) and (13) are identical. Therefore, we can adopt Equation (13) to express the (P⁰) problem approximately.

Equation (13) is a typical optimization problem with an equality constraint. The Lagrange multiplier method can transform such problems into unconstrained optimization problems. The cost function can be constructed as follows,

$$\min L\left(\mathit{s}\right)=\frac{1}{2}{\Vert \mathit{A}\mathit{s}-\mathit{y}\Vert }_2^2+\lambda {\displaystyle \sum }_{i=1}^N{f}_{\sigma }\left(s_i\right)$$

(14)

where λ is the regularization parameter.

2.2. Solution

Generally, the gradient descent method (or the idea of gradient descent) can be used to solve Equation (14). The procedures of gradient descent method are as follows:

The gradient $\nabla L\left(\mathit{s}\right)$ of $L\left(\mathit{s}\right)$ at point $\mathit{s}$ is:

$$\begin{array}{cc}\hfill \nabla L\left(\mathit{s}\right)=\frac{\partial \mathit{L}}{\partial \mathit{s}}& ={\mathit{A}}^{\mathrm{T}}\left(\mathit{A}\mathit{s}-\mathit{y}\right)+\lambda {\left(\frac{\partial f}{\partial s_1},\frac{\partial f}{\partial s_2},\dots ,\frac{\partial f}{\partial s_n}\right)}^{\mathrm{T}}\hfill \\ & ={\mathit{A}}^{\mathrm{T}}\left(\mathit{A}\mathit{s}-\mathit{y}\right)+\lambda {\left(\frac{2\sigma s_1}{{\left(\sigma +s_1^2\right)}^2},\frac{2\sigma s_2}{{\left(\sigma +s_1^2\right)}^2},\dots ,\frac{2\sigma s_n}{{\left(\sigma +s_1^2\right)}^2}\right)}^{\mathrm{T}}\hfill \end{array}$$

(15)

If we want to minimize $L\left(\mathit{s}\right)$, we need to update s to the direction of negative gradient:

$${\mathit{s}}_{\left(k+1\right)}={\mathit{s}}_{\left(k\right)}-h\left[\nabla L\left({\mathit{s}}_{\left(k\right)}\right)\right]$$

(16)

where h is the step size of the iterative process. When the initial value and termination condition are given, the optimal solution of Equation (14) can be obtained by iteration, and that is the standard gradient descent method.

However, the gradient descent method is not intended for non-convex problems, because it is easy to be trapped in local minima. Therefore, we need to adopt strategies to get it out of local minima. The common methods for this include the simulated annealing (SA) algorithm [21], the genetic algorithm (GA) [22], the randomized coordinate descent (RCD) [23] method, and so on. RCD is a special case of stochastic gradient descent (SGD) [24], which is widely used to optimize non-convex problems because of its simple principle and fast convergence speed.

Unlike the traditional gradient descent method, RCD adds random factors to the calculation of gradient. Therefore, even if it falls into the local minima, its gradient may not be zero, so it may jump out of the local minima and continue to iterate [23].

In each iteration process, RCD method does not derive partial derivatives of all variables, but derives partial derivatives of one (or several) variables. The detailed measures are as follows:

Let s_rand be a random element in vector s, and then we can calculate the gradient, as in Equation (17).

$$\begin{array}{cc}\hfill {\nabla }_{s_{\mathrm{rand}}}L\left(\mathit{s}\right)={\frac{\partial \mathit{L}}{\partial \mathit{s}}|}_{s_{\mathrm{rand}}}& ={\mathit{a}}_{\mathrm{rand}}{}^{\mathrm{T}}\left(\mathit{A}\mathit{s}-\mathit{y}\right)+\lambda {\left(0,\dots ,\frac{\partial f}{\partial s_{\mathrm{rand}}},\dots ,0\right)}^{\mathrm{T}}\hfill \\ & ={\mathit{a}}_{\mathrm{rand}}{}^{\mathrm{T}}\left(\mathit{A}\mathit{s}-\mathit{y}\right)+\lambda {\left(0,\dots ,\frac{2\sigma s_{\mathrm{rand}}}{{\left(\sigma +s_{\mathrm{rand}}^2\right)}^2},\dots ,0\right)}^{\mathrm{T}}\hfill \end{array}$$

(17)

where a_rand denotes the column corresponding to s_rand in matrix A. Then Equation (16) can be updated as:

$${\mathit{s}}_{\left(k+1\right)}={\mathit{s}}_{\left(k\right)}-h\left[{\nabla }_{s_{\mathrm{rand}}}L\left({\mathit{s}}_{\left(k\right)}\right)\right]$$

(18)

Now, we discuss how to choose the step size h. In the process of iteration, the value of the function is required to decrease, which is shown as Equation (19).

$$L\left[{\mathit{s}}_{\left(k+1\right)}\right]<L\left[{\mathit{s}}_{\left(k\right)}\right]$$

(19)

That is to say, when the direction of iteration is determined, we need to minimize $L\left[{\mathit{s}}_{\left(k+1\right)}\right]$. This problem can be reduced to Equation (20).

$$\{\begin{array}{c}\min \phi \left(h\right)\\ \mathrm{s}.\mathrm{t}. 0\le h\le h_{\left(\max \right)} \end{array}$$

(20)

where $\phi \left(h\right)={\mathit{s}}_{\left(k\right)}-h\left[{\nabla }_{s_{\mathrm{rand}}}L\left({\mathit{s}}_{\left(k\right)}\right)\right]$, and $h_{\left(\max \right)}=2{\left(\left|{\mathit{s}}_{\left(k\right)}^{\mathrm{T}}{\mathit{s}}_{\left(k\right)}\right|\right)}^{-1}.$ Here we adopt Goldstein method to solve this problem, which can be referred to in the literature [25].

Because our method is to solve the smooth l₀ norm based sparse recovery problem using RCD, we name it L0RCD.

2.3. Algorithm Description

Now, we summarize the sparse recovery algorithm in this paper, as in

Table 1. The description of the smooth l₀ norm based sparse recovery algorithm using RCD (L0RCD).

Algorithm: L0RCD
Input:A, y, λ, h, K, and the termination condition: ${\sigma }_{\min }={10}^{-3}{\sigma }_0$;
Output: The recovery signal s*;
Initialization:s0 = AT(AAT)−1y, ${\sigma }_0$ = 1;
Step 1: Calculate gradient according to Equation 17;
Step 2: Iterate, ${\mathit{s}}_{\left(k+1\right)}={\mathit{s}}_{\left(k\right)}-h\left[{\nabla }_{s_{\mathrm{rand}}}L\left({\mathit{s}}_{\left(k\right)}\right)\right]$;
Step 3: If ${\Vert {\mathit{s}}_{\left(k+1\right)}-{\mathit{s}}_{\left(k\right)}\Vert }_2\le \sigma$, update s and $\sigma$, where ${\sigma }_{\left(k+1\right)}=0.7{\sigma }_{\left(k\right)}$, or else go to Step 1 and continue iterating.
Step 4: Take the updated s and $\sigma$ as initial value, and repeat step 1 to step 3, stop iterating until $\sigma <{\sigma }_{\min }$ and output to s*.

.

In Table 1, ${\sigma }_0$ and s₀ are the initial values of $\sigma$ and s, respectively. A^T(AA^T)⁻¹y is the least square solution of y = As, and λ can be estimated by generalized cross-validation (GCV) method [26]. According to literature [17,18,19,20], in order to ensure that the search scope is more and more refined in the iteration process, we take a set of descending sequences of $\sigma$, where ${\sigma }_{\left(k+1\right)}=0.7{\sigma }_{\left(k\right)}$.

2.4. Computational Cost Analysis

The main computational complexity of L0RCD lies in the product of the vectors, as well as the calculation of partial derivatives in the iterative process. The computational burden of ${\mathit{a}}_{\mathrm{rand}}{}^{\mathrm{T}}\left(\mathit{A}\mathit{s}-\mathit{y}\right)$ is O(m²), while the burden of ${\mathit{A}}^{\mathrm{T}}\left(\mathit{A}\mathit{s}-\mathit{y}\right)$ is O(n²m). In addition, in the iteration process, our method only derives partial derivatives of one variable, while the general gradient descent methods derive partial derivatives of all variables. This means that, in a single iteration, the computational complexity of our method is much less than that of the general gradient descent method. Although adopting RCD will result in more iterations, its computational burden may still be lower [23].

3. Comparative Experiment

In this section, we will test the performance of the sparse recovery algorithms in the field of signal recovery and image denoising. All of the experiments are implemented in MATLAB R2012a on a personal computer with a 2.59 GHz Intel Pentium dual-core processor and 2 GB memory running the Microsoft Windows XP operating system.

3.1. Signal Recovery Experiment

In this section, for matrix A, we adopt a Gaussian random matrix with variance of 1, where the sizes of A are m × n, and m < n. The size of the sparse vector s is n × 1, and the 0 elements distribute randomly in s. The values of non-zero elements in vector s are also random. Vector y can be calculated by $\mathit{y}=\mathit{A}\mathit{s}$. Under the condition that y and A are known, we reconstruct s using the sparse recovery algorithm; the recovery result is denoted as s^*.

Here we use precision to evaluate the performance of the algorithms. The precision is defined as follows:

$$precision=1-\frac{{\Vert \mathit{s}-{\mathit{s}}^{*}\Vert }_2}{{\Vert \mathit{s}\Vert }_2}$$

(21)

3.1.1. Comparison of RCD and the Modified Newton Method

In order to highlight the advantages of RCD, we compared it with the modified Newton (MN) method from the literature [19]. In addition, to give a fair comparison, we constructed the same l₀ norm approximation by Equation (9). That is, we adopted the same cost function but different optimization methods in the signal recovery experiments of this section.

In the first experiment, m = 60 and n = 100 were fixed, and K varied between 5 and 40. The precisions are given in Figure 8.

From Figure 8, we see that when K ≤ 20, both RCD and MN can be used to recover s precisely. However, with the increase of K, RCD performs relatively better than MN. This shows that RCD is more appropriate under these conditions.

In the second experiment, K = 20 and n = 100 were fixed, and m varied from between 30 and 90. The precisions are given in Figure 9.

From Figure 9, we can see that when m ≥ 55, both the two methods perform well. However, with the decrease of m, RCD performs relatively better than MN.

3.1.2. Comparison of L0RCD and the Other Sparse Recovery Algorithms

In this section, we will compare the L0RCD with OMP, MP, basis pursuit(BP), and L0GP. Similarly to the first experiment, m = 60 and n = 100 were fixed, and K varied between 5 and 40. The precision of the five methods are given in Figure 10. From Figure 10, we see that with the increase of K, L0RCD performs better than the other four algorithms in precision.

In the second experiment, K = 20 and n = 100 are fixed, and m varies from between 30 and 90. The precision of the five methods are given in Figure 11. From Figure 11, we see that with the decrease of m, L0RCD performs better than the other four algorithms in precision.

In the rest of this subsection, we give a rough analysis for these sparse recovery algorithms. From Figure 10 and Figure 11, we find that when K ≤ 10 or m ≥ 65, the precision of the five algorithms are all 1 (or very close to 1). For the pursuit algorithms (OMP, MP and BP), they have stricter requirements for sparsity and measurements.

Here we do not intend to claim that the approximate l₀ norm algorithms are “better” than the pursuit algorithms. In fact, approximate l₀ norm algorithms need to appropriate parameters, and this often requires experience or a large number of repeated tests. In comparison, when the sparsity or measurements meet the requirements, the application of pursuit algorithms is more mature. For instance, the MOD (method of optimal directions) and K-SVD (K-singular value decomposition) used for the sparse representation of signals both employ the pursuit algorithms to implement sparse coding [10].

3.2. Image Denoising Experiment

When sparse recovery algorithm is applied to image denoising [27], Equation (4) should be changed to:

$$\{\begin{array}{c}\min {\Vert \mathit{s}\Vert }_0\\ \mathrm{s}.\mathrm{t}. {\Vert \mathit{y}-\mathit{A}\mathit{s}\Vert }_2^2\le ϵ\end{array}$$

(22)

where $ϵ$ is a parameter related to noise level. The solution method of it is the same as that of (P⁰) problem, and ${\Vert \mathit{y}-\mathit{A}\mathit{s}\Vert }_2^2\le ϵ$ can be regard as the relaxation of $\mathit{y}=\mathit{A}\mathit{s}$. According to literature [10], the pursuit algorithms are all applicable to Equation (22).

For approximate l₀ norm method algorithms, the cost function for Equation (22) can be established as:

$$\min \left(\mathit{s}\right)=\frac{1}{2}{\Vert \mathit{A}\mathit{s}-\mathit{y}\Vert }_2^2+\lambda {\Vert \mathit{s}\Vert }_0$$

(23)

In this section, we compared the image denoising performance of these methods under the same experimental conditions. We adopted signal to noise ratio (SNR) and time consumption to evaluate their performance, where SNR was defined as Equation (24).

$$\mathrm{SNR}=10·\lg \frac{{{\displaystyle \sum }}_i|s_i{|}^2}{{{\displaystyle \sum }}_i|s_i^{*}-s_i{|}^2}$$

(24)

where $s_i^{*}$ are the elements in s^*. The higher the SNR, the better the image quality is.

We adopted three images named, respectively Lena, Cameraman, and Boat to test the denoising performance of the algorithms. The sizes of the three images are all 256 × 256 pixels. Gaussian noise with a mean of 0 and a standard deviation of 25 was added to the images.

Because an image is not a sparse signal, first, we used a two dimensional discrete cosine transform (2D-DCT) to represent these images sparsely. For matrix Φ, we adopted a Gauss random matrix with variance of 1, where the size of Φ was 128 × 256.

The intuitive denoising results are shown in Figure 12, Figure 13 and Figure 14. The SNR and time consumption are shown in

Table 2. The quantitative results (signal to noise ratio (SNR)/Time consumption(s)) of the six algorithms.

Algorithm	Lena	Cameraman	Boat
OMP	26.27/43.02	25.98/44.57	25.11/47.54
StOMP	27.01/39.07	27.33/42.16	25.76/45.55
BP	28.34/76.37	28.97/78.39	29.05/84.29
L0GP	32.13/25.14	33.14/27.33	31.06/28.77
L0MN	33.08/18.04	34.59/20.25	33.27/23.05
L0RCD	34.95/17.12	35.65/20.01	34.05/21.86

.

From Figure 12, Figure 13 and Figure 14, we see that the image denoising performance of L0RCD is better than that of other five methods. From Table 2, we see that L0RCD only achieves higher SNR, but also consumes less time.

4. Conclusions

In this paper, a novel method based on RCD and approximate l₀ norm is proposed to solve the optimization problem of sparse recovery. For the local minima defect of approximate l₀ norm, we adopted a solution strategy with RCD. The simulations are carried out on sparse signal recovery and image denoising. Compared with the typical sparse recovery algorithms, our method performs better in precision. However, the success of the algorithm still depends on careful selection of the parameters, and how to overcome this defect is our further research goal.