A Conjugate Gradient Type Method for the Nonnegative Constraints Optimization Problems

Abstract

We are concerned with the nonnegative constraints optimization problems. It is well known that the conjugate gradient methods are efficient methods for solving large-scale unconstrained optimization problems due to their simplicity and low storage. Combining the modified Polak-Ribière-Polyak method proposed by Zhang, Zhou, and Li with the Zoutendijk feasible direction method, we proposed a conjugate gradient type method for solving the nonnegative constraints optimization problems. If the current iteration is a feasible point, the direction generated by the proposed method is always a feasible descent direction at the current iteration. Under appropriate conditions, we show that the proposed method is globally convergent. We also present some numerical results to show the efficiency of the proposed method.

1. Introduction

Due to their simplicity and their low memory requirement, the conjugate gradient methods play a very important role for solving unconstrained optimization problems, especially for the large-scale optimization problems. Over the years, many variants of the conjugate gradient method have been proposed, and some are widely used in practice. The key features of the conjugate gradient methods are that they require no matrix storage and are faster than the steepest descent method.

The linear conjugate gradient method was proposed by Hestenes and Stiefel [1] in the 1950s as an iterative method for solving linear systems 𝐴𝑥=𝑏,𝑥∈𝑅^𝑛, (1) where 𝐴 is an 𝑛×𝑛 symmetric positive definite matrix. Problem (1) can be stated equivalently as the following minimization problem 1min2𝑥^𝑇𝐴𝑥−𝑏^𝑇𝑥,𝑥∈𝑅^𝑛. (2) This equivalence allows us to interpret the linear conjugate gradient method either as an algorithm for solving linear systems or as a technique for minimizing convex quadratic functions. For any 𝑥∈𝑅^𝑛, the sequence {𝑥_𝑘} generated by the linear conjugate gradient method converges to the solution 𝑥^∗ of the linear systems (1) in at most 𝑛 steps.

The first nonlinear conjugate gradient method was introduced by Fletcher and Reeves [2] in the 1960s. It is one of the earliest known techniques for solving large-scale nonlinear optimization problems min𝑓(𝑥),𝑥∈𝑅^𝑛, (3) where 𝑓∶𝑅^𝑛→𝑅 is continuously differentiable. The nonlinear conjugate gradient methods for solving (3) have the following form: 𝑥_𝑘+1=𝑥_𝑘+𝛼_𝑘𝑑_𝑘,𝑑_𝑘=−∇𝑓𝑥_𝑘𝑥,𝑘=0,−∇𝑓_𝑘+𝛽_𝑘𝑑_𝑘−1,𝑘≥1, (4) where 𝛼_𝑘 is a steplength obtained by a line search and 𝛽_𝑘 is a scalar which deternimes the different conjugate gradient methods. If we choose 𝑓 to be a strongly convex quadratic and 𝛼_𝑘 to be the exact minimizer, the nonliear conjugate gradient method reduces to the linear conjugate gradient method. Several famous formulae for 𝛽_𝑘 are the Hestenes-Stiefel (HS) [1], Fletcher-Reeves (FR) [2], Polak-Ribière-Polyak (PRP) [3, 4], Conjugate-Descent (CD) [5], Liu-Storey (LS) [6], and Dai-Yuan (DY) [7] formulae, which are given by 𝛽^HS_𝑘=∇𝑓𝑥_𝑘^⊤𝑦_𝑘−1𝑑^⊤_𝑘−1𝑦_𝑘−1,𝛽^FR_𝑘=∇𝑓𝑥_𝑘²∇𝑓𝑥_𝑘−1²𝛽_PRP𝑘,=∇𝑓𝑥_𝑘^⊤𝑦_𝑘−1∇𝑓𝑥_𝑘−1²,𝛽^CD_𝑘𝑥=−∇𝑓_𝑘²𝑑^⊤_𝑘−1∇𝑓𝑥_𝑘−1𝛽_LS𝑘,𝑥=−∇𝑓_𝑘^⊤𝑦_𝑘−1𝑑^⊤_𝑘−1∇𝑓𝑥_𝑘−1,𝛽^DY_𝑘=∇𝑓𝑥_𝑘²𝑑^⊤_𝑘−1𝑦_𝑘−1, (5)(6)(7) where 𝑦_𝑘−1=∇𝑓(𝑥_𝑘)−∇𝑓(𝑥_𝑘−1) and ‖⋅‖ stands for the Euclidean norm of vectors. In this paper, we focus our attention on the Polak-Ribière-Polyak (PRP) method. The study of the PRP method has received much attention and has made good progress. The global convergence of the PRP method with exact line search has been proved in [3] under strong convexity assumption on 𝑓. However, for general nonlinear function, an example given by Powell [8] shows that the PRP method may fail to be globally convergent even if the exact line search is used. Inspired by Powell’s work, Gilbert and Nocedal [9] conducted an elegant analysis and showed that the PRP method is globally convergent if 𝛽^PRP_𝑘 is restricted to be nonnegative and 𝛼_𝑘 is determined by a line search satisfying the sufficient descent condition 𝑔^⊤_𝑘𝑑_𝑘≤−𝑐‖𝑔_𝑘‖² in addition to the standard Wolfe conditions. Other conjugate gradient methods and their global convergence can be found in [10–15] and so forth.

Recently, Li and Wang [16] extended the modified Fletcher-Reeves (MFR) method proposed by Zhang et al. [17] for solving unconstrained optimization to the nonlinear equations 𝐹(𝑥)=0, (8) where 𝐹∶𝑅^𝑛→𝑅^𝑛 is continuously differentiable, and proposed a descent derivative-free method for solving symmetric nonlinear equations. The direction generated by the method is descent for the residual function. Under appropriate conditions, the method is globally convergent by the use of some backtracking line search technique.

In this paper, we further study the conjugate gradient method. We focus our attention on the modified Polak-Ribière-Polyak (MPRP) method proposed by Zhang et al. [18]. The direction generated by MPRP method is given by 𝑑_𝑘=−𝑔𝑥_𝑘𝑥,𝑘=0,−𝑔_𝑘+𝛽^PRP_𝑘𝑑_𝑘−1−𝜃_𝑘𝑦_𝑘−1,𝑘>0, (9) where 𝑔(𝑥_𝑘)=∇𝑓(𝑥_𝑘),  𝛽^PRP_𝑘=𝑔(𝑥_𝑘)^𝑇𝑦_𝑘−1/‖𝑔(𝑥_𝑘−1)‖²,  𝜃_𝑘=𝑔(𝑥_𝑘)^𝑇𝑑_𝑘−1/‖𝑔(𝑥_𝑘−1)‖²,  and  𝑦_𝑘−1=𝑔(𝑥_𝑘)−𝑔(𝑥_𝑘−1). The MPRP method not only reserves good properties of the PRP method but also possesses another nice property; that it is, always generates descent directions for the objective function. This property is independent of the line search used. Under suitable conditions, the MPRP method with the Armoji-type line search is also globally convergent. The purpose of this paper is to develop an MPRP type method for the nonnegative constraints optimization problems. Combining the Zoutendijk feasible direction method with MPRP method, we propose a conjugate gradient type method for solving the nonnegative constraints optimization problems. If the initial point is feasible, the method generates a feasible point sequence. We also do numerical experiments to test the proposed method and compare the performance of the method with the Zoutendijk feasible direction method. The numerical results show that the method that we propose outperforms the Zoutendijk feasible direction method.

2. Algorithm

Consider the following nonnegative constraints optimization problems: min𝑓(𝑥)s.t.𝑥≥0, (10) where 𝑓∶𝑅^𝑛→𝑅 is continuously differentiable. Let 𝑥_𝑘≥0 be the current iteration. Define the index set 𝐼_𝑘=𝐼𝑥_𝑘=𝑖∣𝑥_𝑘(𝑖)=0,𝐽_𝑘={1,2,…,𝑛}⧵𝐼_𝑘, (11) where 𝑥_𝑘(𝑖) is the  𝑖th component of 𝑥_𝑘. In fact the index set 𝐼_𝑘 is the active set of problem (10) at 𝑥_𝑘.

The purpose of this paper is to develop a conjugate gradient type method for problem (10). Since the iterative sequence is a feasible point sequence, the search directions should be feasible descent directions. Let 𝑥_𝑘≥0 be the current iteration. By the definition of feasible direction, we have that [19] 𝑑∈𝑅^𝑛 is a feasible direction of (10) at 𝑥_𝑘 if and only if 𝑑_𝐼𝑘≥0. Similar to the Zoutendijk feasible direction method, we consider the following problem: min∇𝑓𝑥_𝑘^𝑇𝑑s.t.𝑑_𝐼𝑘≥0,‖𝑑‖≤1. (12) Next, we show that, if 𝑥_𝑘 is not a KKT point of (10), the solution of problem (12) is a feasible descent direction of 𝑓 at 𝑥_𝑘.

Lemma 1. Let 𝑥_𝑘≥0 and let 𝑑 be a solution of problem (12); then ∇𝑓(𝑥_𝑘)^𝑇𝑑≤0. Moreover ∇𝑓(𝑥_𝑘)^𝑇𝑑=0 if and only if 𝑥_𝑘 is a KKT point of problem (10).

Proof. Since 𝑑=0 is a feasible point of problem (12), there must be ∇𝑓(𝑥_𝑘)^𝑇𝑑≤0. Consequently, if ∇𝑓(𝑥_𝑘)^𝑇𝑑≠0, there must be ∇𝑓(𝑥_𝑘)^𝑇𝑑<0. This implies that the direction 𝑑 is a feasible descent direction of 𝑓 at 𝑥_𝑘.
We suppose that ∇𝑓(𝑥_𝑘)^𝑇𝑑=0. Problem (12) is equivalent to the following problem: min∇𝑓𝑥_𝑘^𝑇𝑑s.t.𝑑_𝐼𝑘≥0,‖𝑑‖²≤1. (13) Then there exist 𝜆_𝐼𝑘 and 𝜇 such that the following KKT condition holds: ∇𝑓𝑥_𝑘−𝜆_𝐼𝑘0+2𝜇𝜆𝑑=0,_𝐼𝑘≥0,𝑑_𝐼𝑘≥0,𝜆^𝑇_𝐼𝑘𝑑_𝐼𝑘=0,𝜇≥0,𝑑≤1,𝜇𝑑²−1=0. (14) Multiplying the first of these expressions by 𝑑, we obtain ∇𝑓𝑥_𝑘^𝑇𝑑−𝜆^𝑇𝑑+2𝜇𝑑²=0, (15) where 𝜆=^𝜆𝐼𝑘₀. By combining the assumption ∇𝑓(𝑥_𝑘)^𝑇𝑑=0 with the second and the third expressions of (14), we find that 𝜇=0. Substituting it into the first expressions of (14), we obtain that ∇𝑓_𝐼𝑘𝑥_𝑘−𝜆_𝐼𝑘=0,∇𝑓_𝐽𝑘𝑥_𝑘=0. (16) Let 𝜆_𝑖=0,  𝑖∈𝐽_𝑘; then 𝜆_𝑖≥0,  𝑖∈𝐼_𝑘∪𝐽_𝑘. Moreover, we have ∇𝑓𝑥_𝑘−𝜆_𝐼𝑘𝜆_𝐽𝑘𝜆=0,_𝑖≥0,𝑥_𝑘(𝑖)≥0,𝜆_𝑖𝑥_𝑘(𝑖)=0,𝑖∈𝐼_𝑘∪𝐽_𝑘. (17)
This implies that 𝑥_𝑘 is a KKT point of problem (10).
On the other hand, we suppose that 𝑥_𝑘 is a KKT point of problem (10). Then there exist 𝜆_𝑖,𝑖∈𝐼_𝑘∪𝐽_𝑘, such that the following KKT condition holds: ∇𝑓𝑥_𝑘−𝜆_𝐼𝑘𝜆_𝐽𝑘𝜆=0,_𝑖≥0,𝑥_𝑘(𝑖)≥0,𝜆_𝑖𝑥_𝑘(𝑖)=0,𝑖∈𝐼_𝑘∪𝐽_𝑘. (18) From the second of these expressions, we get 𝜆_𝐽𝑘=0. Substituting it into the first of these expressions, we have ∇𝑓_𝐼𝑘(𝑥_𝑘)=𝜆_𝐼𝑘≥0 and ∇𝑓_𝐽𝑘(𝑥_𝑘)=0, so that ∇𝑓(𝑥_𝑘)^𝑇𝑑=∇𝑓_𝐼𝑘(𝑥_𝑘)^𝑇𝑑_𝐼𝑘=𝜆^𝑇_𝐼𝑘𝑑_𝐼𝑘≥0. However, we had shown that ∇𝑓(𝑥_𝑘)^𝑇𝑑≤0, so ∇𝑓(𝑥_𝑘)^𝑇𝑑=0.
By the proof of Lemma 1 we find that ∇𝑓_𝐼𝑘(𝑥_𝑘)≥0 and ∇𝑓_𝐽𝑘(𝑥_𝑘)=0 are necessary conditions of the fact that 𝑥_𝑘 is a KKT point of problem (10). We summarize these observation results as the following result.

Lemma 2. Let 𝑥_𝑘≥0; then 𝑥_𝑘 is a KKT point of problem (10) if and only if ∇𝑓_𝐼𝑘(𝑥_𝑘)≥0 and ∇𝑓_𝐽𝑘(𝑥_𝑘)=0.

Proof. Firstly, we suppose that 𝑥_𝑘 is a KKT point of problem (10). Similar to the proof of Lemma 1, it is easy to get that ∇𝑓_𝐼𝑘(𝑥_𝑘)≥0 and ∇𝑓_𝐽𝑘(𝑥_𝑘)=0.
Secondly, we suppose that ∇𝑓_𝐼𝑘(𝑥_𝑘)≥0 and ∇𝑓_𝐽𝑘(𝑥_𝑘)=0. Let 𝜆_𝐼𝑘=∇𝑓_𝐼𝑘(𝑥_𝑘)≥0,  𝜆_𝐽𝑘=0; then the KKT condition (18) holds, so that 𝑥_𝑘 is a KKT point of problem (10).
Based on the above discussion, we propose a conjugate gradient type method for solving problem (10) as follows. Let feasible point 𝑥_𝑘 be current iteration. For the boundary of the feasible region 𝑥_𝑘𝐼𝑘=0, we take 𝑑_𝑘𝑖=0,𝑔_𝑖𝑥_𝑘>0,−𝑔_𝑖𝑥_𝑘,𝑔_𝑖𝑥_𝑘≤0,∀𝑖∈𝐼_𝑘, (19) where 𝑔_𝑖(𝑥_𝑘)=∇𝑓_𝑖(𝑥_𝑘). For the interior of the feasible region 𝑥_𝑘𝐽𝑘>0, similar to the direction 𝑑_𝑘 in the MPRP method, we define 𝑑_𝑘𝐽𝑘 by the following formula: 𝑑^MPRP_𝑘𝐽𝑘=−𝑔_𝐽𝑘𝑥_𝑘,𝑘=0,−𝑔_𝐽𝑘𝑥_𝑘+𝛽^PRP_𝑘𝑑_{𝑘−1𝐽𝑘}−𝜃^MPRP_𝑘𝑦_𝑘−1,𝑘>0, (20) where 𝑔_𝐽𝑘(𝑥_𝑘)=∇𝑓_𝐽𝑘(𝑥_𝑘),  𝛽^PRP_𝑘=𝑔_𝐽𝑘(𝑥_𝑘)^𝑇𝑦_𝑘−1/‖𝑔(𝑥_𝑘−1)‖²,  𝜃^MPRP_𝑘=𝑔_𝐽𝑘(𝑥_𝑘)^𝑇𝑑_{𝑘−1𝐽𝑘}/‖𝑔(𝑥_𝑘−1)‖²,  and  𝑦_𝑘−1=𝑔_𝐽𝑘(𝑥_𝑘)−𝑔_𝐽𝑘(𝑥_𝑘−1).
It is easy to see from (19) and (20) that −𝑔_𝐼𝑘𝑥_𝑘²≤𝑔_𝐼𝑘𝑥_𝑘^𝑇𝑑_𝑘𝐼𝑘𝑔≤0,_𝐽𝑘𝑥_𝑘^𝑇𝑑_𝑘𝐽𝑘=−𝑔_𝐽𝑘(𝑥_𝑘)². (21) The above relations indicate that 𝑔𝑥_𝑘^𝑇𝑑_𝑘=𝑔_𝐼𝑘𝑥_𝑘^𝑇𝑑_𝑘𝐼𝑘+𝑔_𝐽𝑘𝑥_𝑘^𝑇𝑑_𝑘𝐽𝑘≤−𝑔_𝐽𝑘𝑥_𝑘²,𝑔𝑥_𝑘^𝑇𝑑_𝑘≥−𝑔_𝐼𝑘(𝑥_𝑘)²−𝑔_𝐽𝑘(𝑥_𝑘)²=𝑔(𝑥_𝑘)², (22)(23) where 𝑔(𝑥_𝑘)=∇𝑓(𝑥_𝑘).

Theorem 3. Let 𝑥_𝑘≥0, 𝑑_𝑘 be defined by (19) and (20) then 𝑔𝑥_𝑘^𝑇𝑑_𝑘≤0. (24) Moreover, 𝑥_𝑘 is a KKT point of problem (10) if and only if 𝑔(𝑥_𝑘)^𝑇𝑑_𝑘=0.

Proof. Clearly, inequality (22) implies that 𝑔𝑥_𝑘^𝑇𝑑_𝑘≤0. (25)
If 𝑥_𝑘 is a KKT point of problem (10), similar to the proof of Lemma 1, we also get that 𝑔(𝑥_𝑘)^𝑇𝑑_𝑘=0.
If 𝑔(𝑥_𝑘)^𝑇𝑑_𝑘=0, by (22), we can get that 𝑔_𝐼𝑘𝑥_𝑘^𝑇𝑑_𝑘𝐼𝑘𝑔=0,_𝐽𝑘𝑥_𝑘^𝑇𝑑_𝑘𝐽𝑘=−𝑔_𝐽𝑘(𝑥_𝑘)²=0. (26) The equality 𝑔_𝐼𝑘(𝑥_𝑘)^𝑇𝑑_𝑘𝐼𝑘=0 and the definition of 𝑑_𝑘𝐼𝑘 (19) imply that 𝑔_𝐼𝑘𝑥_𝑘≥0. (27) Let 𝜆_𝐼𝑘=𝑔_𝐼𝑘(𝑥_𝑘)≥0;  𝜆_𝐽𝑘=0, then the KKT condition (18) also holds, so that 𝑥_𝑘 is a KKT point of problem (10).
By combining (22) with Theorem 3, we conclude that 𝑑_𝑘 defined by (19) and (20) provides a feasible descent direction of 𝑓 at 𝑥_𝑘, if 𝑥_𝑘 is not a KKT point of problem (10).

Based on the above process, we propose an MPRP type method for solving (10) as follows.

Algorithm 4 (MPRP type method).   
Step  0. Given constants 𝜌∈(0,1),  𝛿>0,  𝜖>0. Choose the initial point 𝑥₀≥0; Let 𝑘∶=0.
Step  1. Compute 𝑑_𝑘=(𝑑_𝑘𝐼𝑘,𝑑_𝑘𝐽𝑘) by (19) and (20). If |𝑔(𝑥_𝑘)^𝑇𝑑_𝑘|≤𝜖, then stop. Otherwise, go to the next step.
Step  2. Determine 𝛼_𝑘=max{𝜌^𝑗,𝑗=0,1,2,…} satisfying 𝑥_𝑘+𝛼_𝑘𝑑_𝑘≥0 and 𝑓𝑥_𝑘+𝛼_𝑘𝑑_𝑘≤𝑓𝑥_𝑘−𝛿𝛼²_𝑘𝑑_𝑘². (28)
Step  3. Let the next iteration be 𝑥_𝑘+1=𝑥_𝑘+𝛼_𝑘𝑑_𝑘.
Step  4. Let 𝑘∶=𝑘+1 and go to Step  1.

It is easy to see that the sequence {𝑥_𝑘} generated by Algorithm 4 is a feasible point sequence. Moreover, it follows from (28) that the function value sequence {𝑓(𝑥_𝑘)} is decreasing. In addition if 𝑓(𝑥) is bounded from below, we have from (28) that ^∞_𝑘=0𝛼²_𝑘𝑑_𝑘²<∞. (29) In particular we have lim_𝑘→∞𝛼_𝑘𝑑_𝑘=0. (30)

Next, we prove the global convergence of Algorithm 4 under the following assumptions.

Assumption A. (1) The level set 𝜔={𝑥∈𝑅^𝑛∣𝑓(𝑥)≤𝑓(𝑥₀)} is bound.
(2) In some neighborhood 𝑁 of 𝜔, 𝑓 is continuously differentiable, and its gradient is the Lipschitz continuous; namely, there exists a constant 𝐿>0 such that ∇𝑓(𝑥)−∇𝑓𝑦≤𝐿𝑥−𝑦,∀𝑥,𝑦∈𝑁. (31)

Clearly, Assumption A implies that there exists a constant 𝛾₁ such that ∇𝑓(𝑥)≤𝛾₁,∀𝑥∈𝑁. (32)

Lemma 5. Suppose that the conditions in Assumption A hold; {𝑥_𝑘} and {𝑑_𝑘} are the iterative sequence and the direction sequence generated by Algorithm 4. If there exists a constant 𝜖>0 such that 𝑔𝑥_𝑘≥𝜖,∀𝑘, (33) then there exists a constant 𝑀>0 such that 𝑑_𝑘≤𝑀,∀𝑘. (34)

Proof. By combining (19), (20), and (33) with Assumption A, we deduce that 𝑑_𝑘≤𝑑_𝑘𝐼𝑘+𝑑^MPRP_𝑘𝐽𝑘≤𝛾₁+𝑔_𝐽𝑘𝑥_𝑘𝑔+2^𝐽𝑘𝑥_𝑘𝑦_𝑘−1𝑑^MPRP_{𝑘−1𝐽𝑘}𝑔𝑥_𝑘−1²≤2𝛾₁+2𝛾₁𝐿𝛼_𝑘−1𝑑^MPRP_{𝑘−1𝐽𝑘}𝜖²𝑑^MPRP_{𝑘−1𝐽𝑘}. (35) By (30), there exists a constant 𝛾∈(0,1) and an iteger 𝑘₀ such that the following inequality holds for all 𝑘≥𝑘₀: 2𝐿𝛾₁𝜖²𝛼_𝑘−1𝑑^MPRP_{𝑘−1𝐽𝑘}≤𝛾. (36) Hence, we have for any 𝑘≥𝑘₀𝑑_𝑘≤2𝛾₁+𝛾𝑑_𝑘−1≤2𝛾₁1+𝛾+𝛾²+⋯+𝛾^{𝑘−𝑘0−1}+𝛾^𝑘−𝑘0𝑑_𝑘0≤2𝛾₁+𝑑1−𝛾^𝑘0. (37)
Let 𝑀=max𝑑₁,𝑑₂,…,𝑑_𝑘0,2𝛾₁+𝑑1−𝛾^𝑘0. (38) Then 𝑑_𝑘≤𝑀,∀𝑘. (39)

Theorem 6. Suppose that the conditions in Assumption A hold. Let {𝑥_𝑘} and {𝑑_𝑘} be the iterative sequence and the direction sequence generated by Algorithm 4. Then liminf_𝑘→∞𝑔𝑥_𝑘^𝑇𝑑_𝑘=0. (40)

Proof. We prove the result of this theorem by contradiction. Assume that the theorem is not true; then there exists a constant 𝜀>0 such that 𝑔𝑥_𝑘^𝑇𝑑_𝑘≥𝜖,∀𝑘. (41) So by combining (41) with (23), it is easy to see that (33) holds. (1) If liminf_𝑘→∞𝛼_𝑘>0, we get from (30) that 𝑑_𝑘→0, so that lim_𝑘→∞|𝑔(𝑥_𝑘)^𝑇𝑑_𝑘|=0. This contradicts assumption (41). (2) If liminf_𝑘→∞𝛼_𝑘=0, there is an infinite index set 𝐾 such that lim_{𝑘∈𝐾,𝑘→∞}𝛼_𝑘=0. (42)It follows from Step  2 of Algorithm 4, that when 𝑘∈𝐾 is sufficiently large, 𝜌⁻¹𝛼_𝑘 does not satify 𝑓(𝑥_𝑘+𝛼_𝑘𝑑_𝑘)≤𝑓(𝑥_𝑘)−𝛿𝛼²_𝑘‖𝑑_𝑘‖²; that is 𝑓𝑥_𝑘+𝜌⁻¹𝛼_𝑘𝑑_𝑘−𝑓𝑥_𝑘>−𝛿𝜌⁻²𝛼²_𝑘𝑑_𝑘². (43) By the mean-value theorem, Lemma 1, and Assumption A, there is ℎ_𝑘∈(0,1) such that 𝑓𝑥_𝑘+𝜌⁻¹𝛼_𝑘𝑑_𝑘−𝑓𝑥_𝑘=𝜌⁻¹𝛼_𝑘𝑔𝑥_𝑘+ℎ_𝑘𝜌⁻¹𝛼_𝑘𝑑_𝑘^𝑇𝑑_𝑘=𝜌⁻¹𝛼_𝑘𝑔𝑥_𝑘^𝑇𝑑_𝑘+𝜌⁻¹𝛼_𝑘𝑔𝑥_𝑘+ℎ_𝑘𝜌⁻¹𝛼_𝑘𝑑_𝑘−𝑔𝑥_𝑘^𝑇𝑑_𝑘≤𝜌⁻¹𝛼_𝑘𝑔𝑥_𝑘^𝑇𝑑_𝑘+𝐿𝜌⁻²𝛼²_𝑘𝑑_𝑘². (44) Substituting the last inequality into (43), we get for all 𝑘∈𝐾 sufficiently large 0≤−𝑔𝑥_𝑘^𝑇𝑑_𝑘≤𝜌⁻¹(𝐿+𝛿)𝛼_𝑘𝑑_𝑘². (45) Taking the limit on both sides of the equation, then by combining ‖𝑑_𝑘‖≤𝑀 and recalling lim_{𝑘∈𝐾,𝑘→∞}𝛼_𝑘=0, we obtain that lim_{𝑘∈𝐾,𝑘→∞}|𝑔(𝑥_𝑘)^𝑇𝑑_𝑘|=0. This also yields a contradiction.

3. Numerical Experiments

In this section, we report some numerical experiments. We test the performance of Algorithm 4 and compare it with the Zoutendijk method.

The code was written in Matlab, and the program was run on a PC with 2.20 GHz CPU and 1.00 GB memory. The parameters in the method are specified as follows. We set 𝜌=1/2,  𝛿=1/10. We stop the iteration if |∇𝑓(𝑥_𝑘)^𝑇𝑑_𝑘|≤0.0001 or the iteration number exceeds 10000.

We first test Algorithm 4 on small and medium size problems and compared them with the Zoutendijk method in the total number of iterations and the CPU time used. The test problems are from the CUTE library [20]. The numerical results of Algorithm 4 and the Zoutendijk method are listed in Table 1. The columns have the following meanings.

Table 1 

The numerical results.

𝑃(𝑖) is the number of the test problem, Dim is the dimension of the test problem, Iter is the number of  iterations, and Time is CPU time in seconds.

We can see from Table 1 that Algorithm 4 has successfully solved 12 test problems, and the Zoutendijk method has successfully solved 8 test problems. From the number of iterations, Algorithm 4 has 12 test results better than Zoutendijk method. From the computation time, Algorithm 4 performs much better than the Zoutendijk method did. We then test Algorithm 4 and the Zoutendijk method on two problems with a larger dimension. The problem of VARDIM comes from [20], and the following problem comes from [16]. The results are listed in Tables 2 and 3.

Table 2 

Test results for VARDIM with various dimensions.

Table 3 

Test results for Problem 1 with various dimensions.

Problem 1. The nonnegative constraints optimization problem min𝑓(𝑥)s.t.𝑥≥0, (46) with Engval function 𝑓∶𝑅^𝑛→𝑅 is defined by 𝑓(𝑥)=^𝑛_𝑖=2𝑥²_𝑖−1+𝑥²_𝑖²−4𝑥_𝑖−1+3. (47)

We can see from Table 2 that Algorithm 4 has successfully solved the problem of VARDIM whose scale varies from 1000 dimensions to 5000 dimensions. However, the Zoutendijk method fails to solve the problem of VARDIM with larger dimension. From Table 3, although the number of iterations of Algorithm 4 is more than the Zoutendijk method, the computation time of Algorithm 4 is less than the Zoutendijk method, and this feature becomes more evident as increase of the dimension of the test problem.

In summary, the results from Tables 1–3 show that Algorithm 4 is more efficient than the Zoutendijk method and provides an efficient method for solving nonnegative constraints optimization problems.

Acknowledgment

This research is supported by the NSF (11161020) of China.