# A Modified Conjugate Gradient Method for Solving Large-Scale Nonlinear Equations

Authors: Hongbo Guan; Sheng Wang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, we propose a modified Polak–Ribière–Polyak (PRP) conjugate gradient method for solving large-scale nonlinear equations. Under weaker conditions, we show that the proposed method is globally convergent. We also carry out some numerical experiments to test the proposed method. The results show that the proposed method is efficient and stable.


---

## Body

## 1. Introduction

Solving nonlinear equations is an important problem which appears in various models of science and engineering such as computer vision, computational geometry, signal processing, computational chemistry, and robotics. More specifically, the subproblems in the generalized proximal algorithms with Bergman distances is a monotone nonlinear equations [1], and l1-norm regularized optimization problems can be reformulated as monotone nonlinear equations [2]. Due to its wide applications, the studies in the numerical methods for solving the monotone nonlinear equations have received much attention [3–10]. In this paper, we are interested in the numerical methods for solving monotone nonlinear equations with convex constraints:(1)Fx=0,x∈S,where F:ℝn⟶ℝn is a continuous function and S is a nonempty, closed, and convex set. The monotonicity of the mapping F means that(2)Fx−FyTx−y≥0,∀x,y∈ℝn.The methods for solving monotone nonlinear equations (1) are closely relevant to the methods for solving the following optimization problems:(3)minx∈ℝnfx.Notice that iffx is strictly convex, then ∇fx is strictly monotone which means ∇fx−∇fyTx−y≥0. It is well known that the strictly convex function exists a unique solution x∗, satisfying ∇fx∗=0. To sum up, if there is a convex function fx satisfying ∇fx=Fx, then solving the optimization problems (3) is equivalent to solving monotone nonlinear equations (1). So, a natural idea to solve monotone nonlinear equations (1) is to use the existing efficient methods for solving optimization problems (3). There are many methods for solving optimization problems (3), such as the Newton method, quasi-Newton method, trust region method, and conjugate gradient method. Among these methods, the conjugate gradient method is a very effective method for solving optimization problems (3) due to their simplicity and low storage. A conjugate gradient method generates a sequence of iterates:(4)xk+1=xk+αkdk,k=0,1,…,where αk is the step length and direction dk is defined by(5)dk=−gk,ifk=0,−gk+βkdk−1,ifk≥1,where βk is a parameter and gk is the gradient of the objective function fx. The choice of βk determines different conjugate gradient methods [11–17]. We are interested in the PRP conjugate gradient method in which the parameter βk is defined by(6)βkPRP=gkTyk−1gk−12,where yk−1=gk−gk−1. Based on the idea of [18, 19], Zhang et al. [20] proposed a new modified nonlinear PRP method in which βk is defined by(7)βkNPRP=gkTyk−1zk−1+ηyk−12zk−12gkTdk−1,where zk−1=maxtdk−1,gk−12, t>0 and η>1/4 are two constants. There is a mistake about the definition of βkNPRP. By this definition, we will not be able to prove Lemma 1 in [20]. It should be(8)βkNPRP=gkTyk−1zk−1−ηyk−12zk−12gkTdk−1.There are many conjugate gradient methods for solving nonlinear equations (1). Zhang and Zhou [4] proposed a spectral gradient method by combining the modified spectral gradient and projection method, which can be applied to solve nonsmooth equations. Xiao and Zhou [10] extended the CG-DESCENT method to solve large-scale nonlinear monotone equations and extended this method to decode a sparse signal in compressive sensing. Dai and Zhu [21] proposed a derivative-free method for solving large-scale nonlinear monotone equations and proposed a new line search for the derivative-free method. Other related works can be found [3, 5–8, 10, 22–30]. In this paper, we combined the projection method [3], the modified nonlinear PRP conjugate gradient method for unconstrained optimization [20] and the iterative method [10] and proposed a modified nonlinear conjugate gradient method for solving large-scale nonlinear monotone equations with convex constrains.This paper is organized as follows. In Section2, we propose a modified nonlinear PRP method for solving monotone nonlinear equations with convex constraints. Under reasonable conditions, we prove its global convergence. In Section 3, we make some improvement to the proposed method and give the convergence theorem of the improved method. In Section 4, we do some numerical experiments to test the proposed methods. The results show that our methods are efficient and promising. Furthermore, we use the proposed methods to solve practical problems in compressed sensing.


## 2. A Modified Nonlinear PRP Method

In this section, we develop a modified nonlinear PRP method for solving the nonlinear equations with convex constraints. Based on the modified nonlinear PRP method [20], we now introduce our method for solving (1). Inspired by (8), we define dk as(9)dk=−Fk,ifk=0,−Fk+βkNPRPdk−1,ifk≥1,where yk−1=Fk−Fk−1. The parameter βkNPRP is computed as(10)βkNPRP=FkTyk−1zk−1−ηyk−12zk−12FkTdk−1,where zk−1=maxtdk−1,Fk−12,t>0 and η>1/4 are two constants.The lemma below shows a good property ofdk. The steps of the method are given in Algorithm 1.Lemma 1.
Letdk be generated by Algorithm 1. If zk−1≠0, then there exists a constant c>0 such that(11)FkTdk≤−cFk2.Proof.
Fork=0, we have(12)F0Td0≤−F02.
Fork≥1, we obtain(13)FkTdk=FkT−Fk+βkNPRPdk−1,=−Fk2+FkTyk−1zk−1−ηyk−12zk−12FkTdk−1FkTdk−1,=−Fk2+zk−1FkTdk−1FkTyk−1−ηyk−12FkTdk−12zk−12.
Denote(14)uk=12ηzk−1Fk,vk=2ηFkTdk−1yk−1.
Then, we obtain(15)FkTdk=−Fk2+ukTvk−1/2vk2zk−12=−1−14ηFk2+ukTvk−1/2uk2+vk2zk−12≤−1−14ηFk2.
Letc=−1−1/4η; then, inequality (11) is satisfied.
Next, we establish the global convergence of the proposed method. Without specification, we always suppose that the solution set of equation (1) is nonempty and the following assumption holds.Assumption 1.
(i)
The mappingF is Lipchitz continuous, and it means that the mapping F satisfies(16)Fx−Fy≤Lx−y,∀x,y∈S.(ii)
The projection operatorPS⋅ is nonexpansive, i.e.,(17)PSx−PSy≤x−y,∀x,y∈ℝn.Lemma 2.
Suppose that Assumption1 holds and x∗ is a solution of (1), and the sequence xk and x¯k+1 are generated by Algorithm 1. Then, the sequence xk, x¯k+1, and Fk are bounded.Algorithm 1: Modified NPRP method.
Initial. Given a small constant ε>0 and constants t>0,ξ>0,η>1/4,σ,ρ∈0,1. Choose an initial point x0∈S. Let k=0.Step 1. Stop if Fxk≤ε.Step 2. Compute (9) to get dk.Step 3. Let tk=maxξρi:i=0,1,2… satisfying−Fxk+tkdkTdk≥σtkdk2,denotex¯k+1=xk+tkdk.Step 4. Compute xk+1=PSxk−αkFx¯k+1,whereαk=Fx¯k+1Txk−x¯k+1/Fx¯k+12,andPS⋅ is a projection operator, defined by PSx=argminy−x:y∈S,∀x∈Rn.Step 5. Let k=k+1 and go to Step 1.Proof.
We first show thatxk is bounded. From the monotonicity of function F, we have(18)Fx¯k+1Txk−x∗=Fx¯k+1Txk−x¯k+1+Fx¯k+1−Fx∗Tx¯k+1−x∗,≥Fx¯k+1Txk−x¯k+1.
It is easy to see that(19)xk+1−x∗2=PSxk−αkFx¯k+1−PSx∗2≤xk−αkFx¯k+1−x∗2≤xk−x∗2−2αkFx¯k+1Txk−x¯k+1+αk2Fx¯k+12=xk−x∗2−Fx¯k+1Txk−x¯k+12Fx¯k+12≤xk−x∗2.
The last inequality implies(20)xk+1−x∗2≤xk−x∗2≤⋯x0−x∗2.
It obviously that the sequencesxk is bounded, i.e., there is a constant M>0 such that(21)xk≤M,∀k>0.
Next, we show thatFk is bounded. Since Fx is Lipchitz continuous, we obtain(22)Fxk=Fxk−Fx∗≤Lxk−x∗≤Lx0−x∗.
LetA=Lx0−x∗; then, Fkis bounded. That is,(23)Fk≤A.
At last, we prove that the sequencex¯k+1 is bounded. From (2) and (23) and Algorithm 1, we obtain(24)Fx¯k+1Txk−x¯k+1=−tkFx¯k+1Tdk≥σtk2dk2=σxk−x¯k+12,Fx¯k+1Txk−x¯k+1=Fx¯k+1−FxkTxk−x¯k+1+FxkTxk−x¯k+1,≤Fkxk−x¯k+1≤Axk−x¯k+1.
So, the following inequality holds:(25)xk−x¯k+1≤Aσ.
It implies that the sequencex¯k+1 is bounded.Lemma 3.
Suppose that Assumption1 holds, and the sequence xk and Fk is generated by Algorithm 1. Then, there exists a constant γ>1 such that(26)cFk≤dk≤γFk,∀k≥0.Proof.
We first prove the right side of inequality (26). For k=0, from (9), we have(27)dk=Fk.
Fork≥1, by the definition of βkNPRP and (21), we obtain(28)βkNPRP=FkTyk−1zk−1−ηyk−12zk−12FkTdk−1≤Fkyk−1tdk−1+ηyk−12t2dk−12Fkdk−1≤2LMFktdk−1+η2LM2Fkt2dk−1=2LMt+2ηLMt2Fkdk−1.
By the definition ofdk (9) and the last inequality, we obtain(29)dk≤Fk+βkNPRPdk−1≤1+2LMt+2ηLMt2Fk.
Letγ=1+2LMt+2ηLM/t2; then, we have dk≤γFk.
Now, we turn to prove the left side of the inequality. It follows from (11) that(30)cFk2≤−FkTdk≤Fkdk.
Therefore, we have(31)dk≥cFk.Lemma 4.
Suppose Assumption1 holds; then, the step length tk satisfies(32)tk≥minξ,ρcL+σγ2.Proof.
If the Algorithm1 terminates in a finite number of steps, then there is a k∈ℝ such that xk is a solution of equation (1) and Fk=0. From now on, we assume that Fk≠0, for any k. It is easy to see that dk≠0 from (11).
Iftk≠ξ, by the line search process, we know that ρ−1tk does not satisfy Algorithm 1, that is,(33)−Fxk+ρ−1tkdkTdk≤σρ−1tkdk2.
It follows from (11) and Assumption 1 that(34)cFk2≤−FkTdk,=Fxk+ρ−1tkdk−FxkTdk+Fxk+ρ−1tkdkTdk≤Lρ−1tkdk2+σρ−1tkdk2=ρ−1L+σtkdk2.
From the last inequality and Lemma3, we obtain(35)tk≥ρcL+σFk2dk2≥ρcγ2L+σ.
Hence, it holds that(36)tk≥minξ,ρcL+σγ2.Lemma 5.
Suppose that Assumption1 holds, and the sequence xk and dk are generated by Algorithm 1. Then, we have the following:(37)limk⟶∞tkdk=0.Proof.
Noticed that(38)xk+1−x∗2≤xk−x∗2−Fx¯k+1Txk−x¯k+12Fx¯k+12,≤xk−x∗2−σ2xk−x¯k+14Fx¯k+12.
Since the functionFx is continuous, and the sequence x¯k+1 is bounded, so the sequence Fx¯k+1 is bounded. That is, for all k≥0, there exists a positive constant B>0, such that Fx¯k+1≤B. Then, we obtain(39)∑k=0∞xk−x¯k+14≤∑k=0∞B2σ2xk−x∗2−xk+1−x∗2<+∞.
So, we have(40)limk⟶∞tkdk=limk⟶∞xk−x¯k+1=0.Theorem 1.
Suppose that Assumption1 holds. The sequence xk is generated by Algorithm 1. Then, we have(41)limk⟶∞infFk=0.Proof.
Suppose that (41) does not hold; then, there exists ϵ>0 such that, for any k≥0,(42)Fk≥ε.
From (26) and the last inequality, it is easy to see(43)cε≤dk≤γA.
From (41) and (42), we obtain(44)tkdk≥minξ,ρcγ2L+σdk≥minξcε,ρc2εγ2L+σ.
The last inequality yields a contradiction with (37), so (41) is satisfied.


## 3. An Improvement

In this section, we make some improvement to the modified nonlinear PRP method proposed in Section2. In Algorithm 1, we take the step length αk=Fx¯k+1Txk−x¯k+1/Fx¯k+12. Is there a better choice for αk? This is our purpose to improve Algorithm 1. Under the condition of ensuring the convergence of the algorithm and the related good properties and results, we improve Algorithm 1 in order to get better numerical results.From Algorithm1, to make the inequality xk+1−x∗≤xk−x∗ hold, we only need to satisfy(45)ϕα=−2αkFx¯k+1Txk−x¯k+1+αk2Fx¯k+12≤0.By solving the last inequality, we have(46)0≤αk≤2Fx¯k+1Txk−x¯k+1Fx¯k+12.It is easy to see thatαk=Fx¯k+1Txk−x¯k+1/Fx¯k+12 is the minimum point of the function ϕα. This is the reason why Algorithm 1 takes αk=Fx¯k+1Txk−x¯k+1/Fx¯k+12. Under reasonable conditions, we hope to get a large step length than Algorithm 1. So, we obtain(47)Fx¯k+1Txk−x¯k+1Fx¯k+12≤αk≤2Fx¯k+1Txk−x¯k+1Fx¯k+12.Based on the above arguments, we propose an improved algorithm of Algorithm1. In the improved algorithm, we make the step length:(48)αk=θFx¯k+1Txk−x¯k+1Fx¯k+12,θ∈1,2.Similar to the proof of Theorem1, we have the following results.Theorem 2.
Suppose that Assumption1 holds. The sequence xk is generated by Algorithm 2; then, we have(49)limk⟶∞infFk=0.
The iterative process of the improved method is stated as follows.Algorithm 2: The improved modified nonlinear PRP method.
Initial. Given a small constant ε>0 and constants t>0,ξ>0,η>1/4,σ,ρ∈0,1. θ∈1,2.Choose an initial point x0∈S. Let k=0.Step 1. Stop if Fxk≤ε.Step 2. Compute (9) to get dk.Step 3. Let tk=maxξρi:i=0,1,2,… satisfying−Fxk+tkdkTdk≥σtkdk2,denotex¯k+1=xk+tkdk.Step 4. Computexk+1=PSxk−αkFx¯k+1,whereαk=θFx¯k+1Txk−x¯k+1/Fx¯k+12,θ∈1,2.Step 5. Let k=k+1 and go to Step 1.


## 4. Numerical Results

In this section, we do some numerical experiments to test the performance of the proposed methods. We implemented our methods in MATLAB R2020b and run the codes on a personal computer with 2.3 GHz CPU and 16 GB RAM.We first solve Problems1 and 2.Problem 1.
(see [8]). The mapping F is taken as Fx=f1x,f2x,…,fnxT, where(50)fix=exi−1,i=1,2,…,n,S=ℝ+n=x∈ℝn|xi≥0,i=1,2,…,n.Problem 2.
(see [4]). The mapping F is taken as Fx=f1x,f2x,…,fnxT, where(51)fix=2xi−sinxi,i=1,2,…,n,and S=ℝ+n.
We setξ=1,t=1,ρ=0.95,η=1,σ=0.0001,andθ=2. The stopping criterion of the algorithm is set to Fxk<10−5 or the number of iteration reach to 500. The latter case means that the method is a failure for the test problems. We test both problems with the dimensions of variables n=1000,10000,100000,1000000, and 100000. Start from different initial points, and we list all results in Tables 1 and 2. We compared the performance of the proposed methods with the classical Newton method and an efficient algorithm CGD [10] in the total number of iterations as well as the computational time. The meaning of each column is given below.‘Init’: the initial pointn’: the dimension of the problem‘Iter’: the total number of iterations‘Time’: the CPU time (in seconds) used for the method‘NM’: the Newton method‘CGD’: the conjugate gradient method in [10]‘MNPRP’: the modified nonlinear PRP method‘IMNPRP’: the improved modified nonlinear PRP method
The results in Tables1 and 2 show that our methods performs very well both in the number of iterations and CPU time. The IMNPRP performs best among these methods. It is worth noting that the number of iterations does not increase significantly as n increases. Hence, the proposed method is very suitable for solving large-scale problems. Because of the lack of memory, the dimension of the problems solved by the Newton method is no more than 100,000.
The following example is a signal reconstruction problem from compressed sensing.Table 1
Comparison on Problem1.
InitnNMCGDMNPRPIMNPRPIterTimeIterTimeIterTimeIterTime1,1,…,1T100050.0345100.000640.000310.00021000052.0438110.003040.001510.0009100000Out of memory110.022940.011910.00731000000Out of memory120.400250.159510.092110000000Out of memory134.111251.746810.97062,2,…,2T100060.0325120.000550.000410.00031000062.4640130.004150.003010.0023100000Out of memory130.029960.029910.01751000000Out of memory140.454760.283810.197010000000Out of memory155.140073.403112.08761/2,1/2,…,1/2T100040.0226100.000450.000210.00011000041.6051110.003150.001210.0006100000Out of memory120.024750.008810.00401000000Out of memory120.375760.136810.053110000000Out of memory134.118461.510710.58422,2,…,2T100050.0335100.000440.000210.00011000052.0122110.002740.001510.0010100000Out of memory110.018740.011810.00711000000Out of memory120.379750.157510.088210000000Out of memory134.221751.763811.0510Table 2
Comparison on Problem2.
InitnNMCGDMNPRPIMNPRPIterTimeIterTimeIterTimeIterTime1,1,…,1T100030.0187110.000640.000210.00011000031.1675120.003340.000910.0004100000Out of memory120.024540.005810.00281000000Out of memory130.415650.103610.035010000000Out of memory144.359551.137910.36502,2,…,2T100040.0207110.000550.000310.00021000041.6028120.003650.001610.0009100000Out of memory120.031260.012410.00611000000Out of memory130.415360.162810.077110000000Out of memory144.527561.794010.83001/2,1/2,…,1/2T100030.0226110.000940.000210.00011000031.1711110.003040.000710.0002100000Out of memory120.024840.004910.00151000000Out of memory130.399850.093610.023610000000Out of memory144.380150.991110.26972,2,…,2T100030.0188110.000540.000210.00011000031.1775110.002940.000810.0003100000Out of memory120.023440.004410.00231000000Out of memory130.396440.082810.032010000000Out of memory144.269151.149110.3430Problem 3.
(see [10]). Consider a typical compressive sensing scenario, where we aim to reconstruct a length-n sparse signal form m observations m≪n. In this test, the measurement b contains noise:(52)b=Ax¯+ω,where ω is the Gaussian noise distributed as N0,σ2I with σ2=10−4. The random A is the Gaussian matrix which is generated by command randn (m, n) in Matlab. The merit function(53)fx=τx1+12Ax−b22,where the value τ is forced to decrease as the measure in [31]. The iterative process starts at the measurement image, i.e., x0=ATb, and terminates when the relative change between successive iterates falls below 10−4, i.e.,(54)Tol=fk−fk−1fk−1<10−4,where fk denotes the function value at iterationxk. By the discussion in [10], we know that the l1-norm problem can transformed a monotone nonlinear equation. Hence, it can be solved by Algorithms 1 and 2.
Due to the storage limitations of the PC, we test a small size signal withn=211andm=29, and the original contains 26 randomly nonzero elements. The quality of restoration is measured by the mean of squared error (MSE) to the original signal x¯, that is,(55)MSE=1nx¯−x2,where x is the restored signal. We take the parameters ξ=10,σ=10−4, and ρ=0.5 in CGD, MNPRP, and IMNPRP.
In order to test the effectiveness of the proposed methods, we compare the proposed methods with the CGD method [10] and the solver SGCS which is specially designed to solve monotone equations for recovering a large sparse signal in compressive sensing. The results are listed in Figures 1 and 2.
It can be seen from Figures1 and 2 that all methods have recovered the original sparse signal almost exactly. Among these methods, the IMNPRP method performs best.Figure 1
The original signal (top), the measurement (bottom), and the reconstructed signals by SGCS, CGD, MNPRP, and IMNPRP.Figure 2
Comparison results on Problem3. The y-axes represent the MSE (a, b) and the function values (c, d).
(a)(b)(c)(d)


## 5. Conclusions

In this paper, a modified conjugate gradient method and its improved method are proposed for solving the large-scale nonlinear equations. Under some assumptions, global convergence of the proposed methods are established. Numerical results show that the proposed methods are very efficient and competitive.



---

*Source: 9919595-2021-06-24.xml*