A quasi-Newton algorithm for large-scale nonlinear equations.

In this paper, the algorithm for large-scale nonlinear equations is designed by the following steps: (i) a conjugate gradient (CG) algorithm is designed as a sub-algorithm to obtain the initial points of the main algorithm, where the sub-algorithm's initial point does not have any restrictions; (ii) a quasi-Newton algorithm with the initial points given by sub-algorithm is defined as main algorithm, where a new nonmonotone line search technique is presented to get the step length [Formula: see text]. The given nonmonotone line search technique can avoid computing the Jacobian matrix. The global convergence and the [Formula: see text]-order convergent rate of the main algorithm are established under suitable conditions. Numerical results show that the proposed method is competitive with a similar method for large-scale problems.

Introduction

Consider the following nonlinear equations: Here is continuously differentiable and n denotes large-scale dimension. The large-scale nonlinear equations are difficult to solve since the relations of the variables x are complex and the dimension is larger. Problem (1.1) can model many real-life problems, such as engineering problems, dimensions of mechanical linkages, concentrations of chemical species, cross-sectional properties of structural elements, etc. If the Jacobian of e is symmetric, then problem (1.1) is called a system of symmetric nonlinear equations. Let p be the norm function with , where is the Euclidean norm. Then (1.1) is equivalent to the following global optimization models: In fact, there are many actual problems that can convert to the above problems (1.2) (see [1-9] etc.) and have similar models (see [10-27] etc.). The iterative formula for (1.1) is where is a step length and is a search direction. Now let us review some methods for and , respectively:

(i) Li and Fukashima [28] proposed an approximately monotone technique for : where , are positive constants, is the smallest nonnegative integer i satisfying (1.3) and such that (ii) Gu et al. [29] presented a descent line search technique: (iii) Brown and Saad [30] gave the following technique to obtain : where and .

(iv) Based on this technique, Zhu [31] proposed a nonmonotone technique: , , and , , and M is a nonnegative integer.

(v) Yuan and Lu [32] gave a new technique: and some convergence results are obtained.

Next we present some techniques for the calculation of . At present, there exist many well-known methods for , such as the Newton method, the trust region method, and the quasi-Newton method, etc.

(i) The Newton method has the following form to get : This method is regarded as one of the most effective methods. However, its efficiency largely depends on the possibility to efficiently solve (1.9) at each iteration. Moreover, the exact solution of the system (1.9) could be too burdensome when the iterative point is far from the exact solution [33]. In order to overcome this drawback, inexact quasi-Newton methods are often used.

(ii) The quasi-Newton method is of the form where is generated by a quasi-Newton update formula, where the BFGS (Broyden-Fletcher-Goldfarb-Shanno) update formula is one of the well-known quasi-Newton formulas with where , and . Set to be the inverse of , then the inverse formula of (1.11) has the following form: There exist many quasi-Newton methods (see [31, 32, 34–39]) representing the basic approach underlying most of the Newton-type large-scale algorithms.

The earliest nonmonotone line search framework was developed by Grippo, Lampariello, and Lucidi in [40] for Newton’s methods. Many subsequent papers have exploited nonmonotone line search techniques of this nature (see [41-44] etc.), which shows that the nonmonotone technique works well in many cases. Considering these points, Zhu [31] proposed the nonmonotone line search (1.7). From (1.7), we can see that the Jacobian matrix must be computed at every iteration. Computing the Jacobian matrix may be expensive if n is large and for any n at every iteration. Thus, one might prefer to remove the matrix, leading to a new nonmonotone technique.

Inspired by the above observations, we make a study of inexact quasi-Newton methods with a new nonmonotone technique for solving smooth nonlinear equations. In the kth iteration of our algorithm, the following new nonmonotone technique is used to obtain : where is a constant and is a solution of (1.10). Comparing with (1.7), the new technique (1.13) does not compute the Jacobian matrix . Then the storage and workload can be saved in theory. In Section 3, we will state the technique (1.13) is well defined.

It is well known that the initial point plays an important role in an algorithm. For example, the local superlinear convergence needs the iteration point x lies in the neighborhood of the optimal solution , if the choice of the point x is correct then the Newton method can get the optimal solution just need one step, moreover, the correct initial point can speed up the efficiency of an algorithm. The nonlinear conjugate gradient method is one of the most effective line search methods for unconstrained optimization problems due to its simplicity and low memory requirement, especially for large-scale problems. Many scholars have made many studies and obtained lots of achievements on the CG methods or other similar new methods (see [45-49] etc.), where the results of [46] are especially interesting. It has been proved that the numerical performance of the CG methods is very interesting for large-scale problems in different application fields. These considerations prompt us to design a CG algorithm (sub-algorithm) for solving large-scale nonlinear equations, where the terminated iteration point of the CG algorithm was used as the initial point of the given algorithm (main algorithm). Then there exist two advantages from this process: one is that we can use the CG’s characteristic to get a better initial point and another is that the good convergent results of the main algorithm can be preserved. The main attributes of this paper are stated as follows:

A sub-algorithm is designed to get the initial point of the main algorithm.

A new nonmonotone line search technique is presented, moreover, the Jacobian matrix must not be computed at every iteration.

The given method possesses the sufficient descent property for the normal function .

The global convergence and the -order convergent rate of the new method are established under suitable conditions.

Numerical results show that this method is more effective than other similar methods.

We organize the paper as follows. In Section 2, the algorithms are stated. Convergent results are established in Section 3 and numerical results are reported in Section 4. In the last section, our conclusion is given. Throughout this paper, we use these notations: is the Euclidean norm, and are replaced by and , respectively.

Algorithm

In this section, we will design a sub-algorithm and the main algorithm, respectively. These two algorithms are listed as follows.

Initial point algorithm (sub-algorithm)

Given any , , let .

If , stop. Otherwise let and go to next step.

Choose satisfies (1.4) and let until (1.3) holds.

Let .

If , stop.

Compute , set and go to Step 2.

Remark

(i) of Step 5 is a scalar and different will determine different CG methods.

(ii) From Step 2 and [28], it is easy to deduce that there exists such that (1.3). Thus, this sub-algorithm is well defined.

In the following, we will state the main algorithm. First, assume that the terminated point of sub-algorithm is , then the given algorithm is defined as follows.

Algorithm 1

Main algorithm

Step 0: Choose as the initial point, an initial symmetric positive definite matrix , and constants , a positive integer , , let ;

Stop if . Otherwise solve (1.10) to get .

Let until (1.13) holds.

Let the next iterative be .

Update by quasi-Newton update formula and ensure the update matrix is positive definite.

Let . Go to Step 1.

Step 4 of Algorithm 1 can ensure that is always positive definite. This means that (1.10) has a unique solution . By positive definiteness of , it is easy to obtain . In the following sections, we only concentrate to the convergence of the main algorithm.

Convergence analysis

Let Ω be the level set with Similar to [31, 32, 50], the following assumptions are needed to prove the global convergence of Algorithm 1.

Assumption A

(i) e is continuously differentiable on an open convex set containing Ω.

(ii) The Jacobian of e is symmetric, bounded, and positive definite on , i.e., there exist positive constants such that and

Assumption B

is a good approximation to , i.e., where is a small quantity.

Considering Assumption B and using the von Neumann lemma, we deduce that is also bounded (see [31]).

Lemma 3.1

Let Assumption B hold. Then is a descent direction of at , i.e.,

Proof

By using (1.10), we get Thus, we have It follows from (3.4) that The proof is complete. □

The following lemma shows that the line search technique (1.13) is reasonable, then Algorithm 1 is well defined.

Lemma 3.2

Let Assumptions A and B hold. Then Algorithm 1 will produce an iteration in a finite number of backtracking steps.

From Lemma 3.5 in [32] we have in a finite number of backtracking steps from which, in view of the definition of , we obtain (1.13). Thus we conclude the result of this lemma. The proof is complete. □

Now we establish the global convergence theorem of Algorithm 1.

Theorem 3.1

Let Assumptions A and B hold, and be generated by Algorithm 1. Then

By the acceptance rule (1.13), we have Using and , we obtain This means that the sequence is decreasing for all k. Then is convergent. Based on Assumptions A and B, similar to Lemma 3.4 in [32], it is not difficult to deduce that there exist constants such that By (1.13) and (3.10), for all , we get Since is convergent, from the above inequality, we have This implies that either or If (3.12) holds, following [40], by induction we can prove that and for any positive integer j. As and M is a positive constant, by and (3.14), it can be derived that According to (3.10) and the rule for accepting the step , This means which implies that or If equation (3.18) holds, since is bounded, then holds. The conclusion of this lemma holds. If (3.17) holds. Then acceptance rule (1.13) means that, for all large enough k, such that Since Using (3.19) and (3.20) in [32], we have where is a constant and . So we get Note that and , we have from dividing (3.21) by By (3.10), we have Consider and the bounded again, we complete the proof. □

Lemma 3.3

see Lemma 4.1 in [31]

Let e be continuously differentiable, and be nonsingular at which satisfies . Let If sufficiently small, then the inequality holds.

Theorem 3.2

Let the assumptions in Lemma 3.3 hold. Assume that there exists a sufficiently small such that for each k. Then the sequence converges to superlinearly for . Moreover, if e is q-order smooth at and there is a neighborhood U of satisfying for any , then with order at least , where η is a constant.

Since g is continuously differentiable and is nonsingular at , there exists a constant and a neighborhood U of satisfying where is nonsingular for any . Consider the following equality when : the second term and the third term are . By the von Neumann lemma, and considering that is nonsingular, is also nonsingular. For any and being nonsingular and , then we obtain from Lemma 3.3 this means that the sequence converges to superlinearly for .

If e is q-order smooth at , then we get Consider the second term of (3.27) as , and use (3.26), we can deduce that the second term of (3.27) is also . Therefore, we have The proof is complete. □

In this section, we report results of some numerical experiments with the proposed method. The test functions have the following form: where these functions have the associated initial guess . These functions are stated as follows.

Function 1

Exponential function 2 Initial guess: .

Function 2

Trigonometric function Initial guess: .

Function 3

Logarithmic function Initial guess: .

Function 4

Broyden tridiagonal function [[51], pp. 471-472] Initial guess: .

Function 5

Trigexp function [[51], p. 473] Initial guess: .

Function 6

Strictly convex function 1 [[52], p. 29]. is the gradient of . Initial guess: .

Function 7

Strictly convex function 2 [[52], p. 30]

is the gradient of . Initial guess: .

Function 8

Variable dimensioned function Initial guess: .

Function 9

Discrete boundary value problem [53]. Initial guess: .

Function 10

The discretized two-point boundary value problem similar to the problem in [53] when A is the tridiagonal matrix given by and with , and .

In the experiments, all codes were written in MATLAB r2009a and run on a PC with G1620T@2.40 GHz CPU processor and 4.0 GB memory and Windows XP operation system. In order to compare the performance the given algorithm with CG’s initial points (called new method with CG), we also do the experiment with only the main algorithm with initial points (called the normal method). Aslam Noor et al. [54] presented a variational iteration technique for nonlinear equations, where the so-called VIM1 method has the better numerical performance. The VIM1 method has the following iteration form: where for . In their paper, only low dimension problems (two variables) are tested. In this experiment, we also give the numerical results of this method for large-scale nonlinear equations to compare with our proposed algorithm.

The parameters were chosen as , , and . In order to ensure the positive definiteness of , in Step 4 of the main algorithm: if , update by (1.11), otherwise let . This program will also be stopped if the iteration number of main algorithm is larger than 200. Since the line search cannot always ensure these descent conditions and , an uphill search direction may occur in numerical experiments. In this case, the line search rule maybe fails. In order to avoid this case, the stepsize will be accepted if the searching time is larger than six in the inner circle for the test problems.

In the sub-algorithm, the CG formula is used by the following Polak-Ribière-Polyak (PRP) method [55, 56] For the line search technique, (1.3) is used and the largest search number of times is ten, where , and (NI is the iteration number). The sub-algorithm will also stopped if the iteration number is larger than 150. The iteration number, the function evaluations, and the CPU time of the sub-algorithm are added to the main algorithm for new method with CG. The meaning of the items of the columns of Table 1 is:

Table 1

Numerical results

		New method with CG				Normal method (only main algorithm)
P	Dim	NI/NG	GF	GD	cpu time	NI/NG	GF	GD	cpu time
1	1000	1/1	6.676674e−006	1.335335e−005	0.000000e+000	0/2	6.676674e−006	1.335335e−005	0.000000e+000
	2000	1/1	3.335834e−006	6.671668e−006	0.000000e+000	0/2	3.335834e−006	6.671668e−006	0.000000e+000
	3000	1/1	2.223334e−006	4.446667e−006	1.560010e−002	0/2	2.223334e−006	4.446667e−006	3.120020e−002
2	1000	12/17	1.570352e−007	1.214954e−007	1.544410e+000	199/2879	1.624268e−004	1.228551e+000	1.338801e+002
	2000	200/2927	8.144022e−005	3.138945e−001	8.647135e+002	199/2928	8.144022e−005	3.138945e−001	8.680832e+002
	3000	200/2326	5.434381e−005	2.481121e−003	1.614626e+003	199/2327	5.434381e−005	2.481121e−003	1.622785e+003
3	1000	8/8	4.194859e−006	8.389718e−006	1.560010e−002	115/1009	5.838535e−008	1.171251e−007	7.996611e+001
	2000	8/8	7.775106e−006	1.555021e−005	7.800050e−002	117/1040	1.161670e−007	2.328056e−007	5.662368e+002
	3000	9/9	1.614597e−012	3.231630e−012	1.553770e+001	137/1362	1.739498e−007	3.484891e−007	2.141410e+003
4	1000	92/165	5.632576e−006	5.669442e−006	2.215214e+000	199/285	3.703283e+001	1.196051e+001	1.356897e+002
	2000	87/156	6.245922e−006	6.085043e−006	1.502290e+001	199/230	3.637504e+000	9.570779e+000	9.591097e+002
	3000	94/169	6.678153e−006	6.437585e−006	4.731510e+001	199/234	2.639260e+002	1.904865e+001	3.096277e+003
5	1000	22/51	8.288299e−006	6.268946e−007	2.106014e+000	199/2570	3.195300e+004	4.649782e+006	1.779971e+001
	2000	21/50	4.114462e−006	1.943351e−006	5.179233e+000	199/2652	6.395300e+004	1.354972e+005	8.327333e+001
	3000	21/51	9.843373e−006	8.504719e−006	1.597450e+001	199/2853	9.595300e+004	1.135356e+005	1.314776e+002
6	1000	9/11	5.984185e−012	1.197596e−011	7.176046e−001	6/9	6.069722e−007	1.118437e−006	4.149627e+000
	2000	9/11	1.505191e−006	3.010383e−006	7.800050e−002	6/9	1.210931e−006	2.231765e−006	2.898499e+001
	3000	9/11	2.251571e−006	4.503142e−006	1.404009e−001	6/9	1.814891e−006	3.345093e−006	9.300780e+001
7	1000	200/602	1.208240e−003	4.697005e+000	3.424222e+001	199/573	3.156137e+005	3.893380e+004	1.378113e+002
	2000	200/760	1.612671e+001	9.034319e+000	2.420200e+002	199/644	1.014481e+006	3.131576e+005	9.872367e+002
	3000	200/693	5.570227e−003	9.501149e+001	7.698181e+002	199/743	9.357473e+007	2.087488e+007	3.087119e+003
8	1000	2/2	0.000000e+000	0.000000e+000	0.000000e+000	1/3	0.000000e+000	0.000000e+000	6.552042e−001
	2000	2/2	0.000000e+000	0.000000e+000	0.000000e+000	1/3	0.000000e+000	0.000000e+000	4.820431e+000
	3000	2/2	0.000000e+000	0.000000e+000	6.240040e−002	1/3	0.000000e+000	0.000000e+000	1.538170e+001
9	1000	67/118	7.138941e−006	1.820053e−005	3.010819e+000	2/5	2.358640e−006	4.611203e−006	1.404009e+000
	2000	70/124	6.342724e−006	1.607326e−005	2.062333e+001	2/5	5.917002e−007	1.169969e−006	9.703262e+000
	3000	74/131	7.447187e−006	1.799920e−005	6.450641e+001	2/5	2.632811e−007	5.225655e−007	3.084140e+001
10	1000	26/49	2.044717e−008	3.900140e−008	2.359983e+002	121/125	7.382123e−006	1.467673e−005	4.987196e+002
	2000	24/47	9.030382e−006	2.717060e−006	1.847286e+003	121/125	7.454090e−006	1.481981e−005	3.852538e+003
	3000	27/51	6.468831e−009	1.138377e−008	6.632227e+003	121/125	7.523322e−006	1.495745e−005	1.299774e+004

Dim: the dimension.

NI: the number of iterations.

NG: the number of function evaluations.

cpu time: the cpu time in seconds.

GF: the final norm function evaluations when the program is stopped.

GD: the final norm evaluations of search direction .

fails: fails to find the final values of when the program is stopped.

From Tables 1-2, it is easy to see that the number of iterations and the number of function evaluations of the new method with CG are less than those of the normal method for these test problems. Moreover, the cpu time and the final function norm evaluations of the new method with CG are more competitive than those of the normal method. For the VIM1 method, the results of Problems 1-7 are very interesting, but it fails for Problems 8-10. Moreover, it is not difficult to find that more CUP time is needed for this method. The main reason maybe lies in the computation of the Jacobian matrix at every iteration.

Table 2

Numerical results of VIM1 method

P	Dim	NI/NG	GF	cpu time	P	Dim	NI/NG	GF	cpu time
1	1000	1/1	6.676674e−006	1.560010e−002	6	1000	5/5	4.591162e−011	9.656462e+000
	2000	1/1	3.335834e−006	0.000000e+000		2000	5/5	9.140464e−011	7.439688e+001
	3000	1/1	2.223334e−006	3.120020e−002		3000	5/5	1.368978e−010	2.484628e+002
2	1000	18/18	2.840705e−007	5.494355e+001	7	1000	5/5	4.058902e−006	9.656462e+000
	2000	27/27	2.532474e−006	6.315544e+002		2000	6/6	1.983880e−017	8.993458e+001
	3000	22/22	9.781547e−007	1.669476e+003		3000	6/6	6.708054e−017	3.007543e+002
3	1000	5/5	5.430592e−007	9.578461e+000	8	1000	fails
	2000	5/5	5.619751e−007	7.435008e+001		2000	fails
	3000	5/5	5.870798e−007	2.484160e+002		3000	fails
4	1000	4/4	4.559227e−009	1.243328e+001	9	1000	fails
	2000	4/4	9.082664e−009	1.026487e+002		2000	fails
	3000	4/4	1.360090e−008	3.708768e+002		3000	fails
5	1000	9/9	2.648764e−006	3.196460e+001	10	1000	fails
	2000	9/9	2.649263e−006	2.529244e+002		2000	fails
	3000	9/9	2.649430e−006	8.258849e+002		3000	fails

The tool of Dolan and Moré [57] is used to analyze the efficiency of these three algorithms.

Figures 1-3 show that the performance of these methods are relative to NI, NG, and cpu time of Tables 1-2, respectively. The numerical results indicate that the proposed method performs best among these three methods. To this end, we think that the enhancement of this proposed method is noticeable.

Figure 1

Performance profiles of these three methods (NI).

Figure 3

Performance profiles of these three methods (cpu time).

Performance profiles of these three methods (NG).

Conclusion

In this paper, we focus on two algorithms solved a class of large-scale nonlinear equations. At the first step, a CG algorithm, called a sub-algorithm, was used as the initial points of the main algorithm. Then a quasi-Newton algorithm with the initial points done by a CG sub-algorithm was defined as the main algorithm. In order to avoid computing the Jacobian matrix, a nonmonotone line search technique was used in the algorithms. The convergence results are established and numerical results are reported.

According to the numerical performance, it is clear that the CG technique is very effective for large-scale nonlinear equations. This observation inspires us to design the CG methods to directly solve nonlinear equations in the future.