A Two-Step Spectral Gradient Projection Method for System of Nonlinear Monotone Equations and Image Deblurring Problems

: In this paper, we propose a two-step iterative algorithm based on projection technique for solving system of monotone nonlinear equations with convex constraints. The proposed two-step algorithm uses two search directions which are deﬁned using the well-known Barzilai and Borwein (BB) spectral parameters.The BB spectral parameters can be viewed as the approximations of Jacobians with scalar multiple of identity matrices. If the Jacobians are close to symmetric matrices with clustered eigenvalues then the BB parameters are expected to behave nicely. We present a new line search technique for generating the separating hyperplane projection step of Solodov and Svaiter (1998) that generalizes the one used in most of the existing literature. We establish the convergence result of the algorithm under some suitable assumptions. Preliminary numerical experiments demonstrate the efﬁciency and computational advantage of the algorithm over some existing algorithms designed for solving similar problems. Finally, we apply the proposed algorithm to solve image deblurring problem.


Introduction
Many problems arising from various applications such as optimization, differential equations, variational inequalities problems and so on, can be converted into nonlinear system of equations. Hence the study of iterative algorithms for solving nonlinear equations is of paramount importance especially when analytical method is not feasible or difficult to implement.
Let F : R n → R n be a monotone mapping and Λ be a subset of R n . We wish to find a point x * ∈ Λ, such that F(x * ) = 0.
The feasible set Λ is assumed to nonempty closed and convex. We call problem (1) system of nonlinear monotone equations with convex constraints. This problem appears as a subproblem in generalized proximal algorithms with Bregman distance [1]. In addition, some monotone variational inequality problems of finding y ∈ C for which x − y, F(y) ≥ 0, ∀ x ∈ C can be converted into systems of monotone equations [2]. Furthermore, l 1 − norm regularized optimization problems can be reformulated as monotone nonlinear equations [3]. Consider the following unconstrained optimization problem min x∈R n f (x), (2) where f : R n → R is assumed to be continuous, bounded below and its gradient, denoted by F, exists. Fermat's extremum theorem suggests that if a point x * is the local minimizer of the unconstrained optimization problem (2) then problem (1) holds. In addition, suppose x * is the minimizer of problem (2), then problem (1) is the first order necessary condition for the unconstrained optimization problem (2). This also underlines the importance of problem (1). Starting from a given initial point x 0 ∈ R n , popular iterative methods, such as Newton's method, quasi-Newton method, conjugate gradient method, for solving (2) use an updating rule defined as follows x k+1 = x k + α k d k , k = 0, 1, 2, ..., where α k and d k denote stepsize and search direction respectively. The search direction in (3) is usually defined as d k = −B −1 k F(x k ) where B k is either the exact Hessian matrix ∇ 2 f (x k ) in the case of Newton's method or the approximation of the Hessian matrix in the case of quasi-Newton method. The approximation of the Hessian matrix, B k , is required to satisfy the following secant equation s k−1 = x k − x k−1 and y k−1 = F(x k ) − F(x k−1 ). The Quasi-Newton method was developed to overcome one of the major shortcomings associated with the famous Newton's method which is the need to compute second derivative of the objective function in every iteration. However, it inherits the problem of storing n × n matrices throughout the iteration process which makes it unsuitable for large scale problems. One of the crucial approaches developed to overcome the storage problem of the quasi Newton method is the matrix-free method proposed by Barzilai and Borwein (BB) [4]. The BB method uses (3) to generate the next iterate with the search direction given by d k = −F(x k ) and the stepsize taken as diagonal matrix D k = τ k I which is supposed to satisfy the secant Equation (4). However, since τ k I produces diagonal matrices with identical diagonal elements, it is usually very difficult to find τ k for which D −1 k = τ −1 k I satisfies (4) when the dimension is greater than one. Consequently, Barzilai and Borwein required D −1 k to approximately satisfies (4) by finding τ k ∈ R which minimizes the following least square problems min τ τs k−1 − y k−1 2 , and min τ s k−1 − τy k−1 2 .
By Cauchy Schwarz inequality, we see that the stepsize produced by τ BB1 k is always greater than or equal to the one produced by τ BB2 k whenever y k−1 , s k−1 > 0. Barzilai and Borwein proved that the iterative scheme (3) with α k = τ BB1 k and d k = −F(x k ) converges with R-superlinear rate for two-dimensional strictly convex quadratic problems.
One disadvantage of the BB method, however, is that the stepsizes τ BB1 k and τ BB2 k may become negative if the objective function is not convex. Thus, Dai et al. [5] proposed and analyzed the following positive stepsize The stepsize (8) is the geometric mean of τ BB1 k and τ BB2 k . They showed that the iterative scheme (3) with α k = τ k has the same rate of convergence with the stepsize τ BB1 k under certain conditions for two-dimensional strictly convex quadratic functions. Recently, Dai et al. [6] proposed a family of gradient methods whose stepsize is a convex combination of τ BB1 k and τ BB2 k . The stepsize is obtained by solving the following problem It was shown that if 0 ≤ ξ ≤ 1 and y k−1 , s k−1 > 0, then . They proved that their method is R-superlinearly and R-linearly convergent for two-dimensional strictly convex quadratics and any finite dimensional case respectively. Convergence analysis of the BB stepsizes has been explored and interested reader may refer to the following References [7][8][9][10][11][12].
On the other hand, the BB method with the stepsize τ BB1 k has been extended to solve unconstrained nonlinear equations by La Cruz and Raydan [13]. Their algorithm is built on the strategy of nonmonotone line search technique which guarantees the global convergence of the method. Numerical experiments presented reveal their method competes with some well-established existing methods. However, their algorithm requires descent directions with respect to the squared norm of the residual. This means computation of a directional derivative, or its good approximation is needed at every iteration. Consequently, La Cruz et al. [14] proposed another BB method with a different nonmonotone line-search technique for solving unconstrained nonlinear equations. Their approach has advantage because unlike the former, the computations of directional derivatives are completely avoided. Based on the projection technique of Solodov and Svaiter [15], Zhang and Zhou [16] proposed an interesting projection spectral method which can be viewed as a modification of the method given in References [13,14]. They proposed a new line search strategy which does not require any merit function and takes the monotonicity of F into account. They established the global convergence of the method under some suitable assumptions and present some numerical experiments to demonstrate its computational advantage. In Reference [17], Yu et al. extended the method given by Zhang and Zhou [16] to solve monotone system of nonlinear equations with convex constraints. Their method is globally convergent under some conditions and preliminary numerical results show that the method works well and is more suitable compared to the projection method in Reference [18]. Recently, Mohammad and Abubakar [19] proposed a positive spectral method for unconstrained monotone nonlinear equations based on the projection technique in Reference [15]. The spectral parameter proposed is a convex combination of modified τ BB1 k and τ k . Their method works well and was extended to solve monotone nonlinear equations with convex constraints in Reference [20] as well as signal and image restoration in Reference [21].
Inspired by above contributions, we propose a two step iterative scheme based on the projection technique for solving system of monotone nonlinear equations with convex constraints. We define two search directions using Barzilai and Borwein (BB1 and BB2) spectral parameters with modifications. In addition, we investigate the efficiency of the propose algorithm in restoring blurred images. The symbols ·, · and · denote inner product and Euclidean norm respectively. The remaining part of this paper is organized as follows. In Section 2, we describe the proposed method and its global convergence. We report numerical experiments to show the efficiency of the algorithm in Section 3.
We describe the application of the proposed algorithm in Section 4 and give some conclusions as well as possible future research perspective in Section 5.

Two Step Iterative Scheme and Its Convergence Analysis
We begin this section with the following definition. Definition 1. Let x, y ∈ R n , a mapping F : R n → R n is said to be (ii) Lipschitzian continuous if there exists L > 0 such that From the discussions in the preceding section, we observe that all the methods use the one-step formula (3) to update their respective sequence of iterates. Let I be an identity map in R n , if we set d := (F − I), then formula (3) is closely related to the well-known Mann iterative scheme [22] where 0 ≤ α k < 1. Mann iteration has been applied to solve different kind of nonlinear problems successfully. However, its convergence speed is relatively slow. Different studies have shown that the famous two-step Ishikawa iterative scheme [23] v where α k , β k ∈ [0, 1), converges faster than the one-step Mann iteration.
Letd k = (F − I)u k andd k = F(v k ) − u k , then Ishikawa iterative scheme can be rewritten as follows v k = u k + α kdk , Based on the fact that the two step Ishikawa iterative scheme has faster convergence speed than the one-step Mann iterative scheme, in this paper, we propose a new two-step iterative scheme incorporating nonnegative BB parameters with projection strategy to solve monotone nonlinear equation with convex constraints. Given a starting point x 0 ∈ Λ and α k , β k ∈ (0, 1], we define the updating formula for the proposed two-step scheme as follows where z k = x k + β k d I I k , P Λ (·) is a projection operator defined below and For simplicity we denote λ I (x k ) := λ I k and λ I I (w k ) := λ I I k . The parameters λ I k and λ I I k are modifications of the BB parameters (7) where Assumption 1. Throughout this paper, we assume the following (i) The solution set of problem (1) is nonempty.
The following Lemma shows that the spectral parameters (17) are well-defined and bounded.
Proof of Lemma 1. The monotonicity of F gives F(x k+1 ) − F(x k ), x k+1 − x k ≥ 0. Therefore, by the definition of y I k and y I I k we have On the other hand, by (11) and Cauchy Schwarz inequality, we have Also since t > L, from (23) we can have Therefore, by (19) and (21) (20), (22) and (24) we have Remark 1. We give the following remarks (i) From Lemma 1, it is not difficult to see that the two search directions d I k and d I I k satisfy the descent condition. That is, The two search directions d I k and d I I k satisfy the following inequalities Next, we describe the projection operator in (15) which is usually used in iterative algorithms for solving problems such as fixed point problem, variational inequality problem, and so on. Let x ∈ R n and define an operator P Λ : R n → Λ by P Λ (x) = argmin{ x − y : y ∈ Λ}. The operator P Λ is called a projection onto the feasible set Λ and it enjoys the nonexpansive property, that is, If y ∈ Λ, then P Λ (y) = y and therefore, we have We now state the steps of the proposed algorithm which we call two-step spectral gradient method.

Remark 2.
We quickly note the following remarks (i) We claim that there exists a step-size β k satisfying the line search (1) for any k ≥ 0. Suppose on the contrary that there exists some k 0 such that for any i = 0, 1, 2, ..., the line search (1) is not satisfied, that is Since F is continuous and λ I I k is bounded for all k, letting i → ∞ yields It is clear that the inequality (29) cannot hold. Hence the line search (1) is well-defined. (ii) The line search defined by (1) is more general than that of Reference [24]. (iii) It follows from (15) and Assumption 1 that lim The next Lemma is very crucial to the convergence of Algorithm 1.
Step 1: Compute d I k using (16). if d I k = 0, then x k is a solution and the iteration process stops. end Step 2: Compute w k = x k + α k d I k .
Step 3: Compute d I I k using (16).

Lemma 2.
Let the Assumption 1 holds, then the sequences {w k }, {z k } and {x k } generated by Algorithm 1 are bounded. In addition, there exist some positive constants m 1 , m 2 and m 3 such that Furthermore, lim and lim Proof of Lemma 2. Let x * be a solution of problem (1), then by monotonicity of F, we have By the definition of x k+1 , and (33) we have This implies that x k − x * ≤ x 0 − x * for all k, and therefore the sequence {x k } is bounded and lim It follows from (26), that d I k ≤ µm 1 and d I I k ≤ γm 1 . It further follows from (15) that {w k } is bounded. By Lipschitzian continuity of F, there exists m 2 > 0 such that F(w k ) ≤ m 2 .
Since {d I I k } is bounded, it follows from the definition z k that {z k } is also bounded. By Lipschitzian continuity of F, there exists some constant m 3 for which Since the stepsize β k in Step 4 of Algorithm 1 satisfies β k ≤ 1, ∀ k, then from (1), we have Combining with (34) gives By (36) and (37), we have Taking limit gives Hence, it holds that lim This together with the definition of z k in Step 5 of Algorithm 1 yields By the property of projection (27), we have Theorem 1. Let {x k } be the sequence generated by Algorithm 1. Suppose that Assumption 1 holds, then the sequence {x k } converges to a point x * which satisfies F(x * ) = 0.
Proof of Theorem 1. We begin by proving that Suppose on the contrary that (42) does not hold, then there exists q > 0 for which If β k = κ, since Algorithm 1 uses a backtracking process to compute β k starting from κ, then −1 β k does not satisfy (1), that is, Consequently, we have from Remark 1 (i), where F(x k + −1 β k d I I k ) is bounded above by a positive constant m 5 . This means .
Taking limit on both sides as k → ∞, we have This contradicts (39). Hence (42) must hold. Now, since F is continuous and the sequence {x k } is bounded, then there is some accumulation From the proof of Lemma 2, we know that lim k→∞ x k − x * exists. Therefore, we can conclude that lim k→∞ x k − x * = 0 and the proof is complete.

Numerical Results and Comparison
Attention is now turn to numerical experiments. The experiment is divided into two parts. The first experiment aims to explore the role of the parameter c in the definition of the line search (1). On the other hand, the second experiment discusses the computational advantage of the proposed method in comparison with two existing methods. The two existing methods are: (i) Spectral gradient projection method for monotone nonlinear equations with convex constraints proposed by Yu et al. [17]. (ii) Two spectral gradient projection methods for constrained equations and their linear convergence rate proposed by Liu and Duan [25]. This method has two algorithms i.e., Algorithm 2.1 and Algorithm 2.2. We only compare our proposed method with Algorithm 2.1 since Algorithm 2.2 is similar with that Yu et al. [17].
These two methods were chosen because their search directions are defined based on the BB parameters. For convenience, we respectively denote the two methods by SGPM and TSGP. Algorithm 1 TSSP is implemented using the following parameters κ = 1, σ = 0.01, = 0.5, r = t = 0.01, and α k = 1 (k+1) 2 . The parameters used for the SGPM and TSGP methods were taken respectively from References [17] and [25]. The metrics used for the comparison are: number of iterations (ITER), number of function evaluations (FVAL) and CPU time (TIME). In the course of the experiments, we solved six benchmark test problems using six (6) different starting points (see Table 1) by varying the number of dimension. The test problems are denoted by Pi, i = 1, 2, 3, 4, 5, 6. Since the proposed algorithm is derivative-free, the test problems include two nonsmooth problems. The three solvers were coded in MATLAB R2017a and run on a PC with intel Core(TM) i5-8250u processor with 4 GB of RAM and CPU 1.60 GHZ. The MATLAB code for the TSSP algorithm is available in https://github.com/aliyumagsu/ TSSP_Algorithm. The iteration process is terminated whenever the inequality F(x k ) ≤ 10 −6 or F(z k ) ≤ 10 −6 is satisfied and failure is declared whenever the number of iterations exceeds 1000 and the terminating criterion mentioned above has not been satisfied.
First experiment. This experiment discusses the role of the parameter c in the definition of the line search (1) with regards to the performance of the TSSP algorithm. We solved all the test Problems 1-6 with dimension n = 1000, using all the given initial guesses in Table 1 Table 2. CPU time results are omitted in Table 2 because virtually all are less than 1 s. The results obtained reveal that the parameter c slightly affected the performance of TSSP algorithm when solving Problems 2 and 6. For problem 2, Algorithm 1 TSSP recorded least ITER and FVAL when c = 4 and 5 while different ITER and FVAL values recorded for different values of c may be associated with the random starting points chosen independently by MATLAB. However, extensive numerical experiment is needed to investigate the role of the parameter c in the performance of the TSSP algorithm.
Second experiment. This experiment presents the computational advantage of the proposed method in comparison with the two existing methods mentioned above based on ITER, FVAL and TIME. All the test problems 1 − 6 were solved using the starting points in Table 1 with three (3) different dimensions n = 1000, 50,000 and 100,000. In this experiment, we take c = 2. The results obtained by each solver are reported in Tables 3 and 4. The NORM results presented in Tables 3 and 4 show that each solver successfully obtained solutions of all the test Problems 1-6. However, it is clear that the TSSP algorithm obtained the solutions of virtually all the test problems with least ITER, FVAL and TIME. These information are summarized in Figures 1-3 based on the Dolan and Moré performance profile [26]. This performance profile tells the percentage win by each solver. In all the experiments, we see from Figures 1-3 that the proposed TSSP algorithm performs better with higher percentage win based on ITER, FVAL and TIME for solving all the test problems. In fact, the TSSP algorithm recorded 100 percent least FVAL for all the test problems.

Applications in Image Deblurring
In this section, we apply the proposed Algorithm 1 to solve problems arising from compressive sensing, particularly image deblurring. Consider the following least square problem with 1 − regularization term where x ∈ R n is the underlying images, y ∈ R k is the observed images, E ∈ R k×n (k << n), linear operator, is an m × n blurring matrix, and the parameter µ > 0. Problem (46) is of great importance because it appears in many areas of applications arising from compressive sensing. Recently, problem (46) has been investigated by many researchers and different kinds of iterative algorithms have been proposed in the literature [3,[32][33][34][35]. Many algorithms for solving (46) fall into two categories namely: algorithms that required differentiability assumption and algorithms that are derivative free. Since 1 −norm is a nonsmooth function, algorithms that require the assumption of differentiabilty are not suitable for problem (46) in its original form. Consequently, either x 1 is approximated with some smooth function or problem (46) is reformulated into an equivalent problem. For instance, Figueiredo et al. [3] translate problem (46) into convex quadratic program as follows.
For any x ∈ R n , we can find some vectors, say u, v ∈ R n such that where u i = max{0, x i }, v i = max{0, −x i } for all i = 1, 2, ..., n. Thus, we can write x 1 = e T n (u + v), where e n is an n−dimensional vector with all elements one. Therefore, we can rewrite problem (46) as Furthermore, if we let q = [u v] T , then from Reference [3], we can write (47) as the following min q 1 2 q T Gq + c T q, where . It is not difficult to see that G is a positive semi-definite matrix.
In Reference [36], the resulted constrained minimization problem (48) is further translated into the following linear variable inequality problem Find q ∈ R n such that Since the feasible region of q is R n , problem (49) is equivalent to the following linear complementary problem Find q ∈ R n such that We can see that the point q ∈ R n is a solution of the above linear complementary problem (50) if and only if it satisfies the following system of nonlinear equation The mapping F is a vector-valued and the "min" operator denotes the componentwise minimum of two vectors. Interestingly, Lemma 3 of Reference [37] and Lemma 2.2 of Reference [36] show that the mapping F satisfies Assumption 1 (ii) i.e., is Lipschitzian continuity and monotonicity. Therefore our proposed TSSP algorithm can be applied to solve it.

Image Deblurring Experiment
We tested the performance of the two-step TSSP algorithm in restoring some blurred images in comparison with the one-step spectral gradient method for 1 problems in compressed sensing (SGCS) [36]. The images used for the experiment are the well-known gray test images namely: Lena, House, Pepper, Camera man and Barbara where the size of each image is given in Table 5. The following metrics are used to assess the performance and quality of restoration by each algorithm tested: number of iteration (ITER), CPU time in seconds (TIME), signal-to-noise-ratio (SNR) which is defined as SNR = 20 × log 10 x x −x , and the structural similarity (SSIM) index that measure the similarity between the original image and the restored image [38] for each of the two experiments. The MATLAB implementation of the SSIM index can be obtained at http://www.cns.nyu.edu/~lcv/ssim/. To achieve fairness in comparison, each code was run from same initial point x 0 = E T y and terminate when where f k is the merit function evaluation at x k , with f (x k ) = 1 2 y − Ex k 2 2 + µ x k 1 . The parameters used for both TSSP and SGCS in this experiment come from Reference [36] except for c = 1 in the line search (1) and α k = (0.999 k ) × (10 5 + F(x 0 ) 2 ).
The original, blurred and restored images by each algorithm are given in Figure 4. The two tested algorithms restored the blurred images successfully with different speed and level of quality. The results of the restoration by each algorithm are reported in Table 5. We see from Table 5 that TSSP restored all the five images with less ITER. Taking TIME into consideration, we see that though the SGCS is faster in restoring two of the images (i.e., Camera man and Barbara), TSSP is faster in restoring the remaining three images (i.e., Lena, House and Pepper). In addition, the SNR and SSIM values recorded by each algorithm revealed that TSSP restored the five blurred images with slightly better quality than SGCS except for Camera man. Taking everything together, this experiment shows that the two-step TSSP can deal with the 1 regularization problems effectively and can be a favourable alternative for image reconstruction.  Figure 4. The original images (first row), the blurred images (second row), the restored images by methods TSSP (third row) and SGCS (last row).

Conclusions
This paper presents an efficient derivative-free iterative algorithm called TSSP for nonlinear monotone equations. It utilizes a two-step approach that incorporates the well-known BB parameters with a projection strategy. We showed that the TSSP converges globally under the Lipschitzian and monotonicity assumptions. Also, we proposed a new line search that is more general than the one proposed by Cheng in Reference [24], able to include the line search by Cheng as a special case. Preliminary numerical results reported in Table 2 shows that the parameter c introduced in the new line search defined by (1) may have some effect on the performance of the proposed algorithm. Numerical results presented revealed that the proposed TSSP algorithm has computational advantage and performs better than the two existing algorithms in References [17,25]. These results indicate that the two-step BB-like algorithm is superior to the existing one-step BB-like algorithms, especially on solving nonlinear equations. It is worth emphasizing that the TSSP algorithm improves existing results on monotone nonlinear equations. The results obtained in Section 3 show that the TSSP algorithm possessed excellent numerical performances with evidence of efficiently solving all the test problems considered with minimal number of iterations and function evaluations. The numerical results reported from the experiments of deblurring two-dimensional images from their limited measurements have shown that the two-step algorithm TSSP competes favorably with the one-step SGCS algorithm. Future work includes the extension of TSSP algorithm on different forms of optimization frameworks and applications such as nonlinear least-squares problems [39,40], neural networks [41,42] and machine learning [43,44].