# A Bat-Inspired Sparse Recovery Algorithm for Compressed Sensing

Authors: Wanning Bao; Haiqiang Liu; Dongbo Huang; Qianqian Hua; Gang Hua
Journal: Computational Intelligence and Neuroscience
 (2018)

---

## Abstract

Compressed sensing (CS) is an important research area of signal sampling and compression, and the essence of signal recovery in CS is an optimization problem of solving the underdetermined system of equations. Greedy pursuit algorithms are widely used to solve this problem. They have low computational complexity; however, their recovery performance is limited. In this paper, an intelligence recovery algorithm is proposed by combining the Bat Algorithm (BA) and the pruning technique in subspace pursuit. Experimental results illustrate that the proposed algorithm has better recovery performance than greedy pursuit algorithms. Moreover, applied to the microseismic monitoring system, the BA can recover the signal well.


---

## Body

## 1. Introduction

In recent years, the collected data quantity by the Internet of things (IoT) has increased dramatically. To facilitate storage and transmission, signal sampling based on the traditional Shannon–Nyquist sampling theory works as follows: massive data are collected at the sampling stage and most are discarded at the compression stage. This process is extremely wasteful. CS, which was proposed by D. Donoho, E. Candès, and T. Tao, brings a revolutionary breakthrough in signal sampling [1–3]. If only K nonzero elements exist in a signal, the signal is called K-sparse signal with sparsity K. CS theory indicates that sparse signals can be recovered more accurately using less measurements than the Nyquist sampling principle.The signal recovery problem in CS is a nonconvex combinatorial optimization problem. It can be solved by three types of algorithms. The first type is greedy pursuit algorithms, including Matching Pursuit (MP), Orthogonal Matching Pursuit (OMP) [4], Generalized Orthogonal Matching Pursuit (GOMP) [5], and Subspace Pursuit (SP) [6]. They are two-stage algorithms. In the first stage, they seek support iteratively. In the second stage, the signal can be recovered by using the least square method. This kind of algorithm is effective; however, they have weak recovery performance. The second kind of algorithm is a convex optimization algorithm, which recovers signals by converting nonconvex problems to convex problems. The most common convex optimization algorithms are Basic Pursuit (BP) and Gradient Projection for Sparse Reconstruction (GPSR) [7]. Although they have high accuracy, their computational complexity is high. The third kind of recovery algorithm is Bayesian algorithms [8], which can be divided into the following two types: Maximum A Posteriori (MAP) Estimation and Hierarchical Bayesian. The former underlines a signal distribution, and the latter introduces one or two variables, which control the sparse signal. Bayesian algorithms achieve a balance of high accuracy and short recovery time [9, 10]. In recent years, some scholars have applied Swarm Intelligence Algorithms to signal recovery in CS, such as Particle Swarm Optimization (PSO) [11] and the Grey Wolf Optimizer Algorithm [12]. These methods have good global search ability. However, they have a slow convergence velocity and are easy to be trapped into local optimum.Recently, Yang et al. simulated bat echolocation and proposed a Bat Algorithm (BA) [13] based on stochastic optimization, which has simple structure, fast convergence, strong robustness, and parallel computing. The BA combines the advantages of PSO, Harmony Search (HS), Simulated Annealing (SA), and other algorithms. Many scholars applied it to optimization problems, such as numerical optimization, economic dispatch with random wind power [14], multirobot formation reconfiguration [15], microgrid operation management [16], and image processing [17]. In these areas, the BA shows better performance than do other intelligence optimizer algorithms.In this paper, inspired by the BA and the greedy pursuit algorithm, an intelligence recovery algorithm for CS is proposed. OMP is used to initialize the swarm. The BA and the pruning technique are combined to update populations. The proposed algorithm outperforms the traditional greedy pursuit algorithms. Moreover, it runs faster than other swarm intelligence algorithms such as PSO.The remainder of this paper is organized as follows: Section2 briefly introduces the basic theory of CS and the BA. In Section 3, the proposed algorithm is introduced in detail. In Section 4, experiments are done to evaluate the algorithm performance. Finally, conclusions are drawn in Section 5.The notations used in the rest of this paper include the following:⋅ denotes the absolute value of a real number or the cardinality of a set. ⋅pp=0,1,2,… represents lp−norm. For the vector x∈Rn, x^ is the estimate of x. In a matrix A∈Rm×n, aj is the j-th column of A. E=1≤i≤nxi≠0, which denotes the positions of nonzero elements in x, and it is called support set. S=1,2,⋯,n is the universal set of support. The matrix AE=ajj∈E⊆S is submatrix of A. The pseudoinverse of a matrix A is defined as A+=ATA−1AT, where T denotes the transpose operator.


## 2. Basic Theory

### 2.1. Compressed Sensing

Supposing ann-dimensional signal y=y1,y2,…,ynT can be represented on a basis Ψ∈Rn×n as y=Ψx, where x=x1,x2,…,xnT, if there are only K nonzero elements in x=x1,x2,…,xnT, y is called K−sparse signal with sparsity K, and the positions of nonzero elements in x are called the support set, which can be denoted as E=1≤i≤nxi≠0. The signal y can be sampled by projection onto the measurement matrix Φ∈Rm×nm≪n as(1)b=Φy=ΦΨx=Ax,where b=b1,b2,…,bmT is the measurement signal, y is the original signal, x is the sparse signal, Φ is the measurement matrix, and Ψ is the sparse basis. A∈Rm×n is the sensing matrix.It can be seen that the data quantity is greatly reduced becauseb=Φy is an underdetermined system of equations. There are infinitely many solutions to recover y from b. However, if y is sparse, it can be recovered from b by solving the following l0 minimization problem:(2)minxx0,s.t.Ax=b,where ⋅0 represents l0−norm.


### 2.2. Orthogonal Matching Pursuit (OMP)

In the proposed algorithm, we use OMP to determine the initial solution.The OMP algorithm can be stated in Algorithm1 [4].Algorithm 1: Orthogonal Matching Pursuit (OMP).
Input: Sensing matrix A, measurement signal b, sparsity K.
Output: Estimated signal x^, estimated support set E^
Initialization: Residual signal r0=b, estimated index set E^0=ϕ, estimated atomic set A0=ϕ, iteration number t=1;
The t-th1≤t<Kiteration:
Step 1: Find λt=argmaxj=1,2,…,Nrt−1,aj
Step 2: Update the estimated support set and the estimated atomic set E^t=E^t−1∪λt,At=At−1∪aλt;
Step 3: Estimate the signal x^t=AtTAt−1AtTb;
Step 4: Update the residual rt=b−Atx^t, update iteration counter t=t+1;
Judgement: If t≤K, continue iterating; if t>K, stop the iteration, and output the estimated support set E^ and the estimated signal x^.


### 2.3. Bat Algorithm

In 2010, Yang et al. simulated the characteristics of bats’ predatory behaviors and proposed BA [13, 18], which is a new Swarm Intelligence Optimization Algorithm.The first step is randomly initializing the bat positionsxi0 and the velocities vi0. Then, the BA updates the velocities and positions according to certain strategies. The velocities and positions of the i-th bat at t+1 iteration are updated as follows:(3)fi=fmin+fmax+fminβ,(4)vit+1=vit+xit+x∗⋅fi,(5)xit+1=xit+vit+1,where fi is the pulse frequency of the i-th bat, fi∈fmin,fmax, β∈0,1 is a random number. xit and xit+1 represent the positions of the i-th bat at the t and the t+1 iteration, respectively; vit and vit+1 are the velocities of the i-th bat at the t and the t+1 iteration, respectively; and x∗ denotes the global best position of the swarm.For each bat, the algorithm generates a random numberrand1∈0,1. If rand1>rit, the position of this bat is updated according to formula (6), where rit is the rate of pulse emission. The random walk is regarded as a procedure of the local search [13].(6)xit+1=xj∗+εLt,where xj∗ is a random solution in the best solution set composed of the best solution of each bat, ε∈−1,1 is a random number, and Lt is the mean of the pulse loudness in the t-th iteration.In addition, with the iterative number increasing, the pulse loudnessLi and the rate of pulse emission ri need to be updated. During prey searching, the pulse loudness Li will gradually decline, and the rate of pulse emission ri will also gradually increase. The updating formulas of Li and ri are found in formulas (7) and (8), respectively:(7)Lit+1=αLit,(8)rit+1=ri01−exp−γt,where α and γ are constants, 0<α<1, and 0<γ.The BA algorithm can be described in Algorithm2 [13].Algorithm 2: Bat Algorithm (BA).
Fitness function:
f
x
Initialization:
Randomly initialize the positions xi and velocities vi of the bat population.
Initialize the pulse frequency fi of each bat according to formula (3).
Initialize the rate of pulse emission ri0 and the pulse loudness Li0.
Calculate the fitness value of each bat and find the best solution in the population.
Set the stopping criterion σ.
Termination condition:
Iteration counter reaches tmax or the fitness value is smaller than σ.
Iteration:
Step 1: Update the bat velocity vit and position xit according to formulas (3)∼(5);
Step 2: If (rand1>rit), select a random solution in the best solutions set and update position of this bat by the random walk around the selected best solution as shown in formula (5);
Step 3: Calculate the fitness value of each bat’s new position;
Step 4: If (rand2>Litandfxit<fx∗), accept the new solution and adjust rit, Lit according to formulas (7) and (8);
Step 5: Find the new best solution x∗ and update the best solutions set;
Step 6: If termination condition is satisfied, stop the iteration; otherwise, continue.


## 2.1. Compressed Sensing

Supposing ann-dimensional signal y=y1,y2,…,ynT can be represented on a basis Ψ∈Rn×n as y=Ψx, where x=x1,x2,…,xnT, if there are only K nonzero elements in x=x1,x2,…,xnT, y is called K−sparse signal with sparsity K, and the positions of nonzero elements in x are called the support set, which can be denoted as E=1≤i≤nxi≠0. The signal y can be sampled by projection onto the measurement matrix Φ∈Rm×nm≪n as(1)b=Φy=ΦΨx=Ax,where b=b1,b2,…,bmT is the measurement signal, y is the original signal, x is the sparse signal, Φ is the measurement matrix, and Ψ is the sparse basis. A∈Rm×n is the sensing matrix.It can be seen that the data quantity is greatly reduced becauseb=Φy is an underdetermined system of equations. There are infinitely many solutions to recover y from b. However, if y is sparse, it can be recovered from b by solving the following l0 minimization problem:(2)minxx0,s.t.Ax=b,where ⋅0 represents l0−norm.


## 2.2. Orthogonal Matching Pursuit (OMP)

In the proposed algorithm, we use OMP to determine the initial solution.The OMP algorithm can be stated in Algorithm1 [4].Algorithm 1: Orthogonal Matching Pursuit (OMP).
Input: Sensing matrix A, measurement signal b, sparsity K.
Output: Estimated signal x^, estimated support set E^
Initialization: Residual signal r0=b, estimated index set E^0=ϕ, estimated atomic set A0=ϕ, iteration number t=1;
The t-th1≤t<Kiteration:
Step 1: Find λt=argmaxj=1,2,…,Nrt−1,aj
Step 2: Update the estimated support set and the estimated atomic set E^t=E^t−1∪λt,At=At−1∪aλt;
Step 3: Estimate the signal x^t=AtTAt−1AtTb;
Step 4: Update the residual rt=b−Atx^t, update iteration counter t=t+1;
Judgement: If t≤K, continue iterating; if t>K, stop the iteration, and output the estimated support set E^ and the estimated signal x^.


## 2.3. Bat Algorithm

In 2010, Yang et al. simulated the characteristics of bats’ predatory behaviors and proposed BA [13, 18], which is a new Swarm Intelligence Optimization Algorithm.The first step is randomly initializing the bat positionsxi0 and the velocities vi0. Then, the BA updates the velocities and positions according to certain strategies. The velocities and positions of the i-th bat at t+1 iteration are updated as follows:(3)fi=fmin+fmax+fminβ,(4)vit+1=vit+xit+x∗⋅fi,(5)xit+1=xit+vit+1,where fi is the pulse frequency of the i-th bat, fi∈fmin,fmax, β∈0,1 is a random number. xit and xit+1 represent the positions of the i-th bat at the t and the t+1 iteration, respectively; vit and vit+1 are the velocities of the i-th bat at the t and the t+1 iteration, respectively; and x∗ denotes the global best position of the swarm.For each bat, the algorithm generates a random numberrand1∈0,1. If rand1>rit, the position of this bat is updated according to formula (6), where rit is the rate of pulse emission. The random walk is regarded as a procedure of the local search [13].(6)xit+1=xj∗+εLt,where xj∗ is a random solution in the best solution set composed of the best solution of each bat, ε∈−1,1 is a random number, and Lt is the mean of the pulse loudness in the t-th iteration.In addition, with the iterative number increasing, the pulse loudnessLi and the rate of pulse emission ri need to be updated. During prey searching, the pulse loudness Li will gradually decline, and the rate of pulse emission ri will also gradually increase. The updating formulas of Li and ri are found in formulas (7) and (8), respectively:(7)Lit+1=αLit,(8)rit+1=ri01−exp−γt,where α and γ are constants, 0<α<1, and 0<γ.The BA algorithm can be described in Algorithm2 [13].Algorithm 2: Bat Algorithm (BA).
Fitness function:
f
x
Initialization:
Randomly initialize the positions xi and velocities vi of the bat population.
Initialize the pulse frequency fi of each bat according to formula (3).
Initialize the rate of pulse emission ri0 and the pulse loudness Li0.
Calculate the fitness value of each bat and find the best solution in the population.
Set the stopping criterion σ.
Termination condition:
Iteration counter reaches tmax or the fitness value is smaller than σ.
Iteration:
Step 1: Update the bat velocity vit and position xit according to formulas (3)∼(5);
Step 2: If (rand1>rit), select a random solution in the best solutions set and update position of this bat by the random walk around the selected best solution as shown in formula (5);
Step 3: Calculate the fitness value of each bat’s new position;
Step 4: If (rand2>Litandfxit<fx∗), accept the new solution and adjust rit, Lit according to formulas (7) and (8);
Step 5: Find the new best solution x∗ and update the best solutions set;
Step 6: If termination condition is satisfied, stop the iteration; otherwise, continue.


## 3. A Bat-Inspired Sparse Recovery Algorithm

### 3.1. Fitness Function

If the sparsity levelK satisfying K<sparkA is known as a priori, problem (2) can be approximated as follows [11, 12]:(9)argminxAx−b2,s.t.x0≤K,where b is the measurement signal, and A∈Rm×n is the sensing matrix.From greedy pursuit algorithms, we can see that the signal can be recovered by using a two-step strategy. First, estimate the support set and then estimate the signal using the least square method.(10)xE^=AE^+b,xS−E^=0,where E^ denotes the estimated support set and S=1,2,⋯,n. The vector xE consists of the entries of x∈Rn indexed by i∈E⊆S. The matrix AE=ajj∈E⊆S is the submatrix of A, and aj is the j-th column of A. A+=ATA−1AT is the pseudoinverse of a matrix A, where T denotes the transpose operator.Once the support set is estimated accurately, the following equation must be satisfied.(11)AE^AE^+b−b2=0,where ⋅2 represents l2−norm.BecauseAE^AE^+b−b2≥0, we define the fitness function as(12)fE^=AE^AE^+b−b2.


### 3.2. Initialization

Suppose there areI bats, the position of the first bat is initialized as the estimated support set of the OMP. For the i-th (2≤i≤I) bat, its position is initialized as follows.Randomly chooseq=0.8×m elements from the set S=1,2,…,n and form the set Vi0; here, Vit is regarded as the velocity of the i-th bat at the t iteration. Then, the position of the i-th (2≤i≤I) bat is a set formed by K indices corresponding to the maximum absolute values in AVi0+⋅b. Other related parameters are initialized as follows: the pulsing frequency fi of each bat is initialized according to formula (3), where fmin=0.8andfmax=1. The initial emission rate ri0∈0,1. The pulse loudness Li0∈0,1, β is a random number between 0 and 1, and ε is a random number between −1 and 1, 0<α<1, and 0<γ.


### 3.3. Update Strategy

The main innovation of our algorithm is the update strategy. The iterative process combines the thought of BA and greedy pursuit algorithm. The update strategy can be divided into three parts.The first part is the update of velocity and position: a setQit is formed with β⋅fi⋅K elements selected randomly from the best solution E∗; we define Uit=Eit∪Qit. If Uit≤q, a set Git consists of q−Uit elements selected randomly from the set S−Uit, update Vit+1=Uit∪Git. If Uit>q, the set Git consists of q−K elements selected randomly from the set Qit−Eit, update Vit+1=Eit∪Git. The position of the i-th bat Eit+1 is a set composed of the indices corresponding to the K maximum absolute values in AVit+1+⋅b.The second part is the local search: generate a random numberrand1∈0,1; if rand1>rit, select a random solution Ej∗ from the best solutions set and update the position of this bat by randomly replacing ε⋅Lit⋅K elements of this solution with other elements in S−Ej∗. After updating all positions, calculate the new fitness fEit+1 of each bat. Specifically, in order to accelerate the convergence speed of this algorithm, all the new solutions will be accepted in the algorithm.The third part is the adjustment of the emission raterit and the pulse loudness Lit as formulas (7) and (8). If fEit+1<fEi∗, update the best solutions set Ei∗=Eit+1. If fEit+1<fE∗, update the global best solution E∗=Eit+1.


### 3.4. Stopping Criterion

IffE^<σ, the iteration is terminated; here, σ is the termination threshold. Moreover, in order to avoid too many iterations, set a maximum allowed iteration number N_max. If the iteration number is larger than N_max, the iteration is also terminated.The algorithmic process in Algorithm3.Algorithm 3: A Bat-Inspired Sparse Recovery Algorithm for CS.
Input: Sensing matrix A, measurement vector b, sparsity level K, termination threshold σ, maximum allowed iterative number N_max.
Output: The estimated signal x^
Fitness function:
f
E
^
=
A
E
^
A
E
^
+
b
−
b
2
Initialization: Initialize the positions and velocities of bats (Section 3.2).
Termination condition: If fE∗<σ or t≥N_max, terminated the iteration.
Iteration:
Step 1: A set Qit is formed with β⋅fi⋅K elements selected randomly from the best solution E∗, define Uit=Eit∪Qit.
Step 2: If Uit≤q, a set Git consists of q−Uit elements selected randomly from the set S−Uit, update Vit+1=Uit∪Git. If Uit>q, the set Git consists of q−K elements selected randomly from the set Qit−Eit, update Vit+1=Eit∪Git.
Step 3: The position of the i-th bat Eit+1 is a set composed of the indices corresponding to the K maximum absolute values in AVit+1+⋅b.
Step 4: If rand1>rit, select an optimal solution from the best solutions set and then replace ε⋅Lit⋅K elements of this solution with other elements in S−Ej∗.
Step 5: Calculate the new fitness fEit+1.
Step 6: If rand2>Lit, adjust rit and Lit.
Step 7: If fEit+1<fEi∗, update the best solutions set Ei∗=Eit+1. If fEit+1<fE∗, update the best solution E∗=Eit+1.
Step 8: If termination condition is satisfied, stop iteration; otherwise, continue.


### 3.5. Efficiency Analysis

The computational complexity of the algorithm is mainly dependent on the initialization and iteration.The algorithm initializes solutions using OMP, which has the complexityOKmn [19].In each iteration, the complexity is mainly dependent on three operations. The first is the multiplication of matrix, which has the computational complexityOmn. The second is sorting of an n-dimensional vector, which has the computational complexity Onlogn. The third is the least square, which has the computational complexity Om3. In summary, the computational complexity in the iteration phase is OIm3+Imn+Inlogn, where I is the population size.In general, the computational complexity of the proposed algorithm is similar to OMP when the signal sparsity is small. The complexity depends on the number of iterations when the sparsity is larger or when the measurement number is small. Moreover, this algorithm uses the idea of the BA, which combines the advantages of PSO, HS, SA, and other intelligence algorithms to optimize the convergence speed and search accuracy. It usually runs faster than other intelligence algorithms. OMP can generate an initial solution closed to the target, which accelerates the convergence of the BA. In addition, as swarm intelligence algorithm, it can run in parallel to reduce the running time [12, 20].


## 3.1. Fitness Function

If the sparsity levelK satisfying K<sparkA is known as a priori, problem (2) can be approximated as follows [11, 12]:(9)argminxAx−b2,s.t.x0≤K,where b is the measurement signal, and A∈Rm×n is the sensing matrix.From greedy pursuit algorithms, we can see that the signal can be recovered by using a two-step strategy. First, estimate the support set and then estimate the signal using the least square method.(10)xE^=AE^+b,xS−E^=0,where E^ denotes the estimated support set and S=1,2,⋯,n. The vector xE consists of the entries of x∈Rn indexed by i∈E⊆S. The matrix AE=ajj∈E⊆S is the submatrix of A, and aj is the j-th column of A. A+=ATA−1AT is the pseudoinverse of a matrix A, where T denotes the transpose operator.Once the support set is estimated accurately, the following equation must be satisfied.(11)AE^AE^+b−b2=0,where ⋅2 represents l2−norm.BecauseAE^AE^+b−b2≥0, we define the fitness function as(12)fE^=AE^AE^+b−b2.


## 3.2. Initialization

Suppose there areI bats, the position of the first bat is initialized as the estimated support set of the OMP. For the i-th (2≤i≤I) bat, its position is initialized as follows.Randomly chooseq=0.8×m elements from the set S=1,2,…,n and form the set Vi0; here, Vit is regarded as the velocity of the i-th bat at the t iteration. Then, the position of the i-th (2≤i≤I) bat is a set formed by K indices corresponding to the maximum absolute values in AVi0+⋅b. Other related parameters are initialized as follows: the pulsing frequency fi of each bat is initialized according to formula (3), where fmin=0.8andfmax=1. The initial emission rate ri0∈0,1. The pulse loudness Li0∈0,1, β is a random number between 0 and 1, and ε is a random number between −1 and 1, 0<α<1, and 0<γ.


## 3.3. Update Strategy

The main innovation of our algorithm is the update strategy. The iterative process combines the thought of BA and greedy pursuit algorithm. The update strategy can be divided into three parts.The first part is the update of velocity and position: a setQit is formed with β⋅fi⋅K elements selected randomly from the best solution E∗; we define Uit=Eit∪Qit. If Uit≤q, a set Git consists of q−Uit elements selected randomly from the set S−Uit, update Vit+1=Uit∪Git. If Uit>q, the set Git consists of q−K elements selected randomly from the set Qit−Eit, update Vit+1=Eit∪Git. The position of the i-th bat Eit+1 is a set composed of the indices corresponding to the K maximum absolute values in AVit+1+⋅b.The second part is the local search: generate a random numberrand1∈0,1; if rand1>rit, select a random solution Ej∗ from the best solutions set and update the position of this bat by randomly replacing ε⋅Lit⋅K elements of this solution with other elements in S−Ej∗. After updating all positions, calculate the new fitness fEit+1 of each bat. Specifically, in order to accelerate the convergence speed of this algorithm, all the new solutions will be accepted in the algorithm.The third part is the adjustment of the emission raterit and the pulse loudness Lit as formulas (7) and (8). If fEit+1<fEi∗, update the best solutions set Ei∗=Eit+1. If fEit+1<fE∗, update the global best solution E∗=Eit+1.


## 3.4. Stopping Criterion

IffE^<σ, the iteration is terminated; here, σ is the termination threshold. Moreover, in order to avoid too many iterations, set a maximum allowed iteration number N_max. If the iteration number is larger than N_max, the iteration is also terminated.The algorithmic process in Algorithm3.Algorithm 3: A Bat-Inspired Sparse Recovery Algorithm for CS.
Input: Sensing matrix A, measurement vector b, sparsity level K, termination threshold σ, maximum allowed iterative number N_max.
Output: The estimated signal x^
Fitness function:
f
E
^
=
A
E
^
A
E
^
+
b
−
b
2
Initialization: Initialize the positions and velocities of bats (Section 3.2).
Termination condition: If fE∗<σ or t≥N_max, terminated the iteration.
Iteration:
Step 1: A set Qit is formed with β⋅fi⋅K elements selected randomly from the best solution E∗, define Uit=Eit∪Qit.
Step 2: If Uit≤q, a set Git consists of q−Uit elements selected randomly from the set S−Uit, update Vit+1=Uit∪Git. If Uit>q, the set Git consists of q−K elements selected randomly from the set Qit−Eit, update Vit+1=Eit∪Git.
Step 3: The position of the i-th bat Eit+1 is a set composed of the indices corresponding to the K maximum absolute values in AVit+1+⋅b.
Step 4: If rand1>rit, select an optimal solution from the best solutions set and then replace ε⋅Lit⋅K elements of this solution with other elements in S−Ej∗.
Step 5: Calculate the new fitness fEit+1.
Step 6: If rand2>Lit, adjust rit and Lit.
Step 7: If fEit+1<fEi∗, update the best solutions set Ei∗=Eit+1. If fEit+1<fE∗, update the best solution E∗=Eit+1.
Step 8: If termination condition is satisfied, stop iteration; otherwise, continue.


## 3.5. Efficiency Analysis

The computational complexity of the algorithm is mainly dependent on the initialization and iteration.The algorithm initializes solutions using OMP, which has the complexityOKmn [19].In each iteration, the complexity is mainly dependent on three operations. The first is the multiplication of matrix, which has the computational complexityOmn. The second is sorting of an n-dimensional vector, which has the computational complexity Onlogn. The third is the least square, which has the computational complexity Om3. In summary, the computational complexity in the iteration phase is OIm3+Imn+Inlogn, where I is the population size.In general, the computational complexity of the proposed algorithm is similar to OMP when the signal sparsity is small. The complexity depends on the number of iterations when the sparsity is larger or when the measurement number is small. Moreover, this algorithm uses the idea of the BA, which combines the advantages of PSO, HS, SA, and other intelligence algorithms to optimize the convergence speed and search accuracy. It usually runs faster than other intelligence algorithms. OMP can generate an initial solution closed to the target, which accelerates the convergence of the BA. In addition, as swarm intelligence algorithm, it can run in parallel to reduce the running time [12, 20].


## 4. Experiments

In this section, we compare this algorithm with OMP, GOMP, SP, Fast Laplace, and PSO, where Fast Laplace is a Bayesian-based algorithm, and PSO is a swarm intelligence algorithm for CS. Then, we apply this algorithm to recover the real microseismic signals.Candes and Tao proved that the Gaussian random matrix [2, 21, 22] with independent identically distribution can be a universal measurement matrix. Therefore, the measurement matrices in these experiments are the Gaussian random matrix with the size of m×n. The number of iterations of the OMP, GOMP, and SP algorithms is K. The maximum number of iterations of PSO and the proposed algorithm are Nmax=100. The algorithm will stop when fE∗<σ or the number of iterations reaches N_max.The software and hardware environment of the experiment is as follows:Processor: Intel (R) Pentium (R) CPU G3220 @ 3.00 GHz 3.0Memory: 4.00 GOperating system: 64-bit Windows7Simulation software: MATLAB R2014a


### 4.1. Comparison with Other Algorithms

In the first experiment, we compare the recovery performance of different algorithms against the change of sparsity. We use a 128 × 256 random Gaussian matrix to sample a 256-length sparse signal. The sparsity level is set asK=35,40,45,50,55,60,64,70. The parameters are set as n=256, m=128, σ=1e−8, fmin=0.8,andfmax=1; the rate of pulse emission ri0 and the pulse loudness Li0 are set between 0 and 1; β=0.2, ε=0.2,α=0.9, and γ=0.9. Here, we compare the algorithm performance with the following two metrics: recovery error and exact recovery rate.Recovery error is a metric to evaluate the error between the original sparse signal and the recovered sparse signal. The recovery error is defined as follows (13):(13)Recovery  error=1CNT∑j=1CNTxj−x^j22xj22,where x is the original sparse signal, x^ is the recovered sparse signal, and CNT=100 is the cycle number to obtain the average performance. The simulation result of recovery error is shown in Figure 1(a). The exact recovery rate is obtained by calculating the number of times that the estimated values are exact. If x−x^2<10−6, this experiment is considered successful. In 100 experiments, if there are s successful experiments, the exact recovery rate is s/100. The exact recovery rates of different algorithms are shown in Figure 1(b).Figure 1
Comparison of the algorithm recovery performance against the change of sparsity withn=256, m=128, and the sparsity K=35,40,45,50,55,60,64,70. (a) Recovery error. (b) Exact recovery rate.
(a)
(b)It can be seen that when the sparsity is large, the recovery errors of three commonly used greedy pursuit algorithms such as OMP, GOMP, and SP increase rapidly, and the exact recovery rate decreases quickly. The recovery error of the Fast Laplace algorithm is lower than that of traditional greedy pursuit algorithms. However, experiments demonstrate that the exact recovery rates of the Bayesian algorithm are 0 because there exists no experiment satisfyingx−x^2<10−6. The swarm intelligence algorithms also perform well. The recovery error of the proposed algorithm is lower than 0.4, and the exact recovery rate is higher than other five algorithms when the sparsity is very large.In the second experiment, we compare the recovery performances against the change of measurement number. We use am×256 Gaussian matrix to sample a 256-length signal with sparsity 30. The measurement number is set as m=53,58,63,68,78,88,98,108,118, and other parameters are the same to the first experiment.Figures2(a) and 2(b) illustrate the recovery error and the exact recovery rate, respectively, of the six algorithms under different measurement numbers. As shown in the figure, the recovery of the proposed algorithm performs better than those commonly used greedy pursuit algorithms such as OMP, GOMP, and SP, when the number of measurement vectors is small. It also has more advantages compared with the PSO algorithm and the Fast Laplace algorithm.Figure 2
Comparison of the signal recovery performance against the change of measurement number withn=256, m=53,58,63,68,78,88,98,108,118 and sparsity K=30. (a) Recovery error. (b) Exact recovery rate.
(a)
(b)In the third experiment, we compare the recovery efficiency of different algorithms. Figure3(a) is the running time of the six algorithms with different sparsity. The sparsity level is set as K=1,10,20,30,40,50,55,60,64. Figure 3(b) is the running time of the six algorithms with different measurement number. The measurement number is taken as m=58,63,68,78,88,98,108,118,128. Other parameters are the same as those of the first experiment, and the average time of each algorithm in recovering the signal is counted. As shown in Figure 3, when the sparsity is small or the measurement number is large, the recovery efficiency of this algorithm is similar to the commonly used greedy pursuit algorithms such as OMP. Although the running time of the algorithm increases when the sparsity is large or the measurement number is small, the algorithm is more efficient than other swarm intelligence algorithms such as PSO. However, as a swarm intelligence algorithm, it can run in parallel to reduce the running time [12, 20].Figure 3
Comparison of average running times. (a) Running time of the six algorithms against the change of sparsity withn=256, m=128, and K=1,10,20,30,40,50,55,60,64. (b) Running time of the six algorithms against the change of measurement number with n=256, K=30, and m=53,58,63,68,78,88,98,108,118.
(a)
(b)


### 4.2. Application in Microseismic Signals

Microseismic monitoring is an effective technology for accident prediction in coal mines. Under the background of coal mine IoT, the collected data quantity by numerous microseismic nodes in the monitoring area are massive [23]. The effective compression of microseismic data can reduce the storage space and energy consumption. It plays an important role in improving the data transmission rate. Participants of the project collected three microseismic signals from a coal mine in Anhui Province, China. In this section, we use these signals as the original signal and apply the proposed algorithm to recover signal.We perform experiments on three different microseismic signals. The length of the signal isn=512, as shown in Figure 4. It can be seen that the microseismic signals are one-dimensional time-frequency signals. Microseismic signals are not sparse in the time domain; however, their representations under a certain basis are sparse or compressible.Figure 4
Original waveforms of the three microseismic signals: (a) signal 1, (b) signal 2, and (c) signal 3.
(a)
(b)
(c)In this experiment, the 2-layer db2 wavelet transform [24] is used as sparse basis. Wavelet transform uses a cluster of functions to approximate a signal which can be seen as a decomposition of the signal. It divides the signal into a low-frequency signal and multiple high-frequency signals. The low-frequency signal is the stationary part of the signal. The high-frequency signal indicates the details of the signal. Figure 5 shows the wavelet decomposition coefficients of three groups of microseismic signals.Figure 5
Wavelet decomposition coefficients of three microseismic signals: (a) signal 1, (b) signal 2, and (c) signal 3.
(a)
(b)
(c)The measurement matrix used in this experiment is a 256 × 512 random Gaussian matrix. The sparsity level isK=128 and σ=1e−5. Figure 6 shows the three groups of recovered signals. The recovery errors in the three groups of experiments are 0.1862, 0.1250, and 0.1512.Figure 6
Recovered waveforms of the three microseismic signals: (a) signal 1, (b) signal 2, and (c) signal 3.
(a)
(b)
(c)


## 4.1. Comparison with Other Algorithms

In the first experiment, we compare the recovery performance of different algorithms against the change of sparsity. We use a 128 × 256 random Gaussian matrix to sample a 256-length sparse signal. The sparsity level is set asK=35,40,45,50,55,60,64,70. The parameters are set as n=256, m=128, σ=1e−8, fmin=0.8,andfmax=1; the rate of pulse emission ri0 and the pulse loudness Li0 are set between 0 and 1; β=0.2, ε=0.2,α=0.9, and γ=0.9. Here, we compare the algorithm performance with the following two metrics: recovery error and exact recovery rate.Recovery error is a metric to evaluate the error between the original sparse signal and the recovered sparse signal. The recovery error is defined as follows (13):(13)Recovery  error=1CNT∑j=1CNTxj−x^j22xj22,where x is the original sparse signal, x^ is the recovered sparse signal, and CNT=100 is the cycle number to obtain the average performance. The simulation result of recovery error is shown in Figure 1(a). The exact recovery rate is obtained by calculating the number of times that the estimated values are exact. If x−x^2<10−6, this experiment is considered successful. In 100 experiments, if there are s successful experiments, the exact recovery rate is s/100. The exact recovery rates of different algorithms are shown in Figure 1(b).Figure 1
Comparison of the algorithm recovery performance against the change of sparsity withn=256, m=128, and the sparsity K=35,40,45,50,55,60,64,70. (a) Recovery error. (b) Exact recovery rate.
(a)
(b)It can be seen that when the sparsity is large, the recovery errors of three commonly used greedy pursuit algorithms such as OMP, GOMP, and SP increase rapidly, and the exact recovery rate decreases quickly. The recovery error of the Fast Laplace algorithm is lower than that of traditional greedy pursuit algorithms. However, experiments demonstrate that the exact recovery rates of the Bayesian algorithm are 0 because there exists no experiment satisfyingx−x^2<10−6. The swarm intelligence algorithms also perform well. The recovery error of the proposed algorithm is lower than 0.4, and the exact recovery rate is higher than other five algorithms when the sparsity is very large.In the second experiment, we compare the recovery performances against the change of measurement number. We use am×256 Gaussian matrix to sample a 256-length signal with sparsity 30. The measurement number is set as m=53,58,63,68,78,88,98,108,118, and other parameters are the same to the first experiment.Figures2(a) and 2(b) illustrate the recovery error and the exact recovery rate, respectively, of the six algorithms under different measurement numbers. As shown in the figure, the recovery of the proposed algorithm performs better than those commonly used greedy pursuit algorithms such as OMP, GOMP, and SP, when the number of measurement vectors is small. It also has more advantages compared with the PSO algorithm and the Fast Laplace algorithm.Figure 2
Comparison of the signal recovery performance against the change of measurement number withn=256, m=53,58,63,68,78,88,98,108,118 and sparsity K=30. (a) Recovery error. (b) Exact recovery rate.
(a)
(b)In the third experiment, we compare the recovery efficiency of different algorithms. Figure3(a) is the running time of the six algorithms with different sparsity. The sparsity level is set as K=1,10,20,30,40,50,55,60,64. Figure 3(b) is the running time of the six algorithms with different measurement number. The measurement number is taken as m=58,63,68,78,88,98,108,118,128. Other parameters are the same as those of the first experiment, and the average time of each algorithm in recovering the signal is counted. As shown in Figure 3, when the sparsity is small or the measurement number is large, the recovery efficiency of this algorithm is similar to the commonly used greedy pursuit algorithms such as OMP. Although the running time of the algorithm increases when the sparsity is large or the measurement number is small, the algorithm is more efficient than other swarm intelligence algorithms such as PSO. However, as a swarm intelligence algorithm, it can run in parallel to reduce the running time [12, 20].Figure 3
Comparison of average running times. (a) Running time of the six algorithms against the change of sparsity withn=256, m=128, and K=1,10,20,30,40,50,55,60,64. (b) Running time of the six algorithms against the change of measurement number with n=256, K=30, and m=53,58,63,68,78,88,98,108,118.
(a)
(b)


## 4.2. Application in Microseismic Signals

Microseismic monitoring is an effective technology for accident prediction in coal mines. Under the background of coal mine IoT, the collected data quantity by numerous microseismic nodes in the monitoring area are massive [23]. The effective compression of microseismic data can reduce the storage space and energy consumption. It plays an important role in improving the data transmission rate. Participants of the project collected three microseismic signals from a coal mine in Anhui Province, China. In this section, we use these signals as the original signal and apply the proposed algorithm to recover signal.We perform experiments on three different microseismic signals. The length of the signal isn=512, as shown in Figure 4. It can be seen that the microseismic signals are one-dimensional time-frequency signals. Microseismic signals are not sparse in the time domain; however, their representations under a certain basis are sparse or compressible.Figure 4
Original waveforms of the three microseismic signals: (a) signal 1, (b) signal 2, and (c) signal 3.
(a)
(b)
(c)In this experiment, the 2-layer db2 wavelet transform [24] is used as sparse basis. Wavelet transform uses a cluster of functions to approximate a signal which can be seen as a decomposition of the signal. It divides the signal into a low-frequency signal and multiple high-frequency signals. The low-frequency signal is the stationary part of the signal. The high-frequency signal indicates the details of the signal. Figure 5 shows the wavelet decomposition coefficients of three groups of microseismic signals.Figure 5
Wavelet decomposition coefficients of three microseismic signals: (a) signal 1, (b) signal 2, and (c) signal 3.
(a)
(b)
(c)The measurement matrix used in this experiment is a 256 × 512 random Gaussian matrix. The sparsity level isK=128 and σ=1e−5. Figure 6 shows the three groups of recovered signals. The recovery errors in the three groups of experiments are 0.1862, 0.1250, and 0.1512.Figure 6
Recovered waveforms of the three microseismic signals: (a) signal 1, (b) signal 2, and (c) signal 3.
(a)
(b)
(c)


## 5. Conclusions

In this paper, an Intelligence Sparse Recovery Algorithm is proposed for CS inspired by the BA and the pruning technique of SP. The algorithm inherits the global search ability and the local search ability of the BA, and we use OMP to obtain a better initial solution to accelerate the convergence. Compared with the traditional pursuit algorithm, the algorithm in this paper has better recovery performance under large sparsity and small measurement. At the same time, it runs faster than the intelligent optimization algorithm such as PSO. The proposed algorithm can be applied to the microseismic monitoring system, which can recover the signal well. In future research, performance and efficiency of the algorithm will be more rigorously analyzed.



---

*Source: 1365747-2018-10-29.xml*