A family of global convergent inexact secant methods for nonconvex constrained optimization

Introduction

We assume that the nonlinear equality constrained problem to be solved is stated as

min x ∈ R n f ( x ) s . t . c ( x ) = 0

(1)

The reduced Hessian methods in successive quadratic programming (SQP)^1,2 and secant methods (two-step algorithms) as defined in Fontecilla ³ are two of the most successful methods for solving problem (1). Compared with the widely reduced Hessian methods in which the orthonormal basis for the tangent space of the constraints at the current point x_k changes continuously with k,^4,5 the secant methods have a main advantage which rests in the use of an orthogonal projection operator which is continuous.

Two basic approaches, namely the line search and trust region, have been developed in order to ensure global convergence towards local minima. In Byrd et al., ⁶ an algorithm, which is based on a characterization of inexact sequential quadratic programming (SQP) steps that can ensure global convergence, has been presented for large-scale equality constrained optimization. An exact penalty function was used to determine if a given inexact step makes sufficient progress toward a solution of the nonlinear program. The inexact Newton method was initially proposed by Dembo et al. ⁷ to solve large systems of nonlinear equations. Dembo and Steihaug ⁸ used that procedure to solve unconstrained minimization problems. Currently, Byrd et al. ⁹ proposed an inexact Newton method for nonconvex equality constrained optimization. The method in Byrd et al. ⁹ was allowed for the presence of negative curvature without requiring information about the inertia of the primal-dual iteration matrix for nonconvex problems. Then, an inexact Newton method with a filter line search algorithm is proposed in Wang et al. ¹⁰ for the same problem. In Gu and Zhu ¹¹ and Wang et al., ¹² inexact secant algorithms are proposed for solving the large-scale nonlinear systems of equalities and inequalities and nonlinear constrained convex optimization, respectively.

In this paper, we propose a class of new inexact secant methods with Armijo line search for nonconvex problems. The methods are globalized by line search technique and adopting l₁ penalty function as a merit function. The paper is organized as follows. In the next section, these algorithms are developed. Then, we analyze the global convergence properties. The results of numerical experience with these methods are discussed in Numerical results section. Final remarks are provided in the Conclusions section.

The proposed algorithms

Throughout the paper, we denote the Euclidean vector or matrix norm by | | · | | , l₁ norm by | | · | | 1 , the gradient ∇ f ( x ) by g(x). l is the Lagrangian function defined for x ∈ R n and λ ∈ R m by

l ( x , λ ) = f ( x ) + λ T c ( x )

The gradient of the Lagrangian function is denoted ∇ x l ( x , λ ) . We will be using the BFGS secant update defined by

W k + 1 = W k − W k d k d k T W k d k T W k d k + y k y k T y k T d k

The following is the basic idea about secant methods: At the k − th each iteration

λ k + 1 = U ( x k , λ k , W k )

(2a)

W k u ^ k = − ∇ x l ( x k , λ k + 1 )

(2b)

u k = P k u ^ k

(2c)

v k = − A k † c k

(2d)

y k = ∇ x l ( x k + u k , λ k + 1 ) − ∇ x l ( x k , λ k + 1 )

(2e)

W k + 1 = B F G S ( u k , y k , W k )

(2f)

x k + 1 = x k + u k + v k

(2g)

P ( x ) = I − A ( x ) [ A ( x ) T A ( x ) ] − 1 A ( x ) T

(3a)

P ( x ) = I − W − 1 A ( x ) [ A ( x ) T W − 1 A ( x ) ] − 1 A ( x ) T

(3b)

The multiplier updates U ( x k , λ k , W k ) in equation (2a) can be chosen from one of the following updates

Projection update:

λ k + 1 = − ( A k T A k ) − 1 A k T g k

(4a)

Null-space update:

λ k + 1 = − ( A k T W k − 1 A k ) − 1 A k T W k − 1 g k

(4b)

Newton update:

λ k + 1 = ( A k T W k − 1 A k ) − 1 ( c k − A k T W k − 1 g k )

(4c)

A k † = A k ( A k T A k ) − 1

(5a)

A k † = W k − 1 A k ( A k T W k − 1 A k ) − 1

(5b)

While choosing different projection operators in equation (3), different multiplier updates in equation (4), and different pseudo-inverse in equation (5), we can obtain six algorithms which will be listed in Table 1.

OPEN IN VIEWER

Table 1.

Algorithms.

Algorithm	λ k + 1	P(x)	A k †	Algorithm	λ k + 1	P(x)	A k †
1	(4b)	I	(5b)	2	(4a)	(3b)	(5b)
3	(4c)	(3a)	(5a)	4	(4b)	I	(5a)
5	(4a)	(3b)	(5a)	6	(4c)	(3a)	(5b)

Next, we outline the algorithm and globalization strategy. An integral part of the approach is the mechanism used to determine if a trial step d is acceptable during a given iteration. We commonly use a merit function to determine whether a step is acceptable. The nondifferentiable l₁ merit function and the Fletcher’s exact and differentiable function are representative of most of merit function used in practice. The evaluation of the l₁ merit function, when compared to the Fletcher’s augmented Lagrangian merit function, is very inexpensive. One of the former’s potential drawbacks is that it can suffer from the Maratos effect. Several strategies have been used to remove the damaging effects of the Maratos effect. Therefore, we adopt the l₁ penalty function as the merit function

ϕ ( x ; ω ) : = f ( x ) + ω | | c ( x ) | | 1

In this paper, we perform a backtracking line search to compute a step length coefficient α_k satisfying the Armijo condition

ϕ ( x k + α k d k ; ω k ) ≤ ϕ ( x k ; ω k ) + β α k D ϕ ( d k ; ω k )

m k ( d ; ω ) = f k + g k T d + ω | | c k + A k T d | | 1

Hence, we can estimate the reduction in the merit function by evaluating

Δ m k ( d k ; ω k ) = m k ( 0 ; ω k ) − m k ( d k ; ω k ) = − g k T d k + ω k | | c k | | 1

(6)

Once an acceptable step is obtained, we must ensure that a positive α_k can be calculated to satisfy the Armijo condition. We will consider this issue in Lemma 1.

For nonconvex equality constrained optimization, that is

d k T W k d k < 0 , A k T d k = 0 for some d k

(7)

using the method proposed in Byrd et al., ⁹ we have the reduction condition

Δ m k ( d k ; ω k ) ≥ max ⁡ { 1 2 d k T W k d k , θ | | u k | | 2 } + τ ω k | | c k | |

(8)

| | u k | | 2 ≤ ϒ k ≤ | | d k | | 2

In this paper, we will compute the search direction inexactly, which means that equation (2b) will be substituted by

W k u ^ k = − ∇ x l ( x k , λ k + 1 ) + r k

(9)

The model reduction condition in this work is

Δ m k ( d k ; ω k ) ≥ max ⁡ { 1 2 d k T W k d k , θ ϒ k } + τ ω k | | c k | |

If the model reduction condition is satisfied for the most recent value of the penalty parameter, we have the option of increasing the penalty parameter in order to satisfy the model reduction condition. This, however, should only be done in two circumstances. If 1 2 d k T W k d k ≥ θ ϒ k , it is safe to assume that the problem is sufficiently convex and so d_k should be accepted. If 1 2 d k T W k d k ≤ θ ϒ k , we only consider increasing ω if v_k is a significant component of d_k. We can express this condition as ϒ k ≤ ψ ν k for some ψ > 0 , where 0 ≤ ν k ≤ ‖ v k | | 2 is any lower bound for the squared norm of the normal step component. If none of the above conditions can be satisfied for the model reduction condition, we modify W_k. (It is assumed that after a finite number of modifications we have W k ≻ 2 θ I .) One technique is to modify W_k to increase some or all of its eigenvalues so that the resulting matrix is closer to being positive definite.

We now formally state the global inexact secant methods for solving equation (1).

Algorithms

s1. Initialize. Choose a starting point x₀. Set constants β ∈ ( 0 , 1 ) , τ ∈ ( 0 , 1 ) , ς ∈ ( 0 , 1 ) and θ , ϵ , ε > 0 . Give λ₀ and ω − 1 > 0 , and let k : = 0 .

s2. If ‖ ( ∇ x l ( x k , λ k ) c k ) ‖ ≤ ε , then stop.

s3. Evaluate f_k, g_k, c_k, A_k and W_k from equation (2f). Then compute the multiplier λ k + 1 from (4).

s4. Find an inexact step d_k from

W k u ^ k = − ∇ x l ( x k , λ k + 1 ) + r k

(10a)

u k = P k u ^ k

(10b)

v k = − A k † c k

(10c)

d k = u k + v k

(10d)

Condition I

Δ m k ( d k ; ω k ) ≥ max ⁡ { 1 2 d k T W k d k , θ ϒ k } + τ ω k | | c k | | | | r k | | ≤ η k | | ∇ x l ( x k , λ k + 1 ) | |

(11)

Condition II

| | r k | | ≤ t 2 | | c k | | , and 1 2 d k T W k d k ≥ θ ϒ k or ϒ k ≤ ψ ν k

(12)

s5. If Condition II is satisfied and

ω k ( 1 − ϱ ) | | c k | | 1 < g k T d k + max ⁡ { 1 2 d k T W k d k , θ ϒ k }

for ω k = ω k − 1 , where τ = ϱ ( 1 − ϵ ) , then set

ω k : = g k T d k + max ⁡ { 1 2 d k T W k d k , θ ϒ k } ( 1 − ϱ ) | | c k | | 1 + 10 − 4

(13)

ϕ ( x k + α k d k ; ω k ) ≤ ϕ ( x k ; ω k ) + β α k D ϕ ( d k ; ω k )

s7. Set x k + 1 = x k + α k d k and k : = k + 1 . Go to s2.

Lemma 1. The directional derivative of the merit function ϕ ( x k ; ω k ) along a step d_k satisfies

D ϕ ( d k ; ω k ) = g k T d k − ω k | | c k | | 1

Proof. From equation (3) and equation (10), we know that c ( x k ) + A ( x k ) T d k = 0 . Using Taylor’s theorem, we have that

ϕ ( x k + α d k ; ω k ) − ϕ ( x k ; ω k ) = f ( x k + α d k ) − f ( x k ) + ω k [ ‖ c ( x k + α d k ) ‖ 1 − ‖ c ( x k ) ‖ 1 ] = α g k T d k + O ( α 2 ‖ d k ‖ 2 ) + ω k [ ‖ c ( x k ) + α A ( x k ) d k ‖ 1 − ‖ c ( x k ) | ‖ 1 ] = α g k T d k + O ( α 2 ‖ d k ‖ 2 ) + ω k [ ‖ ( 1 − α ) c ( x k ) ‖ 1 − ‖ c ( x k ) ‖ 1 ] = α [ g k T d k − ω k ‖ c ( x k ) ‖ 1 ] + O ( α 2 ‖ d k ‖ 2 )

Dividing both sides by α and taking the limit as α → 0 yields the result.

It is easy to prove that, if d_k satisfies the following model reduction condition

Δ m k ( d k ; ω k ) ≥ max ⁡ { 1 2 d k T W k d k , θ ϒ k } + τ ω k | | c k | |

Global convergence

In this section, we will establish the global convergence of the proposed algorithms under the following assumptions.

Assumptions G. The sequences x k , λ k generated by Algorithms are contained in a convex set Ω and the following properties hold:

G1. f and c are bounded and twice continuously differentiable on Ω.

G2. The sequences { W k } and { W k − 1 } are bounded.

G3. The sequence { λ k } is bounded.

G4. A ( x k ) have full row rank and their smallest singular values are bounded below by some positive constant.

Lemma 2. The sequence v_k given by equation (2d) is bounded and satisfies that

| | v k | | ≤ κ 1 | | c k | |

(14)

Proof. From equations (2d), (5a) and (5b), we have that

| | v k | | = | | A k ( A k T A k ) − 1 c k | | ≤ | | A k ( A k T A k ) − 1 | | | | c k | |

(15)

| | v k | | ≤ | | W k − 1 A k ( A k T W k − 1 A k ) − 1 | | | | c k | |

(16)

Assumptions G1, G2 and G4 imply that | | A k ( A k T A k ) − 1 | | and | | W k − 1 A k ( A k T W k − 1 A k ) − 1 | | are bounded, then from equations (15) and (16), the result follows.     □

Lemma 3. Suppose Assumptions G hold, the sequence of parameters ω_k is bounded above and for some k ^ , ω k = ω k ^ when k ≥ k ^ .

Proof. In our algorithm, ω_k is increased only if Condition II is satisfied. It means that

| | r k | | ≤ t 2 | | c k | |

From equation (6), we can obtain that

Δ m k ( d k ; ω k ) − max ⁡ { 1 2 d k T W k d k , θ ϒ k } = ω k | | c k | | + { − g k T d k − 1 2 d k T W k d k , if 1 2 d k T W k d k ≥ θ ϒ k − g k T d k − θ ϒ k , otherwise

(17)

For Algorithms 1 and 4, we have that

u k = u ^ k = − W k − 1 g k + W k − 1 A k λ k + 1 + W k − 1 r k = − W k − 1 g k + W k − 1 A k ( A k T W k − 1 A k ) − 1 A k T W k − 1 g k + W k − 1 r k

(18)

Hence

u k T W k u k = u k T [ − g k + A k ( A k T W k − 1 A k ) − 1 A k T W k − 1 g k + r k ] = − u k T g k + g k T W k − 1 A k ( A k T W k − 1 A k ) − 1 A k T W k − 1 r k + u k T r k

(19)

If 1 2 d k T W k d k < θ ϒ k , then equation (12) and 0 ≤ ν k ≤ | | v k | | 2 yield | | u k | | 2 ≤ ϒ k ≤ ψ ν k ≤ ψ | | v k | | 2 . This implies | | d k | | ≤ 1 + ψ | | v k | | is bounded. From Assumption G1, Lemma 2, there exist γ 1 , γ ¯ 1 > 0 such that

− g k T d k − θ ϒ k ≥ − | | g k | | | | d k | | − θ | | d k | | 2 ≥ − γ 1 ( 1 + ψ ) | | v k | | ≥ − γ ¯ 1 | | c k | |

(20)

If 1 2 d k T W k d k ≥ θ ϒ k , from Assumption G2, there exists some constant γ 2 > 0 such that

| | W k | | γ 2 + | | W k | | 2 γ 2 ≤ θ 2

(21)

If | | u k | | 2 < γ 2 | | v k | | 2 , then as above, we can obtain | | d k | | is bounded, and there exist γ 3 , γ ¯ 3 > 0 such that

− g k T d k − 1 2 d k T W k d k ≥ − γ 3 | | d k | | ≥ − γ 3 ( 1 + γ 2 ) | | v k | | ≥ − γ ¯ 3 | | c k | |

(22)

Otherwise, | | u k | | 2 ≥ γ 2 | | v k | | 2 . By equation (21) and 1 2 d k T W k d k ≥ θ ϒ k , we have

1 2 u k T W k u k ≥ θ ϒ k − u k T W k v k − 1 2 v k T W k v k ≥ ( θ − | | W k | | γ 2 − | | W k | | 2 γ 2 ) | | u k | | 2 ≥ θ 2 | | u k | | 2

(23)

By equations (18) and (23), we have

θ | | u k | | 2 ≤ u k T W k u k = u k T ( − g k + A k ( A k T W k − 1 A k ) − 1 A k T W k − 1 g k + r k ) ≤ ( | | g k | | + | | r k | | + | | A k ( A k T W k − 1 A k ) − 1 A k T W k − 1 g k | | ) | | u k | |

− g k T d k − 1 2 d k T W k d k = − 1 2 v k T W k v k − g k T v k + 1 2 u k T W k u k − u k T W k v k − g k T u k − u k T W k u k = − 1 2 v k T W k v k − g k T v k + 1 2 u k T W k u k − u k T W k v k − [ g k T W k − 1 A k ( A k T W k − 1 A k ) − 1 A k T W k − 1 + u k T ] r k ≥ − 1 2 | | v k | | 2 | | W k | | − | | g k | | | | v k | | − | | u k | | | | W k | | | | v k | | − | | g k T W k − 1 A k × ( A k T W k − 1 A k ) − 1 A k T W k − 1 + u k T | | | | r k | | ≥ − γ 4 | | c k | |

(24)

Δ m k ( d k ; ω k ) − max ⁡ { 1 2 d k T W k d k , θ ϒ k } ≥ ( ω k − max ⁡ { γ ¯ 1 , γ ¯ 3 , γ 4 } ) | | c k | |

Thus, if ω k ^ satisfies ω k ^ ≥ max ⁡ { γ ¯ 1 , γ ¯ 3 , γ 4 } / ( 1 − ϱ ) , then ω k = ω k ^ for all k ≥ k ^ . The Lemma is true.

For Algorithms 2 and 5, we have that

u ^ k = − W k − 1 g k + W k − 1 A k ( A k T A k ) − 1 A k T g k + W k − 1 r k

Then

u k = P k u ^ k = − W k − 1 g k + W k − 1 A k ( A k T W k − 1 A k ) − 1 A k T W k − 1 g k + P k W k − 1 r k

u k T W k u k = − g k T u k + u k T W k P k W k − 1 r k

(25)

If 1 2 d k T W k d k < θ ϒ k , similarly, there exist γ 5 > 0 such that

− g k T d k − θ ϒ k ≥ − γ 5 | | c k | |

(26)

If 1 2 d k T W k d k ≥ θ ϒ k , similarly, if | | u k | | 2 < γ 2 | | v k | | 2 , then there exist γ 6 > 0 such that

− g k T d k − 1 2 d k T W k d k ≥ − γ 6 | | c k | |

(27)

Otherwise, | | u k | | 2 ≥ γ 2 | | v k | | 2 . By equations (21) and 1 2 d k T W k d k ≥ θ ϒ k , we have

1 2 u k T W k u k ≥ θ 2 | | u k | | 2

(28)

By equations (25) and (28), we have

θ | | u k | | 2 ≤ u k T W k u k ≤ ( | | g k | | + | | W k P k W k − 1 | | | | r k | | ) | | u k | |

− g k T d k − 1 2 d k T W k d k = − 1 2 v k T W k v k − g k T v k + 1 2 u k T W k u k − u k T W k v k − g k T u k − u k T W k u k = − 1 2 v k T W k v k − g k T v k + 1 2 u k T W k u k − u k T W k v k − u k T W k P k W k − 1 r k ≥ − 1 2 | | v k | | 2 | | W k | | − | | g k | | | | v k | | − | | u k | | | | W k | | | | v k | | − | | u k T W k P k W k − 1 | | | | r k | | ≥ − γ 7 | | c k | |

(29)

Δ m k ( d k ; ω k ) − max ⁡ { 1 2 d k T W k d k , θ ϒ k } ≥ ( ω k − max ⁡ { γ 5 , γ 6 , γ 7 } ) | | c k | |

Thus, if ω k ^ satisfies ω k ^ ≥ max ⁡ { γ 5 , γ 6 , γ 7 } / ( 1 − ϱ ) , then ω k = ω k ^ for all k ≥ k ^ . The Lemma is true.

For Algorithms 3 and 6, we have that

u ^ k = − W k − 1 g k + W k − 1 A k ( A k T W k − 1 A k ) − 1 A k T W k − 1 g k − W k − 1 A k ( A k T W k − 1 A k ) − 1 c k + W k − 1 r k

Then

u k = P k u ^ k = ( I − A k [ A k T A k − 1 A k T ) u ^ k = − W k − 1 g k + W k − 1 A k ( A k T W k − 1 A k ) − 1 A k T W k − 1 g k − W k − 1 A k ( A k T W k − 1 A k ) − 1 c k + A k ( A k T A k ) − 1 c k + P k W k − 1 r k = u ^ k − W k − 1 r k − v k + P k W k − 1 r k

u k T W k u k = u k T W k u ^ k − u k T r k − u k T W k v k + u k T W k P k W k − 1 r k = − u k T g k − u k T W k v k + u k T W k P k W k − 1 r k

(30)

If 1 2 d k T W k d k < θ ϒ k , similarly, there exist γ 8 > 0 such that

− g k T d k − θ ϒ k ≥ − γ 8 | | c k | |

(31)

If 1 2 d k T W k d k ≥ θ ϒ k , similarly, if | | u k | | 2 < γ 2 | | v k | | 2 , then there exist γ 9 > 0 such that

− g k T d k − 1 2 d k T W k d k ≥ − γ 9 | | c k | |

(32)

Otherwise, | | u k | | 2 ≥ γ 2 | | v k | | 2 . By equations (21) and 1 2 d k T W k d k ≥ θ ϒ k , we have

1 2 u k T W k u k ≥ θ 2 | | u k | | 2

(33)

By equations (30) and (33), we have

θ | | u k | | 2 ≤ u k T W k u k ≤ ( | | g k | | + | | W k v k | | + | | W k P k W k − 1 | | | | r k | | ) | | u k | |

− g k T d k − 1 2 d k T W k d k = − 1 2 v k T W k v k − g k T v k + 1 2 u k T W k u k − u k T W k v k − g k T u k − u k T W k u k = − 1 2 v k T W k v k − g k T v k + 1 2 u k T W k u k − u k T W k P k W k − 1 r k ≥ − 1 2 | | v k | | 2 | | W k | | − | | g k | | | | v k | | − | | u k | | | | W k | | | | W k − 1 | | | | r k | | ≥ − γ 10 | | c k | |

(34)

By equations (17), (31), (32) and (34), we have

Δ m k ( d k ; ω k ) − max ⁡ { 1 2 d k T W k d k , θ ϒ k } ≥ ( ω k − max ⁡ { γ 8 , γ 9 , γ 10 } ) | | c k | |

Thus, if ω k ^ satisfies ω k ^ ≥ max ⁡ { γ 8 , γ 9 , γ 10 } / ( 1 − ϱ ) , then ω k = ω k ^ for all k ≥ k ^ . The Lemma is true.   □

Lemma 4. Suppose Assumptions G hold. There exist constants κ 5 > 0 and κ 6 > 0 such that

| | d k | | 2 + | | c k | | ≤ 2 κ 5 ( ϒ k + | | c k | | )

(35)

D ϕ ( d k ; ω k ) ≤ − κ 6 ( ϒ k + | | c k | | )

Proof. Since | | d k | | 2 = | | u k | | 2 + | | v k | | 2 , from Lemma 2 and the fact that ϒ k is an upper bound for | | u k | | 2 , we have

| | d k | | 2 ≤ κ 5 ( ϒ k + | | c k | | )

(36)

Then the inequality equation (35) follows from equation (36).

From equations (6), (35) and (9), we have

D ϕ ( d k ; ω k ) ≤ − max ⁡ { 1 2 d k T W k d k , θ ϒ k } − τ ω k | | c k | | ≤ − θ ϒ k − τ ω k 2 | | c k | |

Since ω k ≥ ω − 1 , the second inequality holds for κ 6 = min ⁡ { θ , σ ω − 1 2 } > 0 .             □

From the above lemma, we can easily prove that the sequence { α k } is bounded below and away from zero (see Lemma 5 in Byrd et al. ⁹ ). We now state the main result as follows.

Theorem 1. Let { x k } be a sequence generated by the proposed algorithm. Suppose Assumptions G hold, then

lim ⁡ k → ∞ | | c k | | = 0 and lim ⁡ k → ∞ | | ∇ x l ( x k , λ k + 1 N ) | | = 0 ,

Proof. By equation (13) and Lemma 4, there exists a k ^ > 0 such that for all k > k ^

ϕ ( x k ^ ; ω k ^ ) − ϕ ( x k ; ω k ^ ) = ∑ j = k ^ k − 1 ( ϕ ( x j ; ω k ^ ) − ϕ ( x j + 1 ; ω k ^ ) )

≥ − ∑ j = k ^ k − 1 β α j D ϕ ( d j ; ω j ) ≥ ∑ j = k ^ k − 1 β α j κ 6 ( ϒ k + | | c k | | ) ≥ β κ 6 2 κ 5 ∑ j = k ^ k − 1 α j ( | | d k | | 2 + | | c k | | )

This, along with the fact that { α k } is bounded below and Assumption G1, implies

lim ⁡ k → ∞ | | d k | | = 0 , ana lim ⁡ k → ∞ | | c k | | = 0

(37)

Using a way similar to Lemma 4.1 and Lemma 4.2 in Gu and Zhu, ¹¹ we can get

W k d k = − ∇ k l ( x k , λ k + 1 N ) + χ k W k [ W k − 1 A k ( A k T W k − 1 A k ) − 1 − A k ( A k T A k ) − 1 ] c k + W k P k W k − 1 r k

Numerical results

In this section, we report some numerical experiments of the proposed algorithms, which were implemented in MATLAB 7.0. The numerical results have been obtained by running our algorithms on the set of 44 equality constrained problems, which can be found on the web page:http://www.gamsworld.org/performance/princetonlib/htm/group5stat.htm.

In each case, the starting point supplied with the problem was used. All attempts to solve the test problems were limited to a maximum of 1000 iterations or half an hour of CPU time. We set the parameters in the algorithm as follows: β = 0.5 , τ = 0.1 , ς = 0.5 , ω − 1 = 1 , θ = 10 − 8 max ⁡ { | | W k | | 1 , 1 } , ϵ = 10 − 3 and ε = 10 − 5 , which seemed to work reasonably well for a broad class of problems. The Jacobian matrices and the Hessian matrices were approximated by finite differences.