Backtracking-Based Iterative Regularization Method for Image Compressive Sensing Recovery

This paper presents a variant of the iterative shrinkage-thresholding (IST) algorithm, called backtracking-based adaptive IST (BAIST), for image compressive sensing (CS) reconstruction. For increasing iterations, IST usually yields a smoothing of the solution and runs into prematurity. To add back more details, the BAIST method backtracks to the previous noisy image using L2 norm minimization, i.e., minimizing the Euclidean distance between the current solution and the previous ones. Through this modification, the BAIST method achieves superior performance while maintaining the low complexity of IST-type methods. Also, BAIST takes a nonlocal regularization with an adaptive regularizor to automatically detect the sparsity level of an image. Experimental results show that our algorithm outperforms the original IST method and several excellent CS techniques.


Introduction
Exploiting a prior knowledge of the original signals is critical to the success of compressive sensing recovery.For compressive sensing (CS) imaging applications, methods that only seek a domain (e.g., discrete cosine transform (DCT), wavelet, and contourlet) to represent and recover images often fail, since they cannot achieve a high enough degree of sparsity.Researchers have therefore explored the use of more elaborate structures, such as minimal total variation [1], block sparsity [2], wavelet trees [3], group sparsity [4], Markov mixture models [5] and nonlocal sparsity [6][7][8][9].
Among these structures, nonlocal sparsity [6][7][8], which exploits the self-repetitive structure exhibited often in natural images, has shown great potential.In NLR-CS (CS via nonlocal low-rank regularization) [6], a patch-based low-rank regularization model is built to enforce the low-rank property over the sets of nonlocal similar patches.In contrast, the D-AMP (denoising-based approximate message passing) algorithm [7] exploits the image self-similarity through the use of nonlocal-based denoisers, such as the BM3D (block-matching and 3D filtering) denoiser [10], which performs hard or soft thresholding with a 3D orthogonal dictionary (3D filtering) on 3D image blocks built by stacking similar patches together (block-matching).Closely relating to D-AMP, BM3D-CS [8] iteratively adds Gaussian noise to approximate the missing part of the spectra and then applies BM3D to the result.
On the other hand, the iterative shrinkage-thresholding (IST) algorithm [11] is frequently adopted to solve sparse models, e.g., in [1,3,6,12], because of the low computational cost.Each IST iteration takes a step in the direction of the negative gradient of the data-fidelity term, followed by the application of the shrinkage function.The two-step IST (TwIST) algorithm [13], the fast IST algorithm (FISTA) [14], and the sparse reconstruction by separable approximation (SpaRSA) algorithm [15] are all accelerated variants of IST.
In this paper, a modified IST method, named backtracking-based adaptive IST (BAIST), is proposed for image CS reconstruction.Similar to IST-type methods, it recovers images with forward-backward iterations.However, unlike IST-type methods that generally yield an over-smoothed solution, BAIST incorporates a backtracking constraint to add back-filtered details and enforce a noisier image through the minimization of the distance between the current solution and the previous ones.Furthermore, the nonlocal regularization with an adaptive regularization parameter set according to the current residual is to adaptively estimate the sparsity level of the image.Experimental results show that our method impressively outperforms previous methods in terms of reconstruction accuracy, and also achieves a competitive advantage in computation complexity.

Compressive Sensing and IST Algorithm
The CS recovery method recovers a signal x 0 ∈ C n from its randomized linear measurements y = Ax 0 , y ∈ C m , where A ∈ C m×n is the measurement matrix.Since m < n, it is an underdetermined problem.More than one solution can yield the same measurements, but if x is sparse, we can recover the original signal by minimizing the following L1 problem: where λ is a regularization parameter.
The above problem can be efficiently solved by the IST algorithm [11].Specifically, the general step of IST is summarized in Algorithm 1. IST recovers the original signal with forward-backward iterations alternating between a gradient descent step (step (a)) and a proximal denoising correction (step (b)).In step (a), the data-fidelity term is f is its gradient at the point x k ; A * is the conjugate transpose of A; and c is a step size.A = SF is the undersampled Fourier matrix used in the paper, where F is the 2D Fourier transform and S is a selection matrix containing m rows of the identity matrix.Step (b) is a proximity operator, which has a close form solution: Algorithm 1 IST [11] Input: y, A, λ, c, x 0 = 0

New Objective Function
There exist abundant similar image patches in natural images, which can be utilized through nonlocal regularizations for high-order sparsity in image CS recovery.In the paper, the image self-similarity is characterized by means of the last item in the following equation: where the last item is derived from BM3D [10], which is used to seek patch correlation in image denoising.To obtain Γ 3D x i , we first stack similar image patches to a 3D group, where x i is the i-th group and G is the total number of groups; then 3D wavelet transform (WT) Γ 3D is conducted on x i , which is achieved by conducting 2D WT on each patch and then conducting 1D WT along the third axis.
Problem (3) can be solved by IST; however, as the iteration increases, the recovered image tends to be over-smoothed.As a result, the algorithm runs into prematurity.To avoid this problem, we add a backtracking constraint to add back-filtered details and force a noisier image: The last item in Equation (4) enforces the current solution which is similar to the previous noisier ones, and the similarity degrees between them are controlled by parameters b 1 , b 2 , . . ., b B .Intuitively, a larger b j means more noise added back.We have observed, in a large number of experiments, that the algorithm always achieves improvements and converges for a wide range of choices of parameters b 1 , b 2 , . . ., b B .

BAIST Algorithm for Solving the New Objective Function
Firstly, we offer the procedure for solving Equation (3) with the IST algorithm outlined in Algorithm 2, named ISTANR (IST with adaptive nonlocal regularization).In step (a), the fidelity term is optimized by gradient descent to get a solution x g , which can be regarded as a noise image.
Step (b) is a proximity denoising correction step implemented through the use of nonlocal thresholding.

Algorithm 2 IST with Adaptive Nonlocal Regularization
Secondly, we solve Equation (4) with ISTANR by letting g and substituting f (x) for g(x).The proposed method named BAIST (backtracking-based adaptive IST) is summarized in Algorithm 3, where ∇g( which obtains the noisy observation of x (i.e., x g ) using the fidelity term, BAIST super-adds the difference between the current solution and the previous ones, and thus obtains a noisier x g for avoiding prematurity.It can be confirmed in Figure 1 that BAIST yields the intermediate result of x g with more details.

Algorithm 3 Backtracking-Based Adaptive IST
The key to success in implementing both BAIST and ISTANR is to compute the proximity denoising correction.We use a lemma in [15], namely the minimization problem: x − a 2 2 + λ x 0 (5) which has a closed form solution, expressed as: Since the orthogonal transform Γ 3D has the property of energy conservation, step (b) can be written as: Its solution is: where x i and (x g ) i are the i-th similar patches 3D group of x and x g , respectively.Since a patch may be assigned to several 3D groups, the patches in the last item of Equation ( 7) are repetitively computed.Let N be a scale factor for balancing the increasing energy caused by the repetitive computation; however, we can directly set the threshold of the hard thresholding operator in Equation ( 8) δ = 6µcN, rather than computing N and selecting an appropriate regularization parameter µ.
Apparently, the threshold δ should be reduced with the increasing iteration number to gain a more and more clear image, which is similar to the continuation strategy [15] in IST-type methods.A natural idea is that δ is proportional to the observation residual in CS reconstruction: where s is a scale factor and can be set empirically.As a result, we obtain an adaptive thresholding method.Note that when more noise is added back to the current solution x g by BAIST, the threshold δ is simultaneously increased, and thus the solution would not get noisier and the algorithm can approach the noise-free image.
Since the orthogonal transform 3D Γ has the property of energy conservation, step (b) can be written as: Its solution is: where i x and ( ) x are the i-th similar patches 3D group of x and xg, respectively.Since a patch may be assigned to several 3D groups, the patches in the last item of Equation ( 7) are repetitively computed.Let N be a scale factor for balancing the increasing energy caused by the repetitive computation; however, we can directly set the threshold of the hard thresholding operator in Equation (8) = 6 cN δ μ , rather than computing N and selecting an appropriate regularization parameter μ .Apparently, the threshold δ should be reduced with the increasing iteration number to gain a more and more clear image, which is similar to the continuation strategy [15] in IST-type methods.A natural idea is that δ is proportional to the observation residual in CS reconstruction: where s is a scale factor and can be set empirically.As a result, we obtain an adaptive thresholding method.Note that when more noise is added back to the current solution xg by BAIST, the threshold δ is simultaneously increased, and thus the solution would not get noisier and the algorithm can approach the noise-free image.

Experiments
To evaluate the reconstruction performances of the proposed algorithm in items of average peak signal to noise ratio (PSNR), runtime and visual quality, we compare the BAIST algorithm with six excellent algorithms, including the ISTANR algorithm; two other algorithms based on the IST method FCSA [1]; WaTMRI [3]; two Bayesian CS methods, Turbo-AMP [5], D-AMP [7]; and the nonparametric method BM3D-CS [8].Three experiments are carried out on 20 standard test images with a size of 256 × 256 (shown in Figure 2) at five sampling ratios (16%, 18%, 20%, 22%, and 24%).

Experiments
To evaluate the reconstruction performances of the proposed algorithm in items of average peak signal to noise ratio (PSNR), runtime and visual quality, we compare the BAIST algorithm with six excellent algorithms, including the ISTANR algorithm; two other algorithms based on the IST method FCSA [1]; WaTMRI [3]; two Bayesian CS methods, Turbo-AMP [5], D-AMP [7]; and the nonparametric method BM3D-CS [8].Three experiments are carried out on 20 standard test images with a size of 256 × 256 (shown in Figure 2) at five sampling ratios (16%, 18%, 20%, 22%, and 24%).All experiments are conducted on a desktop with 3.80 GHz AMD A10-5800K CPU.The Matlab version is R2014a.
All experiments are conducted on a desktop with 3.80 GHz AMD A10-5800K CPU.The Matlab version is R2014a.Several parameters need to be set for running BAIST and ISTANR.We set the step size c = 1, and the scale factor of the threshold in Equation ( 9) s = 10 in the first 10 iterations, and s = 3 in the rest of the iterations because of the presence of more noise in the early iterations.The 3D wavelet transform 3D Γ is composed of 2D bior1.5 and 1D Haar.To construct 3D groups by stacking similar patches, we need to set the following parameters: the size of each patch is 8 × 8; the size of the window for searching matched patches is 25 × 25; the number of best-matched patches is 16; and the sliding step to process every next reference patch is set to six.For BAIST, a wide range of choices of parameters b1, b2, …, bB is available to achieve improvements; in our experiments, we set b1 = 0.9, b2 = 0.7, b3 = 0.4, b4 = 0.3, and backtracking is adopted after 10 iterations.We set the maximum iterations K = 10 for Turbo-AMP, K = 100 for BM3D-CS, and K = 50 for the rest of the algorithms.For the algorithms except for BAIST and ISTANR, the default settings in their codes are used.

Average PSNR Evaluation
In the experiments, we generated CS measurements by randomly sampling the Fourier transform coefficients of input images by following the previous works [1,3].The average PSNR results from running each image five times are shown in Figure 3. From Figure 3, one can see that (1) four nonlocal sparsity-based methods, BAIST, ISTANR, D-AMP, and BM3D-CS, performed better than the others; (2) the highest PSNR results were achieved by the proposed algorithm BAIST, the average PSNR gains of which, compared to the next-best algorithm D-AMP, can be as much as 0.71 dB; (3) through the use of backtracking, BAIST outperformed its original IST version ISTANR by 0.83 dB on average; (4) and the standard deviations of all algorithms at low sampling ratios were much greater than the ones at high sampling ratios because the randomness is higher at low sampling ratios.These results validate the superiority of BAIST in objective quality and the effectiveness of backtracking.Several parameters need to be set for running BAIST and ISTANR.We set the step size c = 1, and the scale factor of the threshold in Equation ( 9) s = 10 in the first 10 iterations, and s = 3 in the rest of the iterations because of the presence of more noise in the early iterations.The 3D wavelet transform Γ 3D is composed of 2D bior1.5 and 1D Haar.To construct 3D groups by stacking similar patches, we need to set the following parameters: the size of each patch is 8 × 8; the size of the window for searching matched patches is 25 × 25; the number of best-matched patches is 16; and the sliding step to process every next reference patch is set to six.For BAIST, a wide range of choices of parameters b 1 , b 2 , . . ., b B is available to achieve improvements; in our experiments, we set b 1 = 0.9, b 2 = 0.7, b 3 = 0.4, b 4 = 0.3, and backtracking is adopted after 10 iterations.We set the maximum iterations K = 10 for Turbo-AMP, K = 100 for BM3D-CS, and K = 50 for the rest of the algorithms.For the algorithms except for BAIST and ISTANR, the default settings in their codes are used.

Average PSNR Evaluation
In the experiments, we generated CS measurements by randomly sampling the Fourier transform coefficients of input images by following the previous works [1,3].The average PSNR results from running each image five times are shown in Figure 3. From Figure 3, one can see that (1) four nonlocal sparsity-based methods, BAIST, ISTANR, D-AMP, and BM3D-CS, performed better than the others; (2) the highest PSNR results were achieved by the proposed algorithm BAIST, the average PSNR gains of which, compared to the next-best algorithm D-AMP, can be as much as 0.71 dB; (3) through the use of backtracking, BAIST outperformed its original IST version ISTANR by 0.83 dB on average; (4) and the standard deviations of all algorithms at low sampling ratios were much greater than the ones at high sampling ratios because the randomness is higher at low sampling ratios.These results validate the superiority of BAIST in objective quality and the effectiveness of backtracking.Several parameters need to be set for running BAIST and ISTANR.We set the step size c = 1, and the scale factor of the threshold in Equation ( 9) s = 10 in the first 10 iterations, and s = 3 in the rest of the iterations because of the presence of more noise in the early iterations.The 3D wavelet transform 3D Γ is composed of 2D bior1.5 and 1D Haar.To construct 3D groups by stacking similar patches, we need to set the following parameters: the size of each patch is 8 × 8; the size of the window for searching matched patches is 25 × 25; the number of best-matched patches is 16; and the sliding step to process every next reference patch is set to six.For BAIST, a wide range of choices of parameters b1, b2, …, bB is available to achieve improvements; in our experiments, we set b1 = 0.9, b2 = 0.7, b3 = 0.4, b4 = 0.3, and backtracking is adopted after 10 iterations.We set the maximum iterations K = 10 for Turbo-AMP, K = 100 for BM3D-CS, and K = 50 for the rest of the algorithms.For the algorithms except for BAIST and ISTANR, the default settings in their codes are used.

Average PSNR Evaluation
In the experiments, we generated CS measurements by randomly sampling the Fourier transform coefficients of input images by following the previous works [1,3].The average PSNR results from running each image five times are shown in Figure 3. From Figure 3, one can see that (1) four nonlocal sparsity-based methods, BAIST, ISTANR, D-AMP, and BM3D-CS, performed better than the others; (2) the highest PSNR results were achieved by the proposed algorithm BAIST, the average PSNR gains of which, compared to the next-best algorithm D-AMP, can be as much as 0.71 dB; (3) through the use of backtracking, BAIST outperformed its original IST version ISTANR by 0.83 dB on average; (4) and the standard deviations of all algorithms at low sampling ratios were much greater than the ones at high sampling ratios because the randomness is higher at low sampling ratios.These results validate the superiority of BAIST in objective quality and the effectiveness of backtracking.

Visual Quality and Runtime Evaluation on Boat Image
For better quality perception, experiments performed on the boat image presented similar PSNR results, shown in Figure 4.The CPU time vs. PSNR and the iteration number vs. PSNR of different algorithms carried out on the boat image with 20% sampling are presented in Figures 5 and 6, respectively.The corresponding recovered images are given in Figure 7. From Figure 7, we can clearly see that four nonlocal sparsity-based methods are still better than the others.Among them, BAIST enjoyed great advantages over ISTANR in producing a clearer image, e.g., on the area of the prow, due to the backtracking strategy.The reconstructed boat image of BAIST and D-AMP was close to some degree at 20% sample rates, but BAIST was still clearer than D-AMP with a zoomed-in version, also achieving higher PSNR results after 19 s as shown in Figure 5, and after 34 iterations as shown in Figure 6.Besides, for BM3D-CS with more iterations, its results are presented for every two iterations in Figure 6.The superiority of the proposed algorithm in visual quality and runtime could be demonstrated by these results.

Visual Quality and Runtime Evaluation on Boat Image
For better quality perception, experiments performed on the boat image presented similar PSNR results, shown in Figure 4.The CPU time vs. PSNR and the iteration number vs. PSNR of different algorithms carried out on the boat image with 20% sampling are presented in Figures 5  and 6, respectively.The corresponding recovered images are given in Figure 7. From Figure 7, we can clearly see that four nonlocal sparsity-based methods are still better than the others.Among them, BAIST enjoyed great advantages over ISTANR in producing a clearer image, e.g., on the area of the prow, due to the backtracking strategy.The reconstructed boat image of BAIST and D-AMP was close to some degree at 20% sample rates, but BAIST was still clearer than D-AMP with a zoomed-in version, also achieving higher PSNR results after 19 s as shown in Figure 5, and after 34 iterations as shown in Figure 6.Besides, for BM3D-CS with more iterations, its results are presented for every two iterations in Figure 6.The superiority of the proposed algorithm in visual quality and runtime could be demonstrated by these results.

Visual Quality and Runtime Evaluation on Boat Image
For better quality perception, experiments performed on the boat image presented similar PSNR results, shown in Figure 4.The CPU time vs. PSNR and the iteration number vs. PSNR of different algorithms carried out on the boat image with 20% sampling are presented in Figures 5  and 6, respectively.The corresponding recovered images are given in Figure 7. From Figure 7, we can clearly see that four nonlocal sparsity-based methods are still better than the others.Among them, BAIST enjoyed great advantages over ISTANR in producing a clearer image, e.g., on the area of the prow, due to the backtracking strategy.The reconstructed boat image of BAIST and D-AMP was close to some degree at 20% sample rates, but BAIST was still clearer than D-AMP with a zoomed-in version, also achieving higher PSNR results after 19 s as shown in Figure 5, and after 34 iterations as shown in Figure 6.Besides, for BM3D-CS with more iterations, its results are presented for every two iterations in Figure 6.The superiority of the proposed algorithm in visual quality and runtime could be demonstrated by these results.

Conclusions
A new iterative regularization method for image compressive sensing recovery is proposed in this paper.By adapting a simple backtracking constraint to the IST method, the proposed method BAIST can flexibly add back some details that are filtered out incorrectly in the previous processing, and hence gives a much better reconstruction performance than the tested IST method ISTANR.Also, it is very appropriate for image CS reconstruction applications since a nonlocal wavelet thresholding method with an adaptive threshold is able to achieve high enough degree of sparsity.In the experiments, our algorithm was shown to outperform six excellent algorithms in reconstruction quality, and also achieved a competitive advantage in runtime.

Conclusions
A new iterative regularization method for image compressive sensing recovery is proposed in this paper.By adapting a simple backtracking constraint to the IST method, the proposed method BAIST can flexibly add back some details that are filtered out incorrectly in the previous processing, and hence gives a much better reconstruction performance than the tested IST method ISTANR.Also, it is very appropriate for image CS reconstruction applications since a nonlocal wavelet thresholding method with an adaptive threshold is able to achieve high enough degree of sparsity.In the experiments, our algorithm was shown to outperform six excellent algorithms in reconstruction quality, and also achieved a competitive advantage in runtime.(No. 61327005, No. 61302120) and the Specialized Research Fund for the Doctoral Program of Higher Education (No. 20130172120045).The authors also gratefully acknowledge the helpful comments and suggestions of the reviewers, which have improved the presentation.

Conclusions
A new iterative regularization method for image compressive sensing recovery is proposed in this paper.By adapting a simple backtracking constraint to the IST method, the proposed method BAIST can flexibly add back some details that are filtered out incorrectly in the previous processing, and hence gives a much better reconstruction performance than the tested IST method ISTANR.Also, it is very appropriate for image CS reconstruction applications since a nonlocal wavelet thresholding method with an adaptive threshold is able to achieve high enough degree of sparsity.In the experiments, our algorithm was shown to outperform six excellent algorithms in reconstruction quality, and also achieved a competitive advantage in runtime.
(No. 20130172120045).The authors also gratefully acknowledge the helpful comments and suggestions of the reviewers, which have improved the presentation.

Figure 2 .
Figure 2. The 20 widely used test images.

Figure 3 .
Figure 3. Average PSNR at different sample rates.

Figure 2 .
Figure 2. The 20 widely used test images.
conducted on a desktop with 3.80 GHz AMD A10-5800K CPU.The Matlab version is R2014a.

Figure 2 .
Figure 2. The 20 widely used test images.

Figure 3 .
Figure 3. Average PSNR at different sample rates.

Figure 3 .
Figure 3. Average PSNR at different sample rates.

Figure 4 .Figure 5 .
Figure 4. Average PSNR on boat image at different sample rates.

Figure 4 .
Figure 4. Average PSNR on boat image at different sample rates.

Figure 4 .Figure 5 .
Figure 4. Average PSNR on boat image at different sample rates.

Figure 5 .
Figure 5. CPU time vs. PSNR on boat image with 20% sampling.Figure 5. CPU time vs. PSNR on boat image with 20% sampling.

Figure 7 .
Figure 7. Visual results on boat image with 20% sampling.