# A Derivative-Free Conjugate Gradient Method and Its Global Convergence for Solving Symmetric Nonlinear Equations

Authors: Mohammed Yusuf Waziri; Jamilu Sabi’u
Journal: International Journal of Mathematics and Mathematical Sciences
 (2015)

---

## Abstract

We suggest a conjugate gradient (CG) method for solving symmetric systems of nonlinear equations without computing Jacobian and gradient via the special structure of the underlying function. This derivative-free feature of the proposed method gives it advantage to solve relatively large-scale problems (500,000 variables) with lower storage requirement compared to some existing methods. Under appropriate conditions, the global convergence of our method is reported. Numerical results on some benchmark test problems show that the proposed method is practically effective.


---

## Body

## 1. Introduction

Let us consider the systems of nonlinear equations:(1)Fx=0,where F:Rn→Rn is a nonlinear mapping. Often, the mapping, F, is assumed to satisfy the following assumptions:(A1)
There existsx∗∈Rn s.t. F(x∗)=0.(A2)
F is a continuously differentiable mapping in a neighborhood of x∗.(A3)
F′(x∗) is invertible.(A4)
The JacobianF′(x) is symmetric.The prominent method for finding the solution of (1) is the classical Newton’s method which generates a sequence of iterates {xk} from a given initial point x0 via(2)xk+1=xk-F′xk-1Fxk,where k=0,1,2,…. The attractive features of this method are rapid convergence and easy implementation. Nevertheless, Newton’s method requires the computation of the Jacobian matrix, which require the first-order derivative of the systems. In practice, computations of some functions derivatives are quite costly and sometime they are not available or could not be done precisely. In this case Newton’s method cannot be applied directly.In this work, we are interested in handling large-scale problems for which the Jacobian either is not available or requires a low amount of storage; the best method is CG approach. It is vital to mention that the conjugate gradient methods are among the popular used methods for unconstrained optimization problems. They are particularly efficient for handling large-scale problems due to their convergence properties, simple implementation, and low storage [1]. Notwithstanding, the study of conjugate gradient methods for large-scale symmetric nonlinear systems of equations is scanty, and this is what motivated us to have this paper.In general, CG methods for solving nonlinear systems of equations generate iterative points{xk} from initial given point x0 using(3)xk+1=xk+αkdk,where αk>0 is attained via line search and direction dk is obtained using(4)dk=-Fxk,ifk=0,-Fxk+βkdk,ifk≥1,where βk is termed as conjugate gradient parameter.These problems which are under study may arise from an unconstrained optimization problem, a saddle point problem, Karush-Kuhn-Tucker (KKT) of equality constrained optimization problem, the discretized two-point boundary value problem, the discretized elliptic boundary value problem, and so forth.Equation (1) is the first-order necessary condition for the unconstrained optimization problem when F is the gradient mapping of some function f:Rn→R,(5)min⁡fx,x∈Rn.For the equality constrained problem,(6)min⁡fx,s.t.hz=0,where h is a vector-valued function, the KKT conditions can be represented as system (1) with x=(z,v), and(7)Fz,v=∇Fz+∇hzv,hz,where v is the vector of Lagrange multipliers. Notice that the Jacobian ∇F(z,v) is symmetric for all (z,v) (see, e.g., [2]).Problem (1) can be converted to the following global optimization problem (5) with our function f defined by(8)fx=12Fx2.A large number of efficient solvers for large-scale symmetric nonlinear equations have been proposed, analyzed, and tested in the last decade. Among them, the most classic one entirely due to Li and Fukushima [3], in which a Gauss-Newton-based BFGS method is developed, and the global and superlinear convergence are also established. Subsequently, its performance is further improved by Gu et al. [4], where norm descent BFGS methods are designed. Since then, norm descent type BFGS methods especially cooperating with trust regions strategy are presented in the literature and showed their moderate effectiveness experimentally [5]. Still the matrix storage and solving of n-linear systems of equations are required in the BFGS type methods presented in the literature. The recent designed nonmonotone spectral gradient algorithm [6] falls within the framework of matrix free methods.The conjugate gradient methods for symmetric nonlinear equations have received a good attention and take an appropriate progress. However, Li and Wang [7] proposed a modified Fletcher-Reeves conjugate gradient method which is based on the work of Zhang et al. [8], and the results illustrate that their proposed conjugate gradient method is promising. In line with this development, further studies on conjugate gradient are inspired for solving large-scale symmetric nonlinear equations. Zhou and Shen [9] extended the descent three-term polak-Ribiere-Polyak of Zhang et al. [10] for solving (1) by combining with the work of Li and Fukushima [3]. Meanwhile the classic polak-Ribiere-Polyak is successfully used to solve symmetric equation (1) by Zhou and Shen [1].Subsequently Xiao et al. [11] proposed a method based on well-known conjugate gradient of Hager and Zhang [12], and the proposed method converges globally. Extensive numerical experiments showed that each of the above-mentioned methods performs quite well. The combination of conjugate gradient algorithms and the Newton method, for the first time, was presented by Andrei [13, 14]. Hence, in this paper, we intended to present a new enhanced CG parameter βk which is matrix- and derivative-free, respectively. This is made possible by combining Birgin and Martínez direction with classical Newton direction.We organized the paper as follows. In the next section, we present the details of the proposed method. Convergence results are presented in Section3. Some numerical results are reported in Section 4. Finally, conclusions are made in Section 5.


## 2. Derivation of the Method

In this section we present a new CG parameterβk, as a result of combining Birgin and Martínez direction with classical Newton direction. Recalling the Birgin and Martínez direction in [15] is defined by(9)dk=-∇fxk,ifk=0,-θ∇fxk+βkdk,ifk≥1,where θk=sTsk/skTyk; see Raydan [16] for detail.In [2] Ortega and Rheinboldt used the term(10)gk=Fxk+αkFk-Fkαk,to approximate the gradient ∇f(xk), which avoids computing exact gradient and αk updated via line search method. It is clear that when Fk is small, then gk≈∇f(xk).Recall, from Newton’s direction,(11)dk+1=-J-1∇fxk+1.Combining (9) and (11), we have(12)-Jxk-1∇fxk+1=-∇fxk+1θk+βkdk.Multiplying both sides of (12) by J(xk) leads to(13)-JxkJxk-1∇fxk+1=-Jxk∇fxk+1θk+Jxkβkdk.After little linear algebra, (13) transforms to(14)-∇fxk+1=-θkJxk∇fxk+1+βkJxkdk.To ensure good approximation, we multiply both sides of (14) by skT to obtain(15)-skT∇fxk+1=-skTθkJxk∇fxk+1+skTβkJxkdk.Equation (15) can be rewritten as(16)-skT∇fxk+1=-θkskTJxk∇fxk+1+βkskTJxkdk.From secant condition, we have(17)Jxksk=yk,sk=Jxk-1yk,(18)skTJxkT=ykT,skT=skTJxk-1T.It is vital to note that, for this work, we claim that J(xk) is symmetric matrix ∀k. Hence, (18) can also be written as(19)skTJxk=ykT.Substituting (19) into (16) yields(20)-skT∇fxk+1=-θkykT∇fxk+1+βkykTdk.After simplification, we obtained our CG parameter (βk) as(21)βk=θkykT∇fxk+1-skT∇fxk+1ykTdk.Motivated by [3, 7] and using (10) we derive our CG parameter(22)βk=θkyk-skTykTdkgk+1.Having derived the CG parameter (βk) in (22) and by using (9), we then present our direction as(23)d0=-gx0,dk+1=-θkgk+1+θkyk-skTgk+1ykTdkdk,k=1,2,…,where θk=sTsk/skTyk.Finally, we present our scheme as(24)xk+1=xk+αkdk.The direction dk given by (23) may not be a descent direction of (8), and then the standard Wolfe and Armijo line searches cannot be used to compute the step size directly. Zhang et al. [8] proved the global convergence of the global PRP method for general nonconvex optimization using some nondescent line search. Motivated by this, we use the nonmonotone line search proposed by Li and Fukushima in [3] to compute our step size αk. Let ω1>0, ω2>0, and r∈(0,1) be constants and let ηk be a given positive sequence such that(25)∑k=0∞ηk<∞.Let αk=max⁡1,rk that satisfy(26)fxk+αkdk-fxk≤-ω1αkFxk2-ω2αkdk2+ηkfxk.Now, we can describe the algorithm for our proposed method as follows.Algorithm 1 (derivative-free CG method (DFCG)).
Consider the following steps.
Step 1. Given x0, α>0, σ∈(0,1), and compute d0=-g0, and set k=0.
Step 2. Compute gk using (10) and test the stopping criterion. If yes, then stop; otherwise continue with Step 3.
Step 3. Compute αk by the line search (26).
Step 4. Compute xk+1=xk+αkdk.
Step 5. Compute the search direction as dk+1=-θkgk+1+(θkyk-sk)Tgk+1/ykTdkdk.
Step 6. Consider k=k+1 and go to Step 2.


## 3. Convergence Result

This section presents global convergence results of the derivative-free conjugate gradient methods. To begin with, let us define the level set(27)Ω=x∣fx≤enfx0.In order to analyze the convergence of our method, we will make the following assumptions on nonlinear systems F.Assumption 2.
Consider the following:(i)
The level setΩ defined by (27) is bounded.(ii)
There existsx∗∈Ω such that F(x∗)=0 and F′(x) is continuous for all x.(iii)
In some neighborhoodN of Ω, the Jacobian is Lipschitz continuous; that is, there exists a positive constant L>0 such that(28)F′x-F′y≤Lx-y,for allx,y∈N.
Properties (i) and (ii) imply that there exist positive constantsM1, M2, and L1 such that(29)Fx≤M1,Jx≤M2,∀x∈N,(30)∇fx-∇fy≤L1x-y,Jx≤M2,∀x,y∈N.Lemma 3 (see [3]).
Let the sequencexk be generated by the algorithms above. Then the sequence Fk converges and xk∈E for all k≥0.Lemma 4.
Let the properties of (1) above hold. Then one has(31)limk→∞⁡αkdk=limk→∞⁡sk=0,(32)limk→∞⁡αkFk=0.Proof.
By (25) and (26) we have, for all k>0,(33)ω2αkdk2≤ω1αkFxk2+ω2αkdk2≤Fk2-Fk+12+ηkFk2by summing the above k inequality; then we obtain(34)ω2∑i=0kαkdk2≤Fk2∑i=0k1-ηi-Fk+12,so from (29) and the fact that ηk satisfies (25) the series ∑i=0kαkdk2 is convergent. This implies (31). By a similar way, we can prove that (32) holds.The following result shows that our derivative-free conjugate gradient methods algorithm is globally convergent.Theorem 5.
Let the properties of (1) above hold. Then the sequence xk which is generated by derivative-free conjugate gradient methods algorithm converges globally; that is,(35)liminfk→∞⁡∇fxk=0.Proof.
We prove this theorem by contradiction. Suppose that (35) is not true, and then there exists a positive constant τ such that(36)∇fxk≥τ,∀k≥0.Since ∇f(xk)=JkFk, (36) implies that there exists a positive constant τ1 satisfying(37)Fk≥τ1,∀k≥0.
Case (i) is as follows: considerlimsupk→∞αk>0. Then by (32), we have liminfk→∞Fk=0. This and Lemma 3 show that limk→∞Fk=0, which contradicts with (36).
Case (ii) is as follows: considerlimsupk→∞αk=0. Since αk≥0, this case implies that(38)limk→∞⁡αk=0,and by definition of gk in (10) and the symmetry of the Jacobian, we have(39)gk-∇fxk=Fxk+αk-1Fk-Fkαk-1-JkTFk=∫01Jxk+tαk-1Fk-JkdtFk≤LM12αk-1,where we use (29) and (30) in the last inequality. Inequalities (25), (26), and (36) show that there exists a constant τ2>0 such that(40)gk≥τ2,∀k≥0.By (10) and (29), we get(41)gk=∫01Jxk+tαk-1FkFkdt≤M1M2,∀k≥0.From (41) and (30), we obtain(42)yk=gk-gk+1≤gk-∇fxk+gk-1-∇fxk-1+∇fxk-∇fxk-1≤LM12αk-1+αk-2+L1sk-1.This together with (38) and (32) shows that limk→∞yk=0. Hence from (41), (42), and (40), we have(43)θk≤skTskskTyk⟶0meaning that there exists a constant λ∈(0,1) such that for sufficiently large k(44)θk≤λ.Again from the definition of our βk∗ we obtain(45)βk∗≤θkyk-skgk+1ykTsk≤M1M2θkyk-skykTsk⟶0which implies that there exists a constant ρ∈(0,1) such that for sufficiently large k(46)βk∗≤ρ.Without loss of generality, we assume that the above inequalities hold for all k≥0. Then we get(47)dk+1≤θkgk+1+βkdk≤λM1M2+ρdkwhich shows that the sequence dk+1 is bounded. Since limk→∞αk=0, then αk′=αk/r does not satisfy (26); namely,(48)fxk+αk′dk>fxk-ω1αk′Fxk2-ω2αk′dk2+ηkfxk,which implies that(49)fxk+αk′dk-fxkαk′>-ω1αk′Fxk2-ω2αk′dk2.By the mean-value theorem, there exists δk∈(0,1) such that(50)fxk+αk′dk-fxkαk′=∇fxk+δkαk′dkTdk.Since xk⊂Ω is bounded, without loss of generality, we assume xk→x∗. By (10) and (23), we have(51)limk→∞⁡dk+1=-limk→∞⁡θk+1gk+1+limk→∞⁡βk∗dk≤-limk→∞⁡gk+1+limk→∞⁡βk∗dk=-∇fx∗,where we use (45) and (26) and the fact that the sequence dk+1 is bounded.
On the other hand, we have(52)limk→∞⁡∇fxk+δkαk′dk=∇fx∗.Hence, from (49)–(52), we obtain -∇f(x∗)T∇f(x∗)≥0, which means ∇f(x∗)=0. This contradicts with (36). The proof is then completed.


## 4. Numerical Results

In this section, we compare the performance of our following methods for solving nonlinear equation (1) with an inexact PRP conjugate gradient method for solving symmetric nonlinear equations [1]:(i)
A derivative-free CG method (DFCG): we setω1=ω2=10-4, α1=0.01, r=0.2, and ηk=1/(k+1)2.(ii)
An inexact PRP (IPRP): we setω1=ω2=10-4, α1=0.01, r=0.2, and ηk=1/(k+1)2.The code for the DFCG method was written in Matlab 7.4 R2010a and run on a personal computer 1.8 GHz CPU processor and 4 GB RAM memory. We stopped the iteration if the total number of iterations exceeds 1000 or Fk≤10-3. We tested the two methods on nine test problems with different initial points and n values. Problems 6–9 are from [9].Problem 1:(53)Fix=1-xi2+xi1+xixn-2xn-1xn-2,i=1,2,…,n.Problem 2:(54)Fix=xi-0.1xi+12,i=1,2,…,n-1,Fnx=xn-0.1x12.Problem 3:(55)F1x=x12-3x1+1+cos⁡x1-x2,Fix=xi2-3xi+1+cos⁡xi-xi-1,i=1,2,…,n.Problem 4:(56)Fix=xicos⁡xi-1n-xnsin⁡xi-1-xi-12-1n∑i=1nxi,i=1,2,…,n.Problem 5:(57)Fix=xi-3xisin⁡xi3-0.66+2,fori=1,2,…,n.Problem 6:(58)F1x=x1x12+x22-1,Fix=xixi-12+2xi2+xi+12,i=2,3,…,n-1,Fnx=xnxn-12+xn2.Problem 7:(59)F3i-2x=x3i-2x3i-1-x3i2-1,F3i-1x=x3i-2x3i-1x3i-x3i-22+x3i-12-2,F3ix=e-x3i-2-e-x3i-1.Problem 8:(60)Fx=2-1-12-1⋱⋱⋱⋱⋱-1-12x+e1x-1,…,enx-1T.Problem 9:(61)Fx=2-102-1⋱⋱⋱⋱⋱-102x+sin⁡x1-1,…,sin⁡xn-1T.In Table1, we listed numerical results, where “Iter” and “Time” stand for the total number of all iterations and the CPU time in seconds, respectively; Fk is the norm of the residual at the stopping point. The numerical results indicate that the proposed method DFCG compared to IPRP has minimum number of iterations and CPU time, respectively. Figures 1 and 2 are performance profile derived by Dolan and Moré [17] which show that our claim is justified, that is, less CPU time and number of iterations for each test problem.Table 1
Test results for the two methods, wheree = ones(n, 1).
Problem (P)
x
0
n
DFCG
IPRP
Iter
Time (s)
‖
F
k
‖
Iter
Time (s)
‖
F
k
‖
1
0.5e
100
7
0.011199
2.03E − 04
15
0.018191
7.56E − 04
500
7
0.019229
4.54E − 04
16
0.039993
9.77E − 04
1000
7
0.025726
6.42E − 04
17
0.058292
7.43E − 04
2000
7
0.042199
9.08E − 04
18
0.099609
6.03E − 04
5000
8
0.138149
1.69E − 05
18
0.263836
9.53E − 04
2
e
10
5
0.005881
7.53E − 06
7
0.007772
9.53E − 05
500
5
0.009256
5.32E − 05
7
0.016192
6.74E − 04
5000
5
0.073385
1.68E − 04
9
0.097816
2.92E − 06
10000
5
0.16378
2.38E − 04
9
0.172704
4.14E − 06
100000
5
2.308894
7.53E − 04
9
1.923969
1.31E − 05
500000
6
8.43021
3.28E − 07
9
10.123595
2.92E − 05
3
0.5e
100
3
0.006565
2.42E − 05
5
0.005039
1.23E − 06
1000
3
0.014514
7.64E − 05
5
0.023714
3.90E − 06
5000
3
0.068582
1.71E − 04
5
0.08417
8.71E − 06
10000
3
0.132546
2.42E − 04
5
0.16412
1.23E − 05
100000
3
1.49153
7.64E − 04
5
2.297944
3.90E − 05
500000
4
9.227381
7.21E − 07
5
9.234416
8.71E − 05
4
0.01e
100
3
0.006919
9.15E − 05
51
0.090634
8.20E − 04
500
3
0.010871
2.06E − 04
57
0.17519
9.40E − 04
1000
3
0.021833
2.92E − 04
61
0.276824
8.50E − 04
5000
3
0.088476
6.53E − 04
67
1.277072
9.73E − 04
20000
4
0.312725
8.75E − 07
73
4.686063
9.97E − 04
5
e
50
4
0.006221
1.67E − 04
11
0.015403
5.79E − 04
500
4
0.014341
5.28E − 04
14
0.040054
5.73E − 04
5000
5
0.090997
2.69E − 06
15
0.236538
2.79E − 04
50000
5
0.725966
8.52E − 06
15
1.627153
8.82E − 04
0.001e
500
6
0.01594
6.07E − 04
17
0.044403
5.63E − 04
5000
7
0.105773
1.22E − 06
19
0.253817
5.54E − 04
50000
7
0.811278
3.86E − 06
20
2.042188
9.98E − 04
6
0.1
50
35
0.038655
8.04E − 04
54
0.06337
7.59E − 04
500
33
0.060248
5.98E − 04
43
0.086961
9.77E − 04
1000
46
0.109835
8.61E − 04
68
0.198439
9.88E − 04
2000
37
0.13327
8.20E − 04
58
0.245378
9.26E − 04
3000
35
0.217145
8.02E − 04
59
0.395438
9.60E − 04
8000
39
0.53251
9.06E − 04
60
1.015941
7.87E − 04
15000
48
1.237265
8.80E − 04
62
1.938341
9.80E − 04
7
e
10
4
0.00539
4.32E − 05
17
0.029569
8.17E − 04
50
4
0.006503
9.98E − 05
21
0.042007
6.78E − 04
500
4
0.011576
3.21E − 04
25
0.08001
7.86E − 04
5000
5
0.099718
5.91E − 07
29
0.442995
8.96E − 04
50000
5
0.731552
1.87E − 06
35
4.100204
6.12E − 04
8
e
10
34
0.366137
8.54E − 04
31
0.403368
8.94E − 04
50
26
0.29999
4.94E − 04
36
0.495837
3.31E − 04
100
28
0.410634
8.53E − 04
29
0.490114
9.47E − 04
500
224
6.894744
9.19E − 04
317
11.140713
7.85E − 04
1000
140
10.284715
7.87E − 04
226
19.377201
7.54E − 04
2000
134
18.935605
4.21E − 04
149
30.958082
9.70E − 04
0.1
10
18
0.202136
9.48E − 04
28
0.373317
9.96E − 04
100
32
0.460852
4.80E − 04
29
0.484755
5.21E − 04
500
32
0.96522
8.97E − 04
37
1.365365
8.93E − 04
1000
23
1.717939
7.32E − 04
32
3.189859
7.99E − 04
2000
27
5.596633
9.00E − 04
34
9.085955
7.93E − 04
5000
27
25.832759
5.92E − 04
28
33.107487
7.68E − 04
9
e
2000
11
0.027851
7.39E − 06
13
0.030663
6.19E − 04
3000
66
25.412837
8.44E − 04
69
34.946683
9.34E − 04
5000
63
72.39168
8.03E − 04
69
83.7963
9.95E − 04Figure 1
Comparison of the performance of DFCG and IPRP methods as the dimension increases (in terms of CPU time).Figure 2
Comparison of the performance of DFCG and IPRP methods as the dimension increases (in terms of number of iterations).Moreover, on average, ourF(xk) is too small which signifies that the solution obtained is true approximation of the exact solution compared to the IPRP.


## 5. Conclusions

Many practical problems possessed symmetrical property. In this paper we present a derivative-free conjugate gradient (DFCG) method for symmetric nonlinear equations and compare its performance with that of an inexact PRP conjugate gradient method [1] by doing some numerical experiments. We also proved the global convergence of our proposed method by using a backtracking type line search, and the numerical results show that our method is efficient.



---

*Source: 961487-2015-09-03.xml*