An Improved Diagonal Jacobian Approximation via a New Quasi-Cauchy Condition for Solving Large-Scale Systems of Nonlinear Equations

Abstract

We present a new diagonal quasi-Newton update with an improved diagonal Jacobian approximation for solving large-scale systems of nonlinear equations. In this approach, the Jacobian approximation is derived based on the quasi-Cauchy condition. The anticipation has been to further improve the performance of diagonal updating, by modifying the quasi-Cauchy relation so as to carry some additional information from the functions. The effectiveness of our proposed scheme is appraised through numerical comparison with some well-known Newton-like methods.

1. Introduction

Let us consider the systems of nonlinear equations 𝐹(𝑥)=0, (1) where 𝐹∶𝑅^𝑛→𝑅^𝑛 is a nonlinear mapping. Often, the mapping 𝐹 is assumed to be satisfying the following assumptions:(A1) there exists an 𝑥^∗∈𝑅^𝑛 s.t 𝐹(𝑥^∗)=0;(A2)𝐹 is a continuously differentiable mapping in a neighborhood of 𝑥^∗;(A3)𝐹^′(𝑥^∗) is invertible.The well-known method for finding the solution to (1) is the classical Newton’s method which generates a sequence of iterates {𝑥_𝑘} from a given initial point 𝑥₀ via 𝑥_𝑘+1=𝑥_𝑘−𝐹^′𝑥_𝑘⁻¹𝐹𝑥_𝑘, (2) where 𝑘=0,1,2…. The attractive features of this method are rapid convergence and being easy to implement. Nevertheless, Newton’s method requires the computation of the matrix entails the first-order derivatives of the systems. In practice, computations of some functions derivatives are quite costly, and sometimes they are not available or could not be done precisely. In this case, Newton’s method cannot be applied directly.

Moreover, some substantial efforts have been made by numerous researchers in order to eliminate the well-known shortcomings of Newton’s method for solving systems of nonlinear equations, particularly large-scale systems (see, e.g., [1, 2]). Notwithstanding, most of these modifications of Newton’s method still have some shortfalls as Newton’s counterpart. For example, Broyden’s method and Chord Newton’s method need to store an 𝑛×𝑛 matrix, and their floating points operations, are 𝑂(𝑛²), respectively.

To tackle these disadvantages, a diagonally Newton’s method has been suggested by Leong et al. [3] and showed that their updating formula is significantly cheaper than Newton’s method and some of its variants. Based on this fact, it is pleasing to present an approach which will improve further the diagonal Jacobian approximation, as well as reducing the computational cost, floating points operations and number of iterations. This is what leads to the idea of this paper. The anticipation has been to further improve the performance of diagonal updating, by modifying the quasi-Cauchy relation so as to carry some additional information from the functions. We organized the paper as follows. In the next section, we present the details of the proposed method. Convergence results are present in Section 3. Some numerical results are reported in Section 4. Finally, conclusions are made in Section 5.

2. Derivation Process

This section presents a new diagonal quasi-Newton-like method for solving large-scale systems of nonlinear equations. The quasi-Newton method is an iterative method that generates a sequence of points {𝑥_𝑘} from a given initial guess 𝑥₀ via the following form:𝑥_𝑘+1=𝑥_𝑘−𝛼_𝑘𝐵_𝑘𝐹𝑥_𝑘𝑘=0,1,2…, (3) where 𝛼_𝑘 is a step length and 𝐵_𝑘 is an approximation to the Jacobian inverse which can be updated at each iteration for 𝑘=0,1,2…; the updated matrix 𝐵_𝑘+1 is chosen in such a way that it satisfies the secant equation, that is, 𝐵_𝑘+1𝑠_𝑘=𝑦_𝑘. (4)

It is clear that the only Jacobian information we have is 𝑦, and this is only approximation information. To this end, we incorporate more information from 𝑠_𝑘 and 𝐹_𝑘 to 𝑦 in order to present a better approximation to the Jacobian matrix. We consider the modification on 𝑦 presented by Li and Fukushima [4]: 𝑦̆=𝑦_𝑘+𝑣_𝑘𝐹_𝑘𝑠_𝑘, (5) where 𝑣_𝑘=1+max{−𝑠^𝑇_𝑘𝑦_𝑘/‖𝑠_𝑘‖²,0}.

Our aim here is to build a square matrix, say 𝐵, using diagonal updating scheme which is an approximation to the Jacobian inverse, and we let 𝐵_𝑘+1 satisfy the quasi-Cauchy equation, that is,𝑦̆_𝑘^𝑇𝑠_𝑘=𝑦̆_𝑘^𝑇𝐵_𝑘+1𝑦̆_𝑘. (6)

In addition, the deviation between 𝐵_𝑘+1 and 𝐵_𝑘 is minimized under some norms; hence, in the following theorem, we state the resulting update formula for 𝐵_𝑘.

Theorem 1. Assume that 𝐵_𝑘+1 be the diagonal update of a diagonal matrix 𝐵_𝑘. Let us denote the deviation between 𝐵_𝑘 and 𝐵_𝑘+1 as Ψ_𝑘=𝐵_𝑘+1−𝐵_𝑘. Suppose that 𝑦̆_𝑘≠0 which is defined by (5). Consider the following problem:1min2Ψ_𝑘²_𝐹s.t𝑦̆_𝑘^𝑇𝐵_𝑘+Ψ_𝑘𝑦̆_𝑘=𝑦̆_𝑘^𝑇𝑠_𝑘, (7) where ‖⋅‖_𝐹 denotes the Frobenius norm. Hence, the optimal solution of (7) is given by Ψ_𝑘=𝑦̆_𝑘^𝑇𝑠_𝑘−𝑦̆_𝑘^𝑇𝐵_𝑘𝑦̆_𝑘tr𝑉²_𝑘𝑉_𝑘, (8) where 𝑉_𝑘=diag((𝑦̆_𝑘⁽¹⁾)²,(𝑦̆_𝑘⁽²⁾)²,…,(𝑦̆_𝑘^(𝑛))²), ∑^𝑛_𝑖=1(𝑦̆_𝑘^(𝑖))⁴=tr(𝑉²_𝑘), and Tr is the trace operation.

Proof. Consider the Lagrangian function of (7): 𝐿Ψ_𝑘=1,𝛼2Ψ_𝑘²_𝐹+𝛼𝑦̆_𝑘^𝑇Ψ_𝑘𝑦̆_𝑘−𝑦̆_𝑘^𝑇𝑠_𝑘+𝑦̆_𝑘^𝑇𝐵_𝑘𝑦̆_𝑘, (9) where 𝛼 is the corresponding Lagrangian multiplier. By differentiating 𝐿 with respect to each Ψ_𝑘⁽¹⁾,Ψ_𝑘⁽²⁾,…,Ψ_𝑘^(𝑛), and setting them all equal to zero, we obtain Ψ_𝑘^(𝑖)=−𝛼𝑦̆_𝑘^(𝑖)2∀𝑖=1,2,…,𝑛. (10)
Multiplying both sides of (10) by (𝑦̆_𝑘^(𝑖))² and summing them all give 𝑦̆^𝑛_𝑖=1𝑘^(𝑖)2Ψ_𝑘^(𝑖)=−𝛼^𝑛_𝑖=1𝑦̆_𝑘^(𝑖)4forevery𝑖=1,2,…,𝑛. (11)
Differentiating 𝐿 with respect to 𝛼, and since 𝑦̆_𝑘^𝑇Ψ_𝑘𝑦̆_𝑘=∑^𝑛_𝑖=1(𝑦̆_𝑘^(𝑖))²Ψ^(𝑖), then we have 𝑦̆^𝑛_𝑖=1𝑘^(𝑖)2Ψ_𝑘^(𝑖)=𝑦̆_𝑘^𝑇𝑠_𝑘−𝑦̆_𝑘^𝑇𝐵_𝑘𝑦̆_𝑘. (12)
Equating (11) and (12) and substituting the relation into (10), finally we have Ψ_𝑘^(𝑖)=𝑦̆_𝑘^𝑇𝑠_𝑘−𝑦̆_𝑘^𝑇𝐵_𝑘𝑦̆_𝑘∑^𝑛_𝑖=1𝑦̆_𝑘^(𝑖)4𝑦̆_𝑘^(𝑖)2∀𝑖=1,2,…,𝑛. (13)
Since 𝐵_𝑘^(𝑖) is a diagonal component of 𝐵_𝑘, 𝑦̆_𝑘^(𝑖) is the 𝑖th component of vector 𝑦̆_𝑘, then 𝑉_𝑘=diag((𝑦̆_𝑘⁽¹⁾)²,(𝑦̆_𝑘⁽²⁾)²,…,(𝑦̆_𝑘^(𝑛))²) and ∑^𝑛_𝑖=1(𝑦_𝑘^(𝑖))⁴=tr(𝑉²_𝑘). We further rewrite (13) as Ψ_𝑘=𝑦̆_𝑘^𝑇𝑠_𝑘−𝑦̆_𝑘^𝑇𝐵_𝑘𝑦̆_𝑘tr𝑉²_𝑘𝑉_𝑘, (14) which completes the proof.

Hence, the best possible updating formula for diagonal matrix 𝐵_𝑘+1 is given by 𝐵_𝑘+1=𝐵_𝑘+𝑦̆_𝑘^𝑇𝑠_𝑘−𝑦̆_𝑘^𝑇𝐵_𝑘𝑦̆_𝑘tr𝑉²_𝑘𝑉_𝑘. (15)

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

Algorithm IDJA

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

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

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

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

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

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

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

Step 8. If ‖Δ𝐹_𝑘‖₂≥𝜖₁ where 𝜖₁=10⁻⁴, compute 𝐵_𝑘+1 as defined by (15); if not, 𝐵_𝑘+1=𝐵_𝑘.

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

3. Convergence Result

This section presents local convergence results of the IDJA methods. To analyze the convergence of these methods, we will make the following assumptions on nonlinear systems 𝐹.

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

We can state the following result on the boundedness of {‖Ψ_𝑘‖_𝐹} by assuming that, without loss of generality, the updating matrix (15) is always used, then we have the following.

Theorem 3. Suppose that {𝑥_𝑘} is generated by Algorithm IDJA where 𝐵_𝑘 is defined by (15). Assume that Assumption 2 holds. There exists 𝛽>0, 𝛿>0, 𝛼>0 and 𝛾>0, such that if 𝑥₀∈𝐸 and 𝐵₀ satisfies ‖𝐼−𝐵₀𝐹^′(𝑥^∗)‖_𝐹<𝛿 for all 𝑥_𝑘∈𝐸 then 𝐼−𝐵_𝑘𝐹^′𝑥^∗_𝐹<𝛿_𝑘, (17) for some constant 𝛿_𝑘>0, 𝑘≥0.

Proof. Since ‖𝐵_𝑘+1‖_𝐹=‖𝐵_𝑘+Ψ_𝑘‖_𝐹, it follows that 𝐵_𝑘+1𝐹≤𝐵_𝑘𝐹+Ψ_𝑘𝐹. (18)
For 𝑘=0 and assuming 𝐵₀=𝐼, we have Ψ₀^(𝑖)=𝑦̆₀^𝑇𝑠₀−𝑦̆₀^𝑇𝐵₀𝑦̆₀tr𝑉²₀𝑦̆₀^(𝑖)2≤𝑦̆₀^𝑇𝑠₀−𝑦̆₀^𝑇𝐵₀𝑦̆₀tr𝑉²₀𝑦̆₀^(max)2, (19) where (𝑦̆₀^(max))² is the largest element among (𝑦̆₀^(𝑖))²,𝑖=1,2,…,𝑛.
After multiplying (19) by (𝑦̆₀^(max))²/(𝑦̆₀^(max))² and substituting tr(𝑉²₀)=∑^𝑛_𝑖=1(𝑦̆₀^(𝑖))⁴, we have Ψ₀^(𝑖)≤𝑦̆₀^𝑇𝑠₀−𝑦̆₀^𝑇𝐵₀𝑦̆₀𝑦̆₀^(max)2∑^𝑛_𝑖=1𝑦̆₀^(𝑖)4𝑦̆₀^(max)4. (20)
Since (𝑦̆₀^(max))⁴/∑^𝑛_𝑖=1(𝑦̆₀^(𝑖))⁴≤1, then (20) turns into Ψ₀^(𝑖)≤𝑦̆₀^𝑇𝐹^′(𝑥)𝑦̆₀−𝑦̆₀^𝑇𝐵₀𝑦̆₀𝑦̆₀^(max)2. (21)
From Assumption 2 and 𝐵₀=𝐼, (21) becomes Ψ₀^(𝑖)≤|𝑐−1|𝑦̆₀^𝑇𝑦̆₀𝑦̆₀^(max)2, (22) where 𝑐=max{|𝑐₁|,|𝑐₂|}.
Since (𝑦̆₀^(𝑖))²≤(𝑦̆₀^(max))² for 𝑖=1,…,𝑛, it follows that Ψ₀^(𝑖)≤𝑛|𝑐−1|𝑦̆₀^(max)2𝑦̆₀^(max)2. (23)
Hence, we obtain Ψ_0𝐹≤𝑛^3/2|𝑐−1|. (24)
Suppose 𝛼=𝑛^3/2|𝑐−1|, then Ψ_0𝐹≤𝛼. (25)
From the fact that ‖𝐵₀‖_𝐹=√𝑛, it follows that 𝐵_1𝐹≤𝛽, (26) where 𝛽=√𝑛+𝛼>0.
Therefore, if we assume that ‖𝐼−𝐵₀𝐹^′(𝑥^∗)‖_𝐹<𝛿, then 𝐼−𝐵₁𝐹^′𝑥^∗_𝐹=𝐼−𝐵₀+Ψ₀𝐹^′𝑥^∗_𝐹≤𝐼−𝐵₀𝐹^′𝑥^∗_𝐹+Ψ₀𝐹^′𝑥^∗_𝐹≤𝐼−𝐵₀𝐹^′𝑥^∗_𝐹+Ψ_0𝐹𝐹^′𝑥^∗_𝐹; (27) therefore, ‖𝐼−𝐵₁𝐹^′(𝑥^∗)‖_𝐹<𝛿+𝛼𝜙=𝛿₁.
Hence, by induction, ‖𝐼−𝐵_𝑘𝐹^′(𝑥^∗)‖_𝐹<𝛿_𝑘 for all 𝑘.

4. Numerical Results

In this section, the performance of IDJA method has been presented, when compared with Broyden’s method (BM), Chord Newton’s method (CN), Newton’s method (NM), and (DQNM) method proposed by [3], respectively. The codes are written in MATLAB 7.4 with a double precision computer; the stopping condition used is 𝑠_𝑘+𝐹𝑥_𝑘≤10⁻⁸. (28)

The identity matrix has been chosen as an initial approximate Jacobian inverse.

We further design the codes to terminates whenever one of the following happens:(i)the number of iteration is at least 200 but no point of 𝑥_𝑘 that satisfies (28) is obtained; (ii)CPU time in seconds reaches 200; (iii)Insufficient memory to initial the run.

The performance of these methods are 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. Spares 1 function of Shin et al. [5]:𝑓_𝑖(𝑥)=𝑥²_𝑖−1𝑖=1,2,…,𝑛,𝑥₀=(5,5,…,5). (29)

Problem 2. Trigonometric function of Spedicato [6] 𝑓_𝑖(𝑥)=𝑛−^𝑛_𝑗=1cos𝑥_𝑗+𝑖1−cos𝑥_𝑖−sin𝑥_𝑖,𝑖=1,…,𝑛,𝑥₀=1𝑛,1𝑛1,…,𝑛. (30)

Problem 3. System of 𝑛 nonlinear equations 𝑓_𝑖(𝑥)=sin1−𝑥_𝑖×^𝑛_𝑖=1𝑥²_𝑖+2𝑥_𝑛−1−3𝑥_𝑛−2−12𝑥_𝑛−4+12𝑥_𝑛−5−𝑥_𝑖ln9+𝑥_𝑖−92exp1−𝑥_𝑛+2𝑖=1,2,…,𝑛,𝑥₀=(0,0,…,0)^𝑇. (31)

Problem 4. System of 𝑛 nonlinear equations 𝑓_𝑖(𝑥)=𝑥²_𝑖−4expsin4−𝑥²_𝑖+sin4−𝑥_𝑖²+𝑖𝑥_𝑛−𝑥_𝑖²+2𝑛−∑^𝑛_𝑖=1𝑥_𝑖cos𝑥_𝑖𝑖=1,…,𝑛,𝑥₀=(2.8,2.8,2.8,…,2.8). (32)

Problem 5. System of 𝑛 nonlinear equations 𝑓_𝑖(𝑥)=^𝑛_𝑖=1𝑥_𝑖𝑥_𝑖−2+cos𝑥_𝑖−2−1𝑖=1,…,𝑛,𝑥₀=(1,1,1,…,1). (33)

Problem 6. System of 𝑛 nonlinear equations 𝑓_𝑖(𝑥)=^𝑛_𝑖=1𝑥²_𝑖−sin𝑥_𝑖−𝑥⁴_𝑖+sin𝑥²_𝑖𝑖=1,…,𝑛,𝑥₀=(.5,.5,.5,…,.5). (34)

Problem 7. System of 𝑛 nonlinear equations 𝑓_𝑗(𝑥)=^𝑛_𝑖=1𝑥²_𝑖−1𝑥_𝑗−1+𝑥_𝑗^𝑛_𝑖=1𝑥_𝑖−1−𝑛+1𝑓_𝑛(𝑥)=^𝑛_𝑖=1𝑥²_𝑖−1𝑥_𝑛−1+cos𝑥_𝑛−1−1𝑗=2,…,𝑛−1,𝑥₀=(.5,.5,.5,…). (35)

Problem 8. System of 𝑛 nonlinear equations 𝑓_𝑖(𝑥)=1−𝑥²_𝑖+𝑥_𝑖+𝑥²_𝑖𝑥_𝑛−2𝑥_𝑛−1𝑥_𝑛−2𝑖=1,2,…,𝑛,𝑥₀=(.5,.5,…,.5). (36)

The numerical results presented in Tables 1, 2, 3, 4, and 5 demonstrate clearly the proposed method (IDJA) shows good improvements, when compared with NM, CN, BM, and DQNM, respectively. In addition, it is worth mentioning, the IDJA method does not require more storage locations than classic diagonal quasi-Newton’s methods. One can observe from the tables that the proposed method (IDJA) is faster than DQNM methods and required little time to solve the problems when compared to the other Newton-like methods and still keeping memory requirement and CPU time in seconds to only 𝑂(𝑛).

Table 1 

Numerical results of NM, CN, BM, DQNM, and IDJA methods.

Table 2 

Numerical Results of NM, CN, BM, DQNM, and IDJA methods.

Table 3 

Numerical Results of NM, CN, BM, DQNM, and IDJA methods.

Table 4 

Numerical results of NM, CN, BM, DQNM, and IDJA methods.

Table 5 

Numerical results of NM, CN, BM, DQNM, and IDJA methods.

5. Conclusions

In this paper, we present an improved diagonal quasi-Newton update via new quasi-Cauchy condition for solving large-scale Systems of nonlinear equations (IDJA). The Jacobian inverse approximation is derived based on the quasi-Cauchy condition. The anticipation has been to further improve the diagonal Jacobian, by modifying the quasi-Cauchy relation so as to carry some additional information from the functions. It is also worth mentioning that the method is capable of significantly reducing the execution time (CPU time), as compared to NM, CN, BM, and DQNM methods while maintaining good accuracy of the numerical solution to some extent. Another fact that makes the IDJA method appealing is that throughout the numerical experiments it never fails to converge. Hence, we can claim that our method (IDJA) is a good alternative to Newton-type methods for solving large-scale systems of nonlinear equations.