Compressive Sensing via Nonlocal Smoothed Rank Function.

Compressive sensing (CS) theory asserts that we can reconstruct signals and images with only a small number of samples or measurements. Recent works exploiting the nonlocal similarity have led to better results in various CS studies. To better exploit the nonlocal similarity, in this paper, we propose a non-convex smoothed rank function based model for CS image reconstruction. We also propose an efficient alternating minimization method to solve the proposed model, which reduces a difficult and coupled problem to two tractable subproblems. Experimental results have shown that the proposed method performs better than several existing state-of-the-art CS methods for image reconstruction.

Introduction

Compressive sensing (CS) [1, 2] allows us to reconstruct high dimensional data with only a small number of samples or measurements, and captures only useful information and has the potential of significantly improving the energy efficiency of sensors in the real-world applications. The key idea behind CS is that the majority of real-world signals including images and videos, can be sparsely represented by given some appropriate basis. Due to the positive theoretical and experimental results, many CS-based imaging methods have been proposed and applied to the various areas such as magnetic resonance imaging [3, 4] in medicine, compressed spectral and hyperspectral imaging [5, 6] in industry, neural network [7] in biotechnology.

CS also makes it possible to well restore corrupted signals at a fast speed and the small memory cost. Conventional CS recovery uses ℓ1 norm to characterize the sparsity of a signal, and the resulting convex optimization problems are tractable. Although there are several methods can be used to efficiently solve the ℓ1 regularization based model for signal recovery [8-10], they only achieve suboptimal recovery performance due to their relaxation of the ℓ0 norm based sparse optimization [11]. More recently, structured or group sparsity based methods [12-14] and nonlocal sparsity based methods [15, 16] have provided better results for CS recovery. Intuitively, the structured or group sparsity can reduce the degrees of freedom in the solution, and the nonlocal sparsity explicitly exploits self-similarities of the signal, thereby obtaining more accurate recovery performance than the common sparsity.

In CS studies, a number of works have suggested that non-convex optimization based approach often yields better results than convex ones though costing higher computational complexity. Therefore, we propose a CS recovery model by considering a non-convex smoothed function to approximate the rank, denoted as smoothed rank function (SRF), as a low-rank regularization. The emphasis behind the proposed model is the utilization of nonlocal sparsity by image patch grouping. Concretely, we group a set of similar image patches to form a matrix X for i-th exemplar image patch, which is extracted from the test image. Then the matrix X is low-rank since similar patches have similar structure. To solve the low-rank optimization problem, we minimize the SRF function to approximate the rank minimization problem, and the resulting problem is a non-convex optimization problem. In order to avoid getting trapped in local solutions, we initialize a rough approximation of the rank, and gradually improve the approximation as the iteration proceeds.

The basic idea of using nonlocal sparsity for image patches has already been used in [16-18] with the very impressive results. There have been an abundant research literatures using the patch-based low-rank as plenty with nonlocal sparsity [19-22]. In [20, 21], the low-rank problem is solved by minimizing the nuclear norm of the low-rank matrix, which leads to a convex minimization problem with many efficient methods available. In [22], a more accurate approximation for rank is proposed by exploiting the logdet function, which can be derived as a weighted nuclear norm. Specifically, compared with previous surrogates for the rank such as nuclear norm and weighted nuclear norm, our surrogate is differentiable and can approximate the rank adaptively.

The outline of this paper is organized as follows. Firstly, we briefly review the background of CS. Secondly, the SRF function is introduced and the model based on SRF for CS image reconstruction is proposed. Thirdly, the optimization algorithm is presented to solve the proposed model. To demonstrate the effectiveness of the proposed method, we show some numerical results on several test images in the following section. In the end, the conclusion and future work are given.

Background

The CS recovery problem aims to find the sparsest solution from the underdetermined linear system y = Φ x, where is the measurements and , M < N is the measurement matrix. It can be formulated as follows: where ||⋅||0 is the ℓ0 norm counting the number of nonzero elements of x. In practical applications, such as signal reconstruction problem, the measurement noise is unavoidable. Then the noisy CS recovery problem is formulated as where ||⋅||2 is the ℓ2 norm and ϵ is the residual error. However, since problem (2) is NP-hard, it is infeasible to solve it directly. Thus it is proposed to replace the ℓ0 norm by the convex ℓ1 norm, namely which is the well-known basis pursuit problem [2]. The problem (3) is the convex relaxation of problem (2) and easy to solve by some methods including iterative shrinkage algorithm [8], alternating direction method of multipliers (ADMM) [23] and Bregman split algorithm [24]. Although problems (2) and (3) are fundamentally different, they return the same solution in many interesting situations [25]. By using an appropriate regularization parameter λ, we can equivalently rewrite the problem (3) as follows:

By replacing the ℓ1 norm by the non-convex ℓ (0 > p > 1) norm, in [26], the non-convex optimization problem based on ℓ norm can achieve the more exact CS reconstruction result than the convex one based on the ℓ1 norm. Moreover, it has been shown in [18, 27] that the nonlocal sparsity exploiting the self-similarity of the natural image leads to the state-of-the-art performance. The nonlocal sparsity is a significant prior for CS recovery. In this work, we will propose a non-convex CS image recovery approach, which exploits the nonlocal sparsity.

Nonlocal smoothed rank function

In this section, we present a procedure of patch grouping that uses image self-similarity and leads to the low-rank problem. Hence one underlying assumption is the image exhibits abundant self-similarities. For a test image, we select some exemplar patches of size , reordered into the column vector lexicographically and denoted as . For each exemplar patch, the underlying assumption makes it possible to find a number of similar patches. Specially, we employ a variant of k-nearest-neighbor search in a local window for each exemplar patch to find its similar patches as follows: where c is a pre-defined threshold and G indexes the positions of corresponding similar patches. After patch grouping, we aline the similar patches as column vectors to form a matrix , for the i-th exemplar patch . Under the above assumption, each matrix X has the low-rank property since the similar patches have the similar structures. In this way, the nonlocal sparsity leads to a rank minimization problem in our image reconstruction approach.

Since the resulting matrix X also contains the noise during patch grouping, in order to obtain a clean and clear image to well match the ground truth, let X = L+W, where L denotes the low-rank matrix and W represents the noise matrix. We can obtain the low-rank matrix L by solving the following problem: where ||⋅|| is the Frobenius norm and ϵ is the residual error. Unfortunately, problem (6) is NP-hard. A popular heuristic method is to adopt the nuclear norm (sum of the singular values) to replace the rank, i.e., where ||⋅||* is the nuclear norm. The problem (7) is a convex surrogate of the rank and can be solved efficiently by the singular value thresholding (SVT) algorithm [28]. However, it is in practice suboptimal due to equally treating each singular value. In [19], Gu et al. have demonstrated that non-convex low-rank approximations adaptively treating the singular values at different scales yield better results than those convex ones.

In [29, 30], they have been shown that for a low-rank matrix , the rank can be approximated by the following SRF function: where σ(X) is the j-th singular value of X, ℓ = min{n, m} and δ is an adjustable parameter. Although the problem of minimizing G(⋅) function using a small δ will lead to many local minima, the solution of minimizing G(⋅) function converges to the minimum rank solution as δ goes to zero [30] as the red curve shown in Fig 1. In order to avoid trapping in local solutions, we initialize a large δ and gradually decrease δ to improve the degree of approximation for the rank. With the decreasing δ, the rank can be better approximated by SRF Eq (8). This technique refines the minimizer of the SRF minimization problem by considering the minimizer of the previous iteration (large δ) as the new initial point of current iteration (small δ), which makes the solution of the SRF minimization problem get closer to the minimum rank solution during the iterations. Fig 1 intuitively illustrates the comparison of SRF, the rank, the nuclear norm and logdet function in the case of a scalar. We observe that SRF can better approximate the rank than the nuclear norm and logdet function. Therefore, we speculate that the problem of minimizing SRF toward rank minimization problem could achieve better performance. Then, we consider the SRF function as the low-rank regularization for CS image recovery.

Fig 1

Performance of SRF = G(x), rank = ||x||0, the nuclear norm = ||x||1 and logdet = log(|x|+ε) in the case of a scalar, where δ = 1 and ε = 1.

For the low-rank matrix L, the rank minimization problem can be replaced to which can be reformulated with an appropriate parameter λ, Obviously, the problem (10) is smoothed and differentiable, and we can adopt the gradient descent method to solve it. For each matrix X, the same method can be used to obtain the corresponding low-rank matrix L.

For CS image reconstruction problem, based on above presented G(L), we propose the global model as follows: where η is a regularization parameter and P x = [P x,P x, ⋯,P x] is the matrix constituted by the set of similar patches for each exemplar patch x. The model (11) exploits the nonlocal sparsity of the image patches and non-convexity of the SRF function G(⋅). Therefore, we conjecture the proposed method can achieve the good performance. To solve the model (11), we develop an efficient alternative minimization method in the next section.

Optimization algorithm

It is difficult to directly solve the global model (11) since x and L is coupled. We use the alternative minimization method to decouple x and L as follows:

For L-subproblem, its optimization problem is

For x-subproblem, its optimization problem is

In the following subsections we will present the optimization algorithm to solve subproblems Eqs (12) and (13) in detail.

Low-rank matrix optimization algorithm

Let . Then the low-rank matrix L can be obtained via the gradient descent method, i.e., where μ( is the step size in the k-th iteration and ∇f(⋅) is the gradient of f(⋅). Before deriving ∇f(⋅), we show the gradient of SRF G(L) at L [30] as follows: where matrixes U and V and singular values σ (i = 1, ⋯, ℓ) come from the singular value decomposition (SVD) of L, which is obtained in previous iteration, Then, it is easy to derive the gradient of the function f(L), According to the gradient descent method, at each iteration, one has In addition, the step size μ in formula (17) should be set in a decreasing order to return a better result. Following the papers [30, 31], we set μ = δ2 to decrease μ proportional to δ2. Then, above formula (17) incorporating the Eq (16) can be rewritten as where ξ = 1 − 2δ2.

Image reconstruction via alternating direction method of multipliers

After obtaining the low-rank matrix L, we reconstruct the image x via solving the problem (13). It is clear that the problem (13) is a quadratic optimization problem with a closed-form solution, However, the inverse of the matrix is very large and difficult to compute. Therefore, we consider to reconstruct the image x in the framework of ADMM, which also leads to the closed-form solutions for each subproblem. The ADMM method is often used to image restoration [32-36].

Based on the definition of ADMM, we present the augmented Lagrangian function of the problem (13) as follows: where z = x is the auxiliary variable, d is the Lagrangian multiplier, and β > 0 is a scalar. With respect to the variables x, z and d, they are decoupled in the framework of ADMM, thus, can be solved separately, leading to the following iterations: Clearly, both subproblems Eqs (21) and (22) are quadratic optimization problem and have the closed-form solutions. For the subproblem Eq (21), its explicit solution is: For the subproblem Eq (22), according to its first-order derivation for x, we can derive the following equation: Considering measurement matrix Φ is a partial Fourier transform matrix, we can transform above problem from image space to Fourier space to efficiently obtain x. Concretely, we let Φ = DF, where D is the down-sampling matrix and F is the fourier transform matrix. It is easy to achieve and then x( can be drawn

We simultaneously obtain the low-rank matrix L and the image x by the alternative minimization method, and the overall procedure is summarized below as Algorithm 1. Empirically we have found that Algorithm 1 is convergent, but in theory the convergence analysis of Algorithm 1 is difficult to give due to the non-convex subproblem. Although there are some papers have proved the convergence of their non-convex optimization problem, these proofs hold in a few unrealistic and rigorous assumptions. Specially, to save computational complexity, we set J = K = 1 in Algorithm 1. Moreover, we estimate an image using the discrete cosine transform (DCT) method as the initial solution for a better initial point, which has been seen in [22]. As iteration increases, the high accuracy results will be achieved.

The complexity of Algorithm 1 is O([T + log N]N) (N is the total number of image pixels and T is the average complexity to compute similar patches per exemplar patch), which is mainly generated by the DCT method in 1 step of initialization and the Fourier transform in 4 (b) step of inner loop. The complexity of SVD in 3 (a) step of inner loop, i.e. O(n × m2), can be ignored due to n, m ≪ N. Therefore, the proposed method is practical feasible and promising.

Algorithm 1 CS via the SRF function

Initialization:

1. Estimate an initial image using the discrete cosine transform (DCT) method;

2. Set regularization parameters λ and η, and parameter β > 0;

3. Let , μ(0) = (δ(0))2, , d(1) = 0;

4. Find a set of similar patches by the method Eq (5) for each exemplar patch using x(1);

Outer loop: for s = 1, 2, …, S do

1. Match similar patches into a matrix X for each exemplar patch using x(;

2. Set ;

3.Inner loop (solving the L-subproblem Eq (12)): for j = 1, 2, …, J do

(a) Compute the SVD of ;

(b) Compute by Eq (16);

(c) Compute ;

(d) Update δ( = cδ(, 0 < c < 1;

(e) Update μ( = (δ()2;

(f) Output until j = J.

End for

4.Inner loop (solving the x-subproblem Eq (13)): for k = 1, 2, …, K do

(a) Update z( by the Eq (24);

(b) Update x( by the Eq (27);

(c) Update d( by the Eq (23);

(d) Output x( = x( until k = K.

End for

5. Output the reconstructed image until s = S.

End for

Numerical results

Here, we present the experimental results of our approach for CS image reconstruction based on the SRF regularization. We generate the CS measurements by random and pseudo-radial sampling of the Fourier transform coefficients of test images respectively. The number of measurements is M = rate*N, where rate is the sampling rate. We use peak signal to noise ratio (PSNR) and Structural SIMilarity (SSIM) index [37] as the quantitative measures in our numerical experiments, and the PSNR is defined as where MAX is the maximum possible pixel value of the image and MSE is the mean squared error, defined as where f and g are the original image and the restored image, respectively.

In all experiments, the main parameters of the Algorithm 1 are set as follows: patch size ; the number of similar patches for each exemplar patch m = 45 (the more similar patches form a more low-rank matrix but lead to high computational complexity); initialize value , where is a constant around two times of the largest singular value of initial L [30]; the parameter c is set as 0.08 experimentally in formula δ( = cδ(; the outer loop iteration number S = 400 (that is selected based on the convergence rate of Algorithm 1). To achieve better performance, we select exemplar patch in each 5 pixels along both horizontal and vertical directions.

In addition, the total variation (TV) method [38], the BM3D based CS method (BM3D-CS) [16], the nuclear norm method (denoted as NLR-CS-baseline) [22] and the logdet function method (denoted as NLR-CS) [22] for CS image recovery are compared with the proposed SRF approach (called as SRF-CS). The TV method just considered the underlying sparsity of an image, and the BM3D-CS method used the nonlocal sparsity with an outstanding performance. The NLR-CS-baseline method exploited the nuclear norm to replace the rank for solving the rank minimization problem. The NLR-CS approach adopted the logdet function to achieve a more accurate surrogate than the nuclear norm for the rank. The source codes of above methods [16, 22, 38] are publicly downloaded from the authors’s websites. These methods are significantly state-of-the-art CS algorithms for image recovery. We have carefully tuned their parameters to obtain the best results for fair comparison. Both noiseless and noisy experiments are presented to demonstrate the performance of the proposed approach for CS recovery. The test images (256 × 256) are exhibited in Fig 2, where Fig 2(a), 2(c), 2(e) and 2(f) can be publicly downloaded from the test image net http://decsai.ugr.es/cvg/dbimagenes/g256.php, Fig 2(b), 2(d) and 2(g) can be publicly downloaded from the author’s homepage of [22] http://see.xidian.edu.cn/faculty/wsdong/NLR_Exps.htm.

Fig 2

The test images: (a) Book; (b) Barbara; (c) Plane; (d) Boat; (e) Number; (f) Moonscape; (g) Head.

Noiseless experiments for CS recovery

We reconstruct the images from less M = 0.05N and M = 0.1N measurements by randomly sampling the Fourier transform coefficients of the test images. Fig 3 shows the restored images by TV [38], NLR-CS-baseline [22], NLR-CS [22] and the proposed methods for CS image recovery. The TV method for CS recovery can’t work well due to less measurements. The close-ups of Number image and Book image obtained by the TV method are indistinct and lose much structured information. The NLR-CS-baseline method gives rise to some artificial shadow like fog as shown in the close-up of Number image. It is well known that the NLR-CS method is a very competitive method for CS recovery due to its remarkable performance. For these test images in Fig 3, however, the performance of the NLR-CS method is undesirable. Obviously, our approach achieves the best visual quality and highest PSNR value among all test methods. Noteworthily, the better performance of our approach compared with the NLR-CS-baseline method verifies that the SRF function Eq (8) is a more accurate surrogate for the rank than the nuclear norm. Moreover, Table 1 displays the PSNR values and the SSIM values obtained by all test methods.

Fig 3

Reconstructed images by randomly sampling.

From top to bottom for each column: Number image, close-up of the red rectangle of Number image, Moonscape image, close-up of the red rectangle of Moonscape image, Book image, close-up of the red rectangle of Book image. From left to right for each row: the original image, the TV method, the NLR-CS-baseline method, the NLR-CS method, the proposed method.

Table 1

The PSNR (dB) and SSIM values of the restoration results achieved by the TV method [38], the NLR-CS-baseline method [22], the NLR-CS method [22] and the proposed method for noiseless images.

Image	Number of measurements	Method	PSNR (dB)	SSIM
Number	M = 0.05N	TV	24.10	0.6545
		NLR-CS-baseline	26.29	0.8970
		NLR-CS	26.70	0.9407
		proposed	29.04	0.9497
Moonscape	M = 0.1N	TV	26.41	0.5450
		NLR-CS-baseline	30.19	0.7037
		NLR-CS	29.92	0.7004
		proposed	30.95	0.7143
Book	M = 0.05N	TV	32.41	0.7304
		NLR-CS-baseline	37.41	0.9371
		NLR-CS	38.44	0.9430
		proposed	39.77	0.9528

Although the proposed method leads to a non-convex optimization problem, we can see that the proposed method for CS image recovery converges in a few iterations in Fig 4. In the following subsection, we will verify that the proposed method is also robust in a noisy situation.

Fig 4

PSNR value vs. Iteration plots for the proposed method and the NLR-CS method.

(a) Book image, measures number M = 0.05N (b) Number image, measurements number M = 0.1N.

Noisy experiments for CS recovery

In the noisy experiments, we add Gaussian noise with 0.05 standard deviation to the original images. Firstly, we recover the Head image from M = 0.15N measurements and the Boat image from M = 0.2N measurements by randomly sampling their Fourier transform coefficients respectively. The restoration results by the TV [38], BM3D-CS [16], NLR-CS-baseline [22], NLR-CS [22] and the proposed methods are presented in Fig 5. Clearly, the restored Head image becomes blocky obtained by the TV method, and contains heavy noise obtained by the NLR-CS-baseline method. Meanwhile, some edges in the close-up of the restored Head image achieved by the NLR-CS method are unclear, and the close-up of Head image generated by the BM3D-CS method is blurring. However, compared with other methods, the reconstructed Head image by the proposed approach returns the clearer edges and less noise. For the Boat image, we can see that the proposed approach also outperforms the other test methods with the better visual quality, the higher PSNR and SSIM values. Thus the proposed method with less measurements is robust in noisy situation.

Fig 5

Reconstructed images by randomly sampling.

From top to bottom for each column: Head image, close-up of the red rectangle of Head image, Boat image, close-up of the red rectangle of Boat image. From left to right for each row: the original image, the TV method, the NLR-CS-baseline method, the NLR-CS method, the BM3D-CS method, the proposed method.

Secondly, we recover the Boy image from M = 0.24N (65 radial lines) measurements and the Barbara image from M = 0.29N (80 radial lines) measurements by pseudo-radial sampling their Fourier transform coefficients. The pseudo-radial sampling way generates streaking artifacts leading to more difficult image reconstruction than randomly sampling way. We illustrate the performance results of all the test methods in Fig 6. Visually, the close-ups of Plane image generated by the TV, NLR-CS-baseline, NLR-CS and BM3D-CS methods are blurring. Nevertheless, the close-up of Plane image obtained by our approach is good with less noise and higher PSNR value. The close-ups of Barbara image achieved by the NLR-CS-baseline method and the NLR-CS method create some artificial strips, which do not exist in the original Barbara image. The bad performance may be caused by the pseudo-radial sampling way or noise. The texture of restored Barbara image by BM3D method is lost. However, the proposed approach removes the artifacts and a lot of noise returning a better result, which demonstrates its robustness. In addition, Table 2 displays the PSNR and SSIM values obtained by all test methods. To sum up, the proposed approach outperforms other competing methods for CS image reconstruction by both randomly sampling and pseudo-radial sampling schemes.

Fig 6

Reconstructed images by pseudo-radial sampling.

From top to bottom for each column: Plane image, close-up of the red rectangle of Plane image, Barbara image, close-up of the red rectangle of Barbara image. From left to right for each row: the original image, the TV method, the NLR-CS-baseline method, the NLR-CS method, the BM3D-CS method, the proposed method.

Table 2

The PSNR (dB) and SSIM values of the restoration results obtained by the TV method [38], the NLR-CS-baseline method [22], the NLR-CS method [22], the BM3D-CS method [16] and the proposed method for noisy images.

Image	Number of measurements	Sampling scheme	Method	PSNR (dB)	SSIM
Head	M = 0.15N	random	TV	23.07	0.5489
			NLR-CS-baseline	27.15	0.6508
			NLR-CS	27.20	0.6515
			BM3D-CS	27.56	0.7203
			proposed	27.76	0.6615
Boat	M = 0.2N	random	TV	22.90	0.5313
			NLR-CS-baseline	27.60	0.7091
			NLR-CS	27.62	0.7093
			BM3D-CS	27.75	0.7281
			proposed	27.98	0.7397
Plane	M = 0.2N	pseudo-radial	TV	25.86	0.6410
			NLR-CS-baseline	28.64	0.8657
			NLR-CS	28.98	0.8866
			BM3D-CS	28.55	0.8203
			proposed	29.38	0.9014
Barbara	M = 0.29N	pseudo-radial	TV	24.09	0.6131
			NLR-CS-baseline	27.79	0.7721
			NLR-CS	28.12	0.7980
			BM3D-CS	28.62	0.8093
			proposed	29.20	0.8231

Conclusion

To better exploit the nonlocal sparsity of similar patches and non-convexity of rank minimization, in this paper, we use the non-convex SRF function surrogating the rank as a low-rank regularization for CS image recovery. This SRF function can better approximate the rank than the standard nuclear norm and the logdet function. We propose an efficient algorithm in the framework of alternative minimization method, which divides this CS problem into the SRF minimization subproblem and the least square subproblem. With respect to the minimization subproblem of the SRF function, we adopt the gradient descent method to solve it since it is differentiable. Simultaneously, the clear image is reconstructed by solving the least square subproblem using the ADMM method. Both noiseless and noisy numerical experiments demonstrate that the proposed approach achieves the better performance and vision quality under the lower sampling rate situation. In the future, we would like to explore better surrogates for the rank to improve performance, and solve other practical problems.