Convergence Gain in Compressive Deconvolution: Application to Medical Ultrasound Imaging

: The compressive deconvolution (CD) problem represents a class of efficient models that is appealing in high-resolution ultrasound image reconstruction. In this paper, we focus on designing an improved CD method based on the framework of a strictly contractive Peaceman–Rechford splitting method (sc-PRSM). By fully excavating the special structure of ultrasound image reconstruction, the improved CD method is easier to implement by partially linearizing the quadratic term of subproblems in the CD problem. The resulting subproblems can obtain closed-form solutions. The convergence of the improved CD method with partial linearization is guaranteed by employing a customized relaxation factor. We establish the global convergence for the new method. The performance of the method is verified via several experiments implemented in realistic synthetic data and in vivo ultrasound images.


Introduction
Ultrasound imaging has become a very important medical imaging scheme, as it is noninvasive, harmless, and cost-effective compared with computed tomography (CT) and nuclear magnetic resonance (MRI).
Moreover, ultrasound requires little energy, which makes ultrasound imaging a good candidate for hand-held applications. An extremely promising application scenario for ultrasound imaging is breast cancer localization, where cancer regions have statistically different acoustic properties compared with benign areas [1].
Compressive sensing (CS) can accelerate the acquisition rate without decreasing the reconstructed signal quality and maintain the image quality with fewer data. Based on CS theory, one can effectively implement the reconstruction under the condition of the restricted isometry property (RIP) [2,3]. The applicability of CS, such as 2D and 3D ultrasound imaging [4] or duplex Doppler [5], has attracted an increasing number of researchers to propose new theories and methods. Among them, several ultrasound imaging devices [6][7][8][9][10] have been proposed in which acquisition and compression are processed simultaneously. The main idea is to combine CS and deconvolution into ultrasound imaging and form the compressive deconvolution problem as follows: where y ∈ R M involves M linear measurements acquired by projecting one RF image Hx ∈ R N onto the CS acquisition matrix Φ ∈ R M×N , with M N. H ∈ R N×N represents a block circulant with a circulant block (BCCB) matrix modeling the 2D convolution between the 2D point spread function (PSF) of the ultrasound system, and the tissue reflectivity function (TRF), ν ∈ R M , represents a zero-mean additive white Gaussian noise.
Some sequential approaches such as YALL1 [11] are prone to first reformulate (1) into an unconstrained optimization model: When the blurred RF imager = Ψa is restored, the TRF x will be inferred by minimizing: The above two-step manipulation is the most intuitive way but inevitably yields larger estimation errors as illustrated in [12]. As such, Chen et al. [13] investigate the following model by solving CS and deconvolution problems simultaneously: where q denotes the shape parameter of the generalized Gaussian distribution. (4) can be solved by alternatively solving three easier sub-problems based on alternating direction method of multipliers (ADMM): 2 , γ is a Lipschitz parameter, and prox represents the proximal operator [14]. For the above three sub-problems, w k+1 has a closed solution and x k+1 has an approximate closed solution, but a k+1 is not easy to compute, even though the effective Sherman-Morrison-Woolbury inversion manipulation can be introduced. For non-orthogonal sparse basis Ψ, approximation algorithms such as Newton's method are quite time-consuming [15]. In terms of saving such computing time, some accelerating strategies and new iterative schemes need to be investigated and explored. Meanwhile, since the compressive deconvolution (CD) method is an inexact ADMM by solving x k+1 approximately, it is extremely important to present a strict and complete convergence analysis.
Motivated by these observations and based on the semi-proximal Peaceman-Rechard splitting method in our previous work [16], in this paper, we propose an improved CD method, referred to as semi-proximal symmetric ADMM. Define υ = Ψa, and the improved compressive deconvolution (ICD) method has the following iterative scheme: Compared with the CD method presented in [13], the contributions of this article can be summarized as follows: 1.
Based on the iterative scheme of the strictly semi-proximal Peaceman-Rechard splitting method, we present an ICD method that will reduce the number of iterations while only involving additional dual update (i.e., λ 1+ 1 2 ) and requiring almost the same computational effort for each iteration.

2.
We prove that the ICD method will converge under mild conditions, while the convergence analysis is not given in the previous CD method. 3.
We introduce some elaborate manipulations that can directly generalize the CD method to more general scenarios with a non-orthogonal sparse basis Ψ.
The rest of this paper is organized as follows. In Section 2, we first present some preliminaries that are useful for subsequent analysis. Then, we illustrate the ICD method to rebuild the sparse coefficients from the measurements of ultrasound imaging, and the convergence analysis is given. In Section 3, we give extensive ultrasound experiments that can be used to evaluate the performance of the proposed reconstruction algorithm in comparison with CD algorithm. Finally, we make some concluding remarks in Section 4. (4) In this section, inspired by He and Yuan's approach [17], we equivalently convert the convex minimization model expressed by Equation (4) to a variational form. It makes sense to perform such reformulation, because convergence analysis becomes more concise under the variational model.

Variational Reformulation of Equation
For succedent analysis of the proposed algorithm, let us denote . Then the ultrasound imaging model can be reformulated as The Lagrangian function and augmented Lagrangian function of Equation (7) can be, respectively, expressed as and where λ ∈ R N represents the Lagrangian multiplier. Then hunting for a saddle point of L(z 1 , That is, for any (z 1 , z 2 , λ), we have Then, resolving Equation (4) is equivalent to seeking w = (z * 1 , z * 2 , λ * ) such that Especially, the mapping F(w) defined in Equation (13) is affine with a skew-symmetric matrix, it is monotone. We express by Ω * the solution set of VI(Ω, F, θ).

Notations
We denote the 2-norm of a vector by · and let z 2 G = z T Gz for z ∈ R N and G ∈ R N×N . For a real symmetric matrix S, S 0 (S 0) represents S, which is positive semidefinite (positive definite). For ease of the analysis, we define the following matrices as and Below we prove three assertions regarding the matrices just defined. These assertions make it possible to present our convergence analysis for the new algorithm compactly with alleviated notation. Lemma 1. Given R 0, let β > 0, µ, ρ ∈ (0, 1), and ∈ (0, 1]. The matrices H, M, and Q defined, respectively, in Equations (15) and (16) have the relationships as follows: and Proof. We consider two cases.
With the matrices H, M, and Q at hand, we easily obtain The second assertion HM = Q is proved. Consequently, we have Using Equations (15) and (16), and the above equation, we have Note that β > 0 and ρ, ∈ (0, 1). Thus, for any w = (z 1 , z 2 , λ) = 0, we have Therefore, the matrix G is positive definite. (II). ρ ∈ (0, 1) and = 1. Note that Thus, it is positive definite, and G is only positive semi-definite. Here, we would emphasize that we do not require the positive definiteness of G. Instead, positive semi-definiteness of G is enough for our algorithmic analysis.

Algorithm
In this section, we will present our new algorithm to solve Equation (4). However, we first present the iterative scheme by using the standard strictly contractive Peaceman-Rechford splitting method with two different relaxation factors: By introducing a customized proximal term, especially for ultrasound imaging, our improved compressive deconvolution (ICD) method has the iterative scheme: is a customized semi-definite matrix. Note that the equivalence of Equations (25b)-(25d) and Equations (6c)-(6e) is evident, while the relationship between the closed-form solution expressed by Equation (25a) and Equations (6a) and (6b) is not evident. We illustrate the latter in Appendix A.

Global Convergence
To make the analysis more elegant, we reformulate ICD Equations (25a)-(25d) into the form Now we analyze the convergence for our proposed ICD method expressed by Equation (26). We prove its global convergence from the contraction perspective. In order to further alleviate the notation in our analysis, we define an auxiliary sequencew k as where (z k+1 1 , z k+1 2 ) is produced by Equations (26a) and (26c), and we immediately have Moreover, we have the following relationship: which can be reformulated as a compact form under the notation of w k andw k : where M is defined in Equation (15). Now we start to prove some properties for the sequence {w k } defined in Equation (27). We are interested in estimating how accurate the pointw k is to a solution point w * of VI (F, Ω, θ). The main result is proved in Theorem 1. Now, we try to find a lower bound in terms of the discrepancy between w − w k+1 2 H and w − w k 2 H for any w ∈ Ω.
Theorem 1. Let {w k } be generated by Equation (26) and let {w k } be defined in Equation (27). Let H and G be defined in Equations (19) and (18), respectively. Then, for any w ∈ Ω, we have Proof. Refer to Appendix B.
The next lemma demonstrates the contraction property of the sequence {w k } generated by Equation (26).

Lemma 2.
Let {w k } be generated by Equation (26) with 0 < ρ < 1 and 0 < < 1, and let {w k } be defined in Equation (27). Let H and G be defined in Equations (18) and (19), respectively. Then, for any w * ∈ Ω * , we have Lemma 3. Let the sequence {w k } be generated by Equation (26) with 0 < ρ < 1 and = 1. Then we have With the above lemmas, we can finally obtain the global convergence theorem of ICD method for solving VI(Ω, F, θ) as follows: Theorem 2. The sequence {w k } generated by Equation (26) converges to some w ∞ that is a solution of VI(Ω, F, θ).

Simulated US Images
Two ultrasound data sets, named group1 and group2, were obtained by 2D convolution between spatially invariant PSFs and the TRFs [15]. The results were quantitatively evaluated in terms of peak signal-to-noise ratio (PSNR), image structural similarity SSIM [18], improvement in SNR (ISNR), and the normalized root mean square error (NRMSE). The metrics are expressed as follows: where x, y,x are respectively the original image, the RF image, and the reconstructed image. L denotes the maximum intensity value in x. σ x and σx are the mean and variance values of x andx; c 1 = (k 1 C) 2 and c 2 = (k 2 C) 2 are two variables stabilizing the division with a weak denominator. C represents the dynamic range of the pixel-values and k 1 , k 2 denote constants. Herein, k 1 = 0.01, k 2 = 0.03, and C = 1.

In Vivo US Images
We consider two real in vivo US images [19]: (a) Mouse bladder: The observed image is with the size 400 × 256 as shown in Figure 1. The number of homogeneous regions was set to K = 3 in this experiment, which is sufficient to represent the anatomical structures of the image. (b) Skin melanoma: The second in vivo image (of size 400 × 298) represents a skin melanoma tumor, as shown in Figure 2. Since the ground truth of the TRF and the label map are not available for in vivo US data, the quality of the deconvolution results is evaluated using the contrast-to-noise ratio (CNR) [20]: where µ 1 and µ 2 are the mean of pixels located in two regions extracted from the image, while and are the standard deviations of the same blocks. All code was written in MATLAB and performed on a ThinkPad computer equipped with Windows 7, 2.60 GHz and 2 GB of memory. Based on the same stopping criterion,

Results
The quantitative results reported in Figures 1-8 confirm that, given the same maximum number of iteration (we set ItMax = 1 × 10 4 ), the ICD method can achieve better reconstruction quality gain than CD method for both simulated and real US images. Moreover, based on Figures 9 and 10, ICD method converges faster than the existing CD method for all CS ratios. We should remark that, for the group2 data set, the metric of the SINR achieved via the CD method is negative. This may be because the CD method is unstable under low CS ratios.

Discussions and Conclusions
In this paper, we proposed an improved compressive deconvolution method by introducing two different parameters in updating the dual variable to improve its convergence rate. We established the relationship between the two parameters under which we proved the global convergence of the algorithm. The ultrasound simulations show that the proposed method can achieve reconstruction US image with better quality gain under the same maximum number of iteration. Moreover, the ICD method has a much faster convergence rate compared with the conventional compressive deconvolution method. It should be noted that parameters such as ρ and are based only on empirical values, and a much greater deconvolution gain can be obtained if the parameters are adaptively optimized. This is the direction in which our research will continue.
Author Contributions: B.G. and S.X. proposed the architecture of ICD method. B.G. and X.L. derived the convergence analysis of ICD method. The simulations were implemented by B.G. and X.L. B.G. and X.L. wrote the manuscript. L.Z. and K.P. read and approved the final manuscript. The above equation is exactly Equation (6b), so the proof is complete.