A Bat-Inspired Sparse Recovery Algorithm for Compressed Sensing

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.


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. is 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][2][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.
e signal recovery problem in CS is a nonconvex combinatorial optimization problem. It can be solved by three types of algorithms. e 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]. ey 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. is kind of algorithm is effective; however, they have weak recovery performance.
e second kind of algorithm is a convex optimization algorithm, which recovers signals by converting nonconvex problems to convex problems. e 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. e 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. e 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]. ese 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. e 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. e BA and the pruning technique are combined to update populations. e proposed algorithm outperforms the traditional greedy pursuit algorithms. Moreover, it runs faster than other swarm intelligence algorithms such as PSO. e remainder of this paper is organized as follows: Section 2 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. e 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. ‖ · ‖ p (p � 0, 1, 2, . . .) represents l p − norm. For the vector x ∈ R n , x is the estimate of x. In a matrix A ∈ R m×n , a j is the j-th column of A. E � 1 ≤ i ≤ n|x i ≠ 0 , which denotes the positions of nonzero elements in x, and it is called support set.
where T denotes the transpose operator.

Compressed Sensing.
Supposing an n-dimensional signal y � (y 1 , y 2 , . . . , y n ) T can be represented on a basis Ψ ∈ R n×n as y � Ψx, where x � (x 1 , x 2 , . . . , x n ) T , if there are only K nonzero elements in x � (x 1 , x 2 , . . . , x n ) T , 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 ≤ n|x i ≠ 0 . e signal y can be sampled by projection onto the measurement matrix x is the sparse signal, Φ is the measurement matrix, and Ψ is the sparse basis. A ∈ R m×n is the sensing matrix. It can be seen that the data quantity is greatly reduced because b � Φy is an underdetermined system of equations.
ere are infinitely many solutions to recover y from b. However, if y is sparse, it can be recovered from b by solving the following l 0 minimization problem: where ‖ · ‖ 0 represents l 0 − norm.

Orthogonal Matching Pursuit (OMP).
In the proposed algorithm, we use OMP to determine the initial solution. e OMP algorithm can be stated in Algorithm 1 [4].

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. e first step is randomly initializing the bat positions x i (0) and the velocities v i (0). en, the BA updates the velocities and positions according to certain strategies. e velocities and positions of the i-th bat at t + 1 iteration are updated as follows: where f i is the pulse frequency of the i-th bat, is a random number. x i (t) and x i (t + 1) represent the positions of the i-th bat at the t and the t + 1 iteration, respectively; v i (t) and v i (t + 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 number rand1 ∈ [0, 1]. If rand1 > r i (t), the position of this bat is updated according to formula (6), where r i (t) is the rate of pulse emission. e random walk is regarded as a procedure of the local search [13].
where x * j is a random solution in the best solution set composed of the best solution of each bat, ε ∈ [−1, 1] is a random number, and L(t) is the mean of the pulse loudness in the t-th iteration.
In addition, with the iterative number increasing, the pulse loudness L i and the rate of pulse emission r i need to be updated. During prey searching, the pulse loudness L i will gradually decline, and the rate of pulse emission r i will also gradually increase. e updating formulas of L i and r i are found in formulas (7) and (8), respectively: where α and c are constants, 0 < α < 1, and 0 < c. e BA algorithm can be described in Algorithm 2 [13].

Fitness Function.
If the sparsity level K satisfying K < spark(A) is known as a priori, problem (2) can be approximated as follows [11,12]: where b is the measurement signal, and A ∈ R m×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.
where E denotes the estimated support set and Once the support set is estimated accurately, the following equation must be satisfied.
where ‖ · ‖ 2 represents l 2 − norm. Because 3.2. Initialization. Suppose there are I 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.
is regarded as the velocity of the i-th bat at the t iteration. en, the position of the i-th (2 ≤ i ≤ I) bat is a set formed by K indices corresponding to the maximum absolute values in Other related parameters are initialized as follows: the pulsing frequency f i of each bat is initialized according to formula (3), where f min � 0.8andf max � 1. e initial emission rate r i (0) ∈ [0, 1]. e pulse loudness L i (0) ∈ [0, 1], β is a random number between 0 and 1, and ε is a random number between −1 and 1, 0 < α < 1, and 0 < c.

Update Strategy.
e main innovation of our algorithm is the update strategy. e iterative process combines the The t-th(1 ≤ t < K)iteration: Step 1: Find λ t � argmax j�1,2,...,N |〈r t−1 , a j 〉| Step 2: Update the estimated support set and the estimated atomic Step 3: Estimate the signal  Iteration counter reaches t max or the fitness value is smaller than σ. Iteration: Step 1: Update the bat velocity v i (t) and position x i (t) according to formulas (3)∼(5); Step 2: If (rand1 > r i (t)), 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: , accept the new solution and adjust r i (t), L i (t) 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.
Computational Intelligence and Neuroscience 3 thought of BA and greedy pursuit algorithm. e update strategy can be divided into three parts. e first part is the update of velocity and position: a set Q i (t) is formed with β · f i · K elements selected randomly from the best solution E * ; we define . e position of the i-th bat E i (t + 1) is a set composed of the indices corresponding to the K maximum absolute values in e second part is the local search: generate a random number rand1 ∈ [0, 1]; if rand1 > r i (t), select a random solution E * j from the best solutions set and update the position of this bat by randomly replacing ε · L i (t) · K elements of this solution with other elements in S − E * j . After updating all positions, calculate the new fitness f(E i (t + 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. e third part is the adjustment of the emission rate r i (t) and the pulse loudness L i (t) as formulas (7) and (8)

Stopping
Criterion. If f(E) < σ, 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. e algorithmic process in Algorithm 3.

Efficiency
Analysis. e computational complexity of the algorithm is mainly dependent on the initialization and iteration. e algorithm initializes solutions using OMP, which has the complexity O(Kmn) [19].
In each iteration, the complexity is mainly dependent on three operations. e first is the multiplication of matrix, which has the computational complexity O(mn). e second is sorting of an n-dimensional vector, which has the computational complexity O(nlogn). e third is the least square, which has the computational complexity O(m 3 ). In summary, the computational complexity in the iteration phase is O(Im 3 + 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. e 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].

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. en, 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. erefore, the measurement matrices in these experiments are the Gaussian random matrix with the size of m × n. e number of iterations of the OMP, GOMP, and SP algorithms is K. e maximum number of iterations of PSO and the proposed algorithm are N max � 100. e algorithm will stop when f(E * ) < σ or the number of iterations reaches N max. e software and hardware environment of the experiment is as follows: Processor: Intel (R) Pentium (R) CPU G3220 @ 3.00 GHz 3.0 Memory: 4.00 G Operating system: 64-bit Windows7 Simulation software: MATLAB R2014a

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.
Recovery error is a metric to evaluate the error between the original sparse signal and the recovered sparse signal. e recovery error is defined as follows (13): 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. e simulation result of recovery error is shown in Figure 1(a). e 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). e exact recovery rates of different algorithms are shown in Figure 1(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. e 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 satisfying x − x 2 < 10 −6 . e swarm intelligence algorithms also perform well. e 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 a m × 256 Gaussian matrix to sample a 256-length signal with sparsity 30. e 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. Figures 2(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. Termination condition: If f(E * ) < σ or t ≥ N max, terminated the iteration. Iteration: Step 1: A set Q i (t) is formed with β · f i · K elements selected randomly from the best solution E * , define Step 2: Step 3: e position of the i-th bat E i (t + 1) is a set composed of the indices corresponding to the K maximum absolute values in Step 4: If rand1 > r i (t), select an optimal solution from the best solutions set and then replace ε · L i (t) · K elements of this solution with other elements in S − E * j .
Step 6: If rand2 > L i (t), adjust r i (t) and L i (t). Step Step 8: If termination condition is satisfied, stop iteration; otherwise, continue.

Computational Intelligence and Neuroscience
In the third experiment, we compare the recovery efficiency of different algorithms. Figure 3(a) is the running time of the six algorithms with different sparsity. e 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. e 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].

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 Computational Intelligence and Neuroscience the monitoring area are massive [23]. e 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. e length of the signal is n � 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.
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 lowfrequency signal and multiple high-frequency signals. e low-frequency signal is the stationary part of the signal. e high-frequency signal indicates the details of the signal. Figure 5 shows the wavelet decomposition coefficients of three groups of microseismic signals. e measurement matrix used in this experiment is a 256 × 512 random Gaussian matrix. e sparsity level is K � 128 and σ � 1e − 5. Figure 6 shows the three groups of recovered signals. e recovery errors in the three groups of experiments are 0.1862, 0.1250, and 0.1512.

Conclusions
In this paper, an Intelligence Sparse Recovery Algorithm is proposed for CS inspired by the BA and the pruning technique of SP. e 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. e 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.
Data Availability e data of microseismic signals used to support the finding of this study have not been made available for the time being because the data are obtained through cooperation between us and the Coal Mine Group, and it must be approved by both parties before it can be disclosed.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.