An Enhanced Matrix-Free Secant Method via Predictor-Corrector Modified Line Search Strategies for Solving Systems of Nonlinear Equations

Abstract

Diagonal updating scheme is among the cheapest Newton-like methods for solving system of nonlinear equations. Nevertheless, the method has some shortcomings. In this paper, we proposed an improved matrix-free secant updating scheme via line search strategies, by using the steps of backtracking in the Armijo-type line search as a step length predictor and Wolfe-Like condition as corrector. Our approach aims at improving the overall performance of diagonal secant updating scheme. Under mild assumptions, the global convergence results have been presented. Numerical experiments verify that the proposed approach is very promising.

1. Introduction

Consider the problem 𝐹(𝑥)=0, (1) with 𝐹∶𝑅^𝑛→𝑅^𝑛.

The mapping 𝐹 is assumed to satisfy the following standard assumptions:(1)𝐹is continuously differentiable in an open convex set Φ;(2)there exists a solution 𝑥^∗ of (1) in Φ such that 𝐹(𝑥^∗)=0 and 𝐹^′(𝑥^∗)≠0;(3)the Jacobian 𝐹^′(𝑥) is local Lipschitz continuous at 𝑥^∗.

The famous scheme for finding the solution to (1) is the Newton method. The method is simple to implement, and it produces an iterative sequence {𝑥_𝑘} from any given initial guess 𝑥₀ in the neighborhood of 𝑥^∗ via 𝑥_𝑘+1=𝑥_𝑘−𝐹^′𝑥_𝑘⁻¹𝐹𝑥_𝑘, (2) where 𝑘=0,1,2…𝑛.

The attractive features of this method are that it is easy to implement and converges rapidly [1]. However, the major difficulty of Newton-type method is matrix storage requirements especially when handling large systems of nonlinear equations [1]. To overcome such difficulty, the simple modification on the Newton method is the fixed Newton method. Fixed Newton method for the determination of solution 𝑥^∗ is given by 𝑥_𝑘+1=𝑥_𝑘−𝐹^′𝑥₀⁻¹𝐹𝑥_𝑘,𝑘=0,1,2,…. (3) This method avoids computation and storing the Jacobian in each iterations (except at 𝑘=0). However, it still requires to solve systems of 𝑛 linear equations and consumes more CPU time as the system’s dimension increases [2, 3].

Quasi-Newton method is another variant of Newton-type methods. It replaces the Jacobian or its inverse with an approximation which can be updated at each iteration and is given as 𝑥_𝑘+1=𝑥_𝑘−𝐵_𝑘⁻¹𝐹𝑥_𝑘, (4) where the matrix 𝐵_𝑘 is the approximation of the Jacobian at 𝑥_𝑘. The main idea behind quasi-Newton’s method is to eliminate the evaluation cost of the Jacobian matrix [4, 5]. However, the most critical part of such quasi-Newton method is storing the full matrix of the approximate Jacobian, which can be a very expensive task as the dimension of systems increases [6]. In contrast, this paper presents an improved variant of quasi-Newton update via the steps of backtracking in the Armijo-type line search and Wolfe–like condition and incorporating restarting strategy whenever the updating matrix is singular or nearly singular. The anticipation has been to improve the overall performance of the diagonal quasi-Newton method.

We organized the rest of this paper as follows. In the next section, we present the proposed method. Section 3 presents convergence results. Numerical experiments are reported in Section 4, and finally, conclusion is given in Section 5.

2. The Improved Secant Diagonal Updating

It is well known that it is not always feasible to compute the full elements of the Jacobian matrix of the given nonlinear function or it may be very expensive. We often have to approximate the Jacobian matrix by some other approach, and the famous method of doing so is quasi-Newton's method [4]. The basic idea underlining this approach has been to reduce the evaluation cost of the Jacobian matrix. This new approach generates a sequence of points {𝑥_𝑘} via 𝑥_𝑘+1=𝑥_𝑘−𝛼_𝑘𝐵_𝑘𝐹𝑥_𝑘, (5) where 𝛼_𝑘 is a step length, and 𝐵_𝑘 is a diagonal approximation of the inverse Jacobian matrix which can be updated in each iterations. To achieve this, we incorporate some new line search strategies, via the steps of backtracking in the Armijo-type line search as predictor and then improved via Wolfe-like condition as a corrector. The Armijo rule is among the inexact line search methods which guarantees a sufficient degree of accuracy to ensure the algorithm convergence. Nevertheless, the scheme needs high floating points operations and function call. In this paper, we present a simple line search strategy, which is less computational cost, floating points operations and CPU time consumptions respectively compared to classical Armijor line search. The approach is as follows.

Given 𝜖∈(0,1) and 𝜎<1, the proposed approach finds the appropriate 𝛼, such that 𝐹𝑥_𝑘+𝛼_𝑘𝑑_𝑘≤𝜎𝐹𝑥_𝑘. (6) In addition, the new strategy is implemented in an iterative way using a fixed initial value of 𝛼 as follows.

Algorithm 1 (Armijo-Like).   

Step  0. Set 𝐾=0, 𝛼₀>0, and 𝜎<1.

Step 1. ‖𝐹(𝑥_𝑘+1)‖≤𝜎‖𝐹(𝑥_𝑘)‖. Choose 𝛼_𝑘 as the step size; stop. Otherwise 𝛼_𝑘+1=𝛼_𝑘/2, 𝑘=𝑘+1.

Waziri et al. [7] have set the step length (𝛼_𝑘=1, for all  𝑘); this approach is mostly used in many Newton-like methods. Here, we continue in the spirit of diagonal updating, using a new line search strategy to obtain a good step length (𝛼_𝑘) in every iteration, anticipating to produce a more accurate approximation of the Jacobian inverse matrix and then employing restating strategy whenever the updating matrix is undefined. To this end, 𝐵_𝑘+1 would be obtained almost similar to the diagonal updating scheme presented in [7] in which 𝑠_𝑘=−𝛼_𝑘𝐵_𝑘𝐹(𝑥_𝑘) instead of 𝑠_𝑘=−𝐵_𝑘𝐹(𝑥_𝑘). Now, let the deviation between 𝐵_𝑘+1 and 𝐵_𝑘 denoted as Φ_𝑘=𝐵_𝑘+1−𝐵_𝑘 be minimized under some norms; the optimal solution is given as Φ_𝑘=𝑦^𝑇_𝑘𝑠_𝑘−𝑦^𝑇_𝑘𝐵_𝑘𝑦_𝑘trΨ²_𝑘Ψ_𝑘, (7) where Ψ_𝑘=diag((𝑦_𝑘⁽¹⁾)²,(𝑦_𝑘⁽²⁾)²,…,(𝑦_𝑘^(𝑛))²), ∑^𝑛_𝑖=1(𝑦_𝑘^(𝑖))⁴=tr(Ψ²_𝑘), 𝑦_𝑘=𝐹(𝑥_𝑘+1)−𝐹(𝑥_𝑘), 𝑠_𝑘=𝑥_𝑘+1−𝑥_𝑘, and tr is the trace operation, respectively. The updated formula for the proposed matrix 𝐵_𝑘 is given as [7] 𝐵_𝑘+1=𝐵_𝑘+𝑦^𝑇_𝑘𝑠_𝑘−𝑦^𝑇_𝑘𝐵_𝑘𝑦_𝑘trΨ²_𝑘Ψ_𝑘. (8) To safeguard on the possibilities of generating undefined 𝐵_𝑘+1, we let  𝐵_𝑘+1=𝐵₀, whenever this situation happens 𝐵_𝑘+1=𝐵_𝑘+𝑦^𝑇_𝑘𝑠_𝑘−𝑦^𝑇_𝑘𝐵_𝑘𝑦_𝑘trΨ²_𝑘Ψ_𝑘,trΨ²_𝑘𝐵≠0,₀,otherwise. (9)

Now, we can describe the algorithm for our proposed approach as follows.

Algorithm 2 (EMFM).  

Step 1. Choose an initial guess 𝑥₀, 𝜎∈(0,1), 𝐵₀=𝐼_𝑛, 𝛼₀>0 and let 𝑘∶=0.

Step 2. Compute  𝐹(𝑥_𝑘), and if ‖𝐹(𝑥_𝑘)‖≤10⁻⁴, stop.

Step 3. Compute 𝑑=−𝐹(𝑥_𝑘)𝐵_𝑘.

Step 4. If ‖𝐹(𝑥_𝑘+𝛼_𝑘𝑑_𝑘)‖≤𝜎‖𝐹(𝑥_𝑘)‖, retain 𝛼_𝑘 and go to 5. Otherwise, set 𝛼_𝑘+1=𝛼_𝑘/2 and repeat 4.

Step 5. If ‖𝐹(𝑥_𝑘+𝛼_𝑘𝑑_𝑘)−𝐹(𝑥_𝑘)‖≥‖𝐹(𝑥_𝑘+𝛼_𝑘𝑑_𝑘)‖−‖𝐹(𝑥_𝑘)‖, retain 𝛼_𝑘 and go to 6. Otherwise, set 𝛼_𝑘+1=𝛼_𝑘×1.1 and repeat 5.

Step 6. Let 𝑥_𝑘+1=𝑥_𝑘+𝛼_𝑘𝑑_𝑘.

Step 7. If ‖𝑥_𝑘+1−𝑥_𝑘‖₂+‖𝐹(𝑥_𝑘)‖₂≤10⁻⁴, stop. Otherwise go to Step 8.

Step 8. If ‖Δ𝐹_𝑘‖₂≥𝜖₁, where 𝜖₁=10⁻⁴, compute 𝐵_𝑘+1, and if not, 𝐵_𝑘+1=𝐵₀.

Step 9. Set 𝑘∶=𝑘+1 and go to 2.

3. Convergence Analysis

We present the convergence result of EMFM method by proving the existence of the step length 𝛼_𝑘>0. We will make the following assumptions on nonlinear system 𝐹.

Assumption 3. (i) 𝐹 is differentiable in an open convex set 𝐸 in ℜ^𝑛.(ii)There exists 𝑥^∗∈𝐸 such that (𝑥^∗)=0, and 𝐹^′(𝑥) is continuous for all 𝑥. (iii)𝐹(𝑥) satisfies Lipschitz condition of order one; that is, there exists a positive constant 𝜇 such that 𝐹(𝑥)−𝐹𝑦≤𝜇𝑥−𝑦, (10) for all 𝑥,𝑦∈ℜ^𝑛.(iv) There exist constants 𝑐₁≤𝑐₂ such that 𝑐₁‖𝜔‖²≤𝜔^𝑇𝐹^′(𝑥)𝜔≤𝑐₂‖𝜔‖² for all 𝑥∈𝐸 and 𝜔∈ℜ^𝑛.

To this end, we proceed by given the following result for the step length generated by the proposed strategies.

Theorem 4. Assume that 𝐹 is a strictly convex function. Suppose that the new strategies are employed with 𝑑_𝑘≠0 and positive 𝛼_𝑘 exist, for all sufficiently large 𝑘. Then, the iterates {𝑥_𝑘} generated by the line search algorithm have the property that lim𝐹𝑥_𝑘+𝛼_𝑘𝑑_𝑘=0,𝑎𝑠𝑘⟶∞. (11)

Proof. From condition (iii) of Assumption 3 and the fact that 𝑠_𝑘=𝑥_𝑘+1−𝑥_𝑘, we have 𝐹𝑥_𝑘−𝐹𝑥_𝑘+𝛼_𝑘𝑑_𝑘𝐹≤𝜇𝑠_𝑘𝐹=𝜇𝛼𝑑_𝑘𝐹≤𝜇|𝛼|𝑑_𝑘𝐹. (12) Recall that 𝐹(𝑥)−𝐹𝑥_𝑘+𝛼_𝑘𝑑_𝑘𝐹≤‖𝐹(𝑥)‖_𝐹−𝐹𝑥_𝑘+𝛼_𝑘𝑑_𝑘𝐹. (13) Equation (13) gives 𝐹(𝑥_𝑘)_𝐹−𝐹𝑥_𝑘+𝛼_𝑘𝑑_𝑘𝐹≤𝛼𝜇𝑑_𝑘𝐹=𝜇𝛼𝑑_𝑘𝐹≤𝜇|𝛼|𝑑_𝑘𝐹. (14) Since ‖𝑑_𝑘‖≠0 and ‖𝐹(𝑥_𝑘+𝛼_𝑘𝑑_𝑘)‖≤𝜎‖𝐹(𝑥_𝑘)‖, hence, it follows that 𝐹(𝑥_𝑘)_𝐹−𝜎𝐹𝑥_𝑘𝐹≤𝛼𝜇𝑑_𝑘𝐹,(1−𝜎)𝐹𝑥_𝑘𝐹≤𝛼𝜇𝑑_𝑘𝐹. (15) After little simplifications, we obtain (𝐹𝑥𝛼≥1−𝜎)_𝑘𝐹𝜇𝑑_𝑘𝐹>0. (16) Hence, 𝛼_𝑘 exist and are positive.
We continue to show the convergence of the iterates {𝑥_𝑘} by recalling that 𝐹𝑥_𝑘+𝛼_𝑘𝑑_𝑘≤𝜎𝐹𝑥_𝑘. (17) Then, we have 𝐹𝑥_𝑘+𝛼_𝑘𝑑_𝑘≤𝜎²𝐹𝑥_𝑘≤𝜎³𝐹𝑥_𝑘 (18) for a finite 𝑘, and (18) yields 𝐹𝑥_𝑘+𝛼_𝑘𝑑_𝑘≤𝜎^𝑘𝐹𝑥_𝑘. (19) Since the algorithm terminates at ‖𝐹(𝑥_𝑘)‖=0, it follows from (19) that 𝐹𝑥_𝑘+𝛼_𝑘𝑑_𝑘𝐹𝑥_𝑘≤𝜎^𝑘. (20) We have 𝐹𝑥lim^𝑘+𝛼_𝑘𝑑_𝑘𝐹𝑥_𝑘≤lim𝜎^𝑘,as𝑘⟶∞. (21) Therefore, due to 𝜎<1, it implies that lim𝜎^𝑘⟶0,as𝑘⟶∞. (22) hence lim𝐹𝑥_𝑘+𝛼_𝑘𝑑_𝑘=0,as𝑘⟶∞. (23)

To show the convergence results of EMFM method, we require to show that the updating matrix 𝐵_𝑘 is bounded above and below by some positive constants. Hence, we can state the following result on the boundedness of {‖Φ_𝑘‖_𝐹} by assuming that, without loss of generality, the updating matrix (8) is always used; then, we have the following.

Theorem 5. Let 𝐹 satisfy Assumption 3, and let ‖𝑦_𝑘‖≠0 for all finite 𝑘. Let {𝐵_𝑘} be the sequence generated by (8). If the given nonsingular 𝐵₀ satisfies 𝜃≤𝐵^𝑖₀≤𝜔,𝑖=1,2,...,𝑛, (24) for some constants 𝜃 and 𝜔, then the sequence {‖𝐵_𝑘‖_𝐹} is bounded for all finite 𝑘.

Proof. Since 𝐵_𝑘+1=𝐵_𝑘+Φ_𝑘, it follows that 𝐵_𝑘+1𝐹≤𝐵_𝑘𝐹+Φ_𝑘𝐹. (25) For 𝑘=0 and assuming 𝐵₀=𝐼, we have Φ₀^(𝑖)=𝑦^𝑇₀𝑠₀−𝑦^𝑇₀𝐵₀𝑦₀trΨ²₀𝑦₀^(𝑖)2≤𝑦^𝑇₀𝑠₀−𝑦^𝑇₀𝐵₀𝑦₀trΨ²₀𝑦₀^(max)2, (26) where (𝑦₀^(max))² is the largest element among (𝑦₀^(𝑖))²,𝑖=1,2,…,𝑛.
After multiplying (26) by (𝑦₀^(max))²/(𝑦₀^(max))² and substituting tr(Ψ²₀)=∑^𝑛_𝑖=1(𝑦₀^(𝑖))⁴, we have Φ₀^(𝑖)≤𝑦^𝑇₀𝑠₀−𝑦^𝑇₀𝐵₀𝑦₀𝑦₀^(max)2∑^𝑛_𝑖=1𝑦₀^(𝑖)4𝑦₀^(max)4. (27) Since (𝑦₀^(max))⁴/∑^𝑛_𝑖=1(𝑦₀^(𝑖))⁴≤1, then (27) turns into Φ₀^(𝑖)≤𝑦^𝑇₀𝐹^′(𝑥)𝑦₀−𝑦^𝑇₀𝐵₀𝑦₀𝑦₀^(max)2. (28) From Assumption 3 and 𝐵₀=𝐼, (28) becomes Φ₀^(𝑖)≤|𝑐−1|𝑦^𝑇₀𝑦₀𝑦₀^(max)2, (29) where 𝑐=max{|𝑐₁|,|𝑐₂|}.
Since (𝑦₀^(𝑖))²≤(𝑦₀^(max))² for 𝑖=1,…,𝑛, it follows that Φ₀^(𝑖)≤𝑛|𝑐−1|𝑦₀^(max)2𝑦₀^(max)2. (30) Hence, we obtain Φ_0𝐹≤𝑛^3/2|𝑐−1|. (31) Suppose that 𝜂=𝑛^3/2|𝑐−1|; then, Φ_0𝐹≤𝜂. (32) From the fact that ‖𝐵₀‖_𝐹=√𝑛, it follows that 𝐵_1𝐹≤𝛽, (33) where 𝛽=√𝑛+𝜂>0.

4. Numerical Results

In this section, we consider some benchmark problems to illustrate the performance of the method proposed in this paper for solving large-scale systems of nonlinear equations when compared to some Newton-like methods. The computations are performed in MATLAB 7.0 using double precision computer, and the stopping rule used is 𝑠_𝑘+𝐹𝑥_𝑘≤10⁻⁴. (34) The identity matrix has been chosen as an initial approximate Jacobian inverse. We further design the codes to terminate whenever one of the following happens:(i)the number of iteration is at least 250, but no point of 𝑥_𝑘 that satisfies (34) is obtained; (ii)CPU time in seconds reaches 250; (iii)There is insufficient memory to initiate the run.

The performances of these methods are be compared in terms of number of iterations and CPU time in seconds. In the following, some details on the benchmarks test problems are presented.

Problem 1. System of 𝑛 nonlinear equations is as follows: 𝑓_𝑖(𝑥)=𝑥_𝑖−3𝑥_𝑖sin𝑥_𝑖3−0.66+2,𝑖=1,2,…,𝑛,𝑥₀=(3,3,…,3). (35)

Problem 2. Extended Trigonometric function of Spedicator [8] is as follows: 𝑓_𝑖(𝑥)=𝑛−^𝑛_𝑗=1cos𝑥_𝑗+𝑖1−cos𝑥_𝑖−sin𝑥_𝑖+exp^{(1−cos𝑥𝑖)}−^𝑛_𝑗=1𝑥_𝑗+1,𝑖=1,…,𝑛,𝑥₀=1𝑛,1𝑛1,…,𝑛. (36)

Problem 3. System of 𝑛 nonlinear equations is as follows: 𝑓_𝑖(𝑥)=𝑥_𝑖−0.1𝑥²_𝑖+1,𝑓_𝑛(𝑥)=𝑥_𝑛−0.1𝑥²₁,𝑖=1,…,𝑛−1,𝑥₀=(7,7,…,7). (37)
(A1) Nonadiabatic Stirred Tank Reactors. The first application is the model of [9]. The model deals with two continuous nonadiabatic stirred tank reactors. The reactors are in a series, at steady state having a recycle component, and with an exothermic first-order irreversible reaction. By eliminating certain variables, the model results into two nonlinear systems of equations as follows: 𝑓₁=(𝐷1−𝜆)101+𝛽₁−𝑥₁exp10𝑥₁1+10𝑥₁/𝛾−𝑥₁,𝑓₂=𝑥₁−1+𝛽₂𝑥₂+()×𝐷1−𝜆10−𝛽₁𝑥₁−1+𝛽₂𝑥₂exp10𝑥₂1+10𝑥₂./𝛾 (38) The dimensionless temperatures of the two reactors are represented by 𝑥₁ and 𝑥₂. The parameters 𝜆, 𝛾, 𝛽₁, 𝛽₂, and 𝐷 are given as 1, 1000, 2, 2, and 22, respectively.
(A2) Navigation by Range Measurements. Consider two beacons determining position by measuring the distances 𝑟1 and 𝑟2, where (𝑝1,𝑝2) is the position of beacon 1, (𝑞1,𝑞2) is the position of beacon 2, and (𝑢,𝑣) is an unknown point (see Figure 1). This can be modeled into a two nonlinear systems of equations as follows: 𝑓₁=𝑝₁−𝑢²+𝑞₁−𝑣²−𝑟₁𝑓=0,₁=𝑝₂−𝑢²+𝑞₂−𝑣²−𝑟₂=0. (39) The parameters (𝑝₁,𝑞₁), (𝑝₂,𝑞₂), 𝑟₁, and 𝑟₂ are given as (10,10), (10,−10), 14, and 16, respectively.

Figure 1 

Table 1 shows the number of iterations and CPU time for these five methods, respectively. In Table 1, the value “Dim” denotes the dimension of the systems, “NI” denotes number of iterations, while “CPU” is the CPU time in seconds respectively. We analyze the performance of each method via execution time, floating points operations, and storage locations. One can observe that EMFM has the smallest number of iterations compared to the classical diagonal updating (I-VDN) proposed by Waziri et al. [7]. This shows that the line search strategies presented in this paper have increased the convergence speed of the classical diagonal updating method.

Table 1 

Numerical comparison of NM, FN, BM, I-VDN, and EMFM methods.

If we compare the performance of all methods, in terms of CPU time, it is clear that EMFM method consumes less CPU time than the others and still keeping memory requirement and CPU time in seconds to only 𝑂(𝑛). All five methods are able to obtain the solution (𝑥^∗=−1.1243,1.5001) of A2, but EMFM method consumes less CPU time in second (0.001) compared to the other 4 methods. Moreover, for A1, still proposed method has shown a promising performance with less storage locations (2 locations) whereas NM, FN, and BM, respectively, required 4 locations for each.

5. Conclusion

From the fact that there is a rapid development in research on solving nonlinear systems, nevertheless, the dimension of the nonlinear system is most of the times so large that it requires a very costly arithmetic operations when using some other Newton-like methods; so we usually employ cheap iterative approach, and the good candidate is classical diagonal updating. Notwithstanding, the classical updating scheme has some lapses. One is that it usually needs high number of iterations. The other is that the method generally has slow convergence behavior. This paper proposes an enhanced diagonal secant updating scheme based on the steps of backtracking in the Armijo-type line search and then improved via Wolfe-like condition. Our approach aims at improving the overall performance of diagonal secant updating scheme. The algorithm only requires to store a row vector while ignoring all the off and low diagonal elements and therefore largely reduces memory locations. In addition, as it uses two line search strategies (predictor and corrector) to obtain a new iterates point, the spectral properties of the diagonal updating scheme is improved, and rapid convergence property is gained.

Computational experiment suggests that it is very vital for diagonal updating scheme to use line search strategy. EMFM method has very good solving speed and the best performance among the Newton-like methods. Finally, it can be concluded that this approach would certainly be quite useful for solving large-scale systems of nonlinear equations.