{'fraction_non_alphanumeric': 0.04541393686718285, 'fraction_numerical': 0.02132966051220965, 'mean_word_length': 4.055513737297704, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fixed Point Algorithm Based on Quasi-Newton Method for Convex Minimization Problem with Application to Image Deblurring

Solving an optimization problem whose objective function is the sum of two convex functions has received considerable interests in the context of image processing recently. In particular, we are interested in the scenario when a non-differentiable convex function such as the total variation (TV) norm is included in the objective function due to many variational models established in image processing have this nature. In this paper, we propose a fast fixed point algorithm based on the quasi-Newton method for solving this class of problem, and apply it in the field of TV-based image deblurring. The novel method is derived from the idea of the quasi-Newton method, and the fixed-point algorithms based on the proximity operator, which were widely investigated very recently. Utilizing the non-expansion property of the proximity operator we further investigate the global convergence of the proposed algorithm. Numerical experiments on image deblurring problem with additive or multiplicative noise are presented to demonstrate that the proposed algorithm is superior to the recently developed fixed-point algorithm in the computational efficiency.


Introduction
The general convex optimization problems that arise in image processing take the form of a sum of two convex functions. Often one function is the data fidelity energy term that is decided by the noise type and one wants to minimize, and the other function is the regularization term to make the solution have certain properties. For instance, the usual l 1 based regularization is used to obtain the sparse solution in the fields such as image restoration [1] and compressed sensing [2,3]. In this paper, we propose an efficient fixed-point algorithm to solve the optimization problem whose objective function is composed of two convex function, i.e., min u f 1 (Bu) + f 2 (u) (1.1) where f 1 : R M → R and f 2 : R N → R are convex function, B : R N → R M is a linear transform, and f 2 is differentiable with a 1/β-Lipschitz continuous gradient, i.e., for any v, w ∈ R N . Despite its simplicity, many variational models in image processing can be formulated in the form of (1.1). For example, the classical total variation (TV) or wavelet sparsity prior based models [1], which are often considered in image restoration under Gaussian noise, have the simple form as follows.
where K ∈ R N×N is a linear blurring operator, b ∈ R N and µ > 0 is a regularization parameter.
Here B denotes the sparse transform such as the gradient operator and the wavelet basis, and the image prior is imposed by using the term Bu 1 , which promotes the sparsity of image under the transform B. With f 1 = µ · 1 and f 2 = 1 2 K · −b 2 2 the problem (1.3) can be seen as a special case of (1.1). It is observed that the main difficulty of solving problem (1.1) is that the function f 1 is non-differentiable.
In the last several years, many optimization algorithms have been developed for efficiently solving the variational models in image processing. The iterative shrinkage/thresholding (IST) algorithm is one of the most successful methods. Consider the general minimization problem min u f (u) + g(u) (1.4) where f and g are convex function, and g is differentiable. The classical IST algorithm for problem (1.4) is given by the following iterative formula u k+1 = T γ f (u k − γ∇g(u k )) = p γ (u k ) (1.5) where γ > 0 is a step parameter.
Here T γ f is called the thresholding operator. It corresponds to the proximity operator prox γ f , which is defined by [5] prox γ f (v) = arg min (1. 6) In different literatures, IST is also called iterative denoising method [4], Landweber iteration [6], proximal forward-back splitting (PFBS) algorithm [7] or fixed-point continuation (FPC) algorithm [8]. In order to further accelerate the convergence speed, many new iterative shrinkage algorithms based on the IST, which include the SpaRSA (Sparse Reconstruction by Separable Approximation) [9], TwIST (Two-step IST) [10], FISTA (Fast iterative shrinkage-thresholding algorithm) [11] were further proposed. Notice that a proximity operator is needed to be computed in each iteration of the iterative shrinkage algorithms. However, the proximity operators prox γ f for the general case of f = f 1 • B often have no closed solutions. For example, if we choose B = ∇ and f = ∇ · 1 , then the minimization problem of (1.6) is just the Rudin-Osher-Fatemi (ROF) denoising problem whose solution cannot be obtained easily. Therefore, inner iterative algorithm is needed for computing the proximity operators in most cases.
In recent years, a class of algorithms based on the splitting methods have been developed and shown to be efficient for computing the proximity operator. For instance, Goldstein and Osher [12] proposed a splitting algorithm based on the Bregman iteration, called the split-Bregman 2 method, to compute the solution of the minimization problem of (1.6) especially for the case of ROF denoising. This algorithm can be successfully applied for solving the general minimization problem (1.1), and theoretically it has been proved to be equivalent to the Douglas-Rachford splitting (DRS) algorithm [13,14] and the alternating direction of multiplier method (ADMM) [15,16]. Although the split-Bregman framework has been shown to be very useful, a sub-minimization problem of solving the system of linear or nonlinear equations is included in each iteration and may time-consuming sometimes. Very recently, alternating direction minimization methods based on the linearized technique [17,18] have been widely investigated to overcome this and further improve the efficiency. Another class of methods is the primal-dual algorithms. Chambolle [19] firstly proposed a dual algorithm for the ROF denoising. Later on, Zhu et al. [29] devised a primal-dual hybrid gradient (PDHG) method, which alternately update the primal and dual variables by the gradient descent scheme and gradient ascend scheme. The theoretical analysis on variants of the PDHG algorithm, and on the connection with the linearized version or variants of ADMM were widely investigated to bridge the gap between different methods. Refer to [17,[21][22][23][24] and the references cited therein for details.
In this paper, we focus our attention on a new class of algorithms that has been developed very recently from the view of fixed-point. In [25], Jia and Zhao proposed a fast algorithm for the ROF denoising by simplify the original split-Bregman framework. Motivated by this idea, Micchelli et al. [26] designed a fixed-point algorithm based on proximity operator (named FP 2 O) for computing prox f 1 •B , which was proved to be more efficient than the splitting methods. Later on, several variants of fixed-point algorithms were proposed for special cases of image restoration. For instance, Micchelli et al. [27,28] further extended the FP 2 O algorithm to solve TV-L1 denoising model where f 2 = 1 2 K · −b 1 in (1.1). Chen et al. [29] proposed a proximity operator based algorithm for solving indicator functions based l 1 -norm minimization problems with application to compressed sample. Krol et al. [30] proposed a preconditioned alternating projection algorithm for emission computed tomography (ECT) restoration, where a diagonal preconditioning matrix is used in the devised fixed-point algorithm. The extension of the FP 2 O algorithm to the more general case of the form of (1.1) has also been investigated very recently [31,32]. Specifically, a primal-dual fixed point algorithm which combines the PFBS algorithm and only one inner iteration of FP 2 O has been proposed in [33].
However, in the previous fixed-point algorithms, we observe that all the iterative formulas are composed of the gradient descent algorithm and the proximity point algorithm. It is well known that the gradient-based algorithms typically have a sub-linear convergence rate, while the Newton method or the quasi-Newton method has been presented with a super-linear convergence rate. This fact motivates us to propose a new fixed-point algorithm which combines the quasi-Newton method and the proximity operator algorithm. Furthermore, the global convergence of the proposed algorithm is investigated under certain assumption.
The rest of this paper is organized in four sections. In section 2 we briefly review the existing fixed-point algorithms based on the proximal operator, and further propose a fixed-point algorithm based on quasi-Newton method. In section 3 the global convergence of the proposed algorithm is further investigated from the point of the view of fixed point theory under certain conditions. The numerical examples on deblurring problem of images contaminated by additive Gaussian noise and multiplicative noise are reported in section 4. The results there demonstrate that FP 2 O κ QN is superior to the recently proposed PDFP 2 O κ in the context of image deblurring.

Existing fixed-point algorithms based on the proximal operator
Motivated by the fast algorithm proposed for the ROF denoising in the literature [25], Micchelli et al. [26] designed a fixed-point algorithm named FP 2 O for the computation of the proximity operator prox f 1 •B (x) for any x ∈ R N . Denote λ max (BB T ) be the largest eigenvalue of BB T . Choose the parameter 0 < λ < 2/λ max , and define the operator Then we can obtain the fixed-point iterative scheme which is just called FP 2 O algorithm as follows where S κ is the κ-averaged operator of S , i.e., S κ = κI + (1 − κ)S for any κ ∈ (0, 1). Calculate the fixed-point v * of the operator S by the formula (2.2), and hence obtain that The key technique served as the foundation of FP 2 O algorithm is the relationship between the proximity operator and the subdifferential of a convex function, as described in (3.2) below. FP 2 O algorithm supplies a simple and efficient method of solving (1.1) with the special case of f 2 (u) = 1 2 u − x 2 2 in the classical framework of fixed-point iteration. In [32], this algorithm has been extended to the more general case that ∇ f 2 (u) is bijective and the inverse can be computed easily. In particular, choose Define the operatorŜ Then the corresponding fixed-point iterative scheme for (2.4) is given by (2.6) and the solution u * of (2.4) can be obtained by the formula where v * is the fixed-point of the operatorŜ κ . In order to deal with the general case of f 2 , the authors in [31] also combined FP 2 O and PFBS algorithms, and proposed a new algorithm named PFBS FP 2 O in which the proximity operator in the PFBS algorithm is calculated by using FP 2 O, i.e., is calculated by FP 2 O. Notice that a inner iteration of solving prox γ f 1 •B is included in PFBS FP 2 O, and it is problematic to set the approximate iteration number to balance the computational time and precision. In order to solve this issue, a primal-dual fixed points algorithm based on proximity operator (PDFP 2 O) [33] was proposed very recently. In this algorithm, instead of implementing FP 2 O for many iteration steps to calculate prox γ f 1 •B (x) in PFBS FP 2 O, only one inner fixed-point iteration is adopted. Suppose κ = 0. Then we can obtain the following iteration scheme (PDFP 2 O): It is obvious that u is the primal variable related to (1.1), and according to the thorough study in [33] we know that the variable v is just the dual variable of the primal-dual form related to (1.1). Therefore, PDFP 2 O also belongs to the class of primal-dual algorithm framework. Similarly to FP 2 O, a relaxation parameter κ ∈ (0, 1) can be introduced to get the algorithm named PDFP 2 O κ . For more details refer to [33].

Proposed fixed-point algorithm based on quasi-Newton method
In the PDFP 2 O algorithm, the iterative formulas consist of the proximity operator and the gradient descent algorithm. Since the Newton-type methods have been shown to have a faster convergence rate compared to the gradient-based methods, a very nature idea is to use the Newtontype methods instead of the gradient descent step in the fixed-point algorithm. Consider the minimization problem (1.1). We use the second-order Taylor expansion of the convex function f 2 (u) at the recent iterative point u k instead of it, i.e., where Q k is a positive definite symmetric matrix to approximate the second derivative ∇ 2 f 2 . Then (1.1) can be reformulated as It is observed that (2.9) corresponds to the minimization problem (2.4) with x = Q k u k − ∇ f 2 (u k ). Therefore, the next iteration scheme u k+1 can be obtained by the fixed point iteration algorithm shown in (2.5)-(2.7). In order to avoid any inner iterations, we use only one inner fixed point iteration in the proposed algorithm. For this, choose κ = 0, and define Using the numerical solution v k for S k as the initial value, and only implementing one iteration of solving the fixed-point of S k+1 (v), we can obtain the following iteration scheme (2.10) It is observed that an intermediate iterative variable u k+1/2 is generated by a quasi-Newton method. Therefore, we called the proposed algorithm a fixed-point algorithm based on quasi-Newton method, and abbreviate it as FP 2 O QN, which is described as Algorithm 1 below. For simplification of convergence analysis below, we set Q k ≡ Q to be unchanged with different k.
Similarly to the literatures [26,32,33], we can introduce a relaxation parameter κ ∈ (0, 1] to obtain Algorithm 2, which is exactly the Picard iterates with the parameter.

Convergence analysis
Let us start with some related notations and conclusions which will serve as the foundation for the proof below.

Definition 3.1. (Nonexpansive operator) A nonlinear operator T
By the application of the Cauchy-Schwarz inequality it is easy to show that a firmly nonexpansive operator is also nonexpansive. [34])For a given initial point v 0 ∈ R M and an operator P :

Definition 3.2. (Picard sequence
For the Picard sequence we have the following conclusion. [34]) Let C be a closed convex set in R M and let P : C → C be a nonexpansive mapping with at least one fixed point. Then for any v 0 ∈ R M and any κ ∈ (0, 1), the Picard sequence of P κ = κI + (1 − κ)P converges to a fixed point of P.
For any convex function f : The following result illustrates the relationship between the proximity operator and the subdifferential of a convex function. This conclusion has appeared in many previous literatures, such as [26,32,33].
In what follows, we establish a fixed-point formulation for the solution of the minimization problem (1.1) based on the conclusion in Proposition 3.4. To this end, we define the operator and the operator T 2 : where λ is a positive parameter. Denote the operator T : is a fixed point of T , then u * is a solution of the minimization problem (1.1).

Proof:
Since u * is one solution of the minimization problem (1.1), by the first-order optimality condition we have Then we obtain that Besides, using Proposition 3.4 we have (3.7) Inserting (3.6) into (3.7) we further obtain that Based on (3.6) and (3.8) we infer that (v * , u * ) is a fixed point of T . Conversely, if there exists w * = (v * , u * ) satisfying w * = T (w * ), we can derive that u * satisfies the first-order optimality condition of (1.1) by the equivalent formulas above. Therefore, we conclude that u * is a solution of (1.1).
According to the formulas in (3.3) and (3.4) we find out that the iterative scheme of FP 2 O QN can be reformulated as which is also equal to w k+1 = T (w k ) with w k = (v k , u k ). This implies that the sequence {v k , u k } generated by FP 2 O QN is just the Picard sequence of the operator T . With the similar discussion we can find that the iterative formulas of FP 2 O κ QN is equal to w k+1 = Tκ(w k ), i.e., the sequence {v k , u k } generated by FP 2 O κ QN is the Picard sequence of the operator T κ . Based on Theorem 3.5 we know that the solution of the minimization problem (1.1) is just equal to the fixed point of the operator T . Therefore, the convergence of FP 2 O κ QN can be guaranteed by verifying the nonexpansion of T according to Proposition 3.3. The proof here is similar to those presented in [33]. However, note that the global convergence of proposed algorithms cannot be directly obtained by applying results in [33], and hence it is included here for completion.
In the following, we give a crucial inequality for showing the nonexpansion of T . Denote Here we assume that 0 < λ ≤ 1/λ max (BQ −1 B T ), and hence M is a symmetric positive semidefinite matrix. Therefore, we can define the semi-norm u M = √ u T Mu, and then define the norm Lemma 3.6. For any two points w 1 = (v 1 , u 1 ) and w 2 = (v 2 , u 2 ) in R M × R N , the following inequality comes into existence.
Proof: According to Lemma 2.4 of [7] we know that I − prox 1 λ f is firmly nonexpansive, and hence obtain Following the definition in (3.4) we also have (3.11) According to (3.10) and (3.11) we further have (3.12) Since f 2 has 1/β-Lipschitz continuous gradient, we get that (3.13) By the definition of M we easily get (3.14) Based on (3.12)-(3.14) we can obtain (3.9) directly. From the results in Lemma 3.6 we know that T is nonexpansive with the norm of · λ,Q . Therefore, we are able to prove the convergence of FP 2 O κ QN according to Proposition 3.3, which is described as follows.

Theorem 3.7. Assume that Q −1
Proof: Note that the solution of (1.1) is just one fixed point of T . From Lemma 3.6 we know that the operator T is nonexpansive, maps the set R M × R N to itself, and has at least one fixed point. According to Opial κ-averaged Theorem, we conclude that, for any w 0 and κ ∈ (0, 1), the Picard sequence of T κ converges to a fixed point of T . With this result we further infer that u k converges to the solution of (1.1).
Finally, we process with the convergence of FP 2 O QN based on the inequality in Lemma 3.6. Proof: Let w * = (v * , u * ) ∈ R M × R N be a fixed point of T . Substitute w 1 and w 2 in (3.9) with w k = (v k , u k ) and w * , we obtain that Summing (3.15) from some k 0 to +∞ we obtain that According to (3.18) we can easily deduce that lim Besides, from (3.6) we know that Q −1 ∇ f 2 (u * ) + λQ −1 B T v * = 0, and hence obtain that Based on (3.21) we immediately get According to (3.15) we know that the sequence w k − w * 2 λ,Q is non-increasing, and hence w k is bounded, which implies that there exists a convergent subsequence of w k j such that for some pointŵ * = (v * ,û * ). Due to the operator T is continuous, we further have lim based on (3.24) and (3.25). Therefore,ŵ * is a fixed point of the operator T . Due to the proof is started with any fixed point w * , we can set w * =ŵ * . In this case, we see that the sequence w k −ŵ * λ,Q is non-increasing. Combining this with the formula (3.25) we infer that lim k→+∞ w k =ŵ * . (3.28) From Theorem 3.5 we know thatû * is a solution of (1.1). This completes the proof.

Numerical examples
In this section, we will compare the proposed FP 2 O κ QN algorithm with PDFP 2 O κ [33] through the experiments of image restoration. Here two cases of additive and multiplicative noise types are considered. One is the additive Gaussian noise which has been extensively investigated over the last decades. In this setting, the data fidelity term can be formulated as where K is the blurring operator, and b is the observed image. The other is the speckle noise which also appears in many real world image processing applications such as laser imaging, synthetic aperture radar (SAR) imaging and ultrasonic imaging. In [35], this speckle noise followed by a Rayleigh distribution is investigated. Under this condition, the observed image can be modeled as corrupted with signal-dependent noise of this form where ε is a zero-mean Gaussian noise with standard deviation σ, i.e., ε ∼ N(0, σ). Based on the model (4.1) and the characteristics of Gaussian distribution, the corresponding fidelity term can be formulated as In the following experiments, we use total variation as the regularization term, and hence choose the function f 1 (Bu) = µ ∇u 1 where ∇ : R N → R 2N is a discrete gradient operator. Here we adopt the isotropic definition of total variation, and the proximity operator prox 1 λ f 1 can be computed easily. For more details refer to [32].

Gaussian image deblurring
In this subsection, we choose three gray-scale images, Cameraman, Barbara (with size of 256 × 256), and Boat (with size of 512 × 512) as the original images, and evaluate FP 2 O κ QN in four typical image blurring scenarios: strong blur with low noise; strong blur with medium noise; mild blur with low noise; mild blur with medium noise, which are summarized in Table 1 (σ and σ a denote the standard deviation). In the following, we discuss the selection of the parameters µ, λ and κ in both fixed point algorithms, the parameter γ in PDFP 2 O κ , and the matrix Q in FP 2 O κ QN. Due to f 2 (u) = 1 2 Ku − b 2 2 , we have ∇ 2 f 2 ≡ K T K. Therefore, we can choose Q = K T K. However, the blurring operator K is ill-posed generally, and (K T K) −1 cannot be used in the proposed algorithm due to the instability. Therefore, we choose Q = K T K + ǫB T B in our experiments. Here ǫ is a small positive number and B is a difference matrix. Notice that the introduction of the term ǫB T B avoids the ill-posed condition, and Q −1 can also be computed efficiently by fast Fourier transforms (FFTs) with periodic boundary conditions.
The regularization parameter µ is decided by the noise level, and the adjustment of the parameters λ and γ does influence the convergence speed and stability of the fixed point algorithms. Through many trials we use the rules of thumb: µ is set to 0.06 and 0.15 for σ = 1.5 and 3.0 respectively; λ is set to 0.125; γ is chosen to be 1.8 for PDFP 2 O κ ; and ǫ = 0.1 for FP 2 O κ QN. Similarly to the literatures [32,33], we find that κ = 0 achieves the best convergence speed compared with other κ ∈ (0, 1), and hence we choose κ = 0 for both algorithms.
The performance of the restored images of the compared algorithms is measured quantitatively by means of the peak signal-to-noise ratio (PSNR), which is defined by PSNR(u,ū) = 10 lg where u and u * denote the original image and the restored image respectively. The stopping criterion for the fixed-point algorithms is defined such that the relative error is below some small constant, i.e., where tol denotes a prescribed tolerance value. In our experiments we choose tol = 5 × 10 −4 .
The PSNR values for the deblurred images, the number of iterations, and the CPU time are listed in Table 2. In this table, the four image blurring scenarios shown in Table 1 are considered, and (·, ·, ·) represents the PSNR values, iteration numbers and CPU time in sequence. From these results we observe that the recovered images obtained by FP 2 O κ QN can achieve better PSNRs than those given by PDFP 2 O κ , and meanwhile, the corresponding iteration number and running time of FP 2 O κ QN is less than those of PDFP 2 O κ . Figures 1-3 show the recovery results of the PDFP 2 O κ and FP 2 O κ QN algorithms. It is observed that the visual qualities of images obtained by both algorithms are more or less the same.  Figure 4 shows the evolution curves of PSNR values obtained by both fixed point algorithms for two cases including in Table 2: one is the Cameraman image blurred by 8 × 8 box average kernel and added with Gaussain noise with σ = 3, the other is the Boat image blurred by 6 × 6 gaussian kernel and added with Gaussain noise with σ = 1.5. From the plots we can implicitly find that FP 2 O κ QN achieves the best solution (with higher PSNRs) much faster than PDFP 2 O κ .

Rayleigh image deblurring
In this subsection, we further discuss the case of images contaminated by Rayleigh noise. The corresponding minimization problem has been introduced above. The two blur kernels shown in Table 1, and Rayleigh noise with σ = 0.5 and 1.0 are considered here.
First of all, we illustrate the setting of the parameters in both fixed point algorithms. Since for any ξ ∈ R N . Notice that the value of ∇ 2 f 2 changes with the iteration number, and the inverse of ∇ 2 f 2 is difficult to be estimated. Therefore, we use Q = βK T K + ǫ to approximate ∇ 2 f 2 in 13    the proposed fixed point algorithm. Here the parameter β is used to replace the unknown b 2 (Kξ) 3 , and the term ǫB T B is included to avoid the ill-posed condition. In the following experiments, we find that β = 0.25 and ǫ = 0.005 are two suitable selection through many trials. Moreover, for parameters in both algorithms we use the following rules of thumb: the regularization parameter µ is chosen to be 0.01 and 0.02 for the noise level of σ = 0.5 and 1.0 respectively; λ is set to 0.125; γ is chosen to be 15.0 for PDFP 2 O κ . We also find out that the selection of κ = 0 is suitable for our experiments here.
In what follows, two images, Pepper (with size of 256×256) and Cameraman, are used for our test. Figures 5-6 show the evolution curves of PSNR (dB) running both fixed point algorithms for the two images. From the plots we observe that the PSNR values obtained by FP 2 O κ QN increase much faster than those by PDFP 2 O κ . This is due to the quasi-Newton method included in FP 2 O κ QN is more efficient than the gradient descent algorithm involved in PDFP 2 O κ . Figure 7 shows the deblurred results of Pepper image convoluted by 8 × 8 box average kernel and contaminated by Rayleigh noise with σ = 1.0. The corresponding PSNR values, the number of iterations, and the CPU time are also included. It is observed that FP 2 O κ QN can obtain higher PSNR values with less iteration number and running time compared to PDFP 2 O κ . The recovery results of Cameraman image blurred by 6 × 6 gaussian kernel and corrupted by Rayleigh noise with σ = 1.0 are also presented in Figure 8. We also observe that FP 2 O κ QN is more efficient than PDFP 2 O κ , especially in the implementation time.

Conclusion
In this article, we propose a fast fixed point algorithm based on the quasi-Newton method, abbreviated as FP 2 O κ QN, for solving the minimization problems with the general form of min u∈R N { f 1 (Bu) + f 2 (u)}. The main distinction between FP 2 O κ QN and previous fixed point algorithms lies in that the quasi-Newton method, rather than the gradient descent algorithm, is included in the algorithm framework. The proposed algorithm framework is applied to solve TV-based image restoration problem. Numerical experiments reported in this paper indicate that FP 2 O κ QN outperform the recently proposed PDFP 2 O κ , especially in the implementation time.

Acknowledgments
The research was supported in part by the National Natural Science Foundation of China under Grant 61271014.