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 αk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\alpha_{k}$\end{document}. The given nonmonotone line search technique can avoid computing the Jacobian matrix. The global convergence and the 1+q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$1+q$\end{document}-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 e : n → n 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 (.) 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 ∇e(x) of e is symmetric, then problem (.) is called a system of symmetric nonlinear equations. Let p be the norm function with p(x) =   e(x)  , where · is the Euclidean norm. Then (.) is equivalent to the following global optimization models: (  .  ) In fact, there are many actual problems that can convert to the above problems (.) (see [-] etc.) and have similar models (see [-] etc.). The iterative formula for (.) is where α k is a step length and d k is a search direction. Now let us review some methods for α k and d k , respectively: (i) Li and Fukashima [] proposed an approximately monotone technique for α k : where e k = e(x k ), δ  > , δ  >  are positive constants, α k = r i k , r ∈ (, ), i k is the smallest nonnegative integer i satisfying (.) and k such that (  .  ) (iii) Brown and Saad [] gave the following technique to obtain α k : where β ∈ (, ) and ∇p(x k ) = ∇e(x k )e(x k ).
and some convergence results are obtained. Next we present some techniques for the calculation of d k . At present, there exist many well-known methods for d k , 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 d k : This method is regarded as one of the most effective methods. However, its efficiency largely depends on the possibility to efficiently solve (.) at each iteration. Moreover, the exact solution of the system (.) could be too burdensome when the iterative point x k is far from the exact solution []. In order to overcome this drawback, inexact quasi-Newton methods are often used.
(ii) The quasi-Newton method is of the form where B k 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 s k = x k+x k , y k = e k+e k , and e k+ = e(x k+ ). Set H k to be the inverse of B k , then the inverse formula of (.) has the following form: There exist many quasi-Newton methods (see [, , -]) 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 [] for Newton's methods. Many subsequent papers have exploited nonmonotone line search techniques of this nature (see [-] etc.), which shows that the nonmonotone technique works well in many cases. Considering these points, Zhu [] proposed the nonmonotone line search (.). From (.), we can see that the Jacobian matrix ∇e(x) must be computed at every iteration. Computing the Jacobian matrix ∇e(x) 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 α k : where σ ∈ (, ) is a constant and d k is a solution of (.). Comparing with (.), the new technique (.) does not compute the Jacobian matrix ∇e(x). Then the storage and workload can be saved in theory. In Section , we will state the technique (.) 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 x * , if the choice of the point x is correct then the Newton method can get the optimal solution x * 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 [-] etc.), where the results of [] 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 ∇e k must not be computed at every iteration. • The given method possesses the sufficient descent property for the normal function p(x).
• The global convergence and the  + q-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 , the algorithms are stated. Convergent results are established in Section  and numerical results are reported in Section . In the last section, our conclusion is given. Throughout this paper, we use these notations: · is the Euclidean norm, e(x k ) and g(x k+ ) are replaced by e k and g k+ , 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)
Step : Given any Step : If e k ≤ , stop. Otherwise let d k = -e k and go to next step.
Step : Let x k+ = x k + α k d k .
Step : Compute d k+ = -e k+ + β k d k , set k = k +  and go to Step .
Step  is a scalar and different β k will determine different CG methods.
(ii) From Step  and [], it is easy to deduce that there exists α k such that (.). 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 x sup , then the given algorithm is defined as follows.

Algorithm  (Main algorithm)
Step : Choose x sup ∈ n as the initial point, an initial symmetric positive definite matrix B  ∈ n×n , and constants r, σ ∈ (, ), main < , a positive integer M > , m(k) = , let k := ; Step : Stop if e k ≤ main . Otherwise solve (.) to get d k .
Step : Let the next iterative be x k+ = x k + α k d k .
Step : Update B k by quasi-Newton update formula and ensure the update matrix B k+ is positive definite.
Step : Let k := k + . Go to Step .

Remark
Step  of Algorithm  can ensure that B k is always positive definite. This means that (.) has a unique solution d k . By positive definiteness of B k , it is easy to obtain e T k d k < . In the following sections, we only concentrate to the convergence of the main algorithm.

Convergence analysis
Let be the level set with Similar to [, , ], the following assumptions are needed to prove the global convergence of Algorithm .
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 M * ≥ m * >  such that Assumption B B k is a good approximation to ∇e k , i.e., where * ∈ (, ) is a small quantity.
Considering Assumption B and using the von Neumann lemma, we deduce that B k is also bounded (see []).
Proof By using (.), we get Thus, we have It follows from (.) that The proof is complete.
The following lemma shows that the line search technique (.) is reasonable, then Algorithm  is well defined.

Lemma . Let Assumptions
A and B hold. Then Algorithm  will produce an iteration x k+ = x k + α k d k in a finite number of backtracking steps.
Proof From Lemma . in [] we have in a finite number of backtracking steps Proof By the acceptance rule (.), we have This means that the sequence {p(x l(k) )} is decreasing for all k.
By (.) and (.), for all k > M, we get Since {p(x l(k) )} is convergent, from the above inequality, we have According to (.) and the rule for accepting the step α k d k , Using (.) and (.) in [], we have where δ * >  is a constant and σ < δ * . So we get If x kx * sufficiently small, then the inequality holds.
Theorem . Let the assumptions in Lemma . hold. Assume that there exists a sufficiently small ε  >  such that B k -∇e(x k ) ≤ ε  for each k. Then the sequence {x k } converges to x * superlinearly for α k = . Moreover, if e is q-order smooth at x * and there is a neighborhood U of x * satisfying for any x k ∈ U, then x k → x * with order at least  + q, where η is a constant.
Proof Since g is continuously differentiable and ∇e(x) is nonsingular at x * , there exists a constant γ >  and a neighborhood U of x * satisfying max ∇e(y) , ∇e(y) - ≤ γ , where ∇e(y) is nonsingular for any y ∈ U. Consider the following equality when α k = : the second term and the third term are o( x kx * ). By the von Neumann lemma, and considering that ∇e(x k ) is nonsingular, B k is also nonsingular. For any y ∈ U and ∇e(y) being nonsingular and max{ ∇e(y) , ∇e(y) - } ≤ γ , then we obtain from Lemma . this means that the sequence {x k } converges to x * superlinearly for α k = . If e is q-order smooth at x * , then we get Consider the second term of (.) as x k → x * , and use (.), we can deduce that the second term of (.) is also O( x kx * q+ ). Therefore, we have The proof is complete.

Numerical results
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 x  . These functions are stated as follows.

Function  Variable dimensioned function
Initial guess:

Function  Discrete boundary value problem []
.

Function  The discretized two-point boundary value problem similar to the problem in []
e(x) = Ax +  (n + )  F(x) = , when A is the n × n tridiagonal matrix given by In the experiments, all codes were written in MATLAB ra and run on a PC with GT@. GHz CPU processor and . 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 x  (called the normal method). Aslam Noor et al.
[] presented a variational iteration technique for nonlinear equations, where the so-called VIM method has the better numerical performance. The VIM method has the following iteration form: where β i ∈ (, ) for i = , , . . . , n. 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 r = ., σ = ., M = , =  - , and main =  - . In order to ensure the positive definiteness of B k , in Step  of the main algorithm: if y T k s k > , update B k by (.), otherwise let B k+ = B k . This program will also be stopped if the iteration number of main algorithm is larger than . Since the line search cannot always ensure these descent conditions d T k e k <  and d T k ∇e(x k )e k < , 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 α k 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 [, ] For the line search technique, (.) is used and the largest search number of times is ten, where δ  = δ  =  - , and k =  NI  (NI is the iteration number). The sub-algorithm will also stopped if the iteration number is larger than . 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   From Tables -, 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 VIM method, the results of Problems - are very interesting, but it fails for Problems -. 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.
The tool of Dolan and Moré [] is used to analyze the efficiency of these three algorithms.
Figures - show that the performance of these methods are relative to NI, NG, and cpu time of Tables -, 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.

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