An efficient two-step trust-region algorithm for exactly determined consistent systems of nonlinear equations

Journal of Computational and Applied Mathematics 367 (2020) 112470

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

An efficient two-step trust-region algorithm for exactly determined consistent systems of nonlinear equations Somayeh Bahrami, Keyvan Amini

∗

Department of Mathematics, Faculty of Science, Razi University, Kermanshah, Iran

article

info

Article history: Received 16 November 2018 Received in revised form 22 June 2019 Keywords: Nonlinear systems Trust-Region algorithm Levenberg–Marquardt algorithm Global and quadratic convergence

a b s t r a c t In this paper, we propose a modified two-step trust-region algorithm for solving nonlinear systems. Two-step trust-region algorithms, at every iteration, use a trust-region step and an approximate step by saving the Jacobian matrix to have fewer computations. We introduce a convex combination to modify the trust-region subproblems along with an additional criterion to verify whether the first step is accepted or not. We establish global and quadratic convergence of the algorithm under some mild assumptions. Numerical results show that the modified algorithm is efficient and promising. © 2019 Elsevier B.V. All rights reserved.

1. Introduction The problem is to solve an exactly determined consistent system such as F (x) = 0, x ∈ Rn n

(1)

n

where F : R → R is a continuously differentiable function. If F (x) has a zero, then every solution of the nonlinear system (1) is also a solution of the following nonlinear least square problem min Φ (x) = x∈Rn

1 2

∥F (x)∥2 ,

(2)

where ∥.∥ refers to Euclidean norm. Conversely, if x∗ is a minimum of (2) and Φ (x∗ ) = 0, then x∗ solves (1). Throughout the paper, we assume that the solution set of (1) denoted by X ∗ is nonempty, where X ∗ = {x : F (x) = 0}. The nonlinear least square problems and nonlinear equations have been studied extensively. Many different iterative methods have been presented to solve them, for example, the Newton method, quasi-Newton method, trust-region method, Levenberg– Marquardt method, etc., see [1–10]. Among these methods, the trust-region (TR) family is an important and efficient algorithm. The global and local convergence properties and good numerical results produced by TR methods have been shown in many kinds of literature, such as [11–14]. A TR method, at each iterate, based on an approximate model of the function F (xk ) at step k, tries to find an area around the current step xk so that the model appropriately fits the function in this region. Usually, the following subproblem is solved to obtain the trial step dk . min ∥Fk + Jk d∥2 s.t . ∥d∥ ≤ ∆k , where Fk = F (xk ), Jk = ∇ F (xk ) is the Jacobian of F at xk , and ∆k is a positive constant called the TR radius. ∗ Corresponding author. E-mail addresses: [email protected] (S. Bahrami), [email protected] (K. Amini). https://doi.org/10.1016/j.cam.2019.112470 0377-0427/© 2019 Elsevier B.V. All rights reserved.

(3)

2

S. Bahrami and K. Amini / Journal of Computational and Applied Mathematics 367 (2020) 112470

To have a proper comparison between the function Fk and the model function ∥Fk + Jk dk ∥, parameter rk is defined as follows: rk =

Aredk Predk

,

(4)

where Aredk denotes the actual reduction of the function Φ between xk and xk + dk , the new trial step, and Predk is the predicted reduction referring to the reduction of the model function. Formally, Aredk = ∥Fk ∥2 − ∥F (xk + dk )∥2 , and Predk = ∥Fk ∥2 − ∥Fk + Jk dk ∥2 . Parameter rk checks the adequacy of the trust radius. If rk is greater than a small positive constant, c0 , it means that the value of the model function decreases similar to the amount of F at xk . In the case, the trial step is accepted, i.e., xk+1 = xk + dk . But if rk is less than c0 , it is predicted that the behavior of the model function does not fit in with the function Φ . Therefore, the trial step is rejected and xk+1 = xk . Meanwhile, when rk is close to 1, it means that there is an appropriate agreement between the function and the model. In this situation, the TR radius can be increased. But if rk is negative or is a small positive, there is no agreement between the function and the model. So, the radius of the area should be reduced in the hope of finding a better step in the smaller region. Based on these subjects, the step of the classical TR algorithm and its radius are updated as follows: xk + d k

{ xk+1 =

xk

if rk ≥ c0 otherwise,

and

∆k+1

⎧ θ1 ∆k ⎪ ⎪ ⎪ ⎨ ∆k = ⎪ ⎪ ⎪ ⎩ θ2 ∆k

if rk < c1 if c1 ≤ rk ≤ c2 otherwise,

where 0 < c0 < c1 < c2 < 1 and 0 < θ1 < 1 < θ2 . According to (3), one of the significant factors that can lead to high cost in many problems is the computation of the Jacobian matrix. The classical TR needs to compute the Jacobin matrix at every step while generally this cost is expensive. It is especially noticeable when the function is complex or the dimension of the problem is large. Based on this fact, some researchers have proposed some ideas to use the current Jacobian more than once. Indeed, they try to decrease the number of Jacobian computations, see for example [15–19]. Fan and Lu in 2015 introduced one of these ideas. They proposed a two-step TR algorithm where the process is in such a way that the calculated Jacobian is used again, in every step, [20]. Briefly, explanation of the two-step TR algorithm, at every iteration, is as follows: Direction d1k is obtained in the first step by solving the subproblem (3), and then the trial step yk is defined by yk = xk + d1k . To obtain a second direction, d2k , the following subproblem is solved. min ∥F (yk ) + Jk d∥2

¯ k. s.t . ∥d∥ ≤ ∆ Note that they replaced J(yk ) with Jk to avoid computing it. Therefore, the step length considered in the process of TR algorithm will be d1k + d2k . They used the following actual reduction. Aredk = ∥Fk ∥2 − ∥F (xk + d1k + d2k )∥2 , and introduced the predicted reduction for their two-step TR algorithm as follows: Predk = ∥Fk ∥2 − ∥Fk + Jk d1k ∥2 + ∥F (yk )∥2 − ∥F (yk ) + Jk d2k ∥2 .

(5)

Additionally, the radii applied in the first step and the second step, respectively, are as follows:

∆k = µk ∥Fk ∥δ , and

¯ k = µk ∥F (yk )∥δ , ∆

(6)

S. Bahrami and K. Amini / Journal of Computational and Applied Mathematics 367 (2020) 112470

3

where δ ∈ (1/2, 1) and µk is updated by

µk+1 =

θ1 µk

⎧ ⎪ ⎪ ⎪ ⎨

if rk < c1

µk

if c1 ≤ rk ≤ c2

⎪ ⎪ ⎪ ⎩

min{θ2 µk , M }

(7)

otherwise,

that c1 , c2 , θ1 and θ2 are constants such that c2 > c1 > 0 and 0 < θ1 < 1 < θ2 . Moreover, M is a sufficiently large positive constant. The algorithm converges globally, and the convergence rate of the algorithm is analyzed. Numerical results illustrate the efficiency of the algorithm for both singular and nonsingular problems, [20]. In the subsequence, some multi-step algorithms have been introduced, for example [15,18,19,21] Another critical factor in the TR family is a determination of the TR radius in each iteration. Some articles have introduced adaptive radii based on the information of the problem in the current iteration, for example [22–29]. Fan, in 2003, presented the following TR radius

∆k = µk ∥Fk ∥.

(8)

The results of this paper show the numerical superiority of the radius in comparison with the similar proposed radii, [30]. Other TR radii for nonlinear systems were proposed, [27,31,32]. In this paper, we try to modify the objective function of the subproblems to improve the two-step TR algorithm performance. In constructing a new subproblem, we introduce two ideas, first of which is a modified linear approximation of the function, and another one is an appropriate convex combination of the function values in the previous step and the trial step. Meanwhile, we take benefit from an additional criterion to determine whether the proposed convex combination is appropriate or not. We investigate the convergence properties of the algorithm, and present some numerical results to show the efficiency of the modified algorithm. The rest of the paper is organized as follows: In Section 2, we introduce a modified two-step TR algorithm based on some related explanations. In Section 3, the convergence property of the modified algorithm is investigated. Section 4 is devoted to verifying the quadratic convergence of the algorithm. Finally, some numerical experiments are reported in Section 5. 2. A modified two-step trust-region algorithm In this section, we propose an effective two-step TR algorithm by introducing two new ideas along with an additional criterion to modify the second step of the two-step TR algorithm proposed by Fan and Lu [20]. Firstly, we define the trial step zk and the next step xk+1 as follows:

{

zk = xk + pd1k , xk+1 = xk + d1k + pd2k ,

(9)

where d1k is the solution of subproblem (3), and the process to introduce a new direction d2k is described below. In the sense, it is easily seen that xk+1 = yk + pd2k ,

(10)

yk = xk + d1k .

(11)

while

Now, by using a linear approximation of function F based on (10), we can consider the following subproblem min ∥F (yk ) + J(yk )(pd)∥2 s. t .

¯ k. ∥d∥ ≤ ∆

Secondly, to improve the numerical experiments, we use Jk instead of J(yk ) to avoid computing a new Jacobian matrix J(yk ), and propose a convex combination of two recent available values of F instead of F (yk ). According to these, we define d2k as the solution of the following subproblem min ∥Dk + pJk d∥2 (12) s.t .

¯ k, ∥d∥ ≤ ∆

where Dk = (1 − θ )Fk + θ F (zk ) 0 ≤ θ ≤ 1, with θ =

1 p

Dk =

and p ≥ 1, i.e. p−1 p

Fk +

1 p

F (zk ) =

1 p

((p − 1)Fk + F (zk )).

(13)

4

S. Bahrami and K. Amini / Journal of Computational and Applied Mathematics 367 (2020) 112470

Remark 2.1. 1. Note that the trial step zk and the next step xk+1 in (9) are similar to the trial step yk and xk+1 mentioned in the paper [19], while we select p ≥ 1, instead of choosing p ∈ (0, 1]. It is noticeable that selection p ∈ (0, 1] leads to the negative multiple of Fk . We think a convex combination is more useful. 2. Similar to (5), we define the two-step predicted reduction as follows: Predk = ∥Fk ∥2 − ∥Fk + Jk d1k ∥2 + ∥Dk ∥2 − ∥Dk + pJk d2k ∥2 .

(14)

3. Given the suitable properties of the radius (8), we set the radii in both steps by

¯ k = ∆k = µk ∥Fk ∥, ∆

(15)

where µk is updated by (7). 4. We add the following criterion to our algorithm to ensure that ∥F (zk )∥ ≤ ∥Fk ∥ before going to the next step. rˆk =

∥Fk ∥2 − ∥F (zk )∥2 . ∥Fk ∥2 − ∥Fk + Jk d1k ∥2

(16)

It is clear that whenever rˆk ≥ 0, the numerator of (16) is also nonnegative. Motivated by the above points, we can present the modified algorithm as follows: Algorithm 2.1 (The Modified Two-step TR Algorithm). The initial point x0 ∈ Rn . Constants 0 < c0 < c1 < c2 < 1, 0 < θ1 < 1 < θ2 , µ0 > 0, and p ≥ 1. The last xk as the solution of the problem. ¯ k = µ0 ∥Fk ∥. Set k = 0, and compute Fk = F (x0 ), Jk = J(x0 ), ∆k = ∆ If ∥JkT Fk ∥ < ε , stop. Compute d1k by solving subproblem (3), and set zk = xk + pd1k . Compute rˆk by (16). If rˆk < c0 , set xk+1 = xk , and rk = rˆk . Compute µk , ∆k by (7), (15), and go to Step 8. Step 4: Obtain d2k by solving the subproblem (12). Step 5: Compute the predicted reduction Predk by (14), the actual reduction Aredk by

Input: Output: Step 0: Step 1: Step 2: Step 3:

Aredk = ∥Fk ∥2 − ∥F (xk + d1k + pd2k )∥2 ,

(17)

and determine rk by (4). Step 6: If rk ≥ c0 , set xk+1 = xk + d1k + pd2k . Otherwise xk+1 = xk . ¯ k by (15). Step 7: Compute Fk , Jk , and calculate ∆k , ∆ Step 8: Set k = k + 1, and go to Step 1. 3. The convergence properties of Algorithm 2.1 In this section, we investigate the global convergence of Algorithm 2.1 under the following standard assumption. Assumption 3.1. Both function F (x) and the Jacobian J(x) are Lipschitz continuous, i.e., there exist positive constants L1 and L2 such that

∥F (y) − F (x)∥ ≤ L1 ∥y − x∥, and

∥J(y) − J(x)∥ ≤ L2 ∥y − x∥, for all x, y ∈ Rn . Using Assumption 3.1 and the mean value theorem, it is clear that

∥J(x)∥ ≤ L1 ,

(18)

∥F (y) − F (x) − J(x)(y − x)∥ ≤ L2 ∥y − x∥2 ,

(19)

and

for all x, y ∈ R . n

The following Lemma shows a property of the predicted reduction applied in the global convergence theorem. Lemma 3.1.

Under Assumption 3.1, the predicted reduction for all k satisfies

{

Predk ≥ ∥JkT Fk ∥ min ∆k ,

{ ∥JkT Fk ∥ } ∥J T D ∥ } T ¯ k /p, k k . + ∥ J D ∥ min ∆ k k ∥JkT Jk ∥ ∥JkT Jk ∥

S. Bahrami and K. Amini / Journal of Computational and Applied Mathematics 367 (2020) 112470

5

Proof. According to Powell’s result, [13], it is known that

{ ∥J T Fk ∥ } ∥Fk ∥2 − ∥Fk + Jk d1k ∥2 ≥ ∥JkT Fk ∥ min ∆k , kT , ∥ Jk Jk ∥ at the first step. Similarly,

{ ∥J T D ∥ } ¯ k /p, k k , ∥Dk ∥2 − ∥Dk + Jk (pd2k )∥2 ≥ ∥JkT Dk ∥ min ∆ ∥JkT Jk ∥ in the second step. So, by considering (14), the proof is completed. To prove the convergence theorem, we need the following lemma. Lemma 3.2.

Under the conditions of Assumption 3.1, we have

|rk − 1| ≤

∥Fk ∥O(∥d1k ∥2 ) Predk

.

Proof. Using (4), (14), and (17), we have

|rk − 1| = |

Aredk − Predk Predk

|≤

|∥Dk ∥2 − ∥Fk + Jk d1k ∥2 | + |∥F (xk + d1k + pd2k )∥2 − ∥Dk + pJk d2k ∥2 | . |Predk |

(20)

Let us define A ≜ ∥Dk ∥2 − ∥Fk + Jk d1k ∥2 ,

B ≜ ∥F (xk + d1k + pd2k )∥2 − ∥Dk + pJk d2k ∥2 .

In the sequel, we will get two upper bounds for A and B

|A| = |∥Dk ∥2 − ∥Fk + Jk d1k ∥2 | = |∥ Dk − Fk − Jk d1k +(Fk + Jk d1k )∥2 − ∥Fk + Jk d1k ∥2 |,    ≜C

so that

|A| ≤ ∥C ∥2 + 2∥C ∥∥Fk + Jk d1k ∥.

(21)

Using (3), we have

∥Fk + Jk d1k ∥ ≤ ∥Fk ∥ ≤ ∥F0 ∥.

(22)

Now, given (9), (13), and (19), we can write

∥C ∥ = ∥Dk − Fk − Jk d1k ∥ = ∥ =

1 p

p−1 p

Fk +

1 p

F (zk ) − Fk − Jk d1k ∥

∥(p − 1)Fk + F (zk ) − pFk − pJk d1k ∥ =

1 p

∥F (zk ) − Fk − pJk d1k ∥ = O(∥d1k ∥2 ).

This relation, (21), and (22) leads to

|A| ≤ O(∥d1k ∥4 ) + ∥Fk ∥O(∥d1k ∥2 ) = ∥Fk ∥O(∥d1k ∥2 ).

(23)

Moreover,

|B| = |∥F (xk + d1k + pd2k )∥2 − ∥Dk + pJk d2k ∥2 | = |∥ F (xk + d1k + pd2k ) − Fk − Jk (d1k + pd2k ) +(Fk + Jk (d1k + pd2k ))∥2    ≜D

− ∥ Dk + pJk d2k − Fk − Jk (d1k + pd2k ) +(Fk + Jk (d1k + pd2k ))∥2 |    ≜E

≤ ∥D∥2 + 2∥D∥∥Fk + Jk (d1k + pd2k )∥ + ∥E ∥2 + 2∥E ∥∥Fk + Jk (d1k + pd2k )∥.

(24)

Note that we have from (18) and (22)

∥Fk + Jk (d1k + pd2k )∥ ≤ ∥Fk + Jk d1k ∥ + p∥Jk ∥∥d2k ∥ ≤ ∥Fk ∥ + O(∥d2k )∥. Now, similar to ∥C ∥ we attain two relations for ∥D∥ and ∥E ∥.

∥D∥ = ∥F (xk + d1k + pd2k ) − Fk − Jk (d1k + pd2k )∥ = O(∥d1k + pd2k ∥2 ),

(25)

6

S. Bahrami and K. Amini / Journal of Computational and Applied Mathematics 367 (2020) 112470

and

∥d1k + pd2k ∥2 ≤ ∥d1k ∥2 + p2 ∥d2k ∥2 + 2p∥d1k ∥∥d2k ∥. (3), (12), and (15), result in

∥d1k ∥ ≤ ∆k ,

(26)

¯ k = ∆k . ∥d2k ∥ ≤ ∆

(27)

Therefore, we can write

∥d1k + pd2k ∥2 ≤ O(∥d1k ∥2 ), and so,

∥D∥ ≤ O(∥d1k ∥2 ).

(28)

Moreover,

∥E ∥ = ∥Dk + pJk d2k − Fk − Jk (d1k + pd2k )∥ = ∥ =

1 p

∥(p − 1)Fk + F (zk ) − pFk − pJk d1k ∥ =

p−1

1 p

p

Fk +

1 p

F (zk ) + pJk d2k − Fk − Jk d1k − pJk d2k ∥

∥F (zk ) − Fk − pJk d1k ∥ = O(∥d1k ∥2 ).

(29)

Combining (24), (25), (28), and (29), it is deduced

|B| ≤ O(∥d1k ∥4 ) + O(∥d1k ∥2 )(∥Fk ∥ + O(∥d2k ∥)) = ∥Fk ∥O(∥d1k ∥2 ).

(30)

Finally, given (20), (23), and (30), result in A

B

      |∥Dk ∥2 − ∥Fk + Jk d1k ∥2 | + |∥F (xk + d1k + pd2k )∥2 − ∥Dk + pJk d2k ∥2 | |rk − 1| ≤ |Predk |

≤

∥Fk ∥O(∥d1k ∥2 ) . |Predk |

Now, we can present the convergence theorem. Theorem 3.3. lim inf ∥ k→∞

Under the conditions of Assumption 3.1, Algorithm 2.1 either terminates in finite iterations or satisfies JkT Fk

∥ = 0.

Proof. By contradiction, there exists a positive constant τ such that

∥JkT Fk ∥ ≥ τ , ∀k.

(31)

By the same way of the proof of Theorem 2.3 [20], we can conclude that lim µk = 0.

(32)

k→∞

Now, we can apply Lemma 3.2 along with (32) to attain a contradiction. It is deduced from Lemmas 3.1, 3.2, (18), (22), (26), and (31) that

|rk − 1| ≤

∥Fk ∥O(∥d1k ∥2 ) Predk

≤

∥F0 ∥O(∥d1k ∥2 ) ∥JkT Fk ∥ min{∆k ,

∥JkT Fk ∥ ∥JkT Jk ∥

≤ }

∥F0 ∥O(∥d1k ∥2 ) ∥F0 ∥O(∥d1k ∥2 ) ≤ . τ min{∆k , Lτ2 } τ min{∥d1k ∥, Lτ2 } 1

(33)

1

(15), (22), (32), and (33) result rk → 1. Using this fact and noticing (7) and θ2 > 1 lead to µk = M for all sufficiently large k. It is a contradiction with (32), and the proof is completed. 4. Quadratic convergence of the modified algorithm In this section, we analyze the convergence rate of Algorithm 2.1. We need the following assumptions and useful lemmas applied in the theorem of convergence rate. Assumption 4.1. (a) Function F (x) and the Jacobian J(x) are Lipschitz continuous on N(x∗ , b), a neighborhood of a solution x∗ , where 0 < b < 1 is a constant.

S. Bahrami and K. Amini / Journal of Computational and Applied Mathematics 367 (2020) 112470

7

(b) ∥F (x)∥ satisfies the following local error bound condition, i.e., there exists a positive constant c > 0 such that

∥F (x)∥ ≥ c dist(x, X ∗ ),

(34)

for all x ∈ N(x , b), where dist(x, A) = infa∈A ∥x − a∥. Suppose x¯ k is a vector in X ∗ satisfying ∗

∥xk − x¯ k ∥ = dist(xk , X ∗ ).

(35)

Similar to (18) and (19), we have

∥J(x)∥ ≤ L1 ,

(36)

∥F (y) − F (x) − J(x)(y − x)∥ ≤ L2 ∥y − x∥2 ,

(37)

and

for all x, y ∈ N(x , b). ∗

Lemma 4.1. Under the conditions of Assumption 4.1, we have Predk ≥ ∥Fk ∥ min{O(∥d1k ∥), O(∥¯xk − xk ∥)}, for adequately large k. Proof. See the proof of Lemma 3.2 [20]. Lemma 4.2. Under the conditions of Assumption 4.1, there exists a positive constant µ ˆ such that 1 c

≤µ ˆ ≤ µk ,

for all sufficiently large k, where c is the constant used in (34). Proof. Lemma 3.2 along with Lemma 4.1 concludes

|rk − 1| ≤

∥Fk ∥O(∥d1k ∥2 ) Predk

≤

O(∥d1k ∥2 ) min{O(∥d1k ∥), O(∥¯xk − xk ∥)}

.

(38)

On the other hand, using Assumption 4.1(a) together with relations (7), (15), (26) and F (x¯ k ) = 0 results in

∥d1k ∥ ≤ ∆k = µk ∥Fk ∥ = µk ∥Fk − F (x¯ k )∥ ≤ ML1 ∥xk − x¯ k ∥.

(39)

Therefore, (38) and (39) result in rk → 1 for adequately large k. So, given (7) and θ2 > 1 conclude that µk = M. It is deduced that there exists a µ ˆ such that µ ˆ ≤ µk , for sufficiently large k, and because M is an adequately large positive constant, we can select µ ˆ ≥ 1c where c is the constant used in (34), for sufficiently large k. Now, we can represent the convergence rate of Algorithm 2.1. Theorem 4.3. Under the condition 4.1, Algorithm 2.1 satisfies the following relation dist(xk+1 , X ∗ ) ≤ dist(xk , X ∗ )2 . Proof. We have from (11), (34), (35), (37), and Assumption 4.1(a) c ∥¯xk+1 − xk+1 ∥ ≤ ∥F (xk+1 )∥ = ∥F (xk + d1k + pd2k )∥ ≤ ∥F (xk + d1k ) + pJ(yk )d2k ∥ + L2 p2 ∥d2k ∥2





xk+1



(40)

   yk

≤ ∥F (yk ) + pJk d2k ∥ + p∥(J(yk ) − Jk )d2k ∥ + L2 p2 ∥d2k ∥2 ≤ ∥F (yk ) + pJk d2k ∥ + pL2 ∥d1k ∥∥d2k ∥ + L2 p2 ∥d2k ∥2 . On the other hand, with using (9) and Taylor expansion instead of F (zk ) and F (yk ), we have Dk − F (yk ) =

1 p

((p − 1)Fk + F (zk ) − pF (yk )) = (p − 1)O(∥d1k ∥2 ).

(41)

Now, we have from (40) and (41) c ∥¯xk+1 − xk+1 ∥ ≤ ∥F (yk ) − Dk ∥ + ∥Dk + pJk d2k ∥ + pL2 ∥d1k ∥∥d2k ∥ + L2 p2 ∥d2k ∥2

≤ (p − 1)O(∥d1k ∥2 ) + ∥Dk + pJk d2k ∥ + pL2 ∥d1k ∥∥d2k ∥ + L2 p2 ∥d2k ∥2 .

(42)

8

S. Bahrami and K. Amini / Journal of Computational and Applied Mathematics 367 (2020) 112470

Similar to (39), (27) results in

∥d2k ∥ ≤ O(∥¯xk − xk ∥).

(43)

Based on (39), (42), and (43), we can conclude that c ∥¯xk+1 − xk+1 ∥ ≤ ∥Dk + pJk d2k ∥ + O(∥¯xk − xk ∥2 ).

(44)

Now, we show ∥Dk + pJk d2k ∥ ≤ O(∥¯xk − xk ∥2 ). Consequently, from (44), we have

∥¯xk+1 − xk+1 ∥ ≤ O(∥¯xk − xk ∥2 ).

(45)

In the sense, we define the following direction Ak ≜

p−1 p2

(x¯ k − xk ) +

1 p2

(z¯k − zk ).

We have from Lemma 4.2 that µ ˆ c ≥ 1. Hence, Lemma 4.2 along with Assumption 4.1 concludes that

∥¯xk − xk ∥ ≤ µ ˆ c ∥¯xk − xk ∥ ≤ µk ∥Fk ∥ = ∆k ,

(46)

and

∥¯zk − zk ∥ ≤ µ ˆ c ∥¯zk − zk ∥ ≤ µk ∥F (zk )∥. Also, using rˆk in step 3 of Algorithm 2.1 results in ∥F (zk )∥ ≤ ∥Fk ∥, hence

∥¯zk − zk ∥ ≤ µk ∥Fk ∥ = ∆k . Therefore,

∥Ak ∥ ≤

1 p

∆k ≤ ∆k .

So, Ak is a feasible direction for subproblem (12). Hence, it can be written

∥Dk + pJk d2k ∥ ≤ ∥Dk + pJk Ak ∥ = ∥ ≤ ≤

p−1 p p−1 p

p−1 p

Fk +

1 p

F (zk ) +

p−1 p

Jk (x¯ k − xk ) +

1 p

Jk (z¯k − zk )∥

1

∥Fk + Jk (x¯ k − xk )∥ + ∥F (zk ) + Jk (z¯k − zk )∥ p 1

1

p

p

∥Fk + Jk (x¯ k − xk )∥ + ∥F (zk ) + J(zk )(z¯k − zk )∥ + ∥(Jk − J(zk ))(z¯k − zk )∥.

Then from the Lipschitz assumption, we can write

∥Dk + pJk d2k ∥ ≤ O(∥¯xk − xk ∥2 ) + O(∥¯zk − zk ∥2 ) + O(∥d1k ∥)O(∥¯zk − zk ∥). On the other hand, from (34), (36), (37), and (46), we get c ∥¯zk − zk ∥ ≤ ∥F (zk )∥ = ∥F (xk + pd1k )∥ ≤ ∥F (xk ) + pJk d1k ∥ + L2 p2 ∥d1k ∥2

≤ ∥F (xk ) + Jk d1k ∥ + (p − 1)∥Jk d1k ∥ + L1 p2 ∥d1k ∥2 ≤ ∥F (xk ) + Jk (x¯ k − xk )∥ + (p − 1)∥Jk ∥∥d1k ∥ + L2 p2 ∥d1k ∥2 . So, c ∥¯zk − zk ∥ ≤ O(∥¯xk − xk ∥2 ) + O(∥¯xk − xk ∥) + O(∥¯xk − xk ∥2 ). Therefore, we have

∥¯zk − zk ∥ ≤ O(∥¯xk − xk ∥), and

∥Dk + pJk d2k ∥ ≤ O(∥¯xk − xk ∥2 ). Therefore, the quadratic convergence of Algorithm 2.1 by (45) is obtained. 5. Numerical experiments In this section, we report some numerical results to show the performance of Algorithm 2.1. We compare Algorithm 2.1 with two following TR algorithms • MNTR, a two-step TR algorithm introduced by Fan and Lu, [20].

S. Bahrami and K. Amini / Journal of Computational and Applied Mathematics 367 (2020) 112470

9

• NTR, the classical TR algorithm while radius in every step is calculated by (6), see [20]. In all algorithms, the following parameter values are used, [11]. c0 = 10−4 , c1 = 0.25, c2 = 0.75, θ1 = 1/4, θ2 = 2, δ = 1, M = 50, µ0 = 0.1∥g0 ∥, where the TR subproblems are solved by nearly exact solution method [9]. The algorithms are terminated when the number of iterations exceeds 10 000 or norm of the Jacobian is less than ϵ , i.e., ∥JkT Fk ∥ ≤ 10−5 . All codes are written in MATLAB R2015.b on a 2.20 GHz core i5 processor laptop with 6 GB RAM. As mentioned in Section 4, we analyzed the rate of convergence under the error bound condition around the solution x∗ , which is weaker than the nonsingularity of the Jacobian. Therefore, it is logical that the numerical results are based on some test problems with singular Jacobian matrix at the point x∗ . Schnabel and Frank constructed some standard singular test problems as follows: Fˆ (x) = F (x) − J(x∗ )A(AT A)−1 AT (x − x∗ ), where F (x) is a standard nonsingular test problem for nonlinear equations introduced by Moré et al. [33], x∗ is a point satisfying F (x∗ ) = 0, and A ∈ Rn×k is a full column rank matrix such that 1 ≤ k ≤ n. It is clear that Fˆ (x∗ ) = 0, and ˆ ∗ )) = n − k where J(x ˆ ∗ ), the Jacobian of Fˆ (x) in x∗ , is computed by rank(J(x

ˆ ∗ ) ≜ J(x∗ )(I − A(AT A)−1 AT ). J(x Similar to [20,21,34], we use the process of [34] to construct some singular test problems, and consider two following selections for matrix A A = [1, 1, . . . , 1]T ∈ Rn×1 , and A=

1, 1, 1, 1, . . . , 1 1, −1, 1, −1, . . . , ±1

[

]T

∈ Rn×2 ,

ˆ ∗ )) = n − 1 and n − 2, respectively. which result in rank(J(x We test Algorithm 2.1 with different amounts of parameter p ∈ [1, 2]. When p = 1, the algorithm is very close to MNTR algorithm, and its numerical results are similar to this algorithm. Also, our numerical experiments indicate that the best results for Algorithm 2.1 are obtained with p = 1.7. We present these results along with the ones of NTR algorithm and MNTR algorithm in Tables 1 and 2. Table 1 demonstrates the results of problems with rank n − 1, and Table 2 exhibits the results of problems with rank n − 2, respectively. We run the test problems with dimension 1000 based on six starting points −100x0 , −10x0 , −x0 , x0 , 10x0 , 100x0 , where x0 is the standard initial point suggested in the paper [33]. In Helical valley problem, we spot dimension 999, because its dimension must be a multiple of 3. Symbols NF and NJ in the tables represent the number of function and Jacobian evaluations, respectively. It is known that the cost of computing every Jacobian matrix in nonlinear equation problems is usually n times of every function computation cost. Hence, we enter the total evaluations NT = NF + n ∗ NJ as another factor in the tables. Moreover, nc is a symbol which shows that the method does not converge to x∗ mentioned in the paper [33], and we use the symbol M to indicate that the number of iterations is more than 10 000. It is noticeable that when three algorithms in one row of the table fail, we eliminate this row from the tables. For example, at extended Powell badly problem, three algorithms are nc for all six starting points. Also, at Discrete boundary value problem, all the algorithms have more than 10 000 iterations for −100x0 , 100x0 . Therefore, we eliminate them from the tables. Also, when the optimal solution is the same with the starting point, the algorithms terminate at the starting step, so we remove these rows from the tables; −x0 in Helical valley problem, and x0 in Discrete boundary value problem. In order to make better comparisons and illustrate the results of the tables, Dolan and Moré in 2002 proposed a plot of the performance profile. This plot illustrates and compares the major characteristics of the algorithms. For ns solvers and np problems, they proposed the performance profile P : R → [0, 1] as follows: Let P and S be the set of problems and the set of solvers, respectively. For each problem p ∈ P and each solver s ∈ S, they set tp,s = the number of computations required (iteration number, the function computations, CPU time, etc.) to solve problem p by solver s. They defined the performance ratio by rp,s =

tp,s mins∈S {tp,s }

,

and the performance profile by P(τ ) =

1 ns

size{p ∈ P|rp,s ≤ τ }

∀τ ∈ R,

where size of a set stands for the number of elements of the set. Note that P(τ ) is the probability for solver s ∈ S that a performance ratio rp,s is within a factor τ ∈ R of the best possible ratio. A parameter rM ≥ rp,s for all p and all s is chosen, and rp,s = rM if and only if algorithm s does not solve problem p, where the selection of rM does not affect the performance evaluations, [35].

10

S. Bahrami and K. Amini / Journal of Computational and Applied Mathematics 367 (2020) 112470

Table 1 Numerical results for singular nonlinear systems with rank(J(x∗ )) = n − 1. Test problem

Extended Rosenbrock

Extended Powell singular

Algorithm

2.1

MNTR

Algorithm

NTR

Algorithm

Starting point

NF /NJ

NT

NF /NJ

NT

NF /NJ

NT

−100x0 −10 x0

22/11 18/9 12/6 29/11 18/9 22/11

11 022 9 018 6 012 11 029 9 018 11 022

26/13 20/10 12/6 18/9 20/10 26/13

13 026 10 020 6 012 9 018 10 020 13 026

17/17 14/14 8/8 12/12 14/14 17/17

17 017 14 014 8 008 12 012 14 014 17 017

18/9 16/8 12/6 12/6 16/8 18/9

9 018 8 016 6 012 6 012 8 016 9 018

24/12 20/10 14/7 14/7 20/10 24/12

12 024 10 020 7 014 7 014 10 020 12 024

17/17 13/13 10/10 10/10 13/13 17/17

17 017 13 013 10 010 10 010 13 013 17 017

26/13 22/11 18/9 20/10 22/11 26/13

13 026 11 022 9 018 10 020 11 022 13 026

30/15 26/13 20/10 22/11 26/13 30/15

15 030 13 026 10 020 11 022 13 026 15 030

21/21 17/17 13/13 14/14 17/17 21/21

21 021 17 017 13 013 14 014 17 017 21 021

20/10 14/7 18/9 24/12 16/8 12/6

10 020 7 014 9 018 12 024 8 016 6 012

24/12 18/9 22/11 26/13 18/9 12/6

12 024 9 018 11 022 13 026 9 018 6 012

16/16 11/11 12/12 14/14 10/10 9/9

16 016 11 011 12 012 14 014 10 010 9 009

27/11 22/10 29/14 36/17 34/13 72/28

11 027 10 022 14 029 17 036 13 034 28 072

28/11 22/9 32/13 36/17 48/19 54/22

11 028 9 022 13 032 17 036 19 048 22 054

20/18 20/18 19/16 18/18 35/27 38/27

18 020 18 020 16 019 18 018 27 035 27 038

20 040 18 036 17 034 17 034 18 036 20 040

50/20 46/23 44/22 42/21 44/22 50/25

25 050 23 046 22 044 21 042 22 044 25 050

35/35 32/32 30/30 29/29 31/31 35/35

35 035 32 032 30 030 29 029 31 031 35 035

-x0 x0 10x0 100x0

−100x0 −10x0 -x0 x0 10x0 100x0

−100x0 −10x0 Extended Wood

Discrete integral equation

-x0 x0 10x0 100x0

−100x0 −10x0 -x0 x0 10x0 100x0

−100x0 −10x0 Trigonometric

-x0 x0 10x0 100x0

−100x0 −10x0 Variably dimensioned

-x0 x0 10x0 100x0

40/20 36/18 34/17 34/17 36/18 40/20

Brown almost linear

-x0 x0

23/7 25/8

7 023 8 025

12/6 12/6

6 012 6 012

−100x0 −10x0

20/10 12/6 41/15 14/7 20/10 28/14

10 020 6 012 15 041 7 014 10 020 14 028

20/10 12/6 52/18 16/8 26/13 34/17

10 020 6 012 18 052 8 016 13 026 17 034

18/9 34/12 22/9 12/6 14/7 18/9

9 018 12 034 9 022 6 012 7 014 9 018

nc 18/9 nc 12/6 18/9 24/12

Broyden banded

-x0 x0 10x0 100x0

−100x0 −10x0 Broyden tridiagonal

Extended Helical valley Discrete boundary value

-x0 x0 10x0 100x0

−100x0 −10x0 x0 100x0

nc nc 64/25 nc

–

4/2 2/1 nc 14/7

−10x0 -x0 10x0

6746/3373 964/482 5360/2680

3 379 746 482 964 2 685 360

13118/6559 1864/932 10424/5212

Notes: NF = the number of function computations, nc = problem does not converge to x∗ , NJ = the number of Jacobian computations, M = problem needs more than 10 000 iterations, NT = NF + n ∗ NJ = the total computations, x0 = the initial point suggested in the paper [33].

– – 25 039

– 9 018 – 6 012 9 018 12 024 2 002 1 001 – 7 007 6 572 118 933 864 5 222 424

8/8 8/8

8 008 8 008

15/15 9/9 25/16 11/11 16/16 23/23

15 015 9 009 16 025 11 011 16 016 23 023

16/16 13/13 15/11 8/8 13/13 16/16

16 016 13 013 11 015 8 008 13 013 16 016

2/2 1/1 nc nc

2 000 1000 – –

M 1861/1861 M

– 1 862 861 –

S. Bahrami and K. Amini / Journal of Computational and Applied Mathematics 367 (2020) 112470

11

Table 2 Numerical results for singular nonlinear systems with rank(J(x∗ )) = n − 2. Algorithm

2.1

MNTR

Algorithm

NTR

Algorithm

Test problem

Starting point

NF /NJ

NT

NF /NJ

NT

NF /NJ

NT

Extended Rosenbrock

x0

48/24

24 048

84/42

42 084

82/82

82 082

−100x0 −10x0 -x0 x0 10x0 100x0

18/9 16/8 12/6 12/6 16/8 18/9

9 018 8 016 6 012 6 012 8 016 9 018

24/12 20/10 14/7 14/7 20/10 24/12

12 024 10 020 7 014 7 014 10 020 12 024

17/17 13/13 10/10 10/10 13/13 17/17

17 017 13 013 10 010 10 010 13 013 17 017

-x0

48/24

24 048

86/43

−100x0 −10x0

742/371 736/368 732/366 662/331 748/374

Extended Powell singular

Extended wood Discrete integral equation

-x0 x0 10x0

−100x0 −10x0

371 742 368 736 366 732 331 662 374 748

43 086

1432/716 1426/713 1420/710 1282/641 1452/726

717 432 714 426 711 420 642 282 727 452

Trigonometric

-x0 x0 10x0 100x0

38/14 24/9 26/13 34/17 49/18 90/35

Brown almost linear

-x0 x0

23/7 25/8

7 023 8 025

12/6 12/6

6 012 6 012

−100x0 −10x0

20/10 12/6 41/15 14/7 20/10 28/14

10 020 6 012 15 041 7 014 10 020 14 028

20/10 12/6 52/19 16/8 26/13 34/17

10 020 6 012 19 052 8 016 13 026 17 034

15/15 9/9 25/16 11/11 17/17 23/23

40/20 36/18 34/17 34/17 36/18 40/20

20 040 18 036 17 034 17 034 18 036 20 040

nc 46/23 44/22 42/21 44/22 50/25

– 23 046 22 044 21 042 22 044 25 050

nc 32/32 30/30 29/29 31/31 nc

18/9 30/12 29/11 12/6 14/7 18/9

9 018 12 030 11 029 6 012 7 014 9 018

nc 36/13 30/10 12/6 18/9 24/12

– 13 036 10 030 6 012 9 018 12 024

23/18 20/15 65/48 8/8 13/13 16/16

Broyden banded

-x0 x0 10x0 100x0

−100x0 −10x0 Variably dimensioned

-x0 x0 10x0 100x0

−100x0 −10x0 Broyden tridiagonal

-x0 x0 10x0 100x0

Extended Helical valley

−100x0 −10x0

Discrete boundary value

−10x0

100x0 -x0 10x0

14 038 9 024 13 026 17 034 18 049 35 090

32/13 22/9 28/12 38/17 62/25 78/31

13 032 9 022 12 028 17 038 25 062 31 078

83/883 1430/1430 1425/1425 1418/1418 1279/1279 1454/1454

nc nc nc

– – –

4/2 2/1 14/7

6746/3373 964/482 5360/2680

3 379 746 482 964 2 685 360

13118/6559 1864/932 10424/5212

2 002 1 001 7 007 6 572 118 933 864 5 222 424

83 083 1 431 430 1 426 425 1 419 418 1 280 279 1 455 454

20/18 20/18 18/15 21/19 31/26 29/19

18 020 18 020 15 018 19 021 26 031 19 029

8/8 8/8

8 008 8 008 15 015 9 009 16 025 11 011 17 017 23 023 – 32 032 30 030 29 029 31 031 – 18 023 15 020 48 065 8 008 13 013 16 016

2/2 1/1 nc

2 000 1 000 –

M 1861/1861 M

– 1 862 861 –

Notes: NF = the number of function computations, nc = problem does not converge to x∗ , NJ = the number of Jacobian computations, M = problem needs more than 10 000 iterations, NT = NF + n ∗ NJ = the total computations, x0 = the initial point suggested in the paper [33].

In this paper, we also plot the performance profile of Dolan and Moré based on the total computations, NT , in Figs. 1 and 2. A glimpse into the tables and figures shows that based on NT, Algorithm 2.1 works better than the other algorithms. In Fig. 1, according to the first point on the left side of the graphs, we observe that Algorithm 2.1 solves 82% of the considered problems in the least number of the total computations NT , while NTR and MNTR algorithms solve only 3% and 24% of the problems in the least number of NT , respectively. On the other hand, the highest points on the right side of graphs in Fig. 1 demonstrate that NTR algorithm and Algorithm 2.1 solve 95% of the considered problems successfully while MNTR algorithm solves 93% of the problems. Also, Fig. 2 displays that Algorithm 2.1 solves 78% of the problems in the minimum number of NT while NTR and MNTR algorithms solve only 6% and 22% of the problems in the minimum number of NT ,

12

S. Bahrami and K. Amini / Journal of Computational and Applied Mathematics 367 (2020) 112470

Fig. 1. Performance profile based on NT of Table 1.

Fig. 2. Performance profile based on NT of Table 2.

respectively. Furthermore, Fig. 2 reveals MNTR algorithm, Algorithm 2.1, and NTR algorithm successfully solve 96%, 93%, and 89% of the considered problems, respectively. Note that Algorithm 2.1 grows up faster than the other two algorithms in both cases, which means even when Algorithm 2.1 is not the best algorithm, its performance index is close to the one of the best algorithm. In another view, we know that Algorithm 2.1 and MNTR algorithm are two-step TR algorithms while NTR algorithm is a one-step TR algorithm. Based on this remark, we can independently compare these two algorithms based on NF , NJ and NT . The data of Table 1 show that Algorithm 2.1 is better, which solves 70%, 63% and 66% of the problems in the least number of NF , NJ and NT , respectively, in comparison with MNTR algorithm. Moreover, based on Table 2, Algorithm 2.1

S. Bahrami and K. Amini / Journal of Computational and Applied Mathematics 367 (2020) 112470

13

also solves 71%, 64% and 66% of the problems in the minimum number of NF , NJ and NT , respectively, in comparison with MNTR algorithm. 6. Conclusion In this study, we introduce a modified two-step TR algorithm for solving determined consistent systems of nonlinear equations. This method tries to present an algorithm by saving the Jacobian computations and adopting a convenient TR radius. We modify the TR subproblem and apply an additional criterion to find out whether the first step is an appropriate step or not. The convergence properties and the quadratic convergence rate of the algorithm are proved. Numerical results demonstrate the superiority of the modified algorithm. References [1] C.G. Broyden, Quasi-Newton methods and their applications to function minimization, Math. Comp. 21 (1967) 577–593. [2] E.M. Gertz, A quasi-Newton trust-region method, Math. Program. 100 (3) (2004) 447–470. [3] C. Kanzow, N. Yamashita, M. Fukushima, Levenberg–Marquardt methods for constrained nonlinear equations with strong local convergence properties, J. Comput. Appl. Math. 172 (2) (2004) 375–397. [4] C.T. Kelley, Iterative methods for linear and nonlinear equations, in: Frontiers in Applied Mathematices, Vol. 16, SIAM, Philadelphia, 1995. [5] C.T. Kelley, Solving nonlinear equations with Newton’s method, in: Fundamentals of Algorithms, SIAM, Philadelphia, 2003. [6] K. Levenberg, A method for the solution of certain nonlinear problems in least squares, Quart. Appl. Math. 2 (1944) 164–168. [7] D.W. Marquardt, An algorithm for least-squares estimation of nonlinear parameters, J. Soc. Ind. Appl. Math. 11 (2) (1963) 431–441. [8] J.J. Moré, D.C. Sorensen, Computing a trust region step, SIAM J. Sci. Stat. Comput. 4 (3) (1983) 553–572. [9] J. Nocedal, J.S. Wright, Numrical Optimization, Springer, New York, 1999. [10] Y.X. Yuan, Trust region algorithms for nonlinear equations, Conference on Scientific Computation, Hong Kong, 17-19 March, 1994. [11] A. Conn, N. Gould, Ph.L. Toint, Trust Region Methods, SIAM, Philadelphia, 2000. [12] M.J.D. Powell, On the global convergence of trust region algorithms for unconstrained optimization, Math. Program. 29 (3) (1984) 297–303. [13] M.J.D. Powell, Convergence properties of a class of minimization algorithms, in: O.L. Mangasarian, R.R. Meyer, S.M. Robinson (Eds.), Nonlinear Programming, Vol. 2, Academic Press, New York, 1975, pp. 1–27. [14] A.G. Shultz, R.B. Schnable, A.R.H. Byrd, A family of trust-region-based algorithm for unconstrained minimization with strong global convergence properties, SIAM J. Numer. Anal. 22 (1) (1985) 47–67. [15] K. Amini, F. Rostami, Three-steps modified Levenberg–Marquardt method with a new line search for systems of nonlinear equations, J. Comput. Appl. Math. 300 (2016) 30–42. [16] J.A. Ezquerro, M.A. Hernández, An optimization of Chebyshev’s method, J. Complexity 25 (4) (2009) 343–361. [17] J.Y. Fan, The modified Levenberg–Marquardt method for nonlinear equations with cubic convergence, Math. Comp. 81 (277) (2012) 447–466. [18] X. Yang, A higher order Levenberg Marquardt method for nonlinear equations, Appl. Math. Comput. 219 (22) (2013) 10682–10694. [19] X. Zhao, J. Fan, On the multi point Levenberg–Marquardt method for singular nonlinear equations, J. Comput. Appl. Math. 36 (1) (2017) 203–223. [20] J.Y. Fan, N. Lu, On the modified trust region algorithm for nonlinear equations, Optim. Methods Softw. 30 (3) (2015) 478–491. [21] K. Amini, F. Rostami, A new modified two steps Levenberg–Marquardt method for nonlinear equations, J. Comput. Appl. Math. 288 (2015) 341–350. [22] M. Ahookhosh, K. Amini, A nonmonotone trust region method with adaptive radius for unconstrained optimization problems, Comput. Math. Appl. 60 (3) (2010) 411–422. [23] A. Sartenaer, Automatic determination of in initial trust region in nonlinear programming, SIAM J. Sci. Comput. 18 (6) (1997) 1788–1803. [24] Z.J. Shi, J.H. Guo, A new trust region methods for unconstrained optimization, J. Comput. Appl. Math. 213 (2) (2008) 509–520. [25] Ph.L. Toint, Nonlinear step size control, trust regions and regularization for unconstrained optimization, Optim. Methods Softw. 28 (1) (2013) 82–95. [26] X.S. Zhang, J.L. Zhang, L.Z. Liao, An adaptive trust region method and its convergence, Sci. Chin. 45 (5) (2002) 620–631. [27] J.L. Zhang, Y. Wang, A new trust region method for nonlinear equations, Math. Methods Oper. Res. 58 (2) (2003) 283–298. [28] Morteza Kimiaei, Nonmonotone self-adaptive Levenberg–Marquardt approach for solving systems of nonlinear equations, Numer. Funct. Anal. Optim. 39 (1) (2018) 47–66. [29] H. Esmaeili, M. Kimiaei, A new adaptive trust-region method for system of nonlinear equations, Appl. Math. Model. 38 (11–12) (2014) 3003–3015. [30] J.Y. Fan, A modified Levenberg–Marquardt algorithm for singular system of nonlinear equations, J. Comput. Math. 21 (5) (2003) 625–636. [31] J.Y. Fan, J.Y. Pan, An improved trust region algorithm for nonlinear equations, Comput. Optim. Appl. 48 (1) (2011) 59–70. [32] J.Y. Fan, Y.X. Yuan, A new trust region algorithm with trust region radius converging to zero, in Proceedings of the 5th Conference on Optimization: Techniques and Applications, Hong Kong, 2001, 786-794. [33] J.J. Moré, B.S. Garbow, K.E. Hillstrom, Testing unconstrained optimization software, ACM Trans. Math. Softw. 7 (1) (1981) 17–41. [34] R.B. Schnabel, P.D. Frank, Tensor methods for nonlinear equations, SIAM J. Numer. Anal. 21 (1984) 815–843. [35] E.D. Dolan, J.J. Moré, Benchmarking optimization software with performance profiles, Math. Program. 91 (2) (2002) 201–213.