Improved Stochastic Gradient Matching Pursuit Algorithm Based on the Soft-Thresholds Selection

1e preliminary atom set exits redundant atoms in the stochastic gradient matching pursuit algorithm, which affects the accuracy of the signal reconstruction and increases the computational complexity. To overcome the problem, an improved method is proposed. Firstly, a limited soft-threshold selection strategy is used to select the new atoms from the preliminary atom set, to reduce the redundancy of the preliminary atom set. Secondly, before finding the least squares solution of the residual, it is determined whether the number of columns of the measurement matrix is smaller than the number of rows. If the condition is satisfied, the least squares solution is calculated; otherwise, the loop is exited. Finally, if the length of the candidate atomic index set is less than the sparsity level, the current candidate atom index set is the support atom set. If the condition is not satisfied, the support atom index set is determined by the least squares solution. Simulation results indicate that the proposed method is better than other methods in terms of the reconstruction probability and shorter running time than the stochastic gradient matching pursuit algorithm.


Introduction
Compressed sensing (CS) [1,2] theory has three core problems: sparse representation of signals, and design of the measurement matrix, and design of reconstruction algorithms.e reconstruction algorithm is directly related to the accuracy of the recovery signal and the convergence rate of the algorithm, which determines whether the theory is feasible.e recovery algorithm is described as the recovery of the high-dimensional original signal from the lowdimensional measurement vectors.At present, many reconstruction algorithms have been proposed to recover the signal, which include convex optimization methods, combinational methods, and greedy pursuit methods.Convex optimization methods approximate the signal by transforming nonconvex problems into convex ones, which include basis pursuit (BP) [3], the gradient projection for sparse reconstruction (GPSR) [4] algorithm, iterative threshold (IT) [5], interior-point method [6], Bergman iteration (BT) [7], and total variation (TV) [8].Although convex optimization algorithms have fewer observations, their higher complexity is not suitable for practical applications.
Combinational methods include Fourier sampling [9,10], chaining pursuit (CP) [11], heavy hitters on steroids pursuit (HHSP) [12].Combinational methods are low in complexity; however, the accuracy of the recovery signal is not as good as convex optimization algorithms.Greedy algorithms have many advantages, such as simple structure, fast convergence, and low complexity, becoming the first choice for the recovery algorithm.
To date, many greedy pursuit algorithms have been proposed, including the first, the matching pursuit (MP) [13] algorithm.Based on this algorithm, the orthogonal matching pursuit (OMP) [14] algorithm was proposed to optimize MP via orthogonalization of the atoms of the support set.However, the OMP algorithm only selects one of the atoms of the support set at each round of iteration, and the efficiency is lower.Successors include the regularized OMP (ROMP) [15], subspace pursuit (SP) [16] algorithm, compressive sampling matching pursuit (CoSaMP) [17,18] algorithm, and stagewise OMP (StOMP) [19] algorithm.e ROMP algorithm realizes the effective selection of the atom via the regularized rule, to improve the speed of the OMP.
e SP and CoSaMP algorithms are similar.Both are proposed with the backtracking strategy.
e difference between the SP and CoSaMP is the number of atoms that are selected from measurement matrix to compose preliminary atomic set in each round of iteration.e SP selects s atoms, while the CoSaMP selects 2s atoms.e StOMP algorithm selects multiple atoms or columns of the measurement matrix in each round of iteration via a threshold parameter.
e greedy algorithms mentioned above all require the sparsity information of the signal, and the choice of the type of atoms is inflexible, which affects the convergence speed, robustness, and reconstruction performance of the algorithms.
Because the traditional greedy pursuit algorithm needs to compute the inverse matrix of the sensing matrix, this process requires a significant amount of computation time and storage space, resulting in lower reconstruction probability.In recent years, some research workers have proposed the fast version of the greedy algorithm, which avoids computation of the inverse matrix [20].e GP algorithm uses the gradient idea of the unconstrained optimization method to replace the computation of the inverse matrix and to reduce the computational complexity and storage space of the traditional greedy algorithms [21].However, the GP algorithm has a slow convergence and lower efficiency.In addition, the GP algorithm cannot solve the problem of large-scale data recovery.To improve those problems, the conjugate gradient pursuit (CGP) [22] algorithm and approximate conjugate gradient pursuit (ACGP) [23] algorithm are proposed, respectively.Although CGP and ACGP algorithms effectively reduce the computational complexity and storage space of traditional greedy algorithms, the reconstruction performance still needs to be improved.Based on the GP algorithm, the stagewise weak selection (SWGP) [24] algorithm was proposed to improve the reconstruction performance and convergence speed of the GP algorithm, via introduction of the conjugate direction and weak selection strategy.Motivated by the stochastic gradient descent methods, the stochastic gradient matching pursuit (StoG-radMP) algorithm [25] was recently proposed for the optimization problem with sparsity constraints.
e atomic selection method of the fixed number in the StoGradMP algorithm (that is, selecting 2s atoms at the preliminary stage of each iteration) will lead to the redundant atoms of the preliminary atomic set.When joined with the candidate atomic set, this will reduce the accuracy of the least squares solution and the inaccuracy of the support atomic set estimation, which affects the precision of the signal reconstruction and increases the computational complexity of the algorithm.In this study, we use the limited softthreshold selection strategy to realize the second selection of the preliminary atom set after the preliminary stage, which will improve the reconstruction accuracy.e combination with reliability verification conditions ensures the correctness and effectiveness of the proposed method.

Compressed Sensing Theory
e recent work in the area of the compressed sensing [26] demonstrated that it was possible to algorithmically recover sparse (and, more generally compressible) signals from incomplete observations.e simplest model is an n-dimensional signal x with a small number of nonzero entries under no noise conditions: where where θ � Ψ T x is the n × 1 column vector of projection coefficients, θ i � 〈x, ψ i 〉 � ψ T i x is the projection coefficient.Obviously, x and θ are the equivalent representations of the same signals in the different domains.at is, x and θ are the signals in the time domain and Ψ domain, respectively.When x is s-sparse, let Ψ � I, I is the unit matrix, and x � θ.We use a measurement matrix Φ ∈ R m×n (m ≤ n) to make a linear measurement of the projection coefficient vector θ, obtaining an observation vector y ∈ R m×1 , expressed as where Equation ( 3) is a linear projection of the original signals x on Φ.Note that the measurement process is nonadaptive; that is, Φ does not depend in any way on the signal x.Clearly, the dimension of y is much lower than the dimension of x. at is, this problem is an underdetermined problem; Equation (3) has infinitely many solutions.It is difficult to recover the projection coefficient vector θ directly from the observation vector y.However, the original signal x is s-sparse and Φ satisfies certain conditions in Equation (3); x can be recovered by solving the l 0 -minimization problem: where ‖.‖ 0 is the l 0 -norm of the vector, representing the number of nonzero entries.Candes and Tao demonstrated that if the s-sparse signal x is to be accurately recovered, the number of measurements m (or the dimension of observation vector y) must satisfy m � O(s ln(n)), and the measurement matrix must satisfy the restricted isometric property (RIP) [27,28].When x is not s-sparse in the time domain, the signal recovery process cannot be directly used in the Equation ( 4).
e signal x can be sparse representations on the sparse basis matrix Ψ. Combining Equations ( 2) and (3), we obtain 2 Journal of Electrical and Computer Engineering where  Φ � ΦΨ T ∈ R m×n is the sensing matrix [29].According to [30], the equivalent condition of the RIP is that the measurement matrix Φ is not correlated with the sparse basis matrix Ψ.Note that if the sensing matrix  Φ also satisfies the RIP, the recovery signal (or projection coefficient) θ can be obtained by solving an optimal l 0 -norm problem similar to Equation (4): From Equation ( 5), we see that Ψ is fixed, so that Φ also satisfies the RIP.
e measurement matrix Φ must meet certain conditions.It is shown in [31] that when the measurement matrix is a Gaussian random matrix, the sensing matrix  Φ can satisfy the RIP condition with a larger probability.
However in most scenarios, the signal contains noise.In this case, the measurement process can be expressed as where the m × n matrix  Φ is the sensing matrix and e is the m-dimensional noise vector.θ is the s-sparse signal in the Ψ domain.In this study, for simplicity, we assume that the signal x is s-sparse, that is, x � θ and Φ �  Φ. erefore, Equation ( 7) can be written as y � Φx + e, and we minimize the following formula to recover x: where s controls the sparsity of the solution to Equation (9).
To analyze (8), we combine Equation ( 2), where θ i are the projection coefficients of signal x. x is considered sparse with respect to Ψ if s is relatively small compared to the ambient dimension n. erefore, we can express the optimization (9) in the form as follows: where f i (x), x ∈ R n , is a smooth function that is nonconvex; ‖x‖ 0,Ψ is defined as the norm that captures the sparsity level of x.In particular, ‖x‖ 0,Ψ is the smallest number of atoms in Ψ such that x can be represented as For sparse signal recovery, the set Ψ consists of n basic vectors, each of size n in the Euclidean space. is problem can be seen as a special case of Equation (10), with f i (x) � (y i − < φ i , x >) 2 and M � m.We decompose the vector y m×1 into nonoverlapping vectors y b i with a size of b and denote Φ b i as Φ b i ×n .Φ b i ×n is the submatrix of measurement matrix Φ.According to Equation ( 9), the objective function is We then denote F(x) as where M � m/b, which is a positive integer.We treat each function and each function f i (x) represents a collection (or block) of observations of size b.Here, we spilt the function F(x) into multiple subfunctions f i (x) or block the measurement matrix Φ into the submatrix Φ b i , which will be beneficial for the calculation of the gradient.

StoGradMP Algorithm
CoSaMP [17] is a popular algorithm for recovering a sparse signal from its linear measurements.e idea of the CoSaMP algorithm is generalized to prove the GradMP algorithm that solves a broader class of sparsity-constrained problems.In this section, we describe the stochastic version of the GradMP, the StoGradMP [25] algorithm, where at each iteration, only the evaluation of the gradient of a function f i is required.e StoGradMP algorithm is described in Algorithm 1, which consists of the following steps in each iteration: Randomize.e measurement matrix Φ is subject to random block processing, to determine the region of the submatrix Φ b i ×n .en, according to Equation (11), calculates the subfunction f i k (x k ).Proxy.Compute the gradient of f i k (x k ).Here, the gradient is an n × 1 column vector.Identify.Select the column indexes of the submatrix Φ b i ×n corresponding to the maximum 2s components in the gradient vector, forming a preliminary index set S k .Merge.Merge the preliminary index set S k and the support index set F k−1 of the previous iteration to form a candidate index set C k . Estimation.
e estimation of signal b k by a suboptimization method is determined, which is the least squares method.Generally, b k is the transition signal.Prune.Select the column indexes of the measurement matrix Φ corresponding to the maximum s components in the signal estimation vector b k that forms a support index set F. Update.Update the signal estimation x � x k .Check.When the residual is less than the tolerance error of the proposed method iteration, stop the iteration.If the loop index k is greater than the maximum number of iterations s, the proposed method ends and the approximation of signal x � x k is output.If the iteration ends, the condition is not satisfied.Otherwise, continue the iteration until the halting condition is met.

Improved StoGradMP Algorithm
e StoGradMP algorithm selects 2s atoms in each iteration, where s is a fixed number.e choice of atoms is inflexible and will increase the redundancy of the preliminary atom set. is affects the reconstruction and computational speed.To solve this problem, we use the limited soft selection strategy to select atoms from the preliminary atom set.
Firstly, according to Equation (11), we randomly block the measurement matrix to obtain a stochastic block matrix (or submatrix) Φ block×n , which can be expressed as where M is the number of the block measurement matrix, M � floor(m/b), and b � min(s, m) is the number of rows of the block matrix.floor(.)and rand(.)generate a maximum integer smaller than m/b and uniformly distributed pseudorandom numbers, respectively.ceil(.)obtains a minimum integer greater than rand * M. block is the row index of the measurement matrix Φ.From Equation ( 12), we know that the area of the block matrix is randomly determined.For simplicity, we express Φ block×n and y block as Φ b i and y b i , respectively.Next, we compute the subfunction f i (x) as follows: where k is the number of iterations, i ∈ [M].Φ b i is the i-th submatrix of the measurement matrix Φ, which is stochastically determined and of size b i × n. e subfunctions f i (x) are also stochastically determined, and Here, f i (x) is the i-th subfunction of F(x).
After determining the area of submatrix Φ b i and the subfunction f i (x), we calculate the gradient of the subfunction, which is expressed as where r k is the gradient of the subfunction f i (x) at the k-th iteration and the gradient is an n × 1 column vector.∇(.) and Φ b i T represent the derivative of the subfunction f i (x) and the transpose matrix of submeasurement matrix Φ b i , respectively.x k−1 is the approximation of sparse signal x at the k − 1-th iteration.
In the atomic preliminary selection stage of the stochastic algorithm, we obtain the gradient vector of the function through the derivation of randomly determined subfunctions, and then select the 2s largest gradient values from the gradient vector r k , thereby determining the column index of the measurement matrix corresponding to the values of the gradient vector, and form a preliminary atom index set S k as where |r k | is the absolute value of the gradient vector r k , k is the number of iterations, and max(|r k |, 2s) is the atomic (or column) index of the measurement matrix corresponding to the maximal value from r k that forms a preliminary atom index set S k .
After the preliminary atomic stage, the preliminary atomic set exits the redundant atoms, which will reduce the accuracy of the least square solutions of signal and the inaccuracy of the support atomic set estimation. is will eventually affect the precision of signal reconstruction, and increase the computational complexity of the algorithm.To improve the reconstruction accuracy, we use the limited soft-threshold selection strategy to realize the second selection of the preliminary atom set.
In the atomic index set S k , the atoms corresponding to the index are expressed as Γ, Γ � [Γ where C k , F k−1 , and S * k represent the candidate atomic index set, the support atomic index set, and the preliminary atomic index set, respectively.Φ C k is the candidate atomic set at the k-th iteration.Φ F k−1 is the support atomic set at the previous iteration.Φ S * k is the new atomic set corresponding to atomic index set S * k at the k-th iteration.Before solving the suboptimization method, we must ensure that the number of rows is greater than the number of columns in the candidate atomic matrix Φ C k ; that is, Φ C k is a full column-rank matrix.erefore, we provide the first reliability verification condition; that is, if then where A k is the candidate atomic set (or matrix) at the k-th iteration.If the condition is not satisfied, namely, the number of rows is smaller than the number of columns, the matrix is not inversed.If this occurs, we exit the loop.Let b k � 0. Next, solving the least squares solution of the sparse approximation signal  x, it can be expressed as where b k is the estimation vector of the sparse signal x at the k-th iteration.is is known as the transition signal.A + C k is the inverse matrix of the candidate atomic set (or matrix), and y is the observing vector.
Since the soft-threshold selection strategy is used to complete the second selection of the preliminary atomic set, the size of the candidate atomic index set C k may be less than s.erefore, to ensure the correctness and effectiveness of the proposed method, we provide the second reliability verification condition, that is, if then where s is the size of sparsity of the signal.length(C k ) is the size of the candidate atomic index set then where |b k | is the absolute value of the transition signal at the k-th iteration.max(|b k |, s) determines the atomic (or column) index of the measurement matrix corresponding to the maximum value from b k and the former support atomic index set F. Next, we update the approximation of sparse signal x, which is expressed as where x k is the approximation of the signal at the k-th iteration, and b kF is the recovery signal corresponding to the support atomic index set F. Finally, we check the iteration stopping condition.If where y − Φx k represents the residual at the k-th iteration and tol is the tolerance error of the algorithm iteration.max Iter is the maximum number of iteration, with special value 500 * M in this study.done � 1 represents the algorithm stops and outputs the signal approximation  x � x k .If the iteration stopping criteria is not satisfied, then the iteration is continued until the condition is satisfied.e entire procedure is shown in Algorithm 2.

Results and Discussion
In this section, we used the signal with s-sparsity as the original signal.e measurement matrix is randomly generated with a Gaussian distribution.All performance is an average after running the simulation 100 times using a computer with a quadcore, 64-bit processor, and 4G memory.
In Figure 1, we compare the reconstruction probability of different sparsity of the proposed algorithm with different measurements in different threshold conditions.We set the sparsity level set as s ∈ [8,12,16,20].
e threshold parameter set is t ∈ [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0].From Figures 1(a)-1(e), we can see that when the size of the threshold parameter is 0.1, 0.2, 0.3, 0.4, and 0.5, the reconstruction probability of the proposed method is very Journal of Electrical and Computer Engineering close, with almost no difference.In Figure 1(f ), when the size of the threshold parameter is 0.6, the proposed method can complete signal reconstruction with less measurements, compared with other thresholds, when the size of sparsity is 16 and 20.In Figures 1(g)-1(j), we see that when the size of the threshold is 0.7, 0.8, 0.9, and 1.0, the proposed method requires more measurements to complete the signal reconstruction under the same sparsity level.
Figure 2 compares the reconstruction probability of different measurements of the proposed method with different thresholds under the same sparsity level.We set the sparsity set as s � [8,12,16,20].e threshold parameter set is consistent with the threshold parameter set in Figure 1.From Figure 2, we see that different threshold parameters have certain influences on the reconstruction probability of the signal for different sparsity levels.In Figures 2(a) and 2(b), we can see that when the sparsity levels are s � 8 and s � 12, the reconstruction probability of all thresholds conditions is very close, except for the thresholds 0.8, 0.9, and 1.0.Meanwhile, from Figures 2(c) and 2(d), we can see that when the sparsity levels are s � 16 or s � 20 and the sizes of the threshold are 0.5, 0.6, and 0.7, the reconstruction performance of the proposed method is better than the reconstruction probability of the other thresholds for different measurements.erefore, from Figure (2), it demonstrates that the soft-threshold strategy is more advantageous for larger sparsity.In particular, when the sparsity level s � 20 and threshold t � 0.6, the reconstruction probability of the proposed method is better than the reconstruction probability of other threshold conditions.Based on the analysis of Figures 1 and 2, we can see that that when the threshold is t � 0.6, the proposed method has better performance.erefore, in the following simulations, we set the threshold as 0.6.
In Figure 3, we compared the average runtime of different thresholds of the proposed algorithm with different measurements.From Figures 1 and 2, we see that the reconstruction probability of the proposed method is 100% when the number of the measurements is greater than or equal to 145.In particular, when the sparsity level is s � 20 and the threshold is 0.6, the reconstruction probability is better than the reconstruction probability of other thresholds.
erefore, we set the range of the number of the measurements as [145 200] in the simulation of Figure 3. From this, we see that the proposed algorithm with t � 0.8 has the shortest runtime, and the next shortest are the proposed algorithm with t � 0.7 and the proposed algorithm with t � 0.6.
e proposed method with t � 0.9 has the longest runtime.
is means that the selection of the threshold parameter t is important to the runtime of the proposed method.
From the aspects of reconstruction probability and runtime of the proposed method, we conclude that when the size of threshold parameter is t � 0.6 and the sparsity level is s � 20, reconstruction performance of the proposed method    Journal of Electrical and Computer Engineering algorithms is 100% when the sparsity is equal to 20 and the number of measurements is greater than or equal to 120.erefore, we set the range of the number of the measurements as [120 250] in the simulation of Figure 7, and the sparsity level is equal to 20.We see that the GradMP algorithm has the lowest runtime, and the next lowest are the StoIHT and proposed algorithms.
e StoGradMP algorithm has the highest running time.From Figures 6  and 7, we see that although the runtime of the proposed method is much more than the GradMP and StoIHT algorithms, its reconstruction probability is much better than the other two algorithms.
is means that the proposed algorithm has lower complexity than other algorithms, except for the GradMP and StoIHT algorithms in different measurements.
Based on the above analysis, the proposed method with threshold t � 0.6 has better reconstruction performance for different measurements than others and lower computational complexity than StoGradMP.

Conclusion
In this study, an improved reconstruction method was proposed.e proposed method utilizes the limited softthreshold selection strategy to select the most relevant atoms from the preliminary atom set, which could reduce the redundancy of the preliminary atoms set and improve the accuracy of the support atom set estimation, thereby improving the reconstruction precision of the signal and reducing the computational complexity of the algorithm.In addition, the combination with reliability verification conditions ensured the correctness and effectiveness of the proposed method.e simulation results proved that the proposed method has better performance than other recovery methods.
1 , Γ 2 , . . .Γ 2s ].We select the atoms that are larger than the threshold in the atomic set Γ, to form a new atomic set.e index of the new atomic set forms the new atomic index set S *

) Signal estimation by the least square method
where sigma is the maximum gradient value at the k-th iteration.‖r k ‖ 2s represents the largest 2s gradient values from the absolute value of the gradient vector, which corresponds to the atomic index set S k .t is the threshold, whose value range is generally [0.1 ∼ 1.0].It should be noted that if the size of the threshold t is greater than 1, the soft-thresholds selection strategy fails.at is to say that the selection of the preliminary atomic index set S * k cannot be completed.Meanwhile, if the value of the threshold t is too small, it will lead to the redundant atoms cannot be effectively eliminated in the preliminary atomic set.t * sigma is the soft-thresholds selection strategy.Find(.)searches the corresponding index that satisfies the soft-threshold condition and forms a new atomic index set S * k .After the soft-threshold selection strategy is completed, we merge the new atomic set Φ S * k and the support set Φ F k−1 to update the candidate atom set Φ C k , which can be expressed as

) Prune to obtain current support index set Reliability verification conditions 2 If
k (s-sparse approximation of signal x)